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) = \cfrac{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

\int_0^t 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~A~F = \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~A~b = \rho~A_0~\cfrac{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~A~b = \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~. $$

Total Form
For a linear elastic material:

{ \sigma(x,t) = E^{\sigma D} \int_0^t 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{align} \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{align} } $$

The  essential boundary conditions are

{ v(x,t) = \bar{v}(x,t) \qquad ~\text{for}~x \in \Gamma_v ~. } $$ The  traction boundary conditions are

{ n~\sigma(x,t) = t_x(x,t) \qquad ~\text{for}~x \in \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~A~b = \rho~A~\cfrac{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$$).

\int_{x_a}^{x_b} w\left[(A~\sigma)_{,x} + \rho~A~b - \rho~A~\cfrac{Dv}{Dt}\right]~dx = 0 ~. $$ Integrate the first term by parts to get

\int_{x_a}^{x_b} w (A~\sigma)_{,x}~dx = \left[w~A~\sigma\right]_{x_a}^{x_b} - \int_{x_a}^{x_b} w_{,x}~(A~\sigma)~dx ~. $$ Plug the above into the weak form and rearrange to get

\int_{x_a}^{x_b} w_{,x}~(A~\sigma)~dx + \int_{x_a}^{x_b} w~\rho~A~\cfrac{Dv}{Dt}~dx = \int_{x_a}^{x_b} w~\rho~A~b~dx + \left[w~A~\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:

{ \int_{x_a}^{x_b} (\delta v)_{,x} (A~\sigma)~dx + \int_{x_a}^{x_b}\delta v~\rho~A~\cfrac{Dv}{Dt}~dx = \int_{x_a}^{x_b} \delta v~\rho~A~b~dx + \left[\delta v~A~\sigma\right]_{x_a}^{x_b} ~. } $$ Recall that

v_{,x} = D(x,t) \qquad\text{and}\qquad t_x = n~\sigma ~. $$ Therefore,

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

\left[\delta v~A~\sigma\right]_{x_a}^{x_b} = \left[\delta v~A~t_x\right]_{\Gamma_t} ~. $$ Hence we can alternatively write the weak form as

{ \int_{x_a}^{x_b}\delta D~A~\sigma~dx + \int_{x_a}^{x_b}\delta v~\rho~A~\cfrac{Dv}{Dt}~dx = \int_{x_a}^{x_b}\delta v~\rho~A~b~dx +\left[\delta v~A~t_x\right]_{\Gamma_t}~. } $$ Comparing with the energy equation, we see that the first term above is the internal virtual power. Using the substitutions $$A~dx = d\Omega$$ and $$Dv/DT = \dot{v}$$, we can also write the weak form as

{ \int_{\Omega}\delta D~\sigma~d\Omega + \int_{\Omega}\delta v~\rho~\dot{v}~d\Omega = \int_{\Omega}\delta v~\rho~b~d\Omega +\left[\delta v~A~t_x\right]_{\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{align} \delta P_{\text{int}} & = \int_{x_a}^{x_b} (\delta v)_{,X}~(A~\sigma)~dx \\ \delta P_{\text{ext}} & = \int_{x_a}^{x_b} \delta v~\rho~A~b~dx + \left[\delta v~A~t_X\right]_{\Gamma_t} \\ \delta P_{\text{kin}} & = \int_{x_a}^{x_b} \delta v~\rho~A~\dot{v}~dx \end{align} } $$

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) = \mathbf{N}(x)~\mathbf{v}(t) } \text{where} \mathbf{N} = \begin{bmatrix} N_1(X) & N_2(X) && N_n(X) \end{bmatrix} ~\text{and}~ \mathbf{v} = \begin{bmatrix} v_1(t) \\ v_2(t) \\ \vdots \\ v_n(t) \end{bmatrix} ~. $$ 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) = \mathbf{N}(x)~\mathbf{a}(t) } \text{where} \mathbf{a} = \begin{bmatrix} a_1(t) \\ a_2(t) \\ \vdots \\ a_n(t) \end{bmatrix} ~. $$ 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 = \varphi^{-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) = \mathbf{N}~\delta \mathbf{v} } \text{where} \delta \mathbf{v} = \begin{bmatrix} \delta v_1 \\ \delta v_2 \\ \vdots \\ \delta v_n \end{bmatrix} ~. $$ 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} = \mathbf{B}~\delta \mathbf{v} } \text{where} \mathbf{B} = \begin{bmatrix} N_{1,x} & N_{2,x} && N_{n,x} \end{bmatrix} ~. $$ 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.

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)\qquad \text{or} { x(\xi, t) = \mathbf{N}(\xi)~\mathbf{x}^e(t) ~. } $$ At $$t = 0$$,

