What numerical methods are used to model deformations in elastic physics?

Question

What numerical methods are used to model deformations in elastic physics? For example, here's an example of a hyperelastic deformation in Ansys:

Perhaps more simply than hyperelasticity, for linear elasticity, we have the equations:

$$ \nabla\cdot\sigma + {F} = \rho\ddot{{u}}\\ {\varepsilon} =\tfrac{1}{2} \left[{\nabla}{u}+({\nabla}{u})^\mathrm{T}\right]\\ {\sigma} = {C}:{\varepsilon} $$

where

• $\sigma$ - Stress tensor
• $\epsilon$ - Strain tensor
• $u$ - Displacement
• $C$ - Stiffness tensor

Say we apply a finite element on a mesh of the domain combined with something like a Runge-Kutta method to handle time, we can solve the above equations and find a solution $u$, which represents a deformation. However, a deformation seems to imply that something needs to move and up until this point we have a static, unmoving mesh of the domain. In this case, what moves? The mesh?

More generally, what's the general class of algorithms used to model motion and deformation of an elastic material similar to what the simulation above shows?

Answer 1

It seems that the type of algorithms differ considerably depending on whether the problem is:

• Quasistatic elastic or
• Hyperelastic

In the quasistatic elastic case, a simple approach is the following: As the first part of each timestep, the displacement field $u$ is computed. Since the displacement $u$ is now known at each node of the mesh, the nodes can be displaced according to $u$ as the second part of each timestep.

A detailed description of such a problem with code examples in C++ is given as one of the tutorial cases of the deal.ii finite element library: link

deal.ii also provides a concrete example for a hyperelastic problem. However, the description of the problem there is so involved that I do not wish to repeat it here in full detail. Just as a short summary, they use a three-field variational principle to derive a set of Euler-Lagrange equations for the displacement $u$ and the stress tensor $\sigma$.

As an in-depth reference they recommend:

A. Holzapfel (2001), Nonlinear Solid Mechanics. A Continuum Approach for Engineering, John Wiley & Sons. ISBN: 0-471-82304-X

Answer 2

If movement is removed there, deformation of the ring under variable load remains. This can be calculated in the model of an elastic body with the addition of Rayleigh damping. The animation shows the deformation of a rubber ring calculated using FEM and Mathematica 12.

Addendum Code (taken from tutorial and modified for this case) to compute animation

Needs["NDSolve`FEM`"]; \[CapitalOmega] = 
 ImplicitRegion[4^2 <= x^2 + y^2 <= 5^2, {x, y}];
mesh = ToElementMesh[\[CapitalOmega], {{-5, 5}, {-5, 5}}, 
   "MaxCellMeasure" -> 0.03];
mesh["Wireframe"]

diffusionCoefficients = 
 "DiffusionCoefficients" -> {{{{-(Y/(1 - \[Nu]^2)), 
       0}, {0, -((Y (1 - \[Nu]))/(2 (1 - \[Nu]^2)))}}, {{0, -((
        Y \[Nu])/(1 - \[Nu]^2))}, {-((Y (1 - \[Nu]))/(
        2 (1 - \[Nu]^2))), 
       0}}}, {{{0, -((Y (1 - \[Nu]))/(2 (1 - \[Nu]^2)))}, {-((
        Y \[Nu])/(1 - \[Nu]^2)), 
       0}}, {{-((Y (1 - \[Nu]))/(2 (1 - \[Nu]^2))), 
       0}, {0, -(Y/(1 - \[Nu]^2))}}}} /. {Y -> 10^2, \[Nu] -> 
    33/100}; massCoefficients = "MassCoefficients" -> {{1, 0}, {0, 1}};
vd = NDSolve`VariableData[{"Time", "DependentVariables", 
     "Space"} -> {t, {u, v}, {x, y}}];
sd = NDSolve`SolutionData[{"Time", "Space"} -> {0., 
     ToNumericalRegion[mesh]}];

methodData = InitializePDEMethodData[vd, sd];

loadCoefficients = "LoadCoefficients" -> {{0}, {0}};

"LoadCoefficients" -> {{0}, {0}};

Subscript[\[CapitalGamma], Nv] = 
 NeumannValue[-3 Sin[Pi t] Sign[y], -1 <= x <= 
   1]; Subscript[\[CapitalGamma], Du] = 
 DirichletCondition[u[x, y] == 0, 
  x == 0]; Subscript[\[CapitalGamma], Nu] = 
 NeumannValue[0, -1 <= x <= 1];
Subscript[\[CapitalGamma], Dv] = 
  DirichletCondition[v[x, y] == 0, y == 0];

