Nonlinear finite elements/Updated Lagrangian approach

Updated Lagrangian Approach

For the total Lagrangian approach, the discrete equations are formulated with respect to the reference configuration. For the updated Lagrangian approach, the discrete equations are formulated in the current configuration, which is assumed to be the new reference configuration.

The independent variables in the total Lagrangian approach are $X$ and $t$. In the updated Lagrangian they are $x$ and $t$ which are with respect to the new reference configuration.

The dependent variable in the total Lagrangian approach is the displacement $u(X,t)$. In the updated Lagrangian approach, the dependent variables are the Cauchy stress $\sigma (x,t)$ and the velocity $v(x,t)$.

Updated Lagrangian Stress and Strain Measures

The stress measure is the Cauchy stress:

$\sigma (x,t)=\frac{T}{A}.$

The strain measure is the rate of deformation (also called velocity strain):

$D(x,t)=\frac{\partial v}{\partial x}=v_{,x}.$

Note that the derivative is a spatial derivative. It can be shown that

$\underset{0}{\overset{t}{\int }}D(x,\tau )d\tau =\ln (F(x,t))$

where $F$ is the deformation gradient. So the integral of the rate of deformation in 1-D is similar to the logarithmic strain (also called natural strain).

Governing Equations in Updated Lagrangian Form

The Updated Lagrangian governing equations are:

Conservation of Mass

$\rho J={\rho }_0.$

For the rod,

$\rho AF={\rho }_0{A}_0.$

These are the same as those for the total Lagrangian approach.

Conservation of Momentum

$\frac{\partial }{\partial x}(A\sigma )+\rho Ab=\rho A_0\frac{Dv}{Dt}.$

In this case $A$ may not be constant at along the length and further simplification cannot be done. In short form,

$(A\sigma {)}_{,x}+\rho Ab=\rho A\dot{v}.$

Conservation of Energy

If we ignore heat flux and heat sources

$\rho \frac{\partial w_{\text{int}}}{\partial t}=\sigma D.$

If we include heat flux and heat sources

$\rho \frac{\partial w_{\text{int}}}{\partial t}=\sigma D-\frac{\partial q}{\partial x}+\rho s.$

In short form:

$\rho {\dot{w}}_{\text{int}}=\sigma D-q_{,x}+\rho s.$

Constitutive Equations

Total Form

For a linear elastic material:

$\sigma (x,t)=E^{\sigma D}\underset{0}{\overset{t}{\int }}D(x,\tau )d\tau =E^{\sigma D}\ln (F(x,t)).$

The superscript $\sigma D$ refers to the fact that this function relates $\sigma$ and $D$.

For small strains:

$E^{\sigma D}=\text{Youngs modulus}.$

Rate Form

For a linear elastic material:

$\dot{\sigma }(x,t)=E^{\sigma D}D(x,t).$

Initial and Boundary Conditions for Updated Lagrangian

For the updated Lagrangian approach, initial conditions are needed for the stress and the velocity. The initial displacement is assumed to be zero.

The initial conditions are:

$\begin{array}{cccc}\sigma (x,0)& ={\sigma }_0(x)& \text{for}& x\in [0,L]\\ v(x,0)& =v_0(x)& \text{for}& x\in [0,L]\\ u(x,0)& =0& \text{for}& x\in [0,L]\end{array}$

The essential boundary conditions are

$v(x,t)=\overline{v}(x,t)\text{for}x\in {\mathrm{\Gamma }}_v.$

The traction boundary conditions are

$n\sigma (x,t)=t_x(x,t)\text{for}x\in {\mathrm{\Gamma }}_t.$

The unit normal to the boundary in the current configuration is $n$. The tractions in the current configuration are related to those in the reference configuration by

$t_{x}A=t_x^0{A}_0.$

The interior continuity or jump condition is

$\lfloor \sigma A\rceil =0.$

Weak Form for Updated Lagrangian

The momentum equation is

