Nonlinear finite elements/Total Lagrangian approach

The strain is

${\displaystyle {\varepsilon (X,t)=F(X,t)-1={\frac {\partial x}{\partial X}}-1={\frac {\partial u}{\partial X}}=u_{,X}~.}}$

For the reference configuration,

${\displaystyle \varepsilon _{0}=\varepsilon (X,0)=F(X,0)-1={\frac {\partial X}{\partial X}}-1=0~.}$

The Cauchy stress (force/current area) is

${\displaystyle {\sigma ={\cfrac {T}{A}}~.}}$

The engineering stress (force/initial area) is

${\displaystyle {P={\cfrac {T}{A_{0}}}~.}}$

The two stresses are related by

${\displaystyle \sigma ={\cfrac {A_{0}}{A}}P~;\qquad P={\cfrac {A}{A_{0}}}\sigma ~.}$

For the reference configuration,

${\displaystyle \sigma =P~.}$

The Total Lagrangian governing equations are:

• Conservation of Mass:

${\displaystyle {\rho ~J=\rho _{0}~J_{0}=\rho _{0}~.}}$

For the axially loaded bar,

${\displaystyle \rho {\cfrac {A}{A_{0}}}F=\rho _{0}\qquad \implies \qquad \rho ~A~F=\rho _{0}~A_{0}~.}$

• Conservation of Momentum:

${\displaystyle {{\frac {\partial }{\partial X}}(A_{0}~P)+\rho _{0}~A_{0}~b=\rho _{0}~A_{0}~{\frac {\partial ^{2}u}{\partial t^{2}}}~.}}$

• For the bar ${\displaystyle A_{0}}$ is constant. Therefore,

${\displaystyle {\frac {\partial P}{\partial X}}+\rho _{0}~b=\rho _{0}~{\frac {\partial ^{2}u}{\partial t^{2}}}~.}$

In short form,

${\displaystyle P_{,X}+\rho _{0}~b=\rho _{0}~{\ddot {u}}~.}$

If ${\displaystyle {\ddot {u}}\approx 0}$, we get the equilibrium equation

${\displaystyle P_{,X}+\rho _{0}~b=0~.}$

• Conservation of Energy:

For the bar:

${\displaystyle {\rho _{0}{\frac {\partial W_{\text{int}}}{\partial t}}={\frac {\partial F}{\partial t}}~P~.}}$

In short form:

${\displaystyle \rho _{0}{\dot {W}}_{\text{int}}={\dot {F}}~P~.}$

• Constitutive Equations:
  • Total Form:

For a linear elastic material:

${\displaystyle {P(X,t)=E^{PF}~\varepsilon (X,t)=E^{PF}~u_{,X}=E^{PF}\left[F(X,t)-1\right]~.}}$

The superscript ${\displaystyle PF}$ refers to the fact that this function relates ${\displaystyle P}$ and ${\displaystyle F}$.

For small strains:

${\displaystyle E^{PF}=~~{\text{Youngs modulus}}~.}$

• • Rate Form:

For a linear elastic material:

${\displaystyle {{\dot {P}}(X,t)=E^{PF}~{\dot {\varepsilon }}(X,t)=E^{PF}{\dot {F}}(X,t)~.}}$

The momentum equation is

${\displaystyle P_{,X}+\rho _{0}~b=\rho _{0}~{\ddot {u}}~.}$

The total form constitutive equation is

${\displaystyle P(X,t)=E^{PF}~\varepsilon (X,t)=E^{PF}~u_{,X}~.}$

Plug constitutive equation into momentum equation to get

${\displaystyle \left(E^{PF}~u_{,X}\right)_{,X}+\rho _{0}~b=\rho _{0}~{\ddot {u}}~.}$

If ${\displaystyle E^{PF}}$ is constant in the bar and the body force is zero,

${\displaystyle {E^{PF}~u_{,XX}=\rho _{0}~{\ddot {u}}\qquad \implies \qquad u_{,XX}={\cfrac {\rho _{0}}{E^{PF}}}~{\ddot {u}}={\cfrac {1}{c_{0}^{2}}}~{\ddot {u}}~.}}$

This is the wave equation (hyperbolic PDE). The wave speed is ${\displaystyle c_{0}}$. If acceleration is zero, the equation becomes elliptic.

The governing equation for the rod is second-order in time. So two initial conditions are needed.

${\displaystyle {\begin{aligned}u(X,0)&=u_{0}(X)&~{\text{for}}~&X\in [0,L]\\{\dot {u}}(X,0)&=v_{0}(X)&~{\text{for}}~&X\in [0,L]\end{aligned}}}$

