Displacement method

The displacement method is the standard formulation of the finite element method (FEM), in which the displacements of the body points are the primary unknowns. In solid mechanics , the displacements include the paths that the body points cover over time and thus the translation, rotation and possibly deformation of a solid . The first application of FEM was the linear treatment of solids and structures (consisting of rods , beams or shells ) and it was based on this that the FEM received its impulses.

One equation on which the displacement method is based is the principle of d'Alembert in the Lagrange version . With this principle, both linear problems, such as the question of natural vibrations, and highly non-linear problems, such as crash tests , can be analyzed. Because they can be used in most problems, isoparametric elements are used in the standard formulation . The Galerkin method is also used.

The displacement method is available in all common finite element programs with which problems in solid mechanics can be calculated, whereby the programs can differ in the strain dimensions used, implemented nonlinearities, material models, time integration methods and / or numerical conversions.

Matrix equations

The principle of d'Alembert in the Lagrange version is a statement equivalent to the momentum balance and reads

${\displaystyle \int _{V}\delta {\vec {u}}\cdot \rho _{0}{\ddot {\vec {u}}}\,{\text{d}}V+\int _{V}{\tilde {\mathbf {T} }}:\delta \mathbf {E} \,{\text{d}}V=\int _{{A}^{\sigma }}\delta {\vec {u}}\cdot {\vec {t}}_{0}\,{\text{d}}A+\int _{V}\delta {\vec {u}}\cdot \rho _{0}{\vec {k}}_{0}\,{\text{d}}V\quad \forall \;\delta {\vec {u}}\in {\mathcal {V}}\,.}$

The first term is the virtual work of the change in momentum on the virtual displacements summed up over the volume V of the body. The factor is the density and the acceleration of the material point. In the second term, the virtual deformation work performed by the stress tensor on the virtual strain tensor, which is formed with the Frobenius scalar product “:” of the tensors, is summed over the volume V of the body. On the right side is the work of the external forces (surface and volume distribution) on the virtual displacements. The set contains all permissible virtual displacements that vanish wherever the displacement boundary conditions are given. The area is that part of the surface of the body on which there are no displacement boundary conditions. ${\displaystyle \rho _{0}{\ddot {\vec {u}}}}$ ${\displaystyle \delta {\vec {u}}\,.}$ ${\displaystyle \rho _{0}}$ ${\displaystyle {\ddot {\vec {u}}}}$ ${\displaystyle {\tilde {\mathbf {T} }}}$ ${\displaystyle \delta \mathbf {E} }$ ${\displaystyle {\mathcal {V}}}$ ${\displaystyle {A}^{\sigma }}$

If the equation above is actually fulfilled for all permitted virtual displacements, then the displacements and the resulting distortions and stresses are in accordance with the momentum balance. The assumed symmetry of the stress tensor also ensures that the angular momentum balance is fulfilled .

Conversion of the principle into a matrix equation

Semicircle (blue) and its FE model (yellow).

In preparation, the body of interest is divided into sub-bodies, the finite elements, see picture. In order for the elements to build up the body without gaps and without overlapping, neighboring elements must be compatible . Points called nodes are located on the boundary surfaces or inside the elements , to which global coordinates for the geometric description of the body and displacements are assigned as node variables to describe the movement.

Using an element, the tensor equation of the principle is converted into a matrix equation:

• The element displacements are represented with a purely time-dependent solution vector of length m with nodal displacements and a purely location-dependent 3 × m shape function matrix N: ${\displaystyle {\hat {u}}}$

${\displaystyle {\vec {u}}=N{\hat {u}}\,.}$

The element accelerations are thus The shape functions are the initial functions for the solution and the restriction to a finite number of them results in the discretization error which is unavoidable in the absence of an analytical solution. ${\displaystyle {\ddot {\vec {u}}}=N{\hat {\ddot {u}}}\,.}$

• In the Galerkin method, the virtual displacements are treated in the same way as the point displacements:

${\displaystyle \delta {\vec {u}}=N\delta {\hat {u}}\,.}$

• The stress tensor is written as a column vector using Voigt's notation :

${\displaystyle {\hat {T}}={\begin{pmatrix}{\tilde {T}}_{xx}&{\tilde {T}}_{yy}&{\tilde {T}}_{zz}&{\tilde {T}}_{xy}&{\tilde {T}}_{yz}&{\tilde {T}}_{xz}\end{pmatrix}}^{\top }\,.}$

The superscript denotes the transposition . ${\displaystyle (\cdot )^{\top }}$

• The virtual distortions are entered into a vector in the same way:

${\displaystyle \delta {\hat {E}}={\begin{pmatrix}\delta E_{xx}&\delta E_{yy}&\delta E_{zz}&2\delta E_{xy}&2\delta E_{yz}&2\delta E_{xz}\end{pmatrix}}^{\top }\,.}$

The factor two at the fourth to sixth position ensures that the scalar product agrees with the matrix product. In addition, the double shear distortion corresponds to the slip γ, so that the components have a clear interpretation. ${\displaystyle {\tilde {\mathbf {T} }}:\delta \mathbf {E} =\delta {\hat {E}}^{\top }{\hat {T}}\,.}$

• The differential virtual distortions are linear in the virtual nodal displacements:

${\displaystyle \delta {\hat {E}}=B\delta {\hat {u}}}$

with the 6 × m distortion displacement matrix B.

These definitions, which are carried out in the section #Element matrices at the integration points below, are entered into the principle:

${\displaystyle {\begin{aligned}\int _{V}\delta {\hat {u}}^{\top }N^{\top }\rho _{0}N{\hat {\ddot {u}}}\,{\text{d}}V+\int _{V}\delta {\hat {E}}^{\top }{\hat {T}}\,{\text{d}}V=\int _{{A}^{\sigma }}\delta {\hat {u}}^{\top }N^{\top }{\vec {t}}_{0}\,{\text{d}}A+\int _{V}\delta {\hat {u}}^{\top }N^{\top }\rho _{0}{\vec {k}}_{0}\,{\text{d}}V&\quad \forall \;\delta {\hat {u}}\in V^{m}\\\rightarrow \delta {\hat {u}}^{\top }{\Biggl [}\underbrace {\int _{V}N^{\top }\rho _{0}N\,{\text{d}}V} _{M}{\hat {\ddot {u}}}+\underbrace {\int _{V}B^{\top }{\hat {T}}\,{\text{d}}V} _{\hat {r}}-\underbrace {\left(\int _{{A}^{\sigma }}N^{\top }{\vec {t}}_{0}\,{\text{d}}A+\int _{V}N^{\top }\rho _{0}{\vec {k}}_{0}\,{\text{d}}V\right)} _{\hat {f}}{\Biggr ]}=0&\quad \forall \;\delta {\hat {u}}\in V^{m}\,.\end{aligned}}}$

