Total Lagrangian SPH

A Total Lagrangian framework is used wherein the governing equations are formulated such that all relevant quantities and operators are measured with respect to the initial configuration (O’Connor & Rogers, 2021, Belytschko et al., 2000).

The governing equations with respect to the initial configuration are given by:

\[\frac{\mathrm{D}\bm{v}}{\mathrm{D}t} = \frac{1}{\rho_0} \nabla_0 \cdot \bm{P} + \bm{g},\]

where the zero subscript denotes a derivative with respect to the initial configuration and $\bm{P}$ is the first Piola-Kirchhoff (PK1) stress tensor.

The discretized version of this equation is given by O’Connor & Rogers (2021):

\[\frac{\mathrm{d}\bm{v}_a}{\mathrm{d}t} = \sum_b m_{0b} \left( \frac{\bm{P}_a \bm{L}_{0a}}{\rho_{0a}^2} + \frac{\bm{P}_b \bm{L}_{0b}}{\rho_{0b}^2} \right) \nabla_{0a} W(\bm{X}_{ab}) + \frac{\bm{f}_a^{PF}}{m_{0a}} + \bm{g},\]

with the correction matrix (see also GradientCorrection)

\[\bm{L}_{0a} := \left( -\sum_{b} \frac{m_{0b}}{\rho_{0b}} \nabla_{0a} W(\bm{X}_{ab}) \bm{X}_{ab}^T \right)^{-1} \in \R^{d \times d}.\]

The subscripts $a$ and $b$ denote quantities of particle $a$ and $b$, respectively. The zero subscript on quantities denotes that the quantity is to be measured in the initial configuration. The difference in the initial coordinates is denoted by $\bm{X}_{ab} = \bm{X}_a - \bm{X}_b$, the difference in the current coordinates is denoted by $\bm{x}_{ab} = \bm{x}_a - \bm{x}_b$.

For the computation of the PK1 stress tensor, the deformation gradient $\bm{F}$ is computed per particle as

\[\bm{F}_a = \sum_b \frac{m_{0b}}{\rho_{0b}} \bm{x}_{ba} (\bm{L}_{0a}\nabla_{0a} W(\bm{X}_{ab}))^T \\ \qquad = -\left(\sum_b \frac{m_{0b}}{\rho_{0b}} \bm{x}_{ab} (\nabla_{0a} W(\bm{X}_{ab}))^T \right) \bm{L}_{0a}^T\]

with $1 \leq i,j \leq d$. From the deformation gradient, the Green-Lagrange strain

\[\bm{E} = \frac{1}{2}(\bm{F}^T\bm{F} - \bm{I})\]

and the second Piola-Kirchhoff stress tensor

\[\bm{S} = \lambda \operatorname{tr}(\bm{E}) \bm{I} + 2\mu \bm{E}\]

are computed to obtain the PK1 stress tensor as

\[\bm{P} = \bm{F}\bm{S}.\]

Here,

\[\mu = \frac{E}{2(1 + \nu)}\]

and

\[\lambda = \frac{E\nu}{(1 + \nu)(1 - 2\nu)}\]

are the Lamé coefficients, where $E$ is the Young's modulus and $\nu$ is the Poisson ratio.

The term $\bm{f}_a^{PF}$ is an optional penalty force. See e.g. PenaltyForceGanzenmueller.

TotalLagrangianSPHSystem(initial_condition,
                         smoothing_kernel, smoothing_length,
                         young_modulus, poisson_ratio;
                         n_fixed_particles=0, boundary_model=nothing,
                         acceleration=ntuple(_ -> 0.0, NDIMS),
                         penalty_force=nothing, source_terms=nothing)

solid = union(beam, fixed_particles)

Penalty Force

In FEM, underintegrated elements can deform without an associated increase of energy. This is caused by the stiffness matrix having zero eigenvalues (so-called hourglass modes). The name "hourglass modes" comes from the fact that elements can deform into an hourglass shape.

Similar effects can occur in SPH as well. Particles can change positions without changing the SPH approximation of the deformation gradient $\bm{F}$, thus, without causing an increase of energy. To ensure regular particle positions, we can apply similar correction forces as are used in FEM.

Ganzenmüller (2015) introduced a so-called hourglass correction force or penalty force $f^{PF}$, which is given by

\[\bm{f}_a^{PF} = \frac{1}{2} \alpha \sum_b \frac{m_{0a} m_{0b} W_{0ab}}{\rho_{0a}\rho_{0b} |\bm{X}_{ab}|^2} \left( E \delta_{ab}^a + E \delta_{ba}^b \right) \frac{\bm{x}_{ab}}{|\bm{x}_{ab}|}\]

The subscripts $a$ and $b$ denote quantities of particle $a$ and $b$, respectively. The zero subscript on quantities denotes that the quantity is to be measured in the initial configuration. The difference in the initial coordinates is denoted by $\bm{X}_{ab} = \bm{X}_a - \bm{X}_b$, the difference in the current coordinates is denoted by $\bm{x}_{ab} = \bm{x}_a - \bm{x}_b$. Note that Ganzenmüller (2015) has a flipped sign here because they define $\bm{x}_{ab}$ the other way around.

This correction force is based on the potential energy density of a Hookean material. Thus, $E$ is the Young's modulus and $\alpha$ is a dimensionless coefficient that controls the amplitude of hourglass correction. The separation vector $\delta_{ab}^a$ indicates the change of distance which the particle separation should attain in order to minimize the error and is given by

\[ \delta_{ab}^a = \frac{\bm{\epsilon}_{ab}^a \cdot \bm{x_{ab}}}{|\bm{x}_{ab}|},\]

where the error vector is defined as

\[ \bm{\epsilon}_{ab}^a = \bm{F}_a \bm{X}_{ab} - \bm{x}_{ab}.\]

PenaltyForceGanzenmueller(; alpha=0.1)