# A 3D cell-centered ADER MOOD Finite Volume method for solving updated Lagrangian hyperelasticity on unstructured grids

Abstract: In this paper, we present a conservative cell-centered Lagrangian Finite Volume scheme for solving the hyperelasticity equations on unstructured multidimensional grids. The starting point of the present approach is the cell-centered FV discretization named EUCCLHYD and introduced in the context of Lagrangian hydrodynamics. Here, it is combined with the a posteriori Multidimensional Optimal Order Detection (MOOD) limiting strategy to ensure robustness and stability at shock waves with piecewise linear spatial reconstruction. The ADER (Arbitrary high order schemes using DERivatives) approach is adopted to obtain second-order of accuracy in time. This strategy has been successfully tested in a hydrodynamics context and the present work aims at extending it to the case of hyperelasticity. Here, the hyperelasticity equations are written in the updated Lagrangian framework and the dedicated Lagrangian numerical scheme is derived in terms of nodal solver, Geometrical Conservation Law (GCL) compliance, subcell forces and compatible discretization. The Lagrangian numerical method is implemented in 3D under MPI parallelization framework allowing to handle genuinely large meshes. A relatively large set of numerical test cases is presented to assess the ability of the method to achieve effective second order of accuracy on smooth flows, maintaining an essentially non-oscillatory behavior and general robustness across discontinuities and ensuring at least physical admissibility of the solution where appropriate. Pure elastic neo-Hookean and non-linear materials are considered for our benchmark test problems in 2D and 3D. These test cases feature material bending, impact, compression, non-linear deformation and further bouncing/detaching motions.

## Summary (9 min read)

### 1 Introduction

- This work is concerned with the accurate multi-dimensional simulation of hyper-elasticity models by updated Lagrangian Finite Volume (FV) numerical scheme.
- Their FV discretization relies on Godunov-type methods.
- Further, the total Lagrangian FV discretization of the physical conservation laws relies on the discretization of the time rate of change of the deformation gradient, its determinant and its co-factor.
- The development of cell-centered FV Lagrangian schemes based on nodal solvers [23,42,22] for gas dynamics has been extended to 2D and 3D updated Lagrangian hyperelasticity [39,14].
- Numerical results involving materials enduring large deformation (bending, twisting, etc.) adopted from [40,36,14] will be presented to assess the ability of this updated Lagrangian numerical scheme to simulate such hyperelastic situations.

### 2 Updated Lagrangian hyperelastic modeling for isotropic materials

- The authors aim at recalling the conservation laws describing the time evolution of isotropic solid materials undergoing large deformations.
- The resulting conservation laws of mass, momentum and total energy shall be written under the updated Lagrangian form.
- The isotropic materials under consideration are characterized by an hyperelastic constitutive law.
- The material frame indifference principle and the thermodynamic consistency are automatically satisfied.
- The interested reader might refer for instance to [35] for further developments on these subtle topics.

### 2.1.1 Lagrange-Euler mapping.

- This smooth function is the Lagrange-Euler mapping which relates the Lagrangian point X to its Eulerian counterpart x.
- By definition, this mapping satisfies Φ(X, 0) = X and its Jacobian, also named the deformation gradient, reads F(X, t) = ∇XΦ(X, t), where the symbol ∇X denotes the gradient operator with respect to the Lagrangian coordinate.
- More precisely, let G(X, t) denotes the Lagrangian representation of a physical quantity.
- In what follows, the same notation shall be employed for both descriptions.

### 2.1.2 Measures of deformation

- Let us consider the infinitesimal material fiber dX in the initial configuration which maps into dx = FdX through the motion.
- The right and the left Cauchy-Green tensors are symmetric positive definite and share the same eigenvalues, refer to [35].
- These tensors are relevant measures of deformation since for a rigid rotation they boil down to the identity tensor.

### 2.1.3 Geometric conservation law (GCL)