The set contains all admissible virtual displacement vectors of dimension m × 1. The nodal variables can be extracted from the integrals because the position dependency is solely due to the shape functions. This is how the principle equation arises: ${\displaystyle V^{m}}$

${\displaystyle \delta {\hat {u}}^{\top }[M{\hat {\ddot {u}}}+{\hat {r}}-{\hat {f}}]=0\quad \forall \;\delta {\hat {u}}\in V^{m}}$		(PvdA)

The m × m matrix M is the constant mass matrix, the m × 1 vector contains nodal reactions due to element stresses and is the nodal force vector of the same size, which represents the external forces. ${\displaystyle {\hat {r}}}$ ${\displaystyle {\hat {f}}}$

This equation was set up for a finite element, which is permissible because the principle applies to every body and every part of it. During the transition from the finite element to the whole body, the element matrices are summed up in global matrices in an assembly step. The global system is also subject to the principle equation (PvdA), only the matrices and vectors are larger.

Partitioning the global system

The set contains all admissible virtual displacement vectors of dimension m × 1. A virtual displacement vector is permissible if its components vanish where boundary conditions are given in the solution vector. The components of the solution vector that are bound to displacement constraints are moved to the end of the solution vector: ${\displaystyle V^{m}}$

${\displaystyle {\hat {u}}={\begin{pmatrix}{\hat {u}}_{u}\\{\hat {u}}_{b}\end{pmatrix}}\quad \rightarrow \quad \delta {\hat {u}}={\begin{pmatrix}\delta {\hat {u}}_{u}\\{\hat {0}}\end{pmatrix}}\,,\quad {\hat {f}}={\begin{pmatrix}{\hat {f}}_{u}\\{\hat {0}}\end{pmatrix}}}$

The vector contains the unknown nodal displacements sought and are specified in boundary conditions. In the vector of external forces , the lower part consists of zeros, because no forces can be specified at locations where displacement boundary conditions exist. The mass matrix is ​​also partitioned: ${\displaystyle {\hat {u}}_{u}}$ ${\displaystyle {\hat {u}}_{b},\,{\hat {f}}_{u}}$ ${\displaystyle {\hat {f}}}$

${\displaystyle M={\begin{pmatrix}M_{uu}&M_{ub}\\M_{bu}&M_{bb}\end{pmatrix}}\,.}$

This results from the equation (PvdA):

${\displaystyle {\begin{aligned}{\begin{pmatrix}\delta {\hat {u}}_{u}^{\top }&{\hat {0}}^{\top }\end{pmatrix}}\left[{\begin{pmatrix}M_{uu}&M_{ub}\\M_{bu}&M_{bb}\end{pmatrix}}{\begin{pmatrix}{\hat {\ddot {u}}}_{u}\\{\hat {\ddot {u}}}_{b}\end{pmatrix}}+{\begin{pmatrix}{\hat {r}}_{u}\\{\hat {r}}_{b}\end{pmatrix}}-{\begin{pmatrix}{\hat {f}}_{u}\\{\hat {0}}\end{pmatrix}}\right]=&\\\delta {\hat {u}}_{u}^{\top }\left(M_{uu}{\hat {\ddot {u}}}_{u}+M_{ub}{\hat {\ddot {u}}}_{b}+{\hat {r}}_{u}-{\hat {f}}_{u}\right)=&0\quad \forall \;\delta {\hat {u}}_{u}\in \mathbb {R} ^{m_{u}}\,,\end{aligned}}}$

if m _{u is} the dimension of the vector . The above equation forces the sum in the parentheses to disappear with the result: ${\displaystyle {\hat {u}}_{u}}$

${\displaystyle M_{uu}{\hat {\ddot {u}}}_{u}={\hat {f}}_{u}-{\hat {r}}_{u}-M_{ub}{\hat {\ddot {u}}}_{b}\,.}$

As soon as the displacements have been completely determined, it is usual - if desired - to calculate the reactions at the bearings from the dropped out lower equation line, for example:

${\displaystyle {\hat {r}}_{b}=-M_{bu}{\hat {\ddot {u}}}_{u}-M_{bb}{\hat {\ddot {u}}}_{b}\,.}$

Simplifications adopted in the text

For the sake of simplicity, only static fixed bearings in which the displacements and accelerations are zero are considered below. Then these shifts and accelerations do not occur in the equation system and a partitioning into unknown and known parts can be dispensed with. The vector with the unknowns still has the dimension m. Then the equation of motion follows directly from the equation (PvdA): ${\displaystyle {\hat {u}}}$ ${\displaystyle \delta {\hat {u}}^{\top }[M{\hat {\ddot {u}}}+{\hat {r}}-{\hat {f}}]=0\quad \forall \;\delta {\hat {u}}\in \mathbb {R} ^{m}}$

${\displaystyle M{\hat {\ddot {u}}}={\hat {f}}-{\hat {r}}\,.}$		(*)

Furthermore, in the following, any dependency on the speeds, as is the case with viscous materials or speed-dependent forces, should be negligible for the sake of simplicity.

Displacements and accelerations

In the matrix equation (*) the nodal accelerations appear , but in the displacement method the displacements are the primary unknowns. The aim here is to clarify how the accelerations are determined by default from the displacements. ${\displaystyle M{\hat {\ddot {u}}}={\hat {f}}-{\hat {r}}}$ ${\displaystyle {\hat {\ddot {u}}}}$

The accelerations are the second derivative of the displacements with respect to time and therefore the accelerations and displacements are not independent of one another, rather one can be calculated from the other via time integration or derivation. Usually, one-step methods are used for time integration , in which the displacements and accelerations are calculated from currently known values. Two time integration methods are particularly widespread: ${\displaystyle t^{n+1}}$ ${\displaystyle t^{n}}$

1. The implicit Newmark beta method and
2. the explicit time integration.

In the Newmark beta method, the shifts are to be used here as primary unknowns, which according to

${\displaystyle {\hat {\ddot {u}}}^{n+1}=p{\hat {u}}^{n+1}+{\hat {p}}^{n}}$		(NbV)

Cause accelerations. The parameter comes from the integration algorithm and the vector is currently known. In the explicit time integration is ${\displaystyle p}$ ${\displaystyle {\hat {p}}^{n}}$ ${\displaystyle t^{n}}$

${\displaystyle {\hat {u}}^{n+1}=q{\hat {\ddot {u}}}^{n}+{\hat {q}}^{n}\,.}$		(EZI)