$(A\sigma {)}_{,x}+\rho Ab=\rho A\frac{Dv}{Dt}.$

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

$\underset{x_a}{\overset{x_b}{\int }}w[(A\sigma {)}_{,x}+\rho Ab-\rho A\frac{Dv}{Dt}]dx=0.$

Integrate the first term by parts to get

$\underset{x_a}{\overset{x_b}{\int }}w(A\sigma {)}_{,x}dx={\left[wA\sigma \right]}_{x_a}^{x_b}-\underset{x_a}{\overset{x_b}{\int }}{w}_{,x}(A\sigma )dx.$

Plug the above into the weak form and rearrange to get

$\underset{x_a}{\overset{x_b}{\int }}{w}_{,x}(A\sigma )dx+\underset{x_a}{\overset{x_b}{\int }}w\rho A\frac{Dv}{Dt}dx=\underset{x_a}{\overset{x_b}{\int }}w\rho Abdx+{\left[wA\sigma \right]}_{x_a}^{x_b}.$

If we think of the weighting function $w$ as a variation of $v$ that satisfies the kinematic admissibility conditions, we get the principle of virtual power:

$\underset{x_a}{\overset{x_b}{\int }}(\delta v{)}_{,x}(A\sigma )dx+\underset{x_a}{\overset{x_b}{\int }}\delta v\rho A\frac{Dv}{Dt}dx=\underset{x_a}{\overset{x_b}{\int }}\delta v\rho Abdx+{\left[\delta vA\sigma \right]}_{x_a}^{x_b}.$

Recall that

$v_{,x}=D(x,t)\text{and}{t}_x=n\sigma .$

Therefore,

$\delta v_{,x}=\delta D(x,t)$

and

${\left[\delta vA\sigma \right]}_{x_a}^{x_b}={\left[\delta vA{t}_x\right]}_{{\mathrm{\Gamma }}_t}.$

Hence we can alternatively write the weak form as

$\underset{x_a}{\overset{x_b}{\int }}\delta DA\sigma dx+\underset{x_a}{\overset{x_b}{\int }}\delta v\rho A\frac{Dv}{Dt}dx=\underset{x_a}{\overset{x_b}{\int }}\delta v\rho Abdx+{\left[\delta vA{t}_x\right]}_{{\mathrm{\Gamma }}_t}.$

Comparing with the energy equation, we see that the first term above is the internal virtual power. Using the substitutions $Adx=d\mathrm{\Omega }$ and $Dv/DT=\dot{v}$, we can also write the weak form as

$\underset{\mathrm{\Omega }}{\int }\delta D\sigma d\mathrm{\Omega }+\underset{\mathrm{\Omega }}{\int }\delta v\rho \dot{v}d\mathrm{\Omega }=\underset{\mathrm{\Omega }}{\int }\delta v\rho bd\mathrm{\Omega }+{\left[\delta vA{t}_x\right]}_{{\mathrm{\Gamma }}_t}.$

The weak form may also be written in terms of the virtual powers as

$\delta P=\delta P_{\text{int}}-\delta P_{\text{ext}}+\delta P_{\text{kin}}=0$

where,

$\begin{array}{cc}\delta P_{\text{int}}& =\underset{x_a}{\overset{x_b}{\int }}(\delta v{)}_{,X}(A\sigma )dx\\ \delta P_{\text{ext}}& =\underset{x_a}{\overset{x_b}{\int }}\delta v\rho Abdx+{\left[\delta vA{t}_X\right]}_{{\mathrm{\Gamma }}_t}\\ \delta P_{\text{kin}}& =\underset{x_a}{\overset{x_b}{\int }}\delta v\rho A\dot{v}dx\end{array}$

Finite Element Discretization for Updated Lagrangian

The trial velocity field is

$v_h(x,t)=\sum _{j=1}^n{N}_j(x)v_j(t)$

In matrix form,