- Time differentiating the deformation gradient, F = ∇XΦ, recalling that the partial time derivative of the mapping is the kinematic velocity v = ∂Φ ∂t , leads to the Geometric Conser- vation Law (GCL) written under total Lagrangian form as ∂F ∂t −∇Xv = 0, (5) where F(X, 0) = Id.
- The GCL is supplemented with the compatibility constraint∇X×F = 0, which ensures that the solution of the foregoing partial differential equation corresponds to the gradient of a mapping.
- The satisfaction of this compatibility constraint at the discrete level is the cornerstone which any proper discretization of the conservation laws written in total Lagrangian form should rely on, refer to [28].
- The notation L = ∇xv represents the velocity gradient tensor with respect to the current configuration.
- Finally, substitution of the GCL (7) into the foregoing equation yields the partial differential equation satisfied by the Jacobian of the deformation gradient dJ dt − Jtr(L) =.

### 2.2 Governing equations

- This section aims at briefly recalling the conservation laws and the constitutive law governing the time evolution of an isotropic material undergoing large deformations.
- The interested reader might consult [35] for further details regarding their derivation.

### 2.2.1 Conservation laws

- The Cauchy stress tensor, T, is symmetric, i.e., T = Tt, which ensures the conservation of angular momentum.
- Let us note that the nabla operator employed in the foregoing equations is expressed in terms of the Eulerian coordinate x.
- This system of conservation laws written under updated Lagrangian representation is supplemented by the trajectory equation already introduced in (1), which is rewritten under the form dx dt = v(x(t), t), x(0) = X. (11).
- To close the system of conservation laws, it remains to provide a constitutive law for expressing the Cauchy stress tensor in terms of the deformation and a thermodynamic variable.
- This will be achieved in the next paragraph introducing the free energy Ψ.

### 2.2.2 Constitutive law for isotropic materials

- The constitutive law is derived invoking the frame indifference principle and the compatibility with thermodynamics.
- Here, the material under consideration is characterized by the free energy expressed in terms of the left Cauchy-Green tensor and the absolute temperature Ψ ≡ Ψ(B, θ).
- Let us point out that the Cauchy stress tensor and the left Cauchy-Green tensor commute, i.e. TB = BT.
- This shows that the system of conservation laws (10) is equipped with a supplementary conservation law which states that entropy is conserved along flow trajectories.
- Thus, constitutive law (14) for isotropic materials is consistent with the second law of thermodynamics.

### 2.2.3 Volumetric shear strain decomposition

- The authors want to study materials that can sustain only limited shear strain but respond elastically to large change in volume.
- To construct this additive decomposition, the authors start by introducing the multiplicative decomposition of the deformation gradient tensor, F, into a volumetric and an isochoric parts.
- This framework provides a constitutive law fulfilling the material frame indifference principle; the thermodynamic consistency with the second law.
- Namely, the time rate of change of the deviatoric stress is expressed in terms of the deviatoric part of the strain rate tensor.

### 2.2.4 Examples of constitutive laws

- Let us point out that the volumetric/shear decomposition allows us to define separately the pressure by introducing an hydrodynamic equation of state characterized by the volumetric free energy Ψv = Ψv(J, θ).
- More generally, one can utilize one’s favorite equation of state regardless of the shearing free energy choice.
- One shall always choose at least a convex equation of state to ensure the hyperbolicity of the hydrodynamic part of the system of conservation laws.
- For the numerical applications, the authors shall consider the particular case a = −1 which corresponds to neo-Hookean materials.
- Finally, material mechanical properties are often described in terms of Young modulus E, Poisson ration ν and shear modulus µ, which also corresponds to the second Lamé coefficient.

### 2.3 Updated Lagrangian hyperelasticity for isotropic materials

- The authors summarize the set of partial differential equations governing the time evolution of the isotropic hyperelastic material under consideration.
- Here, B = J− 23B denotes the isochoric part of the left Cauchy-Green tensor and I1, I2 are respectively its first and second invariants.
- It is remarkable to note that updated Lagrangian isotropic hyperelasticity requires only the knowledge of the left Cauchy-Green tensor.
- The physical admissibility property is defined by a set of so-called admissible states such that the material vector determines a valid state according to the conservation and constitutive laws.

