Continuum Mechanics for CAE Engineers
Deformation Gradient, Strain & Stress Tensors, Nonlinear FEM Formulations
Continuum mechanics is the rigorous mathematical framework underlying all structural and fluid CAE simulations. Linear structural analysis courses often skip over the subtleties of large deformation — but as soon as you simulate a rubber seal, a metal forming process, a crash event, or biological tissue, you encounter multiple stress measures, conjugate strain pairs, and Lagrangian versus Eulerian descriptions. Mastering these concepts is what separates users who can set up nonlinear solvers correctly from those who get wrong answers without knowing it.
1. Kinematics: Motion & Deformation Gradient
I keep seeing the deformation gradient $\mathbf{F}$ mentioned in the Abaqus theory manual. What is it, and why is it more fundamental than just "strain"?
$\mathbf{F}$ is the most complete description of local deformation — it captures stretching, rotation, and shearing all in one tensor. Think of it as the answer to: "How does an infinitesimal material fiber transform when the body deforms?" If a tiny fiber was oriented as $d\mathbf{X}$ in the reference (undeformed) configuration, after deformation it becomes $d\mathbf{x} = \mathbf{F}\,d\mathbf{X}$. That single equation contains everything about local deformation. Strain measures like Green-Lagrange are then derived from $\mathbf{F}$ by stripping out the rotation part — because pure rotation creates no stress in the material, only stretching does.
1.1 Reference and Current Configurations
Let $\mathbf{X}$ denote the position of a material point in the reference (Lagrangian) configuration at time $t=0$, and $\mathbf{x} = \boldsymbol{\chi}(\mathbf{X},t)$ its position in the current (Eulerian) configuration. The displacement vector is:
1.2 Deformation Gradient
Key properties of $\mathbf{F}$:
- $J = \det(\mathbf{F}) > 0$ — the Jacobian determinant equals the volume ratio $dv/dV$
- $J = 1$ for incompressible materials (rubber, soft tissue, plastic deformation)
- $\mathbf{F}$ is not symmetric in general
- In the reference state: $\mathbf{F} = \mathbf{I}$
1.3 Polar Decomposition
$\mathbf{F}$ can be uniquely decomposed as:
where $\mathbf{R}$ is a proper orthogonal rotation tensor ($\mathbf{R}^T\mathbf{R} = \mathbf{I}$, $\det\mathbf{R} = 1$), $\mathbf{U}$ is the right stretch tensor (symmetric, positive definite), and $\mathbf{V}$ is the left stretch tensor. The eigenvalues of $\mathbf{U}$ are the principal stretches $\lambda_1, \lambda_2, \lambda_3$.
Why does it matter that $\mathbf{F} = \mathbf{R}\mathbf{U}$? What does a solver actually use this for?
The polar decomposition is used to define objective (frame-invariant) stress rates and to ensure material models don't spuriously generate stress from rigid body rotations. In corotational formulations — used extensively in crash codes like LS-DYNA — the stress is computed in the rotated frame $\mathbf{R}^T\boldsymbol{\sigma}\mathbf{R}$ where the rotation has been factored out. Without this, a material element that undergoes large rigid body rotation would accumulate artificial shear stresses from the time integration of the constitutive law. The Jaumann, Green-Naghdi, and Truesdell stress rates all involve $\mathbf{R}$ in different ways.
2. Strain Measures
2.1 Green-Lagrange Strain (Material Frame)
The Green-Lagrange strain tensor is defined in the reference configuration and is zero for rigid body motion:
where $\mathbf{C} = \mathbf{F}^T\mathbf{F}$ is the right Cauchy-Green deformation tensor. In component form with displacement gradients $H_{iJ} = \partial u_i/\partial X_J$:
The nonlinear term $H_{kI}H_{kJ}$ is what distinguishes Green-Lagrange from the small-strain engineering strain tensor. For small deformations, $E_{IJ} \approx \varepsilon_{IJ} = \frac{1}{2}(H_{IJ} + H_{JI})$.
2.2 Almansi Strain (Spatial Frame)
The Almansi (Eulerian) strain tensor is defined in the current configuration:
where $\mathbf{B} = \mathbf{F}\mathbf{F}^T$ is the left Cauchy-Green (Finger) tensor. The Almansi strain is used in spatial formulations and is conjugate to the Cauchy stress in the spatial description.
2.3 Logarithmic (Hencky) Strain
The logarithmic or true strain is the most natural measure for large elastic-plastic deformation:
where $\lambda_a$ are principal stretches and $\mathbf{n}_a$ are eigenvectors of $\mathbf{V}$. In 1D: $\varepsilon_H = \ln(L/L_0) = \ln(1 + e)$ where $e = (L-L_0)/L_0$ is engineering strain. The Hencky strain is additive for sequential deformations and is the preferred measure for metal plasticity (used in Abaqus, Ansys, LS-DYNA output as "logarithmic strain").
| Strain Measure | Frame | Formula | Best Used For |
|---|---|---|---|
| Engineering strain $\varepsilon$ | Reference | $(L-L_0)/L_0$ | Small deformations |
| Green-Lagrange $\mathbf{E}$ | Reference (material) | $\frac{1}{2}(\mathbf{C}-\mathbf{I})$ | TL formulation, hyperelasticity |
| Almansi $\mathbf{e}$ | Current (spatial) | $\frac{1}{2}(\mathbf{I}-\mathbf{B}^{-1})$ | UL formulation |
| Logarithmic (Hencky) $\boldsymbol{\varepsilon}_H$ | Current | $\ln\mathbf{V}$ | Metal plasticity, rubber |
3. Stress Measures: Cauchy, PK1, PK2
Abaqus outputs "Cauchy stress" and the theory manual mentions "PK2 stress" in the constitutive law section. Why do we need more than one kind of stress?
Each stress measure tells the same physical story but in a different reference frame, and some are more convenient for certain mathematical operations. Cauchy stress is what you physically feel — force per unit area on the deformed surface. It's what you'd measure with a pressure gauge. But for constitutive laws written in the reference configuration (like hyperelastic models for rubber), you need stress referred to the original undeformed area. That's why PK2 stress exists — it's symmetric, lives in the reference frame, and is the natural conjugate to Green-Lagrange strain for energy calculations. PK1 is intermediate — force in the current config per reference area — and it's non-symmetric, which makes it awkward for constitutive work.
3.1 Cauchy Stress ($\boldsymbol{\sigma}$)
The Cauchy (true) stress tensor is defined in the current configuration. Traction on a surface with unit normal $\mathbf{n}$ in the current configuration:
$\boldsymbol{\sigma}$ is symmetric due to moment equilibrium: $\sigma_{ij} = \sigma_{ji}$. This is what Abaqus outputs as S11, S22, S33, S12, S13, S23.
3.2 First Piola-Kirchhoff Stress (PK1, $\mathbf{P}$)
PK1 maps reference-configuration area normals to current-configuration forces:
PK1 is not symmetric: $P_{iJ} \neq P_{Ji}$. It appears naturally in the weak form for Total Lagrangian formulations but is rarely used directly in constitutive laws.
3.3 Second Piola-Kirchhoff Stress (PK2, $\mathbf{S}$)
PK2 is fully material-frame, symmetric, and energy-conjugate to Green-Lagrange strain:
The internal virtual work per unit reference volume is:
For hyperelastic materials (rubber, biological tissue), the constitutive law is most cleanly expressed as:
where $W(\mathbf{C})$ is the strain energy density function (e.g., Neo-Hookean: $W = \frac{\mu}{2}(I_1 - 3) + \frac{\kappa}{2}(J-1)^2$).
3.4 Stress Transformations Summary
4. Constitutive Frameworks
4.1 Hyperelasticity
A hyperelastic material stores and fully recovers strain energy. The strain energy density $W$ depends only on the deformation state, not on the loading path. Common models:
| Model | $W(I_1, I_2, J)$ | Best For |
|---|---|---|
| Neo-Hookean | $\frac{\mu}{2}(I_1-3) + \frac{\kappa}{2}(J-1)^2$ | Moderate strains, simple rubber |
| Mooney-Rivlin | $C_{10}(I_1-3) + C_{01}(I_2-3)$ | Rubber up to ~150% strain |
| Ogden | $\sum_p \frac{\mu_p}{\alpha_p}(\lambda_1^{\alpha_p}+\lambda_2^{\alpha_p}+\lambda_3^{\alpha_p}-3)$ | Large strains, biological tissue |
| Arruda-Boyce | 8-chain network model | Polymer networks with locking |
4.2 Elastoplasticity at Finite Strains
For metals undergoing large plastic deformation (sheet forming, crash), the multiplicative decomposition is used:
The elastic part $\mathbf{F}^e$ stores energy; the plastic part $\mathbf{F}^p$ is irreversible and volume-preserving ($\det\mathbf{F}^p = 1$ for metals). The Mandel stress $\boldsymbol{\Sigma} = \mathbf{C}^e\mathbf{S}^e$ drives plastic flow in the intermediate configuration.
4.3 Frame Objectivity
A constitutive law must be frame-indifferent (objective): the material response cannot depend on the superimposed rigid body motion of the observer. If the current configuration undergoes a rotation $\mathbf{Q}$, the Cauchy stress must transform as:
This requirement rules out naive rate equations like $\dot{\boldsymbol{\sigma}} = \mathbb{C}:\mathbf{d}$ (which is not objective) and requires use of objective stress rates such as the Jaumann rate $\overset{\nabla}{\boldsymbol{\sigma}} = \dot{\boldsymbol{\sigma}} - \mathbf{W}\boldsymbol{\sigma} + \boldsymbol{\sigma}\mathbf{W}$.
5. Total vs. Updated Lagrangian Formulations
In Abaqus documentation I see "Total Lagrangian" and "Updated Lagrangian." What's the practical difference when I'm setting up a rubber seal compression simulation?
Both formulations give the same final answer for the same problem — they're mathematically equivalent when implemented correctly. The difference is purely in the reference configuration used at each step. Total Lagrangian (TL) always refers quantities back to the original undeformed configuration at $t=0$. Updated Lagrangian (UL) resets the reference to the current deformed configuration at the start of each increment. For your rubber seal, Abaqus uses UL internally for efficiency — the stiffness matrix and residual are computed in the current configuration, which avoids accumulating large displacement gradients in $\mathbf{F}$. But the constitutive law is often written in a material (TL) frame using PK2/Green-Lagrange pairs. The solver handles the conversion.
5.1 Total Lagrangian (TL) Formulation
All integrals are over the reference volume $V_0$. The principle of virtual work:
Advantages: constitutive law always uses the same reference; Green-Lagrange strain is straightforward to compute; suitable for hyperelastic materials. Disadvantage: for very large deformations, the mapping from reference to current is computationally expensive.
5.2 Updated Lagrangian (UL) Formulation
The reference configuration is updated to the current configuration $V_n$ at the start of each increment $n$. The principle of virtual work becomes:
where $\boldsymbol{\tau}$ is the Kirchhoff stress ($\boldsymbol{\tau} = J\boldsymbol{\sigma}$) and $\mathbf{d}$ is the rate-of-deformation tensor. Used in LS-DYNA, Abaqus implicit and explicit for most solid element formulations.
| Aspect | Total Lagrangian | Updated Lagrangian |
|---|---|---|
| Reference config. | Fixed: $t=0$ state | Updated: $t_n$ state each step |
| Stress measure | PK2 ($\mathbf{S}$) | Cauchy ($\boldsymbol{\sigma}$) or Kirchhoff ($\boldsymbol{\tau}$) |
| Strain measure | Green-Lagrange ($\mathbf{E}$) | Almansi ($\mathbf{e}$) or rate-of-deformation ($\mathbf{d}$) |
| Best for | Rubber, hyperelastic, large elastic strains | Metal plasticity, crash, large plastic strains |
6. Nonlinear FEM Implementation
6.1 Sources of Nonlinearity
Three distinct sources of nonlinearity arise in structural FEM:
- Material nonlinearity: Plastic yielding, hyperelasticity, viscoelasticity, damage
- Geometric nonlinearity: Large displacements/rotations, pre-stress (membrane/cable stiffening), buckling
- Contact nonlinearity: Changing boundary conditions, gap elements, friction
6.2 Newton-Raphson Iteration
The nonlinear equilibrium residual $\mathbf{R}(\mathbf{u}) = \mathbf{f}_{\text{ext}} - \mathbf{f}_{\text{int}}(\mathbf{u}) = \mathbf{0}$ is solved iteratively. At iteration $k$:
The tangent stiffness $\mathbf{K}_T$ for the UL formulation contains three contributions:
where $\mathbf{K}_{\text{geo}}$ is the geometric (stress) stiffness matrix — the term responsible for buckling and pre-stress stiffening.
6.3 Convergence Criteria
Standard convergence checks (Abaqus defaults):
- Force residual: $\|\mathbf{R}^{(k)}\| / \|\mathbf{f}_{\text{ref}}\| < 5 \times 10^{-3}$
- Displacement correction: $\|\Delta\mathbf{u}^{(k)}\| / \|\mathbf{u}^{(k)}\| < 10^{-2}$
- Energy criterion: $\Delta\mathbf{u}^{(k)} \cdot \mathbf{R}^{(k)} / \Delta\mathbf{u}^{(0)} \cdot \mathbf{R}^{(0)} < 10^{-5}$
My Abaqus simulation keeps failing to converge at high compressive loads on a rubber gasket. What should I check first?
For rubber compression, check these in order. First, is your hyperelastic model calibrated from real test data? A Neo-Hookean model fitted to only uniaxial tension data often goes unstable under compression because the biaxial/volumetric behavior isn't constrained. Second, check your mesh — highly distorted elements as the rubber is squeezed will cause Jacobian $J \to 0$ issues. Try using hybrid elements (e.g., C3D8H in Abaqus) which handle near-incompressibility better. Third, reduce the initial increment size and enable automatic incrementation. Fourth, check the contact definition — if the gasket is sliding as well as compressing, sticky versus frictionless contact makes a big difference to the convergence behavior.