$v_h(x,t)=𝐍(x)𝐯(t)\text{where}𝐍=\left[\begin{array}{cccc}{N}_1(X)& N_2(X)& & N_n(X)\end{array}\right]\text{and}𝐯=\left[\begin{array}{c}{v}_1(t)\\ v_2(t)\\ ⋮\\ v_n(t)\end{array}\right].$

The resulting acceleration field is the material time derivative of the velocity

$a_h(x,t)={\dot{v}}_h(x,t)=\sum _{j=1}^n{N}_j(x){\dot{v}}_j(t).$

Hence,

$a_h(x,t)=\sum _{j=1}^n{N}_j(x)a_j(t)\text{or}{a}_h(X,t)=𝐍(x)𝐚(t)\text{where}𝐚=\left[\begin{array}{c}{a}_1(t)\\ a_2(t)\\ ⋮\\ a_n(t)\end{array}\right].$

Note that the shape functions are functions of $X$ and not of $x$. We could transform them into function of $x$ using the inverse mapping

$x={\phi }^{-1}(x,t).$

However, in that case the shape functions become functions of time and the procedure becomes more complicated.

The test (weighting) function is

$\delta v_h(x)=\sum _{i=1}^n{N}_i(x)\delta v_i.\text{or}\delta v_h(x)=𝐍\delta 𝐯\text{where}\delta 𝐯=\left[\begin{array}{c}\delta v_1\\ \delta v_2\\ ⋮\\ \delta v_n\end{array}\right].$

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

$(\delta v_h{)}_{,x}=\sum _{i=1}^n{N}_{i,x}\delta v_i.\text{or}(\delta v_h{)}_{,x}=𝐁\delta 𝐯\text{where}𝐁=\left[\begin{array}{cccc}{N}_{1,x}& N_{2,x}& & N_{n,x}\end{array}\right].$

It is convenient to use the isoparametric concept to compute the spatial derivatives. Let us reexamine what this approach involves.

Figure 1 shows a two-noded 1-D element in the reference and current configurations along with the parent element.

File:Config1D.png Figure 1. Reference and Current Configurations of a 1-D element.

The map between the Eulerian coordinates and the parent coordinates is

$x(\xi ,t)=(1-\xi )x_1^e(t)+\xi x_2^e(t)\xi =N_1(\xi )x_1^e(t)+N_2(\xi )x_2^e(t)\text{or}x(\xi ,t)=𝐍(\xi ){𝐱}^e(t).$

At $t=0$,

$x(\xi ,0)=x(\xi )=𝐍(\xi ){𝐗}^e.$

Therefore the displacement in the parent coordinates is

$u(\xi ,t)=x(\xi ,t)-x(\xi )=𝐍(\xi )[{𝐱}^e(t)-{𝐗}^e]\text{or}u(\xi ,t)=𝐍(\xi ){𝐮}^e(t).$

Similarly, the velocity and its variation in the parent coordinates are given by

$v(\xi ,t)=𝐍(\xi ){𝐯}^e(t),\delta v(\xi ,t)=𝐍(\xi )\delta {𝐯}^e(t).$

The acceleration in the parent coordinates is given by

$a(\xi ,t)=𝐍(\xi ){\dot{𝐯}}^e(t).$

The rate of deformation is given by

$D(x,t)=v_{,x}(x,t)={𝐍}_{,x}(X(x,t)){𝐯}^e(t).$

We can evaluate the derivative with respect to $x$ by simply using the map to the parent coordinate instead of mapping back to the reference coordinates. Using the chain rule (for the two noded element)

$N_{1,\xi }=N_{1,x}{x}_{,\xi }\text{and}{N}_{2,\xi }=N_{2,x}{x}_{,\xi }.$

In matrix form,

${𝐍}_{,\xi }={𝐍}_{,x}{x}_{,\xi }⟹{𝐍}_{,x}={𝐍}_{,\xi }{\left(x_{,\xi }\right)}^{-1}=:𝐁(\xi ).$