The rod is initially at rest. Therefore, the initial conditions are

${\displaystyle {\begin{aligned}u(X,0)&=0&~{\text{for}}~&X\in [0,L]\\{\dot {u}}(X,0)&=0&~{\text{for}}~&X\in [0,L]\end{aligned}}}$

Since the problem is one-dimensional, there are two boundary points. Also, the momentum equation is second-order in the displacements. Therefore, at each end, either ${\displaystyle u}$ or ${\displaystyle u_{,X}}$ must be prescribed. In mechanics, instead of ${\displaystyle u_{,X}}$, the traction is prescribed.

Let ${\displaystyle \Gamma _{u}}$ be the part of the boundary where displacements are prescribed. Let ${\displaystyle \Gamma _{t}}$ be the part of the boundary where tractions (force vector per unit area) are prescribed.

Then the displacement boundary conditions are

${\displaystyle {u(X,t)={\bar {u}}(X,t)\qquad ~{\text{for}}~X\in \Gamma _{u}~.}}$

The traction boundary conditions are

${\displaystyle {n^{0}~P(X,t)=t_{X}^{0}(X,t)\qquad ~{\text{for}}~X\in \Gamma _{t}~.}}$

The unit normal to the boundary in the reference configuration is ${\displaystyle n^{0}}$.

For the axially loaded bar, the displacement boundary condition is

${\displaystyle u(0,t)=0\qquad ~{\text{for}}~X=0}$

The unit normal to the boundary is

${\displaystyle {\begin{aligned}n^{0}&=-1&~{\text{at}}~&X=0\\n^{0}&=1&~{\text{at}}~&X=L\end{aligned}}}$

Therefore, the traction boundary condition is

${\displaystyle P(L,t)=T(t)\qquad ~{\text{for}}~X=L~.}$

For shock problems or fracture problems, an addition interior continuity or jump condition may be needed. This condition is written as

${\displaystyle {\lfloor P\rceil =0}}$

where

${\displaystyle \lfloor P(X)\rceil =P(X+\epsilon )-P(X-\epsilon )~,\qquad \epsilon \rightarrow 0~.}$

The first term represents the virtual internal work

${\displaystyle \delta W_{\text{int}}=\int _{X_{a}}^{X_{b}}(\delta u_{h})_{,X}(A_{0}P)~dX~.}$

Plugging in the derivative of the test function, we get

${\displaystyle {\begin{aligned}\delta W_{\text{int}}&=\int _{X_{a}}^{X_{b}}\left[\sum _{i=1}^{n}N_{i,X}\delta u_{i}\right](A_{0}P)~dX\\&=\sum _{i=1}^{n}\delta u_{i}\left[\int _{X_{a}}^{X_{b}}N_{i,X}(A_{0}P)~dX\right]\\&=\sum _{i=1}^{n}\delta u_{i}f_{i}^{\text{int}}~.\end{aligned}}}$

The last substitution is made because the virtual internal work is the work done by internal forces in moving through a virtual displacement. The matrix form of the expression for virtual internal work is

${\displaystyle {\delta W_{\text{int}}=\delta \mathbf {u} ^{T}\mathbf {f} _{\text{int}}}~~{\text{where}}~~\mathbf {f} _{\text{int}}={\begin{bmatrix}f_{1}^{\text{int}}\\f_{2}^{\text{int}}\\\vdots \\f_{n}^{\text{int}}\end{bmatrix}}~.}$

The internal force is

${\displaystyle {f_{i}^{\text{int}}=\int _{X_{a}}^{X_{b}}N_{i,X}(A_{0}P)~dX}~~{\text{or}}~~{\mathbf {f} _{\text{int}}=\int _{X_{a}}^{X_{b}}\mathbf {B} _{0}^{T}~(A_{0}P)~dX=\int _{\Omega _{0}}\mathbf {B} _{0}^{T}~P~d\Omega _{0}~.}}$

The second term represents the virtual kinetic work

${\displaystyle \delta W_{\text{kin}}=\int _{X_{a}}^{X_{b}}\delta u_{h}\rho _{0}A_{0}{\ddot {u}}_{h}~dX~.}$

Plugging in the test function, we get

${\displaystyle {\begin{aligned}\delta W_{\text{kin}}&=\int _{X_{a}}^{X_{b}}\left[\sum _{i=1}^{n}\delta u_{i}N_{i}\right]\rho _{0}A_{0}{\ddot {u}}_{h}~dX\\&=\sum _{i=1}^{n}\delta u_{i}\left[\int _{X_{a}}^{X_{b}}\rho _{0}A_{0}N_{i}{\ddot {u}}_{h}~dX\right]\\&=\sum _{i=1}^{n}\delta u_{i}f_{i}^{\text{kin}}~.\end{aligned}}}$