### 3.1.1 Geometrical entities

- The cell center, xc, corresponds to its centroid, that is xc = 1 |ωc(t)| ∫ ωc(t) x dv. 0 is the initial density distribution, and dV refers to the integral measure over volume.
- (36) This equation allows us to distinguish the fundamental geometrical object apcnpc = ∂|ωc| ∂xp , (37) known as the corner normal vector.
- It has been introduced by several authors in the context of Lagrangian hydrodynamics, refer for instance to [2].
- Here, apc represents a (d−1)-measure (length in 2D, area in 3D) and npc is a unit outward pointing vector.

### 3.1.2 Conservative and constitutive discrete variables

- The time dependent conserved or constitutive variables are the cell-centered approximate mass-averaged values gathered into vector Qc(t) = (τc(t),vc(t), ec(t),Bc(t)).
- From now on the authors implicitly assume the dependence on time and to lighten the notation they omit it.

### 3.2 Discrete velocity gradient operators

- This paragraph aims at presenting a FV approximation of the velocity gradient for the updated Lagrangian representation.
- The interested reader might refer for instance to [22] for finding the demonstration of this result.
- The authors end up presenting an important by-product of the discrete Eulerian velocity gradient operator.
- This is nothing but the discrete cell-centered divergence operator.

### 3.3.1 Conservation laws - GCL, momentum and total energy

- The geometrical conservation law (GCL) is a fundamental consistency property in Lagrangian framework.
- The space discretization of the momentum and the total energy equations is performed introducing the subcell force fpc, which is the traction force exerted on the outer boundary of the subcell ωpc.
- The discretization of this equation, which relies on geometrical considerations, shall be detailed in section 3.4.2.
- The thermodynamic consistency of the semi-discrete scheme; .
- The first requirement allows us to express the subcell force fpc in terms of the node velocity vp.

### 3.3.2 The subcell force expression ensuring thermodynamic consistency

- The constitutive law leads to the definition of the following discrete Cauchy stress tensor: Tc = 2ρc ∂Ψ ∂B (Bc, θc)Bc. Starting from the Gibbs identity (16) let us compute the time evolu- tion of the entropy within cell ωc mcθc dηc dt = −1 2 |ωc|TcB−1c :.
- Therefore, in order to ensure a proper entropy dissipation, the authors propose to design the subcell force as fpc = apcTcnpc + Mpc(vp − vc), (54) where the subcell matrix Mpc is symmetric positive definite.
- The authors easily verify that mcθc dηc dt = ∑ p∈P(c) Mpc(vp − vc) · (vp − vc) ≥ 0, (55) which expresses the consistency with the second law of thermodynamics at the semi-discrete level.
- Now, it remains to determine the subcell matrix Mpc, which genuinely characterizes the numerical scheme.
- Several possibilities have already been explored by different authors in [44,18,52].

### 3.3.3 The nodal solver ensuring the conservation of total energy and momentum

- Ignoring the boundary conditions, the conservation of total energy and momentum at the semi-discrete level is ensured provided that the sum of the subcell forces impinging at node xp is equal to zero, that is ∑ c∈C(p) fpc = 0. (56).
- The interested reader might refer for instance to [46,45] for the justification of this result.
- Substituting the expression of the subcell force (54) in the foregoing condition leads to the linear system Mpvp = ∑ c∈C(p) Mpcvc − ∑ c∈C(p) apcTcnpc, where Mp = ∑ c∈C(p) Mpc. (57) Notice that Mp is symmetric positive definite and, thus, invertible.
- Once the node velocity is determined thanks to (57) then the trajectory equation can be invoked to update the node position.

### 3.4.1 Physical conservation laws