The rate of deformation in parent coordinates is then given by

$D(\xi ,t)=𝐁(\xi ){𝐯}^e(t).$

To derive the finite element system of equations, we follow the usual approach of substituting the trial and test functions into the approximate weak form

$\underset{x_a}{\overset{x_b}{\int }}(\delta v_h{)}_{,x}(A\sigma )dx+\underset{x_a}{\overset{x_b}{\int }}\delta v_h\rho A\frac{D{v}_h}{Dt}dx=\underset{x_a}{\overset{x_b}{\int }}\delta v_h\rho Abdx+{\left[\delta v_{h}A\sigma \right]}_{x_a}^{x_b}.$

Let us proceed term by term.

First LHS Term

The first term represents the virtual internal power

$\delta P_{\text{int}}=\underset{x_a}{\overset{x_b}{\int }}(\delta v_h{)}_{,x}(A\sigma )dx.$

Plugging in the derivative of the test function, we get

$\begin{array}{cc}\delta P_{\text{int}}& =\underset{x_a}{\overset{x_b}{\int }}\left[\sum _{i=1}^n{N}_{i,x}\delta v_i\right](A\sigma )dx\\ & =\sum _{i=1}^n\delta v_i[\underset{x_a}{\overset{x_b}{\int }}{N}_{i,x}(A\sigma )dx]\\ & =\sum _{i=1}^n\delta v_i{f}_i^{\text{int}}.\end{array}$

The matrix form of the expression for virtual internal work is

$\delta P_{\text{int}}=\delta {𝐯}^T{𝐟}_{\text{int}}\text{where}{𝐟}_{\text{int}}=\left[\begin{array}{c}{f}_1^{\text{int}}\\ f_2^{\text{int}}\\ ⋮\\ f_n^{\text{int}}\end{array}\right].$

The internal force is

$f_i^{\text{int}}=\underset{x_a}{\overset{x_b}{\int }}{N}_{i,x}(A\sigma )dx=\underset{-1}{\overset{1}{\int }}{N}_{i,\xi }(x_{,\xi }{)}^{-1}(A\sigma )x_{,\xi }d\xi =\underset{-1}{\overset{1}{\int }}{N}_{i,\xi }(A\sigma )d\xi .$

Note that the above simplification only occurs in 1-D. In matrix form,

${𝐟}_{\text{int}}=\underset{x_a}{\overset{x_b}{\int }}{𝐁}^T(A\sigma )dx=\underset{-1}{\overset{1}{\int }}({𝐍}_{,\xi }{)}^T(A\sigma )d\xi .\text{or}{𝐟}_{\text{int}}=\underset{\mathrm{\Omega }}{\int }{𝐁}^T\sigma d\mathrm{\Omega }.$

Remark:

Note that we can write

$N_{,x}=N_{,X}\frac{dX}{dx}⟹N_{,x}dx=N_{,X}dX.$

If we use the above relation, we get

${𝐟}_{\text{int}}=\underset{x_a}{\overset{x_b}{\int }}({𝐍}_{,x}{)}^T\sigma Adx=\underset{X_a}{\overset{X_b}{\int }}({𝐍}_{,X}{)}^T\sigma AdX.$

Using the relation

$\sigma A=P{A}_0$

we get

${𝐟}_{\text{int}}=\underset{X_a}{\overset{X_b}{\int }}({𝐍}_{,X}{)}^{T}P{A}_{0}dX.$

The above is the same as the expression we had for the internal force in the total Lagrangian formulation.

Second LHS Term

The second term represents the virtual kinetic power

$\delta P_{\text{kin}}=\underset{x_a}{\overset{x_b}{\int }}\delta v_h\rho A\frac{D{v}_h}{Dt}dx.$

Plugging in the test function, we get