The inertial (kinetic) force is

${\displaystyle {f_{i}^{\text{kin}}=\int _{X_{a}}^{X_{b}}\rho _{0}A_{0}N_{i}{\ddot {u}}_{h}~dX~.}}$

Now, plugging in the trial function into the expression for the inertial force, we get

${\displaystyle {\begin{aligned}f_{i}^{\text{kin}}&=\int _{X_{a}}^{X_{b}}\rho _{0}A_{0}N_{i}\left[\sum _{j=1}^{n}{\ddot {u}}_{j}N_{j}\right]~dX\\&=\sum _{j=1}^{n}\left[\int _{X_{a}}^{X_{b}}\rho _{0}A_{0}N_{i}N_{j}~dX\right]{\ddot {u}}_{j}\\&=\sum _{j=1}^{n}M_{ij}a_{j}\end{aligned}}}$

where

${\displaystyle {M_{ij}:=\int _{X_{a}}^{X_{b}}\rho _{0}A_{0}N_{i}N_{j}~dX\qquad \leftarrow \quad {\text{Consistent Mass Matrix Term.}}}}$

and

${\displaystyle a_{j}:={\ddot {u}}_{j}\qquad \leftarrow \quad {\text{Acceleration.}}}$

Note that the mass matrix is independent of time!

In matrix form,

${\displaystyle {\mathbf {f} _{\text{kin}}=\mathbf {M} \mathbf {a} }~~{\text{where}}~~\mathbf {a} ={\begin{bmatrix}a_{1}\\a_{2}\\\vdots \\a_{n}\end{bmatrix}}={\begin{bmatrix}{\ddot {u}}_{1}\\{\ddot {u}}_{2}\\\vdots \\{\ddot {u}}_{n}\end{bmatrix}}~~~{\text{and}}~~\mathbf {f} _{\text{kin}}={\begin{bmatrix}f_{1}^{\text{kin}}\\f_{2}^{\text{kin}}\\\vdots \\f_{n}^{\text{kin}}\end{bmatrix}}~.}$

The consistent mass matrix in matrix notation is

${\displaystyle {\mathbf {M} =\int _{X_{a}}^{X_{b}}\rho _{0}A_{0}\mathbf {N} ^{T}\mathbf {N} ~dX=\int _{\Omega _{0}}\rho _{0}~\mathbf {N} ^{T}~\mathbf {N} ~d\Omega _{0}~.}}$

Plugging the expression for the inertial force into the expression for virtual kinetic work, we get

${\displaystyle {\begin{aligned}\delta W_{\text{kin}}&=\sum _{i=1}^{n}\delta u_{i}\left[\sum _{j=1}^{n}M_{ij}a_{j}\right]\\&=\sum _{i=1}^{n}\sum _{j=1}^{n}\delta u_{i}M_{ij}a_{j}\end{aligned}}}$

In matrix form,

${\displaystyle {\delta W_{\text{kin}}=(\delta \mathbf {u} )^{T}\mathbf {M} \mathbf {a} =(\delta \mathbf {u} )^{T}~\mathbf {f} _{\text{kin}}~.}}$

The right hand side terms represent the virtual external work

${\displaystyle \delta W_{\text{ext}}=\int _{X_{a}}^{X_{b}}\delta u_{h}\rho _{0}A_{0}b~dX+\left[\delta u_{h}A_{0}t_{X}^{0}\right]_{\Gamma _{t}}~.}$

Plugging in the test function into the above expression gives

${\displaystyle {\begin{aligned}\delta W_{\text{ext}}&=\int _{X_{a}}^{X_{b}}\left[\sum _{i=1}^{n}N_{i}\delta u_{i}\right]\rho _{0}A_{0}b~dX+\left[\left(\sum _{i=1}^{n}N_{i}\delta u_{i}\right)A_{0}t_{X}^{0}\right]_{\Gamma _{t}}\\&=\sum _{i=1}^{n}\delta u_{i}\left[\int _{X_{a}}^{X_{b}}N_{i}\rho _{0}A_{0}b~dX\right]+\left[\sum _{i=1}^{n}\delta u_{i}N_{i}A_{0}t_{X}^{0}\right]_{\Gamma _{t}}\\&=\sum _{i=1}^{n}\delta u_{i}\left[\int _{X_{a}}^{X_{b}}N_{i}\rho _{0}A_{0}b~dX+\left[N_{i}A_{0}t_{X}^{0}\right]_{\Gamma _{t}}\right]\\&=\sum _{i=1}^{n}\delta u_{i}~f_{i}^{\text{ext}}~.\end{aligned}}}$