- The first-order time discretization simply considers t∗ = tn and the cell-centered values of the state vector Qnc = (τc,vc, ec,Bc)n.
- This reconstruction is not limited, thus it can be efficiently obtained following the procedure detailed in [13].
- The limiting will be done a posteriori via the MOOD paradigm.
- Once the spatial reconstruction polynomial qh(ξ, 0) is known, the ADER methodology performs a local time evolution of the governing equations (10)-(11)- (32).

### 3.4.2 Space-time discretization of the left Cauchy-Green tensor equation

- The closure of the foregoing discrete system of conservation laws relies on the generic constitutive law (15) which requires the knowledge of the cell-centered left Cauchy-Green tensor Bc within the cell ωc.
- Subtracting the gradient of this velocity field from the time derivative of (73) leads to ∂Fh ∂t (X, t)− ∑ p∈P(c) vp(t)⊗∇Xλp = 0. (74) This shows that the foregoing Finite Element mapping provides a consistent spatial discretization of the kinematics which is fully compatible with the Geometrical Conservation Law (5).
- Let us point out that the foregoing space discretization of the deformation gradient are expressed in terms of the Lagrangian grid.
- Let us emphasize the crucial role played by (78).
- This suggests the alternative time discretization of the Jacobian equation Jn+1c = ( 1 + ∆tI1 + ∆t2 2 I21 ) Jnc . (81) This is a positivity preserving time discretization since the quadratic polynomial between parenthesis is always strictly positive regardless ∆t.

### 3.5 Limiting: a posteriori MOOD loop

- While in the original ADER schemes the limiting relies on a priori limited WENO reconstructions for all variables [26,25], here the authors adopt an a posteriori MOOD paradigm, see [21,13].
- In the present implementation, the MOOD loop simply embraces the main evolution routines of the ADER method and iterates to recompute those cells with invalid values, detected by the admissibility criteria.
- In the worst case scenario all cells in the domain are updated with the parachute scheme, leading to the true first-order accurate and robust numerical solution.
- In any other case, the MOOD loop always converges and produces an acceptable numerical solution, assuming that the parachute scheme does so.
- Any cell which does not belong to Ah or does not fulfill (83) is declared troubled, and, sent back to tn along with its neighbors for their re-computation using the next scheme in the cascade, see [13].

### 3.6 Internal consistency consideration

- The first of such equality is the link between the constant cell mass, the specific volume τ (or the density ρ = 1/τ) and volume |ω|: m = |ω| τ .
- The primary variable is the specific volume: τn+1c computed by (64).
- Therefore, from a consistency point of view, the authors should monitor that the following equality holds ω = ||ωn+1c | − τn+1c mc|. (84) This constraint is nothing but the discrete version of the GCL which is fulfilled by construction for modern cell-centered Lagrangian schemes thanks to the exact integration of the geometry in (67).
- The authors therefore monitor their difference as an internal consistency criteria as B = ∣∣∣√detBn+1c − τ 0cτn+1c ∣∣∣. (85) Here B must be small compared to the accuracy of the scheme, and, in fact, the current scheme ensures that B = O(∆t3) due to the use of Simpson’s rule in (67), see section 5.9 for a numerical validation.
- The difference between these two “measures” was monitored to assess the internal consistency of the scheme.

### 3.7 Time-step monitoring

- Notice that the a posteriori detection allows to ensure the positivity of the cell volume and the internal energy provided that the parachute first-order scheme does.
- As such the time-step control must be suited for the parachute scheme.
- Notice that the a posteriori MOOD loop may also be used to try to exceed the time-step restrictions (86) at the price of creating more troubled cells, for instance by setting CCFL closer to one.

### 3.8 Boundary condition treatments

- The Boundary Conditions (BCs) play a crucial role in the time evolution of the numerical solution.
- In the context of an hyperelasticity model solved by the Lagrangian numerical scheme the authors consider several types of BCs, such as free traction, restricted normal/tangential displacement and contact/symmetry plans.
- To enlarge even further the ability of the code to handle complex situations, the authors have added the possibility for a BCs to change its type during the simulation, for instance transitioning from free-traction to null normal velocity.
- Later, if the medium happens to detach from the target, then the distance function becomes again greater than the threshold value and the original BC is restored, that is BCB → BCA.
- Then they must be applied in a hierarchical manner, taking into account the most important one first, possibly relaxing the fulfillment of the other ones.