$\begin{array}{cc}\delta P_{\text{kin}}& =\underset{x_a}{\overset{x_b}{\int }}\left[\sum _{i=1}^n\delta v_i{N}_i\right]\rho A\frac{D{v}_h}{Dt}dx\\ & =\sum _{i=1}^n\delta v_i\left[\underset{x_a}{\overset{x_b}{\int }}\rho A{N}_i\frac{D{v}_h}{Dt}dx\right]\\ & =\sum _{i=1}^n\delta v_i{f}_i^{\text{kin}}.\end{array}$

The inertial (kinetic) force is

$f_i^{\text{kin}}=\underset{x_a}{\overset{x_b}{\int }}\rho A{N}_i\frac{D{v}_h}{Dt}dx.$

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

$\begin{array}{cc}{f}_i^{\text{kin}}& =\underset{x_a}{\overset{x_b}{\int }}\rho A{N}_i\left[\sum _{j=1}^n\frac{D{v}_j}{Dt}{N}_j\right]dx\\ & =\sum _{j=1}^n\left[\underset{x_a}{\overset{x_b}{\int }}\rho A{N}_i{N}_{j}dx\right]\frac{D{v}_j}{Dt}\\ & =\sum _{j=1}^n{M}_{ij}{a}_j\end{array}$

where

$M_{ij}:=\underset{x_a}{\overset{x_b}{\int }}\rho A{N}_i{N}_{j}dx\leftarrow \text{Consistent Mass Matrix Coefficient.}$

and

$a_j:=\frac{D{v}_j}{Dt}={\dot{v}}_j\leftarrow \text{Acceleration.}$

In matrix form,

${𝐟}_{\text{kin}}=𝐌𝐚\text{where}𝐚=\left[\begin{array}{c}{a}_1\\ a_2\\ ⋮\\ a_n\end{array}\right]=\left[\begin{array}{c}{\dot{v}}_1\\ {\dot{v}}_2\\ ⋮\\ {\dot{v}}_n\end{array}\right]\text{and}{𝐟}_{\text{kin}}=\left[\begin{array}{c}{f}_1^{\text{kin}}\\ f_2^{\text{kin}}\\ ⋮\\ f_n^{\text{kin}}\end{array}\right].$

The consistent mass matrix in matrix notation is

$𝐌=\underset{x_a}{\overset{x_b}{\int }}\rho A{𝐍}^T𝐍dx=\underset{\mathrm{\Omega }}{\int }\rho {𝐍}^T𝐍d\mathrm{\Omega }.$

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

$\begin{array}{cc}\delta P_{\text{kin}}& =\sum _{i=1}^n\delta v_i\left[\sum _{j=1}^n{M}_{ij}{a}_j\right]\\ & =\sum _{i=1}^n\sum _{j=1}^n\delta v_i{M}_{ij}{a}_j\end{array}$

In matrix form,

$\delta P_{\text{kin}}=(\delta 𝐯{)}^T𝐌𝐚=(\delta 𝐯{)}^T{𝐟}_{\text{kin}}.$

Remark:

Note that since the integration domain is a function of time, the mass matrix for the updated Lagrangian formulation is also a function of time. However, from the conservation of mass, we have

${\rho }_0{A}_{0}dX=\rho Adx.$

Therefore, we can write the mass matrix as

$𝐌=\underset{x_a}{\overset{x_b}{\int }}\rho A{𝐍}^T𝐍dx=\underset{X_a}{\overset{X_b}{\int }}{\rho }_0{A}_0{𝐍}^T𝐍dX=\underset{{\mathrm{\Omega }}_0}{\int }{\rho }_0{𝐍}^T𝐍d{\mathrm{\Omega }}_0.$

Hence the mass matrices are the same for both formulations.

RHS Terms

The right hand side terms represent the virtual external power

$\delta P_{\text{ext}}=\underset{x_a}{\overset{x_b}{\int }}\delta v_h\rho Abdx+{\left[\delta v_{h}A{t}_x\right]}_{{\mathrm{\Gamma }}_t}.$

Plugging in the test function into the above expression gives

