5. Construction of the matrix system

We will now present the various ingredients leading to the construction of the matrix system that will solve the problem.

5.1. New notation (Voigt notation)

In order to fully understand the construction of discrete terms in the finite element method, we will use a more compact notation:

\(\langle V\rangle\) is a line vector

\(\left\{V\right\}\) is a column vector

\(\left[A\right]\) is a matrix

Thus geometric interpolation is written according to the three dimensions of space:

(5.1)\[{x}^{e}={x}_{1}^{e}=\langle {x}_{\mathrm{1,}i}^{e}\rangle \mathrm{.}\left\{{N}_{i}^{e}\right\}=\langle {N}_{i}^{e}\rangle \mathrm{.}\left\{{x}_{\mathrm{1,}i}^{e}\right\} {y}^{e}={y}_{1}^{e}=\langle {y}_{\mathrm{1,}i}^{e}\rangle \mathrm{.}\left\{{N}_{i}^{e}\right\}=\langle {N}_{i}^{e}\rangle \mathrm{.}\left\{{y}_{\mathrm{1,}i}^{e}\right\} {z}^{e}={z}_{1}^{e}=\langle {z}_{\mathrm{1,}i}^{e}\rangle \mathrm{.}\left\{{N}_{i}^{e}\right\}=\langle {N}_{i}^{e}\rangle \mathrm{.}\left\{{z}_{\mathrm{1,}i}^{e}\right\}\]

Or even more compactly in vector form:

(5.2)\[\left\{{x}^{e}\right\}=\left[{N}_{i}^{e}\right]\mathrm{.}\left\{{x}_{i}^{e}\right\}=\langle {x}_{i}^{e}\rangle \mathrm{.}{\left[{N}_{i}^{e}\right]}^{T}\]

With the \(N\) matrix of shape functions. By considering an element with two nodes, we obtain in expanded form:

(5.3)\[\begin{split}\left\{\begin{array}{c}{x}^{e}\\ {y}^{e}\\ {z}^{e}\end{array}\right\}=\left[\begin{array}{cccccc}{N}_{1}^{e}& 0& 0& {N}_{2}^{e}& 0& 0\\ 0& {N}_{1}^{e}& 0& 0& {N}_{2}^{e}& 0\\ 0& 0& {N}_{1}^{e}& 0& 0& {N}_{2}^{e}\end{array}\right]\mathrm{.}\left\{\begin{array}{c}{x}_{1}^{e}\\ {y}_{1}^{e}\\ {z}_{1}^{e}\\ {x}_{2}^{e}\\ {y}_{2}^{e}\\ {z}_{2}^{e}\end{array}\right\}\end{split}\]

5.2. Discretized system

In the hyperelastic case, we place ourselves in small deformations, the mechanical problem to be solved is written in a more compact way:

(5.4)\[\textrm{Trouver} u\mathrm{\in }{E}^{h} \textrm{avec} a(u,\tilde{u})+l(\tilde{u})\mathrm{=}0\]

With \(a(u,\tilde{u})\) a bilinear, symmetric form that represents the potential energy of the structure and \(l(\tilde{u})\) the potential [1] volume and surface forces:

(5.5)\[\begin{split}\begin{array}{}a(u,\tilde{u})=\underset{{\Omega }^{h}}{\int }\varepsilon (\tilde{u}):\sigma (u)\mathrm{.}d{\Omega }^{h}\\ l(\tilde{u})=\underset{{\Omega }^{h}}{\int }f\mathrm{.}\tilde{u}\mathrm{.}d{\Omega }^{h}+\underset{{\Gamma }_{N}^{h}}{\int }g\mathrm{.}\tilde{u}\mathrm{.}d{\Gamma }^{h}\end{array}\end{split}\]

Discretization consists in choosing a base from space \({\Omega }^{h}\) and numerically calculating the terms of the matrix \(A\) and the vector \(L\). To do this, we express the bilinear form \(a(\mathrm{.},\mathrm{.})\) and the linear form \(l(\mathrm{.})\) as a sum over elements, defined by division of the base domain:

(5.6)\[\begin{split}\mathrm{\{}\begin{array}{cc}a({u}_{i},{\tilde{u}}_{j})& \text{=}\mathrm{\sum }_{\text{éléments}{\Omega }^{e}}\underset{{\Omega }^{e}}{\mathrm{\int }}{\sigma }_{\mathit{kl}}({u}_{i})\mathrm{.}{\varepsilon }_{\mathit{kl}}({\tilde{u}}_{j})\mathrm{.}d{\Omega }^{e}\\ l({\tilde{u}}_{i})& \text{=}\mathrm{\sum }_{\text{éléments}{\Omega }^{e}}\underset{{\Omega }^{e}}{\mathrm{\int }}{f}_{i}\mathrm{.}{\tilde{u}}_{i}\mathrm{.}d{\Omega }^{e}+\underset{{\Gamma }_{N}^{e}}{\mathrm{\int }}{g}_{i}\mathrm{.}{\tilde{u}}_{i}\mathrm{.}d{\Gamma }_{N}^{e}\end{array}\end{split}\]

The terms \({A}_{\mathrm{ij}}\), which represent the interaction between two degrees of freedom \(i\) and \(j\), are constructed by assembling (the operation noted \(\sum _{\text{éléments}{\Omega }^{e}}(\mathrm{...})\)) the contributions coming from each of the elements that contain the corresponding nodes; we proceed in the same way to build the second member vector \({L}_{i}\). These contributions, called elementary terms, are calculated during a loop over the elements and only depend on the variables of element \({\Omega }^{e}\):

(5.7)\[\begin{split}\{\begin{array}{cc}{a}^{e}& \text{=}\underset{{\Omega }^{e}}{\int }{\sigma }_{\mathrm{kl}}\mathrm{.}{\varepsilon }_{\mathrm{kl}}\mathrm{.}d{\Omega }^{e}\\ {l}^{e}& \text{=}\underset{{\Omega }^{e}}{\int }{f}_{i}\mathrm{.}{w}_{i}\mathrm{.}d{\Omega }^{e}+\underset{{\Gamma }_{N}^{e}}{\int }{g}_{i}\mathrm{.}{w}_{i}\mathrm{.}d{\Gamma }_{N}^{e}\end{array}\end{split}\]

The relationship between the Cauchy stress tensor \(\sigma\) and the displacements \(u\) is given by the behavioral relationship, and is independent of the writing of the variational formulation. In the elastic case, we have:

(5.8)\[{\sigma }_{\mathrm{ij}}({w}_{i})={\Lambda }_{\mathrm{ijkl}}\mathrm{.}{\varepsilon }_{\mathrm{ku}}({w}_{i})\]

\({\Lambda }_{\mathrm{ijkl}}\) is Hooke’s elasticity tensor. This tensor form is not very practical; the Voigt notation is preferentially used, which makes it possible to write:

(5.9)\[\sigma :\varepsilon =\langle \sigma \rangle \mathrm{.}\{\varepsilon \}\]

In Cartesian coordinates, we have:

(5.10)\[\mathrm{\langle }\sigma \mathrm{\rangle }\mathrm{=}\langle \begin{array}{cccccc}{\sigma }_{\mathit{xx}}& {\sigma }_{\mathit{yy}}& {\sigma }_{\mathit{zz}}& {\sigma }_{\mathit{xy}}& {\sigma }_{\mathit{xz}}& {\sigma }_{\mathit{yz}}\end{array}\rangle\]

And a modified form of the deformation components to allow the contracted product to be expressed, i.e.:

(5.11)\[\mathrm{\langle }\varepsilon \mathrm{\rangle }\mathrm{=}\langle \begin{array}{cccccc}{\varepsilon }_{\mathit{xx}}& {\varepsilon }_{\mathit{yy}}& {\varepsilon }_{\mathit{zz}}& 2.{\varepsilon }_{\mathit{xy}}& 2.{\varepsilon }_{\mathit{xz}}& 2.{\varepsilon }_{\mathit{yz}}\end{array}\rangle\]

Important note:

In the integration of laws of behavior, the stress shear and deformation components used by Code_Aster*are:*

\(\mathrm{\langle }\sigma \mathrm{\rangle }\mathrm{=}\langle \begin{array}{cccccc}{\sigma }_{\mathit{xx}}& {\sigma }_{\mathit{yy}}& {\sigma }_{\mathit{zz}}& \sqrt{2}\mathrm{.}{\sigma }_{\mathit{xy}}& \sqrt{2}\mathrm{.}{\sigma }_{\mathit{xz}}& \sqrt{2}\mathrm{.}{\sigma }_{\mathit{yz}}\end{array}\rangle\)

\(\mathrm{\langle }\varepsilon \mathrm{\rangle }\mathrm{=}\langle \begin{array}{cccccc}{\varepsilon }_{\mathit{xx}}& {\varepsilon }_{\mathit{yy}}& {\varepsilon }_{\mathit{zz}}& \sqrt{2}\mathrm{.}{\varepsilon }_{\mathit{xy}}& \sqrt{2}\mathrm{.}{\varepsilon }_{\mathit{xz}}& \sqrt{2}\mathrm{.}{\varepsilon }_{\mathit{yz}}\end{array}\rangle\)

The product of these two vectors gives the same result as the double contracted product ( ) .

With this new notation, we have in elasticity:

(5.12)\[\{\sigma \}=[A]\mathrm{.}\{\varepsilon \}\]

We start from the EF writing of the travel field:

(5.13)\[\left\{{u}^{e}\right\}=\left[{N}_{i}^{e}\right]\mathrm{.}\left\{{u}_{i}^{e}\right\}=\langle {u}_{i}^{e}\rangle \mathrm{.}{\left[{N}_{i}^{e}\right]}^{T}\]

And, similarly, the field of virtual travel:

(5.14)\[\left\{{\tilde{u}}^{e}\right\}=\left[{N}_{i}^{e}\right]\mathrm{.}\left\{{\tilde{u}}_{i}^{e}\right\}=\langle {\tilde{u}}_{i}^{e}\rangle \mathrm{.}{\left[{N}_{i}^{e}\right]}^{T}\]

For the sake of simplifying the notations, the reference to the element will be omitted. It is necessary to express the tensor of the deformations (virtual or real):

(5.15)\[\left\{\varepsilon \right\}=\left[B\right]\mathrm{.}\left\{u\right\}=\langle u\rangle \mathrm{.}{\left[B\right]}^{T}\]

For the matrix relating to the bilinear form, we then obtain:

(5.16)\[\begin{array}{cc}[a]& \text{=}\langle \tilde{u}\rangle \mathrm{.}\underset{{\Omega }^{e}}{\int }{\left[B\right]}^{T}\mathrm{.}[\Lambda ]\mathrm{.}\left[B\right]\mathrm{.}d{\Omega }^{e}\mathrm{.}\left\{u\right\}\end{array}\]

The matrices \(\left[B\right]\) and \([\Lambda ]\) contain the possible non-linearity of the behavior and will depend on the displacements:

1. \(\left[B\right]\) is a function of movements if we are in the situation of large deformations or large transformations (large rotations and/or large displacements).

2. \([\Lambda ]\) is the behavior matrix. It becomes dependent on displacements (and other variables) in the case of non-linear and/or non-elastic behaviors.

In these two cases, the process of solving the equations will involve a specific treatment (most often, a linearization of the Newton-Raphson type). In a similar manner, the elementary form will easily be obtained for the second member.

5.3. Calculation of elementary terms

The basic terms to be calculated are of the form:

(5.17)\[\underset{{\Omega }_{e}}{\int }f(u(x),\frac{\partial u(x)}{\partial x})\mathrm{.}dx\]

Three types of operations are to be performed:

1. the transformation of derivatives with respect to \(x\) into derivatives with respect to \(\xi\);

2. the transition from integration on the real element to integration on the reference element,

3. the numerical realization of this integration which is generally done by a quadrature formula.

5.3.1. Transformation of derivatives

The transformation of derivatives is carried out using the Jacobian matrix \(J\), according to the chain rule of derivation:

(5.18)\[\frac{\partial {u}_{\alpha }}{\partial {x}_{\beta }}=\frac{\partial {\xi }_{\gamma }}{\partial {x}_{\beta }}\mathrm{.}\frac{\partial {u}_{\alpha }}{\partial {\xi }_{\gamma }}={J}^{\text{-1}}\mathrm{.}{\left[\frac{\partial N}{\partial \xi }\right]}^{T}\mathrm{.}{u}_{\alpha }^{\text{nod}}\]

where \({u}_{\alpha }^{\text{nod}}\) is the vector of the nodal values of the \(\alpha\) component of the displacement.

5.3.2. Change of integration domain

The transition to integration on the reference element is carried out by multiplying the integrand by the determinant of the Jacobian matrix, called Jacobian:

(5.19)\[\underset{{\Omega }_{e}}{\int }f(u(x),\frac{\partial u(x)}{\partial x})\mathrm{.}dx=\underset{{\Omega }_{r}}{\int }f(u(\xi ),\frac{\partial u(\xi )}{\partial \xi })\mathrm{.}\text{det}(J(\xi ))\mathrm{.}d\xi\]

The transition from the reference element to the real element implies the bijectivity of the \(\tau\) transformation. We must therefore \(\text{det}(J(\xi ))\ne 0\), which implies that the element must not turn around or degenerate (for example, the quadrangle must not degenerate into a triangle).

5.3.3. Digital integration

5.3.3.1. In some particular cases, integrals can be calculated analytically. For example, for a triangle in two dimensions, the Jacobian is constant on the triangle, and the integrands are reduced to monomials that we know how to integrate exactly thanks to the numerical integration formula known as « Gauss » [2]

(5.20)\[\underset{0}{\overset{1}{\int }}\underset{0}{\overset{1-\xi }{\int }}{\xi }_{1}^{\alpha }\mathrm{.}{\xi }_{2}^{\beta }\mathrm{.}d{\xi }_{1}\mathrm{.}d{\xi }_{2}=\frac{\alpha !\beta !}{(\alpha +\beta +2)!}\]

However, these particular cases are rare, and integrals are preferred to be evaluated numerically using quadrature formulas. These give an approximation of the integral in the form of a weighted sum of the values of the integrand at a certain number of points in the element called integration points:

(5.21)\[\underset{{\Omega }_{r}}{\mathrm{\int }}g(\xi )\mathrm{.}d\xi \mathrm{\approx }\mathrm{\sum }_{g\mathrm{=}1}^{r}{\omega }_{g}\mathrm{.}g({\xi }_{g})\]

Scalars \({\omega }_{g}\) are called integration weights, and coordinates \({\xi }_{g}\) are the coordinates of the \(r\) integration points in the reference element.

In Gauss integration methods, integration points and weights are determined in such a way as to exactly integrate polynomials of a given order. This type of method is used in Code_Aster, the integration points are then called Gauss points.

The number of Gauss points chosen makes it possible to integrate exactly into the reference element. In fact, because of the possible non-linearity of the geometric transformation or the spatial dependence of the coefficients (for example for deformed or second-order elements), the integration is not exact in the real element.

For each element \({\Omega }_{e}\), we were able to calculate the so-called elementary terms: elementary matrix \({\mathrm{A}}^{\mathrm{e}}\) and elementary vector \({L}^{e}\). The matrix \(A\) and the vector \(L\) are obtained by a procedure called the assembly of elementary terms.

If we go back to the elementary form of stiffness:

(5.22)\[\mathrm{[}a\mathrm{]}\text{=}\underset{{\Omega }^{e}}{\mathrm{\int }}\mathrm{\{}\sigma \mathrm{\}}\mathrm{.}\mathrm{\langle }\varepsilon \mathrm{\rangle }\mathrm{.}d{\Omega }^{e}\]

Digital integration involves evaluating the constraints and deformations at the integration points:

(5.23)\[[a]\text{=}\underset{{\Omega }^{e}}{\int }\{\sigma \}\mathrm{.}\langle \varepsilon \rangle \mathrm{.}d{\Omega }^{e}\approx \sum _{g=1}^{r}{\omega }_{g}\mathrm{.}\{\sigma ({\xi }_{g})\}\mathrm{.}\langle \varepsilon ({\xi }_{g})\rangle\]

This means that the constraints and deformations are the most exact (or the least false) at the integration points (fields called « ELGA » in Code_Aster). Simply extrapolating these values to the nodes for display introduces an error. In fact, it is a method for evaluating error, called the Zhu-Zienkiewicz error indicator.

In 2D elasticity, a triangle exhibiting a constant Jacobian, a single Gauss point is sufficient to exactly integrate the terms of the matrix and the second member (if it is constant).

The calculation cost increases with the number of integration points, especially for non-linear laws of behavior. For example, a hexahedron with 27 knots needs 27 Gauss points to integrate quantities. It therefore often happens that we « under-integrate », that is to say that we use fewer integration points than the minimum required, thus making an error that will eventually be compensated for by a finer mesh. In addition to this systematic error, this under-integration must be done carefully because it can produce matrix rank defects and thus make the linear system non-invertible.

[1]

The potential for external efforts does not depend on the movement of the structure, it’s what we call dead or non-tracking a loading. In the case of large deformations, loads of pressure type cannot respect this hypothesis.

[2]

For abuse of language, numerical integration schemes are frequently called « Gauss schemes » although there are several types of such schemes (Hammer for the triangle, Gauss-Radau, Newton-Cotes, etc.).