### 4.1 Algorithm

- In this section the authors recall the main steps of the MOOD loop applied to this cell-centered Lagrangian scheme sketched in figure 3.
- Then a nodal solver and the ADER methodology allow to compute a candidate solution at tn+1, leading to an P1 (unlimited 2nd order accurate) candidate solution.
- This is the purpose of the ’Detection’ box to determine which cells are troubled, and, on the contrary to accept the admissible ones.
- Those troubled cells and their Voronoi neighbors are solely sent back into the Lagrangian FV solver for re-computation.
- This part of the solution which has been recomputed is re-tested against the detection criteria.

### 4.2 Meshing and parallelization

- The 3D Lagrangian simulation algorithm is fully coded in FORTRAN and relies on MPI protocol for the parallelization and the free graph partitioning software METIS [38].
- This primary mesh is then partitioned among the total number of threads NCPU, see figure 4- right for NCPU = 4 in 2D and the coarse mesh in black.
- A local isotropic recursive refinement is further applied.
- Notice that their aim is not to produce a perfectly scaled and massively parallel code, but rather being able to simulate with large enough 3D mesh reasonably fast.
- To split one single tetrahedron the authors insert new vertices at the midpoints of each edge and connect the vertices together and form four new sub-tetrahedra associated to the vertices.

### 5 2D and 3D test problems

- The authors refer to the current Lagrangian ADER-MOOD method with the acronym ’LAM’.
- The first cascade was employed for the hydrodynamics equations in [13], the second one is tested in this work.
- The time-dependent computational domain is addressed with ω(t), while Q(x, t = 0) ≡ Q0(x) = (1/ρ0,v0, p0,B0) denotes the vector of initial conservative variables typically used to setup the test problems.
- The unstructured tetrahedral meshes are obtained by meshing software, such as GMSH [32] which takes a characteristics target length h as input parameter.

### 5.1 Swinging plate

- The 2D swinging plate test problem, see [40,53], is employed to evaluate the numerical order of convergence.
- The velocity and displacement fields are divergence-free, leading to the exact pressure pex =.
- Space-time dependent boundary conditions are prescribed for the normal velocities, according to the exact solution (89).
- Notice that the exact solution is a smooth one and the final time is set to tfinal = π/ω, and the final displacement corresponds to the initial one.
- The unstructured mesh made of triangles is successively refined and the final characteristics length Lc(ω(tfinal)) is measured and further used to compute the numerical order of convergence O as reported in table 2, where one can notice that the numerical scheme is able to retrieve the second-order of convergence on this regular solution.

### 5.2 Elastic vibration of a Beryllium plate

- This test case describes the elastic vibration of a beryllium plate or bar, see [49,14] for instance.
- Free boundary conditions are applied on the plate faces.
- In figure 6 the authors present the numerical results obtained at different output times for the pressure (left panels) and cell orders (right panels).
- Yellow cells are dealt with the unlimited second order scheme (maximal order, prone to oscillation), while the blue ones employ a piecewise reconstruction limited by BJ slope limiter, via the a posteriori MOOD loop.
- As expected the nominally second order scheme is able to follow the barycenter with lower dissipation.

### 5.3 Finite deformation of a cantilever thick beam

- Free boundary conditions are considered apart from the fixed-wall bottom of the bar.
- The results are qualitatively in agreement with the published ones from the literature.
- Moreover the authors observe on the right panels that the yellow cells (unlimited second-order scheme) are massively represented, while only few of them demand dissipation (blue cells).
- For comparison purposes the authors also superimpose in black line the shapes obtained with the simpler cascade P1 → P0 from [13].
- At last the right panel presents the percentage of troubled cells encountered as a function of time.

