Nonlinear finite elements/Total Lagrangian approach

Total Lagrangian Approach

In the total Lagrangian approach, the discrete equations are formulated with respect to the reference configuration. The independent variables are $X$ and $t$. The dependent variable is the displacement $u(X,t)$.

Total Lagrangian Stress and Strain Measures

The strain is

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

For the reference configuration,

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

The Cauchy stress (force/current area) is

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

The engineering stress (force/initial area) is

${ P = \cfrac{T}{A_0} ~. }$

The two stresses are related by

$\sigma = \cfrac{A_0}{A} P ~; \qquad P = \cfrac{A}{A_0} \sigma ~.$

For the reference configuration,

$\sigma = P ~.$

Governing Equations in Total Lagrangian Form

The Total Lagrangian governing equations are:

• Conservation of Mass:

${ \rho~J = \rho_0~J_0 = \rho_0~. }$

For the axially loaded bar,

$\rho\cfrac{A}{A_0} F = \rho_0 \qquad \implies \qquad \rho~A~F = \rho_0~A_0 ~.$

• Conservation of Momentum:

${ \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 $A_0$ is constant. Therefore,

$\frac{\partial P}{\partial X} + \rho_0~b = \rho_0~\frac{\partial^2 u}{\partial t^2} ~.$

In short form,

$P_{,X} + \rho_0~b = \rho_0~\ddot{u} ~.$

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

$P_{,X} + \rho_0~b = 0 ~.$

• Conservation of Energy:

For the bar:

${ \rho_0\frac{\partial W_{\text{int}}}{\partial t} = \frac{\partial F}{\partial t}~P ~. }$

In short form:

$\rho_0\dot{W}_{\text{int}} = \dot{F}~P ~.$

• Constitutive Equations:
  • Total Form:

For a linear elastic material:

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

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

For small strains:

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

• • Rate Form:

For a linear elastic material:

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

The Wave Equation

The momentum equation is

$P_{,X} + \rho_0~b = \rho_0~\ddot{u} ~.$

The total form constitutive equation is

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

Plug constitutive equation into momentum equation to get

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

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

${ 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 $c_0$. If acceleration is zero, the equation becomes elliptic.

Initial and Boundary Conditions

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

${ \begin{align} 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{align} }$

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

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

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 $u$ or $u_{,X}$ must be prescribed. In mechanics, instead of $u_{,X}$, the traction is prescribed.

Let $\Gamma_u$ be the part of the boundary where displacements are prescribed. Let $\Gamma_t$ be the part of the boundary where tractions (force vector per unit area) are prescribed.

Then the displacement boundary conditions are

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

The traction boundary conditions are

${ n^0~P(X,t) = t^0_X(X,t) \qquad ~\text{for}~X \in \Gamma_t ~. }$

The unit normal to the boundary in the reference configuration is $n^0$.

For the axially loaded bar, the displacement boundary condition is

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

The unit normal to the boundary is

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

Therefore, the traction boundary condition is

$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

${ \lfloor P\rceil = 0 }$

where

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

Weak Form for Total Lagrangian

The momentum equation (in its full form) is

$(A_0 P)_{,X} + \rho_0 A_0 b = \rho_0 A_0 \ddot{u} ~.$

To get the weak form over an element, we multiply the equation by a weighting function and integrate over the length of the element (from $X_a$ to $X_b$).

$\int_{X_a}^{X_b} w\left[(A_0 P)_{,X} + \rho_0 A_0 b - \rho_0 A_0 \ddot{u}\right]~dX = 0 ~.$

Integrate the first term by parts to get

$\int_{X_a}^{X_b} w (A_0 P)_{,X}~dX = \left[w A_0 P\right]_{X_a}^{X_b} - \int_{X_a}^{X_b} w_{,X} (A_0 P)~dX ~.$

Plug the above into the weak form and rearrange to get

$\int_{X_a}^{X_b} w_{,X} (A_0 P)~dX + \int_{X_a}^{X_b} w\rho_0 A_0 \ddot{u}~dX = \int_{X_a}^{X_b} w \rho_0 A_0 b~dX +\left[w A_0 P\right]_{X_a}^{X_b} ~.$

If we think of the weighting function $w$ as a variation of $u$ that satisfies the kinematic admissibility conditions, we get

${ \int_{X_a}^{X_b} (\delta u)_{,X} (A_0 P)~dX + \int_{X_a}^{X_b} \delta u\rho_0 A_0 \ddot{u}~dX = \int_{X_a}^{X_b} \delta u \rho_0 A_0 b~dX +\left[\delta u A_0 P\right]_{X_a}^{X_b} ~. }$

Recall that

$u = \varphi(X,t) - X ~.$

Therefore,

$\delta u = \delta \varphi(X,t) = \delta x ~.$

Taking a derivative of this variation, we have

$\delta u_{,X} = \delta x_{,X} = \delta F ~.$

Also,

$\begin{align} \left[\delta u A_0 P\right]_{X_a}^{X_b} &= [\delta u A_0 P]_{X_b} - [\delta u A_0 P]_{X_a} \\ &=[\delta u A_0 P n_0]_{X_b} + [\delta u A_0 P n_0]_{X_a} \\ &=[\delta u A_0 t^0_X]_{X_b} + [\delta u A_0 t^0_X]_{X_a} \\ &= \left[\delta u A_0 t^0_X\right]_{\Gamma_t} ~. \end{align}$

Therefore, we can alternatively write the weak form as

${ \int_{X_a}^{X_b} \delta F A_0 P~dX + \int_{X_a}^{X_b} \delta u\rho_0 A_0 \ddot{u}~dX = \int_{X_a}^{X_b} \delta u \rho_0 A_0 b~dX +\left[\delta u A_0 t^0_X\right]_{\Gamma_t}~. }$

Comparing with the energy equation, we see that the first term above is the internal virtual work and the weak form is the principle of virtual work for the 1-D problem. The weak form may also be written as

${ \delta W = \delta W_{\text{int}} - \delta W_{\text{ext}} + \delta W_{\text{kin}} = 0 }$

where,

${ \begin{align} \delta W_{\text{int}} & = \int_{X_a}^{X_b} (\delta u)_{,X} (A_0 P)~dX \\ \delta W_{\text{ext}} & = \int_{X_a}^{X_b} \delta u \rho_0 A_0 b~dX + \left[\delta u A_0 P\right]_{X_a}^{X_b} \\ \delta W_{\text{kin}} & = \int_{X_a}^{X_b} \delta u\rho_0 A_0 \ddot{u}~dX \end{align} }$

Finite Element Discretization for Total Lagrangian

The trial solution is

${ u_h(X,t) = \sum_{j=1}^n N_j(X) u_j(t) }$

where $n$ is the number of nodes. In matrix form,

${ u_h(X,t) = \mathbf{N}~\mathbf{u} } ~~\text{where}~~ \mathbf{N} = \begin{bmatrix} N_1(X) & N_2(X) && N_n(X) \end{bmatrix} ~\text{and}~ \mathbf{u} = \begin{bmatrix} u_1(t) \\ u_2(t) \\ \vdots \\ u_n(t) \end{bmatrix} ~.$

The test (weighting) function is

${ \delta u_h(X) = \sum_{i=1}^n N_i(X) \delta u_i ~. } ~~\text{or}~~ { \delta u_h(X) = \mathbf{N}~\delta \mathbf{u} } ~~\text{where}~~ \delta \mathbf{u} = \begin{bmatrix} \delta u_1 \\ \delta u_2 \\ \vdots \\ \delta u_n \end{bmatrix} ~.$

The derivatives of the test functions with respect to $X$ are

${ (\delta u_h)_{,X} = \sum_{i=1}^n N_{i,X} \delta u_i ~. } ~~\text{or}~~ { (\delta u_h)_{,X} = \mathbf{B}_0~\delta \mathbf{u} } ~~\text{where}~~ \mathbf{B}_0 = \begin{bmatrix} N_{1,X} & N_{2,X} && N_{n,X} \end{bmatrix} ~.$

We will derive the finite element system of equations after substituting these into the approximate weak form

$\int_{X_a}^{X_b} (\delta u_h)_{,X} (A_0 P)~dX + \int_{X_a}^{X_b} \delta u_h\rho_0 A_0 \ddot{u}_h~dX = \int_{X_a}^{X_b} \delta u_h \rho_0 A_0 b +\left[\delta u_h A_0 P\right]_{X_a}^{X_b} ~.$

Let us proceed term by term.

First LHS Term

The first term represents the virtual internal work

$\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

$\begin{align} \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{align}$

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

${ \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

${ 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~. }$

Second LHS Term

The second term represents the virtual kinetic work

$\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

$\begin{align} \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{align}$

The inertial (kinetic) force is

${ f^{\text{kin}}_i =\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

$\begin{align} f^{\text{kin}}_i & = \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{align}$

where

${ 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

$a_j := \ddot{u}_j \qquad\leftarrow \quad\text{Acceleration.}$

Note that the mass matrix is independent of time!

In matrix form,

${ \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

${ \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

$\begin{align} \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{align}$

In matrix form,

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

RHS Terms

The right hand side terms represent the virtual external work

$\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^0_X\right]_{\Gamma_t}~.$

Plugging in the test function into the above expression gives

$\begin{align} \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^0_X\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^0_X\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^0_X\right]_{\Gamma_t}\right] \\ & = \sum_{i=1}^n \delta u_i~f^{\text{ext}}_i ~. \end{align}$

In matrix notation,

${ \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

${ f_i^{\text{ext}} =\int_{X_a}^{X_b} N_i \rho_0 A_0 b~dX + \left[N_i A_0 t^0_X\right]_{\Gamma_t}~. }$

In matrix notation,

${ \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^0_X\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^0_X\right]_{\Gamma_t}~. }$

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

$\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 $1$. Then, the finite element system of equations becomes

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

or,

${ \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 $1$ is known, the acceleration at node $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 $C^1$ function.

We can take the known acceleration $a_1$ to the RHS to get

$\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

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

In matrix form,

${ \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,

$\begin{align} M^{\text{diagonal}}_{ii} & = \sum_{j=1}^n M^{\text{consistent}}_{ij} \\ & = \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{align}$

Discrete Equations for Total Lagrangian

The finite element equations in total Lagrangian form are:

${ \begin{align} u(X,t) & = \sum_i N_i(X)~u_i^e(t) = \mathbf{N}~\mathbf{u}^e \\ \varepsilon(X,t) & = \sum_i N_{i,X}~u_i^e(t) = \mathbf{B}_0~\mathbf{u}^e \\ F(X,t) & = \sum_i N_{i,X}~x_i^e(t) = \mathbf{B}_0~\mathbf{x}^e \\ \mathbf{f}^e_{\text{int}} & = \int_{\Omega_0^e} \mathbf{B}_0^T P~d\Omega_0 \\ \mathbf{f}^e_{\text{ext}} & = \int_{\Omega_0^e} \rho_0 \mathbf{N}^Tb~d\Omega_0 + \left[\mathbf{N}^T A_0 t_X^0 \right]_{\Gamma_t^e} \\ \mathbf{M}^e & = \int_{\Omega_0^e} \rho_0 \mathbf{N}^T \mathbf{N} ~d\Omega_0 \\ \mathbf{M}~\ddot{\mathbf{u}} & = \mathbf{f}_{\text{ext}} - \mathbf{f}_{\text{int}} ~. \end{align} }$