$\begin{array}{cc}\delta P_{\text{ext}}& =\underset{x_a}{\overset{x_b}{\int }}\left[\sum _{i=1}^n{N}_i\delta v_i\right]\rho Abdx+{\left[\left(\sum _{i=1}^n{N}_i\delta v_i\right)A{t}_x\right]}_{{\mathrm{\Gamma }}_t}\\ & =\sum _{i=1}^n\delta v_i\left[\underset{x_a}{\overset{x_b}{\int }}{N}_i\rho Abdx\right]+{\left[\sum _{i=1}^n\delta v_i{N}_{i}A{t}_X\right]}_{{\mathrm{\Gamma }}_t}\\ & =\sum _{i=1}^n\delta v_i[\underset{x_a}{\overset{x_b}{\int }}{N}_i\rho Abdx+{\left[N_{i}A{t}_x\right]}_{{\mathrm{\Gamma }}_t}]\\ & =\sum _{i=1}^n\delta v_i{f}_i^{\text{ext}}.\end{array}$

In matrix notation,

$\delta P_{\text{ext}}=\delta {𝐯}^T{𝐟}_{\text{ext}}\text{where}{𝐟}_{\text{ext}}=\left[\begin{array}{c}{f}_1^{\text{ext}}\\ f_2^{\text{ext}}\\ ⋮\\ f_n^{\text{ext}}\end{array}\right].$

The external force is given by

$f_i^{\text{ext}}=\underset{x_a}{\overset{x_b}{\int }}{N}_i\rho Abdx+{\left[N_{i}A{t}_x\right]}_{{\mathrm{\Gamma }}_t}.$

In matrix notation,

${𝐟}_{\text{ext}}=\underset{x_a}{\overset{x_b}{\int }}{𝐍}^T\rho Abdx+{\left[{𝐍}^{T}A{t}_x\right]}_{{\mathrm{\Gamma }}_t}\text{or}{𝐟}_{\text{ext}}=\underset{\mathrm{\Omega }}{\int }\rho {𝐍}^{T}bd\mathrm{\Omega }+{\left[{𝐍}^{T}A{t}_x\right]}_{{\mathrm{\Gamma }}_t}.$

Remark:

Using the conservation of mass

${\rho }_0{A}_{0}dx=\rho Adx$

and the relations between the traction sin the reference and the current configurations

$A_0{t}_x^0=A{t}_x$

we can transform the integral in the expression for the external force to one over the reference coordinates as follows:

${𝐟}_{\text{ext}}=\underset{x_a}{\overset{x_b}{\int }}{𝐍}^T{\rho }_0{A}_{0}bdx+{\left[{𝐍}^T{A}_0{t}_x^0\right]}_{{\mathrm{\Gamma }}_t}.$

The above is the same as the expression we had for the external force in the total Lagrangian formulation.

Discrete Equations for Updated Lagrangian

The finite element equations in updated Lagrangian form are:

$\begin{array}{cc}v(x,t)& =\sum _i{N}_i(x)v_i^e(t)=𝐍{𝐯}^e\\ D(x,t)& =\sum _i{N}_{i,x}{v}_i^e(t)=𝐁{𝐯}^e\\ {𝐟}_{\text{int}}^e& =\underset{{\mathrm{\Omega }}^e}{\int }{𝐁}^T\sigma d\mathrm{\Omega }\\ {𝐟}_{\text{ext}}^e& =\underset{{\mathrm{\Omega }}^e}{\int }\rho {𝐍}^{T}bd\mathrm{\Omega }+{\left[{𝐍}^{T}A{t}_x\right]}_{{\mathrm{\Gamma }}_t^e}\\ {𝐌}^e& =\underset{{\mathrm{\Omega }}_0^e}{\int }{\rho }_0{𝐍}^T𝐍d{\mathrm{\Omega }}_0\\ 𝐌\ddot{𝐮}& ={𝐟}_{\text{ext}}-{𝐟}_{\text{int}}.\end{array}$