### 5.5 Twisting column

- A twisting column test case aims at examining the effectiveness of the proposed methodology in highly nonlinear scenarios, see [36] and the reference therein.
- In figure 13 the authors plot the shape of the column colored by the pressure distribution for different output times.
- The main behaviors are reproduced by the numerical simulation.
- Notice that there is no spurious oscillations nor suspicious pressure distribution.
- For a numerical simulation recall that the twisting period does not only depend on the material but also on the numerical dissipation of the scheme.

### 5.6 Rebound of a hollow circular bar

- The bar impacts against a rigid friction-less wall with an initial velocity of v0 = (0, 0,−100)t m/s and the separation distance between the bar and wall is 4 mm.
- Upon impact, the bar undergoes large compressive deformation until t = 150 µs when all the kinetic energy of the bar is converted into internal strain energy.
- At last, on the right panel of figure 16, the authors show the percentage of bad cells detected by the a posteriori limiter and observe that, on average, less than 3% demands limiting at each iteration.

### 5.7 Impact of a jelly-like droplet

- For this test case the authors consider the impact of a jelly-like material onto a flat rigid horizontal surface, inspired by the test in [51].
- They are put in respect to each other for comparison purposes.
- Regardless of the constitutive model, i.e. the value of a, the jelly-like material is compressed after the impact and deforms back and forth due to its elastic behavior.
- The neo-Hookean model produces faster and more pronounced elastic behaviors compared to the non-linear model which retrieves a ratio closer to one faster.
- The experimental results in [41] provide approximate values 2.25 and 2.75, respectively, while their simulations produce 1.8 and 2.5 in accordance to the numerical results in [51].

### 5.8 L-shaped block

- To assess the linear and angular momentum preservation, a specific test case has been designed in [8,33].
- It consists in simulating the motion of a 3D L-shaped block subjected to initial time-dependent impulse traction forces, F1(t) and F2(t), at two of its sides, creating de facto an external torque, see figure 12.
- The results are in agreements with the ones presented for instance in [8,33].
- In the figure the authors clearly see that after t = 5, when the boundary forces vanish, the angular momenta are nicely preserved even if nothing in the method has been designed to ensure so.
- Notice that the linear momentum is preserved up to machine precision and the authors omit these results.

### 5.9 Monitoring of the internal consistency

- Moreover for the two 2D problems the authors compare their results to the ones given by a genuine first order Lagrangian scheme.
- Table 3 reports the maximum error in time through the entire simulations for both 2D and 3D tests, demonstrating that it remains at machine precision and thus that the internal consistency of the update of B is satisfied by their second order scheme.

### 6 Conclusions and perspectives

- This paper considers the second-order accurate cell-centered Lagrangian scheme originally designed for the hydrodynamics system of conservation laws [13], and, extends it to solve the hyperelasticity model for materials in 2D and 3D.
- The authors have focused the first part of the paper on presenting the hyperelasticity model and its consistency in the Lagrangian frame.
- Robustness and stability are gained by the use of an a posteriori MOOD limiting strategy, that is a second-order unlimited candidate solution at tn+1 is tested against appropriate detection criteria to determine troubled cells.
- Moreover evolving boundary conditions have been implemented to allow for impacting and detaching of materials onto walls.

Did you find this useful? Give us your feedback

...read more

##### References

4,113 citations

3,818 citations

2,043 citations

1,619 citations

1,109 citations

##### Related Papers (5)

##### Frequently Asked Questions (2)

###### Q2. What future works have the authors mentioned in the paper "A 3d cell-centered ader mood finite volume method for solving updated lagrangian hyperelasticity on unstructured grids" ?

A plan for future work involves the introduction of plasticity into this hyperelasticity model. The authors also plan to investigate the high-order extension over curvilinear simplicial grids of the present FV discretization. Another direction of evolution would be to add some ArbitraryLagrangian-Eulerian capability and the possibility to let two elastic materials interacting with each other, for instance impacting, deforming and further detaching from each others.