The scalar is also a parameter of the integration algorithm and the vector is currently known. ${\displaystyle q}$ ${\displaystyle {\hat {q}}^{n}}$ ${\displaystyle t^{n}}$

These results are needed when solving the matrix equations.

Applications

Linear case

In the linear case the tensions hanging (but by assumption linearly from the nodal displacements not on the speeds) from: ${\displaystyle {\hat {T}}}$

${\displaystyle {\hat {T}}={\hat {T}}^{0}+C_{L}{\hat {E}}_{L}={\hat {T}}^{0}+C_{L}B_{L}{\hat {u}}\,.}$

The vector contains internal stresses and the material matrix is linear , independent of the distortions and thus of the nodal displacements. The distortions are also linear in the displacements, which also applies to the virtual distortions. The reactions are thus calculated ${\displaystyle {\hat {T}}^{0}}$ ${\displaystyle C_{L}}$ ${\displaystyle ({\hat {E}}_{L}=B_{L}{\hat {u}})\,,}$

${\displaystyle {\hat {r}}=\int _{V}B_{L}^{\top }{\hat {T}}\,{\text{d}}V=\underbrace {\int _{V}B_{L}^{\top }{\hat {T}}^{0}\,{\text{d}}V} _{{\hat {r}}_{L}^{0}}+\underbrace {\int _{V}B_{L}^{\top }C_{L}B_{L}\,{\text{d}}V} _{K_{L}}{\hat {u}}={\hat {r}}_{L}^{0}+K_{L}{\hat {u}}}$

The linear stiffness matrix K _L is independent of the displacements. Insertion into the equation of motion (*) yields the equation for the nodal displacements at a time t ^{n + 1} : ${\displaystyle M{\hat {\ddot {u}}}={\hat {f}}-{\hat {r}}}$

${\displaystyle \rightarrow M{\hat {\ddot {u}}}^{n+1}+K_{L}{\hat {u}}^{n+1}={\hat {f}}^{n+1}-{\hat {r}}_{L}^{0}\,.}$

Free oscillation of a beam clamped on one side

This linear equation with a mass matrix and stiffness matrix independent of the distortions can be transferred into the modal space, with a result like the one in the picture. ${\displaystyle M}$ ${\displaystyle K_{L}}$

By expressing the accelerations with the displacements according to the Newmark beta method , the displacement can be calculated in each time step: ${\displaystyle {\hat {\ddot {u}}}^{n+1}=p{\hat {u}}^{n+1}+{\hat {p}}^{n}}$

${\displaystyle (pM+K_{L}){\hat {u}}^{n+1}={\hat {f}}^{n+1}-{\hat {r}}_{L}^{0}-M{\hat {p}}^{n}\,.}$

If only the equilibrium position is sought, the displacements result from ${\displaystyle {\hat {\ddot {u}}}={\hat {0}}}$

${\displaystyle K_{L}{\hat {u}}^{n+1}={\hat {f}}^{n+1}-{\hat {r}}_{L}^{0}\,.}$

Even if a linear dependence of the external forces on the nodal displacements would simply have to be taken into account in the linear case , this is normally not used. ${\displaystyle {\hat {f}}}$

Implicit solution to nonlinear problems

Geometrically non-linear analysis of a bulging, elastoplastic carrier

In the nonlinear case, the nodal reaction in the equation of motion (*) depends nonlinearly on the nodal displacements . Causes of the non-linearity can be: ${\displaystyle {\hat {r}}}$ ${\displaystyle M{\hat {\ddot {u}}}={\hat {f}}-{\hat {r}}}$ ${\displaystyle {\hat {u}}}$

Material nonlinearity

Plasticity , nonlinear elasticity .

Displacement-dependent boundary conditions

Forces depending on the deformation, body contact .

Geometric non-linearity

Large twists or deformations, kinks , bumps .

The animation shows a geometrically non-linear, implicit analysis of a bulging, elastoplastic beam.

In order to determine the nodal displacements, Newton's method is used as standard , which provides for a linearization of the equation. Linearization of the reactions delivers in the point ${\displaystyle {\hat {r}}}$ ${\displaystyle {\hat {u}}^{i}\,:}$

${\displaystyle {\begin{aligned}{\hat {r}}({\hat {u}}^{i}+\Delta {\hat {u}})=&\int _{V}B^{\top }({\hat {u}}^{i}+\Delta {\hat {u}}){\hat {T}}({\hat {u}}^{i}+\Delta {\hat {u}})\mathrm {d} V\\[-2ex]\approx &\underbrace {\int _{V}B^{i\top }{\hat {T}}^{i}\mathrm {d} V} _{{\hat {r}}^{i}}+\underbrace {\int _{V}{\frac {\mathrm {d} }{\mathrm {d} {\hat {u}}}}\left.(B^{i\top }{\hat {T}}^{i})\right|_{{\hat {T}}^{i}={\text{const.}}}\mathrm {d} V} _{G^{i}}\Delta {\hat {u}}+\underbrace {\int _{V}B^{i\top }\overbrace {\frac {\mathrm {d} {\hat {T}}}{\mathrm {d} {\hat {E}}}} ^{C}\overbrace {\frac {\mathrm {d} {\hat {E}}}{\mathrm {d} {\hat {u}}}} ^{B^{i}}\mathrm {d} V} _{K^{i}}\Delta {\hat {u}}\\=&{\hat {r}}^{i}+(G^{i}+K^{i})\Delta {\hat {u}}\,.\end{aligned}}}$

The sizes marked with the superscript can depend on the displacements . The B matrix only does this in the geometrically non-linear case and only in this case does the geometrical stiffness matrix have to be established. ${\displaystyle (\cdot )^{i}}$ ${\displaystyle {\hat {u}}^{i}}$ ${\displaystyle G^{i}}$

In the case of displacement-dependent external forces, the nodal force vector is also linearized ${\displaystyle {\hat {f}}}$

${\displaystyle {\hat {f}}({\hat {u}}^{i}+\Delta {\hat {u}})\approx {\hat {f}}({\hat {u}}^{i})+{\frac {\mathrm {d} {\hat {f}}}{\mathrm {d} {\hat {u}}}}\Delta {\hat {u}}=:{\hat {f}}^{i}+F^{i}\Delta {\hat {u}}\,,}$

which leads to an m × m matrix F ⁱ . In dynamic systems, the shift increment also causes an acceleration increment: Inserting these results into the equation of motion (*) yields with the abbreviation ${\displaystyle \Delta {\hat {\ddot {u}}}=p\Delta {\hat {u}}\,.}$ ${\displaystyle M{\hat {\ddot {u}}}={\hat {f}}-{\hat {r}}}$ ${\displaystyle K_{NL}=G^{i}+K^{i}-F^{i}\,:}$

