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

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?

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

• Thanks for the deal.ii reference as that was very helpful. As for the book, the Holzapfel reference is helpful in discussing the formulations, but it doesn't appear to have much information about the mesh management or nuances such as updating the stress at each step. Do you know of another book that provides additional information about these details? Alternatively, do you know what these kind of algorithms are called, so that I can search further? – wyer33 Mar 11 '20 at 21:47
• @wyer33 Thank you. Maybe some of the other references in the step-44 tutorial could be more helpful. This seems promising: C. Miehe (1994), Aspects of the formulation and finite element implementation of large strain isotropic elasticity International Journal for Numerical Methods in Engineering 37 , 12, 1981-2004. DOI: 10.1002/nme.1620371202; – 123 Mar 12 '20 at 19:39
• I also realize that my answer is very focused on the Finite Element method and as pointed out by @DanielShapiro mesh-free/particle methods are also an option for larger deformations. Computation of elastic deformations is a very broad topic. – 123 Mar 12 '20 at 19:45

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["NDSolveFEM"]; \[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 = NDSolveVariableData[{"Time", "DependentVariables",
"Space"} -> {t, {u, v}, {x, y}}];
sd = NDSolveSolutionData[{"Time", "Space"} -> {0.,
ToNumericalRegion[mesh]}];

methodData = InitializePDEMethodData[vd, sd];

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,

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},
NDSolveSetSolutionDataComponent[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]
`
• Could you add the Mathematica code that you used to compute this? – user21 May 19 at 13:03
• @user21 Thank you for your question. See Addendum with code you required. – Alex Trounev May 19 at 15:09
• Ah, you made use of Rayleigh Damping, I see thanks! – user21 May 19 at 15:17
• @user21 Yes it is not contact problem like we discussed on this page mathematica.stackexchange.com/questions/124324/… – Alex Trounev May 19 at 15:22

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.

• Thanks for the response. Quick clarification. Say we use a Galerkin finite element method to solve the elastic equations. Coming out this, we have the displacements $u$, the strain $\epsilon$ and the stress $\sigma$. Next, we could move the mesh according to $u$ modulo the difficulties you mentioned. Once this is accomplished, don't we still have to map the stress $\sigma$ to the new mesh in order to capture the accumulated forces? In addition, do we need to monitor the displacements and stop the time integrator as soon as a new collisions occurs? Is there a general name for this process? – wyer33 Mar 11 '20 at 3:13
• That depends. Some formulations of the problem will solve only for the displacement and calculate the strain / stress after the fact, while others will use a mixed formulation for both stress and displacement. In the first case you'd just recalculate the stress on the new mesh, while in the second you'd be solving for it directly. For the contact problem, I imagine you'd need some kind of explicit collision detection (which is the term to search for) but I haven't tried that myself. – Daniel Shapero Mar 12 '20 at 16:19