General materials in rate/incremental form

For materials defined by rate-constitutive equations (e.g., hypoplasticity, plasticity) or incremental formulations, numgeo utilizes the objective integration approach proposed by Hughes & Winget (1980)¹. This method ensures that the stress update is invariant to rigid body rotations.

Stress update

The Cauchy stress \(\boldsymbol{\sigma}^{t+\Delta t}\) at time \(t + \Delta t\) is calculated by rotating the stress from the previous configuration to the current configuration and adding the material stress increment:

\[ \boldsymbol{\sigma}^{t+\Delta t} = \boldsymbol{Q}_{1/2}\boldsymbol{\sigma}^{t}\boldsymbol{Q}_{1/2}^T + \Delta \boldsymbol{\sigma} \]

where \(\boldsymbol{\sigma}^{t}\) is the converged stress from the previous step, \(\boldsymbol{Q}_{1/2}\) is the incremental rotation matrix, and \(\Delta \boldsymbol{\sigma}\) is the constitutive response.

Incremental Rotation

The rotation matrix \(\boldsymbol{\sigma}^{t}\) is determined at the midpoint configuration using the Cayley transform (a Crank-Nicolson approximation of the exponential map). This preserves the orthogonality of the rotation matrix:

\[ \boldsymbol{Q}_{1/2} = \left( \boldsymbol{I} - \frac{1}{2} \Delta t~\boldsymbol{\omega}_{1/2} \right)^{-1} \left( \boldsymbol{I} + \frac{1}{2} \Delta t~\boldsymbol{\omega}_{1/2} \right) \]

where \(\boldsymbol{I}\) is the identity matrix and \(\boldsymbol{\omega}_{1/2}\) is the spin tensor (skew-symmetric part of the velocity gradient):

\[ \boldsymbol{\omega}_{1/2} = \frac{1}{2}\left(\boldsymbol{l}_{1/2}-\boldsymbol{l}_{1/2}^T\right) \]

The velocity gradient \(\boldsymbol{l}_{1/2}\) is approximated using the displacement increment \(\Delta \boldsymbol{u}\) evaluated on the midpoint geometry:

\[ \boldsymbol{l}_{1/2} \approx \frac{\partial (\Delta \boldsymbol{u})}{\partial \boldsymbol{x}_{1/2}} \frac{1}{\Delta t} \]

\[ \boldsymbol{x}_{1/2} = \frac{1}{2}\left(\boldsymbol{x}^{t+\Delta t} + \boldsymbol{x}^{t}\right) \]

Strain Increment

The stress increment \(\Delta \boldsymbol{\sigma}\) is obtained from the constitutive model. The model is called using a "rotation-neutralized" strain increment \(\Delta \boldsymbol{\epsilon}\), often referred to as the logarithmic strain increment in this context:

\[ \Delta \boldsymbol{\epsilon} = \text{sym}\left( \frac{\partial (\Delta \boldsymbol{u})}{\partial \boldsymbol{x}_{1/2}} \right) = \boldsymbol{B}_{1/2} \Delta \boldsymbol{u} \]

Here, \(\boldsymbol{B}_{1/2}\) is the strain-displacement matrix evaluated using the midpoint coordinates \(\boldsymbol{x}_{1/2}\).

Geometric stiffness

To achieve quadratic convergence in the Newton-Raphson scheme, the tangent stiffness matrix must account for geometric changes. The total tangent stiffness consists of a material part and a geometric (initial stress) part:

\[ \boldsymbol{K}^{tan} = \boldsymbol{K}^\sigma + \boldsymbol{K}^{geo} \]

Material Stiffness

The material stiffness arises from the linearization of the constitutive law:

\[ \boldsymbol{K}^\sigma = \int_V \boldsymbol{B}^T \boldsymbol{D} \boldsymbol{B}~\text{d} v \]

where: * \(\boldsymbol{B}\) is the strain-displacement matrix defined on the current geometry \(\boldsymbol{x}\). * \(\boldsymbol{D}\) is the spatial tangent modulus consistent with the objective rate used in the stress update (approximating the Jaumann rate in the Hughes-Winget formulation).

Geometric Stiffness

The geometric stiffness accounts for the effect of existing stresses on the equilibrium equations during deformation. It is defined as:

\[ \boldsymbol{K}^{geo} = \int_V \boldsymbol{G}^T \boldsymbol{\tau} \boldsymbol{G}~\text{d} v \]

where: * \(\boldsymbol{G}\) is the matrix of shape function derivatives with respect to current coordinates (\(\partial \boldsymbol{N}/\partial \boldsymbol{x}\)). * \(\boldsymbol{\tau}\) is the matrix form of the current Cauchy stresses \(\boldsymbol{\sigma}^{t+\Delta t}\).

For a standard continuum element, \(\boldsymbol{\tau}\) is constructed such that \(\boldsymbol{G}^T \boldsymbol{\tau} \boldsymbol{G}\) represents the term \(\nabla \boldsymbol{N} \cdot \boldsymbol{\sigma} \cdot (\nabla \boldsymbol{N})^T\) in the weak form.

---

1. Thomas J. R. Hughes and James Winget. Finite rotation effects in numerical integration of rate constitutive equations arising in large-deformation analysis. International Journal for Numerical Methods in Engineering, 15(12):1862–1867, 1980-12. doi:10.1002/nme.1620151210. ↩