Interpretation of Stiffness Matrix and Mass Matrix in Finite Element Method

I would like to have a general interpretation of the coefficients of the stiffness matrix that appears in FEM. For instance if we are solving a linear elasticity problem and we modelize the relation between a node $i$ and a node $j$ as a spring system, then $K_{i,j}$ (where $K$ is the stiffness matrix of the system) can be seen as the stiffness constant of the virtual spring between the two nodes. But does there exist a more general interpretation? Perhaps in terms of internal work?

Another similar question is: What could be an interpretation of the coefficients $M_{i,j}$ of the mass matrix $M$?

Let us consider that we are dealing with linear elasticity. In that case, the Lagrangian of the system, without considering body forces, is given by

$$L = \int\limits_\Omega \sigma_{ij} \varepsilon_{ij}\, \mathrm{d}\Omega - \frac{1}{2}\int\limits_\Omega \rho \dot{u}_i \dot{u}_i\, \mathrm{d}\Omega\, ,$$

where $$\dot{u}_i$$ are the components of the velocity vector, $$\sigma_{ij}$$ are the components of the stress tensor, $$\varepsilon_{ij}$$ are the components of the strain tensor, $$\rho$$ is the mass density, and we have used Einstein notation.

Using a displacement-based Finite Element formulation we have and interpolation/approximation of the displacements ($$U$$) given by

$$U = \sum_{j} U^{j} N^{j}(x)\, .$$

We also have the following relations

$$\varepsilon_{ij} = \frac{1}{2}(u_{i,j} + u_{j, i})\, ,$$

$$\sigma_{ij} = C_{ijkl} \varepsilon_{kl}\, .$$

When replacing the (approximated) displacements in the Lagrangian, we obtain

$$L_\text{approx} = \frac{1}{2}\{U^T\}[K]\{U\} - \frac{1}{2}\{\dot{U}^T\}[M]\{\dot{U}\}\, .$$

Thus, we can think the components $$K_{ij}$$ as the elastic energy related to the relative deformation between degrees of freedom $$U_i$$ and $$U_j$$. Similarly, we can think the components $$M_{ij}$$ as the kinetic energy related to relative velocities between $$\dot{U}_i$$ and $$\dot{U}_j$$.

As a final coment, this interpretation is for elasticity. Nevertheless, matrices $$K$$ and $$M$$ are commonly termed stifness and mass matrices even for other applications. In other applications they can have different interpretations, for example, in Quantum mechanics, $$K$$ would be the Hamiltonian matrix and $$M$$ the overlap matrix.

I can not see a good geometrical interpretation of the stiffness matrix components, but there is a problem if we try to see them as spring constants between nodes.

The simplest case is an isotropic solid loaded in the range where linear elasticity applies. Linear elasticity means that the stress tensor is proportional to the strain tensor, and for an isotropic solid, there are only 2 elastic constants ($$\lambda$$ and $$\mu$$), relating the components of that tensors:

For any point inside the solid, we can always choose axis x, y and z in order to have only normal stresses, and the relation between that stresses and strain components are:

$$\sigma_{xx}=(\lambda + 2*\mu)*\epsilon_{xx} + \lambda*\epsilon_{yy}+\lambda*\epsilon_{zz}$$ $$\sigma_{yy}=\lambda*\epsilon_{xx} + (\lambda + 2*\mu)*\epsilon_{yy}+\lambda*\epsilon_{zz}$$ $$\sigma_{zz}=\lambda*\epsilon_{xx} + \lambda*\epsilon_{yy}+(\lambda + 2*\mu)*\epsilon_{zz}$$

It is easy to verify that, when only one of that stresses is different from zero, there is a linear relation between the remaining non zero stress and its related strain, (ex: $$\sigma_{xx}$$ and $$\epsilon_{xx}$$), and the constant of proportionality is the Young's modulus given by:

$$E = \mu(3\lambda + 2\mu) / (\lambda + \mu)$$

So, if we model a linear elastic solid as a collection of springs that link nodes, we are taking only that particular case when two of the principal stresses are zero everywhere. It is only true for very simple situations, where finite elements is not necessary.

Modelling the solid with finite elements, and considering only one of that elements, the equation for the force in the x-direction at the node #1 is:

$$\sum_{i=1}^{ne}{(u_{xi}*(\lambda + 2*\mu)*\int_{V}\frac{\partial {N_1}}{\partial x}*\frac{\partial {N_i}}{\partial x}} \ dV) +\sum_{i=1}^{ne}({u_{yi}*(\lambda + \mu)*\int_{V}(\frac{\partial {N_1}}{\partial y}*\frac{\partial {N_i}}{\partial y} + \frac{\partial {N_1}}{\partial x}*\frac{\partial {N_i}}{\partial x}}) \ dV)\sum_{i=1}^{ne}({u_{zi}*(\lambda + \mu)*\int_{V}(\frac{\partial {N_1}}{\partial z}*\frac{\partial {N_i}}{\partial z} + \frac{\partial {N_1}}{\partial x}*\frac{\partial {N_i}}{\partial x}}) \ dV) = F_{1x}$$

where $$u_{xi}$$, $$u_{yi}$$, $$u_{zi}$$ are the displacements of each node for each direction, $${N_i}$$ are interpolation functions for each node, (that depends on the element type), and ne is the number of nodes of the element.

As can be seen, the stiffness matrix components (all the stuff multiplying the displacements) are the product of the elastic constants by an integral involving derivatives of the interpolation functions.