Updated Lagrangian formulation for geometrically non-linear elastic problems

Introduction

This post introduces an implementation of the updated-Lagrangian formulation within the finite element method. The aim is to tackle geometrically non-linear problems involving significant deformations. The provided example revolves around a basic 3D cantilever beam scenario.

Delving into non-linear analysis introduces considerations like divergence issues and heightened problem complexity. Nevertheless, in certain scenarios such as buckling instabilities, resorting to non-linear finite element solvers becomes imperative. This raises inquiries into the disparity between results from linear and non-linear problems, as well as the extent to which reasonable predictions can be derived from the more simplified linearized problem.

In this context, a straightforward approach rooted in the updated-Lagrangian method is unveiled, focusing on the cantilever beam case. This approach effectively highlights the constraints of relying solely on linearization.

Formulations

Kinematics

Imagine the transformation of a solid object from its original shape $\Omega_0$ to its ultimate form $\Omega$, as depicted in the accompanying diagram.

In this context, $\vec{X}$ signifies the material coordinate representing the initial, undeformed arrangement of a point, while $\vec{x}$ denotes the spatial coordinate that represents the current, altered configuration of that point. The shift in position for any given point can be elucidated through the displacement vector $\vec{u} = \vec{x} - \vec{X}$.

A valuable gauge of this transformation is the deformation gradient,

\[\begin{split} \boldsymbol{F} &= \frac{\partial \vec{x}}{\partial \vec{X}}\\ F_{ij} & = \frac{\partial x_i}{\partial X_j} \end{split}\]

a parameter that can also be articulated in relation to the displacement field:

\[\begin{split} \boldsymbol{F} &= \boldsymbol{I} + \boldsymbol{\nabla}_0 \vec{u}\\ F_{ij} &= \delta_{ij} + \frac{\partial u_i}{\partial X_j} \\ \end{split}\]

To anticipate the distribution of displacement, strain, and stress within the domain $\Omega$ under specific loading circumstances, it is imperative to initially lay the foundation for the deformation’s kinematics and the constitutive equations governing the relationship between stress and strain. To embark upon this, we first address the nuances of accurate and linearized strain measurement.

The initial divergence between linear and non-linear models surfaces through the characterization of deformation measurements. An accurate measure of deformation should inherently possess objectivity, remaining invariant irrespective of the chosen coordinate system. Among the commonly employed measures of deformation, the Lagrangian strain emerges as one noteworthy example:

\[\begin{split} \boldsymbol{E} &= \frac{1}{2}\left[\left(\boldsymbol{\nabla}_0 \vec{u}\right)^T + \boldsymbol{\nabla}_0 \vec{u} + \boldsymbol{\nabla}_0\vec{u} \cdot \left(\boldsymbol{\nabla}_0 \vec{u}\right)^T\right]\\ E_{ij} &= \frac{1}{2} \left[\frac{\partial u_i}{\partial X_j} + \frac{\partial u_j}{\partial X_i} + \frac{\partial u_k}{\partial X_i} \frac{\partial u_k}{\partial X_j}\right] \end{split}\]

It can be easily derived that \(\begin{split} \boldsymbol{E} &= \frac{1}{2}\left(\boldsymbol{F}^T \boldsymbol{F} - \boldsymbol{I}\right)\\ E_{ij} &= \frac{1}{2}\left(F_{ki}F_{kj} - \delta_{ij}\right) \end{split}\)

Observe the presence of non-linear terms within this equation. Given the premise of minor deformations, these terms are expected to be significantly smaller in magnitude when contrasted with the remaining components of the equation. Hence, adopting this perspective allows us to omit these non-linear terms. Consequently, we arrive at a linearized strain measurement:

\[\begin{split} \boldsymbol{\varepsilon} &= \frac{1}{2}\left(\boldsymbol{F} + \boldsymbol{F}^T\right) - \boldsymbol{I}\\ \varepsilon_{ij} &= \frac{1}{2}\left(F_{ij} + F_{ji}\right) - \delta_{ij} \end{split}\]

Constitutive model

A simple linear constitutive model (Hook’s law) can be established \(\boldsymbol{S}^{2PK} = \boldsymbol{C} : \boldsymbol{E}\)

which for an isotropic material model can be written in the following form: \(S^{2PK}_{ij} = \frac{E}{1+\nu} \left[E_{ij} + \frac{\nu}{1-2\nu} E_{kk}\delta_{kk}\right]\)

where $E$ is the Young’s modulus and $\nu$ is the Poisson’s ratio. As before $\boldsymbol{E}$ is the Green-Lagrange strain and $\boldsymbol{S}^{2PK}$ is the second Piola-Kirchhoff stress. The relationship between second Piola-Kirchhoff stress and Cauchy stress can be established: $\boldsymbol{\sigma} = J^{-1}\boldsymbol{F}\boldsymbol{S}^{2PK}\boldsymbol{F}^T$, where $J$ is the Jacobian of deformation $J = \det \boldsymbol{F}$.

