A1. General formulation for linear elastic calculations

A1. A general formulation for linear elastic calculations

The two main formulations for linear calculations are displacements-based and stress-based. Most of the software nowadays rely on a displacement formulation because it can easily consider the boundary conditions, unlike a stress-based formulation. Mixed formulations, using displacements and stresses also exist. There are other formulations, but they will not be discussed in this chapter.

In this document, only the displacement formulation (the most commonly used) is introduced: in this context, the basic unknowns are the displacements at the nodes, typically represented by q, from which all the other values are deduced (strains ε, stresses σ  ....)

The displacement at any point of an element, represented by u(x,y,z), are obtained by interpolation between the values at the nodes using interpolation functions (known as shape functions) represented by N(x,y,z), which are low order polynomials (generally of degree 1 to 3):

The strain-displacement relation is obtained by computing the derivative using the differential operator (matrix) D (different according to the situation,) from which the matrix B appears, connecting strains and displacements to the nodes:

To obtain the stresses, it is essential to introduce the constitutive law of the material. Thus, in linear elasticity, the strain-stress relation is written (H representing Hooke’s matrix):

One can observe that all the interesting quantities (displacement at any point, strains, stresses) can be obtained from the displacements at the nodes of the mesh.

Note: as the deformations (and consequently the strains) are obtained by computing the derivative of the displacements at the nodes, there is a loss in accuracy when calculating those values.

For a static calculation in mechanics of structures, identifying the numerical solution of the problem by the Finite Elements Method (FEM) can be summed up in 4 main steps:

Step 1. Determining the elementary matrices and vectors of each element of volume V by applying the following equations (with BT representing the transposed of B):

 elementary stiffness matrix

 elementary vector of equivalent loads 

with f^V the volumetric load vector and f^S the surface load vector. The vector of equivalent loads allows us to “bring back to the nodes” the loads applied to the elements (in the volume or the surface,) according to the chosen interpolation.

Be careful! Usually, it does not correspond to an even distribution at the nodes of the total load applied on an element. 

Step 2: Determining the stiffness matrix K and the force vector F with respect to the global coordinate system by assembling the elementary quantities and introducing the boundary conditions.

Step 3: Determining the vector of nodal displacements q by solving the linear system (large scale):


The vector q gives the node displacements of the mesh, likewise, the vector F contains the values of the equivalent loads. 

Step 4: Resolution of the quantities of interest such as strains and stresses by post-processing the vector of nodal displacements q.

Going back to the elements, the displacement, the strain, and the stress values at any point can be deducted from the displacements at the nodes by using the previous equations. As opposed to the displacements, the stresses and the strains are calculated element by element (and these quantities are not always continuous from one element to another.)

Note: it is possible to make an analogy between solving, using finite elements, a mechanical problem, and a thermal problem in a linear permanent regime: in this case, the temperatures T at the nodes are obtained by solving a linear system: Λ. T = Φ, with Λ representing the conductivity matrix and Φ the flux vector at the nodes.  

The order of the shape functions N(x,y,z), which is deducted for the solid elements (or continuum mechanics – CM) from the number of nodes in the edge of a finite element (2 nodes: linear element, 3 nodes: quadratic element) has an impact on the degree of strains and stresses evaluated using finite elements. For example, if the differential operator D corresponds to the first-order derivative (which is the case in CM,) the strains and stresses will be constant along with the element for linear elements, and linear along with the element for quadratic elements… As a reminder: in Strength of Materials (SM) to compute the strains we make use of the second-order derivatives (bending of beams, plates….)

Shape functions – Triangular element with 3 nodes (linear functions)

Shape functions – Triangular element with 6 nodes (quadratic functions)

Be careful! It is not possible to mix in the same mesh finite element with shape functions of a different order (for example linear and quadratic elements.) Actually, the continuity of strains cannot be ensured at the interface between elements with shape functions of a different order.

In the software, most finite elements are defined with the same parameters (isoparametric,) meaning that the different finite elements of a mesh are built as a geometric transformation of a reference element (which is itself expressed in unitary coordinate system ξ,η,χ.)

Isoparametric element (quadratic): different transformations (2D)

These geometric transformations, which generally rely on the same nodes used for the interpolation of the displacements (hence “isoparametric”, with the same parameters,) are characterized by their Jacobian matrix J (matrix of the geometric transformation’s partial derivatives of the first order.)  

Isoparametric element (quadratic): local and global coordinate systems (2D)

One of the advantages of this transformation is the possibility of approximating the boundaries of a given domain by using polynomial geometries. However, the main problem is that the calculation of elementary matrices and vectors can become more complex. Their evaluation is made by numerical integration using a quadrature scheme: weighted summation of the values at a certain point located inside the elements (weights ω_i, points, ξ_i) particular to the method (Gauss method, Hammer method….) In this case, the calculation of the elementary stiffness matrix and the elementary vector of equivalent loads are written:

When the matrices (by applying the B matrix) are evaluated this way, it is obvious that one must calculate the strains and stresses at the same points (also as a function of B.) The software provides these quantities at the integration points (Gauss or Hammer points, according to the integration scheme.) Moreover, talking about the values at the nodes of a mesh for strains and stresses does not make sense because there are as many values in a node as the number of elements connected to it. Be careful, there is not necessarily a continuity between the elements! Wanting at all costs a nodal value (of strain or stress) of the mesh implies some data processing: extrapolating from the values at the integration points and the average values of different elements sharing the node, which will, at the very least, smoothen the results (nonetheless, talking about values at the nodes in a particular element is justified.) When the quantities are obtained by analytical integration (in the case of beam elements for instance,) they can be expressed at the nodes, at the center of gravity… the previous observations about the continuity between elements being still valid.

Quadratic elements and integration points (Gauss points)

The way to consider the displacements at the boundary conditions (restraint, support…) might lead to different results depending on the adopted numerical technique.  There are two ways of considering the displacements at the boundaries:

• Penalization method: it is simple to implement but it has the limitation of being sensitive to the order of magnitude of the stiffness matrix terms.

• Lagrange multipliers method: it is not affected by the above-mentioned accuracy problems, but slightly increases the size of the system to be solved.    

Note: in dynamics, the displacements at any point x,y,z at any moment t is obtained, as in statics, by interpolating the values at the nodes: ux,y,z=Nx,y,z.q(t). These nodes verify the equilibrium condition, given at any time t by:



The mass matrix can be obtained by assembling the elementary mass matrices M_V (volumetric mass) and C the damping matrix (viscous). To compute the aforementioned matrix, a system of differential equations needs to be solved using adequate numerical methods (see the chapter on dynamics.)