x(\xi, 0) = x(\xi) = \mathbf{N}(\xi)~\mathbf{X}^e ~. $$ Therefore the  displacement in the parent coordinates is

u(\xi, t) = x(\xi, t) - x(\xi) = \mathbf{N}(\xi)~\left[\mathbf{x}^e(t) - \mathbf{X}^e\right] \qquad \text{or} { u(\xi, t) = \mathbf{N}(\xi)~\mathbf{u}^e(t) ~. } $$ Similarly, the  velocity and its variation in the parent coordinates are given by

{ v(\xi, t) = \mathbf{N}(\xi)~\mathbf{v}^e(t) ~,\qquad \delta v(\xi, t) = \mathbf{N}(\xi)~\delta \mathbf{v}^e(t) ~. } $$ The  acceleration in the parent coordinates is given by

{ a(\xi, t) = \mathbf{N}(\xi)~\dot{\mathbf{v}}^e(t) ~. } $$ The  rate of deformation is given by

D(x, t) = v_{,x} (x, t) = \mathbf{N}_{,x}(X(x,t))~\mathbf{v}^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} \qquad\text{and}\qquad N_{2,\xi} = N_{2,x}~x_{,\xi}~. $$ In matrix form,

{ \mathbf{N}_{,\xi} = \mathbf{N}_{,x}~x_{,\xi} } \qquad\implies\qquad { \mathbf{N}_{,x} = \mathbf{N}_{,\xi}~\left(x_{,\xi}\right)^{-1} =: \mathbf{B}(\xi)~. } $$ The  rate of deformation  in parent coordinates is then given by

{ D(\xi, t) = \mathbf{B}(\xi)~\mathbf{v}^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

\int_{x_a}^{x_b} (\delta v_h)_{,x} (A~\sigma)~dx + \int_{x_a}^{x_b}\delta v_h~\rho~A~\cfrac{Dv_h}{Dt}~dx = \int_{x_a}^{x_b} \delta v_h~\rho~A~b~dx + \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}} = \int_{x_a}^{x_b} (\delta v_h)_{,x} (A~\sigma)~dx ~. $$ Plugging in the derivative of the test function, we get

\begin{align} \delta P_{\text{int}} & = \int_{x_a}^{x_b} \left[\sum_{i=1}^n N_{i,x}\delta v_i\right] (A~\sigma)~dx\\ & = \sum_{i=1}^n \delta v_i \left[\int_{x_a}^{x_b} N_{i,x}(A~\sigma)~dx\right]\\ & = \sum_{i=1}^n \delta v_i f_i^{\text{int}} ~. \end{align} $$ The matrix form of the expression for virtual internal work is

{ \delta P_{\text{int}} = \delta \mathbf{v}^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~\sigma)~dx = \int_{-1}^{1} N_{i,\xi}(x_{,\xi})^{-1}(A~\sigma)~x_{,\xi}~d\xi = \int_{-1}^{1} N_{i,\xi}(A~\sigma)~d\xi ~. } $$ Note that the above simplification only occurs in 1-D. In matrix form,

{ \mathbf{f}_{\text{int}} = \int_{x_a}^{x_b} \mathbf{B}^T~(A~\sigma)~dx = \int_{-1}^{1} (\mathbf{N}_{,\xi})^T~(A~\sigma)~d\xi~. }\text{or} { \mathbf{f}_{\text{int}} = \int_{\Omega} \mathbf{B}^T~\sigma~d\Omega ~. } $$

 Remark:

Note that we can write

N_{,x} = N_{,X}~\cfrac{dX}{dx} \qquad \implies \qquad N_{,x}~dx = N_{,X}~dX ~. $$ If we use the above relation, we get

\mathbf{f}_{\text{int}} = \int_{x_a}^{x_b} (\mathbf{N}_{,x})^T~\sigma~A~dx = \int_{X_a}^{X_b} (\mathbf{N}_{,X})^T~\sigma~A~dX ~. $$ Using the relation

\sigma~A = P~A_0 $$ we get