In matrix notation,

${\displaystyle {\delta W_{\text{ext}}=\delta \mathbf {u} ^{T}~\mathbf {f} _{\text{ext}}}~~{\text{where}}~~\mathbf {f} _{\text{ext}}={\begin{bmatrix}f_{1}^{\text{ext}}\\f_{2}^{\text{ext}}\\\vdots \\f_{n}^{\text{ext}}\end{bmatrix}}~.}$

The external force is given by

${\displaystyle {f_{i}^{\text{ext}}=\int _{X_{a}}^{X_{b}}N_{i}\rho _{0}A_{0}b~dX+\left[N_{i}A_{0}t_{X}^{0}\right]_{\Gamma _{t}}~.}}$

In matrix notation,

${\displaystyle {\mathbf {f} _{\text{ext}}=\int _{X_{a}}^{X_{b}}\mathbf {N} ^{T}\rho _{0}A_{0}b~dX+\left[\mathbf {N} ^{T}A_{0}t_{X}^{0}\right]_{\Gamma _{t}}}~~{\text{or}}~~{\mathbf {f} _{\text{ext}}=\int _{\Omega _{0}}\rho _{0}~\mathbf {N} ^{T}~b~d\Omega _{0}+\left[\mathbf {N} ^{T}A_{0}t_{X}^{0}\right]_{\Gamma _{t}}~.}}$

Collecting the terms, we get the finite element system of equations

${\displaystyle \sum _{i=1}^{n}\delta u_{i}\left[f_{i}^{\text{int}}+f_{i}^{\text{kin}}-f_{i}^{\text{ext}}\right]=0}$

Now, the weighting function is arbitrary except at nodes where displacement BCs are prescribed. At these nodes the weighting function is zero.

For the bar, let us assume that a displacement is prescribed at node ${\displaystyle 1}$. Then, the finite element system of equations becomes

${\displaystyle f_{i}^{\text{int}}+f_{i}^{\text{kin}}-f_{i}^{\text{ext}}=0\qquad {\text{for}}~~i=2\dots n}$

or,

${\displaystyle {\sum _{j=1}^{n}M_{ij}~a_{j}+f_{i}^{\text{int}}-f_{i}^{\text{ext}}=0\qquad {\text{for}}~~i=2\dots n~.}}$

Since the displacement at node ${\displaystyle 1}$ is known, the acceleration at node ${\displaystyle 1}$ is also known. Note that we have to differentiate the displacement twice to get the acceleration and hence the prescribed displacement must be a ${\displaystyle C^{1}}$ function.

We can take the known acceleration ${\displaystyle a_{1}}$ to the RHS to get

${\displaystyle \sum _{j=2}^{n}M_{ij}~a_{j}+f_{i}^{\text{int}}-f_{i}^{\text{ext}}=-M_{i1}~a_{1}\qquad {\text{for}}~~i=2\dots n~.}$

The above equation shows that prescribed boundary accelerations (and hence prescribed boundary displacements) contribute to nodes which are not on the boundary.

We can avoid that issue by making the mass matrix diagonal (using ad-hoc procedures that conserve momentum). If the mass matrix is diagonal, we get

${\displaystyle M_{ii}~a_{i}+f_{i}^{\text{int}}-f_{i}^{\text{ext}}=0\qquad {\text{for}}~~i=2\dots n~.}$

In matrix form,

${\displaystyle {\mathbf {M} ~\mathbf {a} =\mathbf {f} _{\text{ext}}-\mathbf {f} _{\text{int}}}~~{\text{or}}~~{\mathbf {f} =\mathbf {M} ~\mathbf {a} .}~.}$

One way of generating a diagonal mass matrix or lumped mass matrix is the row-sum technique. The rows of the mass matrix are summed at lumped at the diagonal of the matrix. Thus,

${\displaystyle {\begin{aligned}M_{ii}^{\text{diagonal}}&=\sum _{j=1}^{n}M_{ij}^{\text{consistent}}\\&=\sum _{j=1}^{n}\int _{X_{a}}^{X_{b}}\rho _{0}N_{i}N_{j}A_{0}~dX\\&=\int _{X_{a}}^{X_{b}}\rho _{0}N_{i}\left[\sum _{j=1}^{n}N_{j}\right]A_{0}~dX\\&=\int _{X_{a}}^{X_{b}}\rho _{0}N_{i}A_{0}~dX\qquad {\text{since}}~\sum _{j=1}^{n}N_{j}=1~.\end{aligned}}}$