${\displaystyle (pM+K_{NL}^{i})\Delta {\hat {u}}={\hat {f}}^{i}-{\hat {r}}^{i}-M{\hat {\ddot {u}}}^{i}=:{\hat {b}}_{I}^{i}\,.}$		(I)

If only the position of equilibrium is sought, this is reduced to ${\displaystyle {\hat {\ddot {u}}}={\hat {0}}}$

${\displaystyle K_{NL}^{i}\Delta {\hat {u}}={\hat {f}}^{i}-{\hat {r}}^{i}=:{\hat {b}}_{II}^{i}\,.}$		(II)

The node displacements at a point in time are calculated from the currently known node displacements and accelerations using the following scheme. The scheme can also be used in the static case. There the accelerations disappear and time only has an ordering function for the successive equilibrium positions. However, this has no influence on the scheme. ${\displaystyle t^{n+1}}$ ${\displaystyle t^{n}}$

1. The solution sought is initialized with the known shifts in the last increment (or the zero vector) and the iteration counter with i = 0. The mass matrix is provided in dynamic systems . In general (not here) the boundary conditions are entered in the solution vector and partitioned into a known and an unknown part. ${\displaystyle {\hat {u}}^{n+1}}$ ${\displaystyle {\hat {u}}^{n+1,0}={\hat {u}}^{n}}$ ${\displaystyle M}$
2. The matrix and the residual are set up with the approximate solution provided . ${\displaystyle K_{NL}^{i}}$ ${\displaystyle {\hat {b}}_{I/II}^{i}}$ ${\displaystyle {\hat {u}}^{n+1,i}}$
3. The increment is calculated using equations (I) or (II) . ${\displaystyle \Delta {\hat {u}}}$
4. Fallen suitable norms of the vectors and below a predetermined limit , is the approximate solution accepted and in the solution vector transmitted, and - if desired - the accelerations calculated. The elements are given the opportunity to update their inner variables (see #Tangent operator C ). The counter is incremented and continued with step 1 or the analysis is ended. ${\displaystyle \Delta {\hat {u}}}$ ${\displaystyle {\hat {b}}_{I/II}^{i}}$ ${\displaystyle {\hat {u}}^{n+1,i}}$ ${\displaystyle {\hat {u}}^{n+1}}$ ${\displaystyle {\hat {\ddot {u}}}^{n+1}}$ ${\displaystyle n}$
5. However, if the norms of the vectors are unacceptable or not, the approximate solution is updated using , the counter is incremented and step 2 is continued. ${\displaystyle \Delta {\hat {u}}}$ ${\displaystyle {\hat {b}}_{I/II}^{i}}$ ${\displaystyle {\hat {u}}^{n+1,i+1}={\hat {u}}^{n+1,i}+\Delta {\hat {u}}}$ ${\displaystyle i}$

Explicit solution to nonlinear problems

Visualization of an FEM simulation of the deformation of a car in an asymmetrical frontal collision

In the case of explicit time integration, the equation of motion (*) is already the determining equation for the only unknown than the vectors and is calculated from quantities known at the time : ${\displaystyle M{\hat {\ddot {u}}}={\hat {f}}-{\hat {r}}}$ ${\displaystyle {\hat {\ddot {u}}}^{n+1}\,,}$ ${\displaystyle {\hat {r}}}$ ${\displaystyle {\hat {f}}}$ ${\displaystyle t^{n}}$

${\displaystyle M{\hat {\ddot {u}}}^{n+1}={\hat {f}}^{n}-{\hat {r}}^{n}\,.}$

The speeds and displacements for calculating the reaction forces for the next increment are determined from the accelerations , see explicit time integration compared to the Newmark beta method . A further simplification is achieved by diagonalizing the mass matrix ("lumped mass matrix") so that ${\displaystyle {\hat {r}}}$ ${\displaystyle M}$

${\displaystyle {\hat {\ddot {u}}}^{n+1}=M^{-1}({\hat {f}}^{n}-{\hat {r}}^{n})}$

can be evaluated particularly quickly. This is also necessary because this method is only stable below a critical time step size, which is measured according to the Courant-Friedrichs-Lewy condition according to how long a signal needs to get from one node to the next. The minimum node distance is calculated for steel components with a modulus of elasticity ( megapascal ) and a density : ${\displaystyle \Delta t_{c}}$ ${\displaystyle l_{c}=10\,\mathrm {mm} \,,}$ ${\displaystyle E=200.000\,\mathrm {MPa} }$ ${\displaystyle \rho =7870\,{\frac {\mathrm {kg} }{\mathrm {m} ^{3}}}}$

${\displaystyle \Delta t_{c}\leq {\frac {l_{c}}{c}}={\frac {l_{c}}{\sqrt {\frac {E}{\rho }}}}\approx 2\cdot 10^{-6}\mathrm {s} }$

where is the wave propagation speed in steel. With these values, which are common in practice, the critical time step size is in the range of microseconds . For movements lasting tenths of a second, tens of thousands of time steps must therefore often be calculated. It is advantageous that non-linearities can be taken into account without linearization, which is why this method is used for non-linear, dynamic and short-term processes such as crash test simulations , see figure. Another advantage is that the effort for calculating the accelerations only increases linearly with the dimension of the solution vector , so that this method is also suitable for very large problems under quasi-static conditions. ${\displaystyle c}$ ${\displaystyle {\hat {u}}}$

Element matrices at the integration points

The integrals, which occur in the principle of d'Alembert in the Lagrangian version, cannot be exactly integrated in the general application. Instead, the volume integrals are calculated using numerical integration methods such as Gaussian quadrature , in which the integral is approximated as the sum of weighted integrands at integration points .

In this section the matrices appearing above, which are to be built up at each integration point, are given.

Shape functions N and their derivatives

An eight-nodular hexahedron element as part of a gearwheel

Each element models a three-dimensional volume of the body spanned by the nodes. The picture shows an eight-nodular hexahedron element as part of the body of a gear. The coordinates of the points in the element are determined with shape functions depending on local coordinates ${\displaystyle {\vec {X}}}$ ${\displaystyle N^{i}}$

${\displaystyle {\vec {\Theta }}={\begin{pmatrix}\xi \\\eta \\\zeta \end{pmatrix}}\in [-1,1]^{3}}$

interpolated between the k nodes of the element:

${\displaystyle {\vec {X}}={\begin{pmatrix}X\\Y\\Z\end{pmatrix}}=\sum _{i=1}^{k}N^{i}({\vec {\Theta }}){\vec {X}}^{i}=N({\vec {\Theta }}){\hat {X}}\,,}$

The 3 k × 1 vector contains all components of the nodal coordinates ${\displaystyle {\hat {X}}}$ ${\displaystyle {\vec {X}}^{1\ldots k}}$

${\displaystyle {\hat {X}}={\begin{pmatrix}X^{1}&Y^{1}&Z^{1}&X^{2}&Y^{2}&Z^{2}&\ldots &X^{k}&Y^{k}&Z^{k}\end{pmatrix}}^{\top }}$

and the 3 × 3 k matrix

${\displaystyle N={\begin{pmatrix}N^{1}&0&0&N^{2}&0&0&\ldots &N^{k}&0&0\\0&N^{1}&0&0&N^{2}&0&\ldots &0&N^{k}&0\\0&0&N^{1}&0&0&N^{2}&\ldots &0&0&N^{k}\end{pmatrix}}}$

The shape functions The argument of the shape functions has been omitted here for the sake of clarity and this will also be done in the following. ${\displaystyle N^{1\ldots k}\,.}$ ${\displaystyle {\vec {\Theta }}}$

The derivative of the shape functions according to the global coordinates is calculated with the Jacobi matrix : ${\displaystyle X,\,Y,\,Z}$ ${\displaystyle J}$

${\displaystyle {\begin{pmatrix}N_{,\xi }^{i}\\N_{,\eta }^{i}\\N_{,\zeta }^{i}\end{pmatrix}}=\underbrace {\begin{pmatrix}X_{,\xi }&Y_{,\xi }&Z_{,\xi }\\X_{,\eta }&Y_{,\eta }&Z_{,\eta }\\X_{,\zeta }&Y_{,\zeta }&Z_{,\zeta }\end{pmatrix}} _{J^{\top }}{\begin{pmatrix}N_{,X}^{i}\\N_{,Y}^{i}\\N_{,Z}^{i}\end{pmatrix}}\quad \Leftrightarrow \quad {\begin{pmatrix}N_{,X}^{i}\\N_{,Y}^{i}\\N_{,Z}^{i}\end{pmatrix}}=\underbrace {\begin{pmatrix}\xi _{,X}&\eta _{,X}&\zeta _{,X}\\\xi _{,Y}&\eta _{,Y}&\zeta _{,Y}\\\xi _{,Z}&\eta _{,Z}&\zeta _{,Z}\end{pmatrix}} _{J^{\top -1}}{\begin{pmatrix}N_{,\xi }^{i}\\N_{,\eta }^{i}\\N_{,\zeta }^{i}\end{pmatrix}}\,.}$

An index after a comma means here, as in the following, a derivation according to the variable The matrix is the transposed inverse Jacobian matrix . Because the inversion of the Jacobian matrix always succeeds in the coordinate transformation and the derivation according to the local coordinates is analytically feasible, the right equation results in the sought-after derivatives according to the global coordinates. The vectorial surface elements and volume shapes required for the integration are converted with the determinant of the Jacobian matrix : ${\displaystyle (\cdot )_{,x}}$ ${\displaystyle x\,.}$ ${\displaystyle J^{\top -1}}$ ${\displaystyle \xi ,\,\eta ,\,\zeta }$

${\displaystyle {\begin{aligned}\mathrm {d} {\vec {A}}=&{\vec {\varphi }}_{,X}\times {\vec {\varphi }}_{,Y}\,\mathrm {d} X\,\mathrm {d} Y=\mathrm {det} (J)J^{\top -1}{\vec {\varphi }}_{,\xi }\times {\vec {\varphi }}_{,\eta }\,\mathrm {d} \xi \,\mathrm {d} \eta \\\mathrm {d} V=&\mathrm {d} X\,\mathrm {d} Y\,\mathrm {d} Z=\mathrm {det} (J)\,\mathrm {d} \xi \,\mathrm {d} \eta \,\mathrm {d} \zeta \,.\end{aligned}}}$

An area that can be described by the coordinates X and Y was assumed here as an example . ${\displaystyle {\vec {\varphi }}(X,Y)}$

Displacements and their gradient H

The displacement vector is interpolated in isoparametric elements analogous to the position vector:

${\displaystyle {\vec {u}}={\begin{pmatrix}u\\v\\w\end{pmatrix}}=\sum _{i=1}^{k}N^{i}{\vec {u}}^{i}=N{\hat {u}}}$

The 3 k × 1 vector contains the displacement components and in the x, y and z directions at the nodes ${\displaystyle {\hat {u}}}$ ${\displaystyle u,\,v}$ ${\displaystyle w}$

${\displaystyle {\hat {u}}={\begin{pmatrix}u^{1}&v^{1}&w^{1}&u^{2}&v^{2}&w^{2}&\ldots &u^{k}&v^{k}&w^{k}\end{pmatrix}}^{\top }}$

The displacement gradient is also created

${\displaystyle \mathbf {H} =\mathrm {GRAD} ({\vec {u}})={\frac {\mathrm {d} {\vec {u}}}{\mathrm {d} {\vec {X}}}}={\begin{pmatrix}{\vec {u}}_{,X}&{\vec {u}}_{,Y}&{\vec {u}}_{,Z}\end{pmatrix}}={\begin{pmatrix}u_{,X}&u_{,Y}&u_{,Z}\\v_{,X}&v_{,Y}&v_{,Z}\\w_{,X}&w_{,Y}&w_{,Z}\end{pmatrix}}\,.}$

The derivatives of the displacement components are calculated using the derived shape functions, e.g. B.

${\displaystyle {\vec {u}}_{,X}={\begin{pmatrix}u_{,X}\\v_{,X}\\w_{,X}\end{pmatrix}}=\sum _{i=1}^{k}N_{,X}^{i}{\vec {u}}^{i}=N_{,X}{\hat {u}}\,.}$

In the Galerkin method, the virtual displacements that occur in d'Alembert's principle are treated in the same way as the nodal displacements:

${\displaystyle \delta {\vec {u}}={\begin{pmatrix}\delta u\\\delta v\\\delta w\end{pmatrix}}=\sum _{i=1}^{k}N^{i}\delta {\vec {u}}^{i}=N\delta {\hat {u}}\,.}$

Distortion Displacement Matrix B

The symmetrical Green-Lagrange strain tensor results from the displacement gradient according to

${\displaystyle \mathbf {E} ={\frac {1}{2}}(\mathbf {H} +\mathbf {H} ^{\top }+\mathbf {H} ^{\top }\cdot \mathbf {H} )\,.}$

Its six independent components are entered in a vector ( Voigt notation ):

${\displaystyle {\hat {E}}:={\begin{pmatrix}E_{xx}&E_{yy}&E_{zz}&2E_{xy}&2E_{yz}&2E_{xz}\end{pmatrix}}^{\top }\,.}$

The distortion displacement matrix is the derivative of the vector according to the nodal displacements: ${\displaystyle B}$ ${\displaystyle {\hat {E}}}$

${\displaystyle B={\frac {\mathrm {d} {\hat {E}}}{\mathrm {d} {\hat {u}}}}\in \mathbb {R} ^{6\times 3k}\quad \leftrightarrow \quad B_{ij}={\frac {{\textrm {d}}E_{i}}{{\textrm {d}}{\hat {u}}_{j}}}\quad {\text{für}}\quad i=1,2,...,6\quad {\text{und}}\quad j=1,2,...,3k\,.}$

The differential virtual distortions then result from the virtual node displacements with the B matrix  : ${\displaystyle \delta {\hat {E}}}$ ${\displaystyle \delta {\hat {u}}}$

${\displaystyle \delta {\hat {E}}=B\delta {\hat {u}}\,.}$

Geometrically linear case

In the geometrically linear case, the distortions are

${\displaystyle \mathbf {E} _{L}:={\frac {1}{2}}(\mathbf {H} +\mathbf {H} ^{\top })\quad \rightarrow \quad {\hat {E}}_{L}={\begin{pmatrix}u_{,X}\\v_{,Y}\\w_{,Z}\\u_{,Y}+v_{,X}\\v_{,Z}+w_{,Y}\\u_{,Z}+w_{,X}\end{pmatrix}}}$

linear in the node displacements, which is why the B-matrix results from "factoring out" the node displacements and the clear form

${\displaystyle B=B_{L}:={\begin{pmatrix}N_{,X}^{1}&0&0&N_{,X}^{2}&0&0&\ldots &N_{,X}^{k}&0&0\\0&N_{,Y}^{1}&0&0&N_{,Y}^{2}&0&\ldots &0&N_{,Y}^{k}&0\\0&0&N_{,Z}^{1}&0&0&N_{,Z}^{2}&\ldots &0&0&N_{,Z}^{k}\\N_{,Y}^{1}&N_{,X}^{1}&0&N_{,Y}^{2}&N_{,X}^{2}&0&\ldots &N_{,Y}^{k}&N_{,X}^{k}&0\\0&N_{,Z}^{1}&N_{,Y}^{1}&0&N_{,Z}^{2}&N_{,Y}^{2}&\ldots &0&N_{,Z}^{k}&N_{,Y}^{k}\\N_{,Z}^{1}&0&N_{,X}^{1}&N_{,Z}^{2}&0&N_{,X}^{2}&\ldots &N_{,Z}^{k}&0&N_{,X}^{k}\end{pmatrix}}}$

owns. Then: ${\displaystyle {\hat {E}}_{L}=B_{L}{\hat {u}}\,.}$

Geometrically nonlinear case

In geometrically nonlinear case has to geometrically linear B-matrix nor an interest in ${\displaystyle B_{L}}$

${\displaystyle \mathbf {E} _{NL}={\frac {1}{2}}\mathbf {H} ^{\top }\cdot \mathbf {H} }$

that are added to the matrix

${\displaystyle B_{NL}={\begin{pmatrix}B_{NL}^{1}&B_{NL}^{2}&\dots &B_{NL}^{k}\end{pmatrix}}}$

with the blocks

${\displaystyle B_{NL}^{i}={\frac {\mathrm {d} {\hat {E}}_{NL}}{\mathrm {d} {\vec {u}}^{i}}}={\begin{pmatrix}u_{,X}N_{,X}^{i}&v_{,X}N_{,X}^{i}&w_{,X}N_{,X}^{i}\\u_{,Y}N_{,Y}^{i}&v_{,Y}N_{,Y}^{i}&w_{,Y}N_{,Y}^{i}\\u_{,Z}N_{,Z}^{i}&v_{,Z}N_{,Z}^{i}&w_{,Z}N_{,Z}^{i}\\u_{,X}N_{,Y}^{i}+u_{,Y}N_{,X}^{i}&v_{,X}N_{,Y}^{i}+v_{,Y}N_{,X}^{i}&w_{,X}N_{,Y}^{i}+w_{,Y}N_{,X}^{i}\\u_{,Y}N_{,Z}^{i}+u_{,Z}N_{,Y}^{i}&v_{,Y}N_{,Z}^{i}+v_{,Z}N_{,Y}^{i}&w_{,Y}N_{,Z}^{i}+w_{,Z}N_{,Y}^{i}\\u_{,X}N_{,Z}^{i}+u_{,Z}N_{,X}^{i}&v_{,X}N_{,Z}^{i}+v_{,Z}N_{,X}^{i}&w_{,X}N_{,Z}^{i}+w_{,Z}N_{,X}^{i}\end{pmatrix}}}$

leads. The resulting B matrix

${\displaystyle B=B_{L}+B_{NL}}$

in this case depends on the node displacements . ${\displaystyle {\hat {u}}}$

Tangent operator C

The stresses how the strains entered into a vector

${\displaystyle {\hat {T}}={\begin{pmatrix}{\tilde {T}}_{xx}&{\tilde {T}}_{yy}&{\tilde {T}}_{zz}&{\tilde {T}}_{xy}&{\tilde {T}}_{yz}&{\tilde {T}}_{xz}\end{pmatrix}}^{\top }\,.}$

At the element level, a material routine must calculate these stresses from stresses in the last time step, a distortion increment and possibly other internal variables of the material model: ${\displaystyle {\mathcal {T}}}$ ${\displaystyle {\hat {T}}^{n}}$ ${\displaystyle t^{n}\,,}$ ${\displaystyle \Delta {\hat {E}}}$ ${\displaystyle V^{n}}$

${\displaystyle {\hat {T}}^{i}={\mathcal {T}}({\hat {T}}^{n},\Delta {\hat {E}},V^{n})\,.}$

With linear elasticity, the stresses are linear in the distortions. The 6 × 6 matrix is then the elasticity matrix independent of the distortions. For the application of Newton's method in the calculation of problems with nonlinear material behavior, the derivative of the stresses according to the distortions must be provided, which leads to the consistent tangent operator, which is also a 6 × 6 matrix. In the case of linear elasticity, the tangent operator results from the derivation of the stresses according to the distortions at the location of the current stresses (and internal variables, the current value of which is obtained when calculating the stresses): ${\displaystyle {\hat {T}}^{i}={\hat {T}}^{n}+C_{L}\Delta {\hat {E}}\,.}$ ${\displaystyle C_{L}}$ ${\displaystyle C}$ ${\displaystyle C=C_{L}\,.}$ ${\displaystyle {\hat {T}}^{i}}$

${\displaystyle C_{NL}H=\lim _{s\rightarrow 0}{\frac {{\mathcal {T}}({\hat {T}}^{i},sH,V^{i})-{\hat {T}}^{i}}{s}}\quad \forall \;H\in \mathbb {R} ^{6}\quad \rightarrow \;C=C_{NL}\,.}$

Then can

${\displaystyle C={\frac {\mathrm {d} {\hat {T}}^{i}}{\mathrm {d} {\hat {E}}}}}$

to be written. The consistency refers to the fact that the tangent operator is calculated from the derivation of the material routine and not in the analytical material model, which is implemented numerically by. ${\displaystyle {\mathcal {T}}}$ ${\displaystyle {\mathcal {T}}}$

Geometric stiffness matrix G

The geometric stiffness matrix

${\displaystyle G^{i}:=\int _{V}{\frac {\mathrm {d} }{\mathrm {d} {\hat {u}}}}\left.(B^{i\top }{\hat {T}}^{i})\right|_{{\hat {T}}^{i}={\text{const.}}}\mathrm {d} V=\int _{V}{\begin{pmatrix}D(g_{11})&\ldots &D(g_{1k})\\\vdots &\ddots &\vdots \\D(g_{k1})&\ldots &D(g_{kk})\end{pmatrix}}\,\mathrm {d} V}$

has a block structure with blocks of diagonal matrices

${\displaystyle D(g_{mn}):={\frac {1}{2}}{\begin{pmatrix}g_{mn}&0&0\\0&g_{mn}&0\\0&0&g_{mn}\end{pmatrix}}}$

and the diagonal links

${\displaystyle {\begin{aligned}g_{mn}=g_{nm}=&{\tilde {T}}_{1}^{i}N_{,X}^{m}N_{,X}^{n}+{\tilde {T}}_{2}^{i}N_{,Y}^{m}N_{,Y}^{n}+{\tilde {T}}_{3}^{i}N_{,Z}^{m}N_{,Z}^{n}+{\tilde {T}}_{4}^{i}(N_{,X}^{m}N_{,Y}^{n}+N_{,Y}^{m}N_{,X}^{n})\\&+{\tilde {T}}_{5}^{i}(N_{,Y}^{m}N_{,Z}^{n}+N_{,Z}^{m}N_{,Y}^{n})+{\tilde {T}}_{6}^{i}(N_{,X}^{m}N_{,Z}^{n}+N_{,Z}^{m}N_{,X}^{n})\end{aligned}}}$

example

Bar with variable cross-section under single force loading.

The elongation of a tension rod clamped on one side with a linearly decreasing cross-sectional area under individual force at the end, as in the picture, should be calculated. For this purpose, a one-dimensional, two-nodal rod element lying in the x-direction is constructed. The x-coordinate of the points in the member, the shape function matrix and node coordinates form the relationship

${\displaystyle {\begin{aligned}X=&N{\hat {X}}={\frac {1}{2}}{\begin{pmatrix}1-\xi &1+\xi \end{pmatrix}}{\begin{pmatrix}X^{1}\\X^{2}\end{pmatrix}}\\=&{\frac {1}{2}}(X^{1}+X^{2})+{\frac {1}{2}}(X^{2}-X^{1})\xi =:{\frac {1}{2}}(R+L\xi )\,,\end{aligned}}}$

with The Jacobi matrix degenerates to a number: ${\displaystyle \xi \in [-1,1]\,.}$

${\displaystyle J=X_{,\xi }={\frac {L}{2}}\,.}$

The derivation of the shape function matrix is ​​calculated with its inverse:

${\displaystyle N_{,X}=J^{-1}N_{,\xi }={\frac {2}{L}}{\frac {1}{2}}{\begin{pmatrix}-1&1\end{pmatrix}}={\frac {1}{L}}{\begin{pmatrix}-1&1\end{pmatrix}}\,.}$

The displacements, the shape function matrix and node displacements form the relationship

${\displaystyle u=N{\hat {U}}={\frac {1}{2}}{\begin{pmatrix}1-\xi &1+\xi \end{pmatrix}}{\begin{pmatrix}U^{1}\\U^{2}\end{pmatrix}}\,.}$

In the geometrically linear domain, the strains and the distortion displacement matrix are

${\displaystyle E_{L}=u_{,X}=N_{,X}{\hat {U}}={\frac {1}{L}}{\begin{pmatrix}-1&1\end{pmatrix}}{\begin{pmatrix}U^{1}\\U^{2}\end{pmatrix}}=B_{L}{\hat {U}}\quad \rightarrow \quad B_{L}={\frac {1}{L}}{\begin{pmatrix}-1&1\end{pmatrix}}\,.}$

Of the strains, there is only one constant component in the direction of the member and the same applies to the stresses:

${\displaystyle T=CE_{L}=CB_{L}{\hat {U}}\,,}$

in which only one component C is required from the elasticity matrix C _L. The cross-sectional area is described with two parameters a and b: A = a-bX. This becomes the volume element ${\displaystyle \mathrm {d} V=A\operatorname {det} (J)\,\mathrm {d} \xi ={\tfrac {L}{2}}(a-bX)\,\mathrm {d} \xi \,.}$

The reactions can now be calculated:

${\displaystyle {\begin{aligned}{\hat {r}}=&\int _{V}B_{L}^{\top }CB_{L}\,\mathrm {d} V{\hat {U}}=\int _{-1}^{1}{\frac {C}{L^{2}}}{\begin{pmatrix}1&-1\\-1&1\end{pmatrix}}{\frac {L}{2}}(a-bX)\,\mathrm {d} \xi {\hat {U}}\\=&\int _{-1}^{1}\left[a-{\frac {b}{2}}(R+L\xi )\right]\,\mathrm {d} \xi {\frac {C}{2L}}{\begin{pmatrix}1&-1\\-1&1\end{pmatrix}}{\hat {U}}\,.\end{aligned}}}$

The integral is calculated as

${\displaystyle \int _{-1}^{1}\left[a-{\frac {b}{2}}(R+L\xi )\right]\,\mathrm {d} \xi =\left[a\xi -{\frac {b}{2}}\left(R\xi +L{\frac {\xi ^{2}}{2}}\right)\right]_{-1}^{1}=2\left(a-b{\frac {R}{2}}\right)=2A_{m}\,,}$

if A _{m is} the cross-sectional area of ​​the rod at the center of the element. The stiffness matrix is ​​thus fixed:

${\displaystyle {\hat {r}}=K_{L}{\hat {U}}\quad \rightarrow \quad K_{L}=C{\frac {A_{m}}{L}}{\begin{pmatrix}1&-1\\-1&1\end{pmatrix}}\,.}$

The external force is calculated from a constant stress σ at the end of the rod at ${\displaystyle \xi =1:}$

${\displaystyle {\hat {f}}=\int _{A^{\sigma }}N^{\top }{\vec {t}}_{0}\,\mathrm {d} A={\begin{pmatrix}0\\A(\xi =1)\sigma \end{pmatrix}}={\begin{pmatrix}0\\F\end{pmatrix}}\,.}$

Two half-length rod elements can be assembled to form a finite element model of the rod. The second element results analogously to the above explanations with the node coordinates X ² and X ³ and the node displacements U ² and U ³ , see figure. Both elements share the degrees of freedom at node 2. Assembling both element contributions results in:

${\displaystyle {\hat {r}}=C{\begin{pmatrix}{\frac {A_{m}^{1}}{L^{1}}}&-{\frac {A_{m}^{1}}{L^{1}}}&0\\-{\frac {A_{m}^{1}}{L^{1}}}&{\frac {A_{m}^{1}}{L^{1}}}+{\frac {A_{m}^{2}}{L^{2}}}&-{\frac {A_{m}^{2}}{L^{2}}}\\0&-{\frac {A_{m}^{2}}{L^{2}}}&{\frac {A_{m}^{2}}{L^{2}}}\end{pmatrix}}{\begin{pmatrix}U^{1}\\U^{2}\\U^{3}\end{pmatrix}}\,.}$

The superscript indices indicate the element number, i.e. no exponentiation . The element parameters result from the node coordinates and the definitions for A _m and L. By setting U ¹ to zero and a force F acting on the third node, the displacements of the second and third node are calculated from the matrix equation ${\displaystyle A_{m}^{1,2},\,L^{1,2}}$

${\displaystyle C{\begin{pmatrix}{\frac {A_{m}^{1}}{L^{1}}}+{\frac {A_{m}^{2}}{L^{2}}}&-{\frac {A_{m}^{2}}{L^{2}}}\\-{\frac {A_{m}^{2}}{L^{2}}}&{\frac {A_{m}^{2}}{L^{2}}}\end{pmatrix}}{\begin{pmatrix}U^{2}\\U^{3}\end{pmatrix}}={\begin{pmatrix}0\\F\end{pmatrix}}\,.}$

With the parameters in the table below, U ² = 0.032 and U ³ = 0.109 are calculated . In a similar way, any number of rod elements can in principle be combined.

Comparison of the FEM solution (red) with the analytical (black).

There is also an analytical solution to the problem. The strain in the rod is obtained - just as in the element - of ε = u _{, X} and the stresses are proportional to: σ = ε C = C u _{, X} . The longitudinal force in the rod is constantly equal to the tensile force, but is distributed over a decreasing cross-sectional area: F = (a-b X) σ C = (a-b X) u _{, X} . This differential equation can be solved uniquely with the boundary condition u (0) = 0:

${\displaystyle u_{,X}={\frac {F}{C(a-bX)}}\quad \rightarrow \quad u(X)={\frac {F}{Cb}}\ln \left({\frac {a}{a-bX}}\right)\,.}$

In the picture on the right, this analytical solution for the displacement at the end of the rod at X = L is plotted with the FEM solution using the parameters and variable number of elements used in the table.

parameter	L.	a	b	C.	F.
unit	mm	mm 2	mm	MPa	N
value	1000	100	0.09	200,000	1000

The convergence of the FEM solution towards a limit value with increasing network refinement has a typical profile here. It results from the fact that the rod element has constant elongation and tension over its length, but the tension in the rod increases continuously due to the decreasing cross-sectional area. The approximation of this smooth, monotonously rising course by a step function in the FE model is always more precise, the smaller the steps and thus the smaller the elements.

Footnotes

1. ↑ ^{a b c d} Comparison of the equation

   ${\displaystyle {\hat {\ddot {u}}}^{n+1}=p{\hat {u}}^{n+1}+{\hat {p}}^{n}}$		(NbV)

   with the regulation for updating the variables

   ${\displaystyle {\ddot {x}}_{n+1}={\frac {2}{\beta \Delta t^{2}}}(x_{n+1}-x_{n})-{\frac {2}{\beta \Delta t}}{\dot {x}}_{n}-{\frac {1-\beta }{\beta }}{\ddot {x}}_{n}}$

   results in: and ${\displaystyle p={\tfrac {2}{\beta \Delta t^{2}}}}$ ${\displaystyle {\hat {p}}^{n}=-{\tfrac {2}{\beta \Delta t^{2}}}x_{n}-{\frac {2}{\beta \Delta t}}{\dot {x}}_{n}-{\frac {1-\beta }{\beta }}{\ddot {x}}_{n}\,.}$

2. ↑ Comparison of the equation

   ${\displaystyle {\hat {u}}^{n+1}=q{\hat {\ddot {u}}}^{n}+{\hat {q}}^{n}\,.}$		(EZI)

   with the regulation for updating the variables

   ${\displaystyle x_{n+1}=x_{n}+\Delta t({\dot {x}}_{n-1/2}+\Delta t{\ddot {x}}_{n})}$

   results in: and ${\displaystyle q=\Delta t^{2}}$ ${\displaystyle {\hat {q}}^{n}=x_{n}+\Delta t{\dot {x}}_{n-1/2}\,.}$

3. ↑ ^{a b c d} The derivation of an n ₁ × 1 vector y after an n ₂ × 1 vector x is the n ₁ × n ₂ matrix Z with the entries

   ${\displaystyle Z_{ij}={\tfrac {\partial y_{i}}{\partial x_{j}}}\,,\quad i=1,...,n_{1},\,j=1,...,n_{2}\,.}$

   Then it is also written. ${\displaystyle Z={\frac {\mathrm {d} y}{\mathrm {d} x}}}$

literature

• Klaus-Jürgen Bathe: Finite element methods: Matrices and linear algebra, the finite element method, solving equilibrium conditions and principle equations . Springer, 1986, ISBN 3-540-15602-X .
• Peter Wriggers: Nonlinear Finite Element Methods . Springer, 2001, ISBN 3-540-67747-X .
• OC Zienkiewicz, RL Taylor, DD Fox: The Finite Element Method for Solid and Structural Mechanics . Elsevier, 2013, ISBN 9781856176347 .