{ \mathbf{f}_{\text{int}} = \int_{X_a}^{X_b} (\mathbf{N}_{,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}} = \int_{x_a}^{x_b}\delta v_h~\rho~A~\cfrac{Dv_h}{Dt}~dx ~. $$ Plugging in the test function, we get

\begin{align} \delta P_{\text{kin}} & = \int_{x_a}^{x_b}\left[\sum_{i=1}^n \delta v_i N_i\right] \rho~A~\cfrac{Dv_h}{Dt}~dx\\ & = \sum_{i=1}^n \delta v_i \left[\int_{x_a}^{x_b} \rho~A~N_i~\cfrac{Dv_h}{Dt}~dx\right]\\ & = \sum_{i=1}^n \delta v_i f_i^{\text{kin}} ~. \end{align} $$ The inertial (kinetic) force is

{ f^{\text{kin}}_i = \int_{x_a}^{x_b}\rho~A~N_i~\cfrac{Dv_h}{Dt}~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~A~N_i \left[\sum_{j=1}^n \cfrac{Dv_j}{Dt} N_j\right]~dx \\ & = \sum_{j=1}^n \left[\int_{x_a}^{x_b}\rho~A~N_i~N_j~dx\right] \cfrac{Dv_j}{Dt} \\ & = \sum_{j=1}^n M_{ij} a_j \end{align} $$ where

{ M_{ij} := \int_{x_a}^{x_b}\rho~A~N_i~N_j~dx \qquad \leftarrow \quad \text{Consistent Mass Matrix Coefficient.} } $$ and

a_j := \cfrac{Dv_j}{Dt} = \dot{v}_j \qquad\leftarrow \quad\text{Acceleration.} $$ 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} \dot{v}_1 \\ \dot{v}_2 \\ \vdots \\ \dot{v}_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~A~\mathbf{N}^T~\mathbf{N}~dx = \int_{\Omega}\rho~\mathbf{N}^T~\mathbf{N}~d\Omega ~. } $$ Plugging the expression for the inertial force into the expression for virtual kinetic power, we get

\begin{align} \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{align} $$ In matrix form,

{ \delta P_{\text{kin}} = (\delta \mathbf{v})^T \mathbf{M} \mathbf{a} = (\delta \mathbf{v})^T~\mathbf{f}_{\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~A~dx ~. $$ Therefore, we can write the mass matrix as

{ \mathbf{M} =\int_{x_a}^{x_b}\rho~A~\mathbf{N}^T~\mathbf{N}~dx = \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 ~. } $$ 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}} = \int_{x_a}^{x_b}\delta v_h~\rho~A~b~dx + \left[\delta v_h~A~t_x\right]_{\Gamma_t}~. $$ Plugging in the test function into the above expression gives

\begin{align} \delta P_{\text{ext}} & = \int_{x_a}^{x_b}\left[\sum_{i=1}^n N_i~\delta v_i\right]\rho~A~b~dx + \left[\left(\sum_{i=1}^n N_i~\delta v_i\right) A~t_x\right]_{\Gamma_t} \\ & = \sum_{i=1}^n \delta v_i\left[\int_{x_a}^{x_b} N_i~\rho~A~b~dx\right] + \left[\sum_{i=1}^n \delta v_i~N_i~A~t_X\right]_{\Gamma_t} \\ & = \sum_{i=1}^n \delta v_i\left[\int_{x_a}^{x_b} N_i~\rho~A~b~dx + \left[N_i~A~t_x\right]_{\Gamma_t}\right] \\ & = \sum_{i=1}^n \delta v_i~f^{\text{ext}}_i ~. \end{align} $$ In matrix notation,

{ \delta P_{\text{ext}} = \delta \mathbf{v}^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~A~b~dx + \left[N_i~A~t_x\right]_{\Gamma_t}~. } $$ In matrix notation,

{ \mathbf{f}_{\text{ext}} = \int_{x_a}^{x_b} \mathbf{N}^T~\rho~A~b~dx + \left[\mathbf{N}^T~A~t_x\right]_{\Gamma_t} }\text{or} { \mathbf{f}_{\text{ext}} = \int_{\Omega} \rho~\mathbf{N}^T~b~d\Omega + \left[\mathbf{N}^T~A~t_x\right]_{\Gamma_t}~. } $$

 Remark:

Using the conservation of mass

\rho_0~A_0~dx = \rho~A~dx $$ 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:

\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}~. $$ 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{align} v(x,t) & = \sum_i N_i(x)~v_i^e(t) = \mathbf{N}~\mathbf{v}^e \\ D(x,t) & = \sum_i N_{i,x}~v_i^e(t) = \mathbf{B}~\mathbf{v}^e \\ \mathbf{f}^e_{\text{int}} & = \int_{\Omega^e} \mathbf{B}^T~\sigma~d\Omega \\ \mathbf{f}^e_{\text{ext}} & = \int_{\Omega^e} \rho~\mathbf{N}^T~b~d\Omega + \left[\mathbf{N}^T~A~t_x \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} } $$