6.7. Nonlinear solid mechanics

The possibility of having a moving mesh, i.e. having the mesh coordinates as degrees of freedom, is obviously a good foundation to account for deformable solids. In this section, we cover the possibilities of pyoomph for nonlinear elasticity of solid bodies. The essential idea in pyoomph’s implementation is that the Lagrangian coordinates represent the undeformed body, whereas the Eulerian coordinates reflect the current deformed configuration.

Note

The implementation of nonlinear elasticity is essentially a one-to-one copy of the implementation in oomph-lib. Oomph-lib’s documentation covers the individual aspects in considerably more detail, and with the exact co- and contravariant definitions, than this documentation here.

As governing equations, we have the principle of virtual displacements, which reads in pyoomph’s notation

\[\left( \rho \partial_t^2 \vec{x} - \vec{f}, \Gamma \vec{\chi} \right)_\xi + \left( \left(\nabla_\xi \vec{x}\right) \boldsymbol{\sigma},\Gamma \nabla_\xi\vec{\chi}\right)_\xi-\langle \vec{T},\vec{\chi}\rangle=0\]

Again, \(\vec{x}(\vec{\xi},t)\) represents the deformed configuration (with test function \(\vec{\chi}\)) and \(\vec{\xi}\) is the undeformed configuration. \(\vec{f}\) is a bulk force density (in the undeformed configuration), \(\rho\) is the mass density in undeformed configuration and \(\Gamma\) is an isotropic growth factor, i.e. the factor the undeformed solid wants to grow in terms of the undeformed length/area/volume (depending on the dimension). The bulk part of the equation is available in the class DeformableSolidEquations, whereas a traction \(\vec{T}\) can be applied on the boundary of the solid by the class SolidTraction. For pressures, tractions in negative normal direction, the class SolidNormalTraction can be used.

The constitutive law is entirely wrapped in the second Piola-Kirchhoff stress tensor \(\boldsymbol{\sigma}\). Following oomph-lib’s GeneralisedHookean class, a reasonable definition is given by the following

\[\begin{split}\begin{aligned} \sigma_{ij}&=E_{ijkl}\gamma_{kl}\\ E_{ijkl}&=\frac{E}{1+\nu}\left(\frac{\nu}{1-2\nu}G^{-1}_{ij}G^{-1}_{kl}+\frac{1}{2}\left(G^{-1}_{ik}G^{-1}_{jl}+G^{-1}_{il}G^{-1}_{jk}\right)\right)\\ \boldsymbol{G}&=\left(\left(\nabla_\xi \vec{x}\right)^\mathrm{t}\left(\nabla_\xi \vec{x}\right)\right)\\ \boldsymbol{\gamma}&=\frac{1}{2}\left(\boldsymbol{G}-\boldsymbol{g}\right)\\ \boldsymbol{g}&=\Gamma^{2/D}\left(\left(\nabla_\xi \vec{\xi}\right)^\mathrm{t}\left(\nabla_\xi \vec{\xi}\right)\right) \end{aligned}\end{split}\]

Please kindly excuse the sloppy notation not respecting the co- and contravariances here correctly. This is done considerably better in oomph-lib’s documentation. However, the way it is written here corresponds to the implementation in pyoomph, where we heavily use matrices and the vector gradients, since they automatically adjust correctly based on the considered coordinate system (e.g. \(\boldsymbol{g}\) is the (scaled) identity matrix in Cartesian coordinates, but in axisymmetric coordinates, it is not). Entering parameters are Young’s modulus \(E\) and the Poisson ratio \(\nu\). For small deformations, it perfectly recovers Hook’s law. \(D\) is the dimension of the solid, which is e.g. 3 for a 2d mesh with axisymmetric coordinate system. This constitutive law is implemented in the class GeneralizedHookeanSolidConstitutiveLaw.

If \(\nu=1/2\), the solid is incompressible and the definition becomes singular. In that case, we have to use

\[\begin{align} \sigma^{ij}&=-pG^{-1}_{ij}+\frac{E}{3}\left(G^{-1}_{ik}G^{-1}_{jl}+G^{-1}_{il}G^{-1}_{jk}\right)\gamma_{kl} \end{align}\]

instead. Here, \(p\) is a pressure field, which enforces that

\[\begin{align} \det \boldsymbol{G}-\det \boldsymbol{g}=0 \end{align}\]

i.e. that the deformation locally conserves the volume. For this law, we have to use the class IncompressibleHookeanSolidConstitutiveLaw.

Note

For compressible constitutive laws, the mass density \(\rho\) is to be understood as mass density in the undeformed configuration. The mass density in the deformed configuration can be accessed by var("deformed_mass_density"), which is calculated by the conservation of mass. In the incompressible case, there is no difference and var("deformed_mass_density") will expand to \(\rho\).

Warning

The nonlinear elasticiy equations cannot be used in combination with azimuthal or Cartesian normal mode stability analysis as outline in section Section 10.3 and Section 10.4. For this approach, we would have to expand the mesh coordinates also by an azimuthal or additional Cartesian mode, which is not implemented (yet).