initCoeffs = 
 InitializePDECoefficients[vd, 
  sd, {diffusionCoefficients, massCoefficients, loadCoefficients}];

initBCs = 
 InitializeBoundaryConditions[vd, 
  sd, {{Subscript[\[CapitalGamma], Du], Subscript[\[CapitalGamma], 
    Nu]}, {Subscript[\[CapitalGamma], Dv], Subscript[\[CapitalGamma], 
    Nv]}}];

sdpde = DiscretizePDE[initCoeffs, methodData, sd, "Stationary"];
sbcs = DiscretizeBoundaryConditions[initBCs, methodData, sd, 
  "Stationary"];

rhs[t_?NumericQ, uv_, duv_] := 
 Module[{l, s, d, m, tdpde, tbcs, rayleighDamping},
  NDSolve`SetSolutionDataComponent[sd, "Time", t];
  {l, s, d, m} = sdpde["SystemMatrices"];
  tdpde = DiscretizePDE[initCoeffs, methodData, sd, "Transient"];
  tbcs = DiscretizeBoundaryConditions[initBCs, methodData, sd, 
    "Transient"];
  {l, s, d, m} += tdpde["SystemMatrices"];
  rayleighDamping = 0.1*m + 0.04*s;
  DeployBoundaryConditions[{l, s, rayleighDamping, m}, tbcs];
  DeployBoundaryConditions[{l, s, rayleighDamping, m}, sbcs];
  l - s . uv - rayleighDamping . duv
  ]

dof = methodData["DegreesOfFreedom"];
init = dinit = ConstantArray[0, {dof, 1}];

mass = sdpde["MassMatrix"];
stiff = sdpde["StiffnessMatrix"];
rd = 0.1*mass + 0.04*stiff;

sparsity = 
  ArrayFlatten[{{mass["PatternArray"], 
     mass["PatternArray"]}, {rd["PatternArray"], rd["PatternArray"]}}];

Dynamic["time: " <> ToString[CForm[currentTime]]]
AbsoluteTiming[
 tfun = NDSolveValue[{
    mass . uv''[ t] == rhs[t, uv[t], uv'[t]]
    , uv[ 0] == init, uv'[ 0] == dinit}, uv, {t, 0, 2}
   , Method -> {"EquationSimplification" -> "Residual"}
   , Jacobian -> {Automatic, Sparse -> sparsity}
   , EvaluationMonitor :> (currentTime = t;)
   ]]

Visualization

split = Span @@@ 
  Transpose[{Most[# + 1], Rest[#]} &[methodData["IncidentOffsets"]]];

graphics = Function[t,
    dmesh = ElementMeshDeformation[mesh, Part[tfun[t], #] & /@ split];
    Show[{
      mesh["Wireframe"["MeshElement" -> "BoundaryElements"]],
      dmesh[
       "Wireframe"[
        "ElementMeshDirective" -> 
         Directive[EdgeForm[Red], FaceForm[]]]]
      }, PlotRange -> {{-6, 6}, {-5.5, 5.5}}]] /@ Range[0, 2, .1];

ListAnimate[graphics, SaveDefinitions -> True]

Answer 3

TL,DR: either the good old Galerkin finite element method, or mesh-free / particle methods.

There are a few things to unpack here. First, the simulation you show includes contact between an elastic body (the ring) and a hard boundary, so the constraints are non-holonomic. Contact problems are much more challenging than, say, elastic deformation under gravity, which only has holonomic constraints.

In the animation you showed above, the only thing that appears to change are the positions of the mesh nodes. This is a reasonable approach to take when the deformations from the reference configuration of the material are small, which is exactly when the elasticity equations are linear. For very large deformations, the transformed mesh in the body configuration could become tangled and as you might imagine this is bad news. There are a few ways around this.

First, you could change the internal workings of the simulation by periodic remeshing. When the triangles become too degenerate, you can flip edges to restore mesh quality. This approach can eliminate the possibility of mesh tangling but it's hard to code up right.

Second, you could change the physics problem that you're simulating in the first place. If the deformations are that large, using a linearized system of equations is no longer a good way to describe the problem. Nonlinear elasticity, for example see Neo-Hookean solids, includes terms that go to infinity when the Jacobian of the reference-to-body transformation becomes singular. This prevents weird deformations at a mathematical level but it also makes the numerical solution of the problem more challenging.

Both of these ideas are very much in the mindset of Galerkin finite element methods on simplicial meshes. A much larger departure from that way of thinking is to use mesh-free methods. This includes things like smoothed particle hydrodynamics, the material point method, etc. These are often more expensive computationally because you need to do neighbor queries on every timestep and this is hard to make cache-friendly. But they can adapt much better to very deformed or moving geometries.