Similarly, for the linearized case, one can establish the following constitutive model: \(\boldsymbol{\sigma} = \boldsymbol{C}:\boldsymbol{\varepsilon}\)

where $\boldsymbol{\sigma}$ is a measure of stress.

Mechanical equilibrium

The fundamental equation at play here is the conservation of linear momentum. For the purpose of this context, we consider a scenario where inertial forces can be disregarded. This leads us to the following simplified expression:

\[\begin{split} \boldsymbol{\nabla} \cdot \boldsymbol{\sigma} + \rho \boldsymbol{b} &= \boldsymbol{0}\\ \frac{\partial \sigma_{ij}}{\partial x_j} + \rho b_i &= 0 \end{split}\]

Updated Lagrangian formulation

One way to handle a non-linear problem is to solve it as a sequence of small linear problems. Here, the approach is to solve the governing equations under small load increment and use the updated domain $\vec{X} + \vec{u}_i$ as the next domain configuration. Thus, the finite element mesh is updated at every time step. We can show that by such linearization the solution to the non-linear problem can be obtained. A schmeatic of a the evolution of the domain from is initial configuration $\Omega_0$ to its final configuration $\Omega$ in a sequence of small loadsteps is shown in the figure below.

One possible strategy to manage a non-linear problem involves approaching it through a series of incremental linear problems. In this methodology, the governing equations are addressed stepwise, with each step corresponding to a small increase in the applied load. The evolving domain, denoted as $\vec{X} + \vec{u}_i$, becomes the updated configuration for subsequent steps, leading to the adjustment of the finite element mesh at each time increment. Notably, this approach’s efficacy lies in its ability to yield a solution for the non-linear problem by virtue of this incremental linearization.

The visual representation below depicts the progression of the domain from its initial configuration $\Omega_0$ to its final state $\Omega$, realized through a sequence of minor load increments. It’s important to highlight that the crux of this linearization technique hinges on expressing dependent variables, like stress $\boldsymbol{\sigma}$, in terms of material coordinates, i.e., $\boldsymbol{\sigma} = \boldsymbol{\sigma}(\vec{X})$, while derivative operations are carried out relative to spatial coordinates. Thus, in the implementation of the finite element method, a fundamental assumption is the relatively modest disparity between the current and preceding configurations.

Below, an illustrative instance involving a cantilever beam is provided.

Solution - cantilever beam

For an illustration purposes, a simple cantilever beam simulation is considered. The beams’ length is of 10 units with a square cross-sectional area measuring 1 by 1 unit. The domain is discretized into $200\times 20\times 20$ hexahedral finite elements. The beam experiences loading by subjecting each node of the finite element mesh to an incremental force $\Delta f_x = 0.125 \cdot 10^{-5}$. Pertinent material properties include a Young’s modulus of 210 GPa and a Poisson’s ratio of 0.3. The non-linear analysis is undertaken over the course of 10 incremental (linear) load steps. The outcomes corresponding to each load step are aggregated and depicted in the subsequent diagram.

In constrast, the results of the linear analysis are presented in the following figure at the same load step intervals.

Noteworthy disparities become evident between the two scenarios. In the non-linear solution, the dominant deformation mechanism involves the bending of the cantilever beam, which leads to a curved trajectory for the endpoint—manifesting as movement along a circular path. In contrast, the linearized solution results in the endpoint remaining at a consistent $z$ coordinate while tracing a linear path along the $x$ direction. This divergence in behavior is a core consequence of linearization.

Additionally, it’s crucial to observe that in the linearized case, the volume of the beam appears to experience significant expansion. This phenomenon is linked to the linearized treatment causing the endpoint’s movement along a line. This, in turn, can give rise to the generation of artificial stresses that should not exist if rotational effects were comprehensively considered. It’s important to note, however, that when dealing with minor deformations, the linear and non-linear solutions exhibit close resemblance.

The subsequent inquiry regarding the optimal number of load steps needed can be addressed through a sensitivity analysis. By conducting simulations across varying quantities of load steps, the outcomes pertaining to the ultimate deformed configuration are compiled and illustrated in the subsequent figure.

Evidently, the results distinctly converge with just 10 load steps, establishing that this quantity is adequate for the analysis.

Conclusions

This post delves into addressing geometrically non-linear problems through a relatively simple numerical technique associated with finite element simulations. It explores the implications of linearization, shedding light on the circumstances where a linear solver suffices and when the adoption of a non-linear solver becomes imperative.