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

TL;DR: In this article, a cell-centered Lagrangian finite volume (FV) discretization 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.

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

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.

Figures

Fig. 1. Sketch of the Lagrangian-Eulerian mapping relating a material Lagrangian point X at t = 0 and a spatial Eulerian one x at t > 0 through Φ, and its Jacobian F(X, t) = ∂Φ

Fig. 14. Twisting column — Time evolution of non-dimensionalised height of the column measured at initial point xT = (0, 0, 6) (left) — Numerical dissipation (88) of second order scheme (center) — Percentage of bad cells detected at each time step (right).

Fig. 9. Cantilever thick beam — Comparison between the MOOD cascades: P1 → P0 (LAM P0-lim) and P1 → Plim1 → P0 (LAM P1-lim) for the computed numerical dissipation (88) as a function of time (left) and percentage of bad cells detected at each time step (right).

Fig. 19. L-shaped block test problem — From top-left to bottom-right: time t = 2.5, 5.0, 7.5, 10.0, 12.5 and 15 s.

Table 3 Maximal value of B for the first 75000 time-steps of the 2D and 3D tests cases in this paper.

Fig. 20. L-shaped block test problem — Angular momentum Ac in x (green), y (red) and z (blue) directions as a function of time.

Fig. 10. Blake’s problem — Computational mesh of the needle domain ω = [r, θ, φ] = [0.9, π/180, π/180] with h = 1/1000 (left) and pressure distribution at the final time tfinal = 1.6 · 10−4 s (right).

Table 2 Numerical errors in L2 norm and convergence rates for the 2D swinging plate test computed with second order of accuracy Lagrange ADER scheme at time tfinal = π/ω. The error norms refer to the variables u (horizontal velocity), B11 (first component of the left Cauchy-Green tensor B) and T11 (first component of the Cauchy stress tensor T).

Fig. 13. Twisting column — Beam shape and pressure distribution at output times t = 0.00375, t = 0.075, t = 0.1125, t = 0.15, t = 0.1875, t = 0.225, t = 0.2625 and t = 0.3 (from top left to bottom right). The shape is compared with respect to the initial configuration (hollow box).

Fig. 4. Left: Split of a tetrahedron into sub-tetrahedra by inserting six new midpoint edge vertices to get four corner sub-tetrahedra (colored ones). After choosing a diagonal (yellow line) to split the remaining central octahedron into fours more sub-tetrahedra, it yields a total of eight sub-tetrahedra — Right: example of 2D partitioning on NCPU = 4 threads (colors), the refinement is performed locally to each thread, only the coarse partition of large black triangles is actually built across the threads.

Fig. 5. Sketch for the elastic vibration of a beryllium plate in section 5.2 (left) and the finite deformation of a cantilever thick beam in section 5.3 (right).

Fig. 6. Elastic vibration of a beryllium plate — Numerical results at output times t = 10−5 (top), t = 2 · 10−5 (middle) and t = 3 · 10−5 (bottom) for pressure (left) and cell order map (right), the cells in yellow are at unlimited order 2, while the blue ones are the BJ limited ones. No first-order updated cell is observed.

Fig. 18. Impact of a jelly droplet — Time evolution of the maximum spreading of the droplet L/L0 in the case neo-Hookean model (a = −1, black line) or non-linear one (a = 0, red line) — The impact velocity is 2 m.s−1 (left) and 3 m.s−1 (right).

Fig. 16. Rebound of a hollow circular bar — Time evolution of vertical displacement of the points on the top plane xT = (1.6, 0, 32.4) · 10−3m and on the bottom plane xB = (1.6, 0, 4) · 10−3m (left) and percentage of bad cells detected at each time step (right).

Fig. 17. Impact of a jelly droplet with impact velocity 3 m.s−1 — Time evolution of the droplet shape at different output times for neo-Hookean model (a = −1, petroleum shade) or non-linear one (a = 0, black shade).

Fig. 11. Blake’s problem— Mesh convergence of the second order solution towards the reference solution (red) for the radial pressure (top row) and radial deviatoric stress (bottom row) at time tfinal = 1.6 · 10−4 (left) and zoom across the shock (right). Remind that LAM P1-lim corresponds to cascade: P1 → Plim1 → P0.

Fig. 7. Elastic vibration of a beryllium plate — Comparison between the MOOD cascades: P1 → P0 (LAM P0-lim) and P1 → Plim1 → P0 (LAM P1-lim) for the vertical displacement at the barycenter of the plate (left) and the computed numerical dissipation (88) as a function of time (right).

Table 1 Strong scaling computed for the L-shaped block test case presented in section 5.8.The computational time needed to perform an element update with second order of accuracy in space and time is addressed with τEU .

Fig. 15. Rebound of a hollow circular bar — Time evolution of the deformation and pressure distribution at output times at times t = 50 µs then 75, 100, 125, 150, 200, 300 and the final time t = 325 µs (from top left to bottom right).

Fig. 12. Sketch for the twisting column in section 5.5 (left) and the rebound of a hollow circular bar from section 5.6 (middle), and, L-shaped block test problem from section 5.8 (right).

Fig. 8. Cantilever thick beam test case — Pressure distribution with deformed shapes (left column) and cell order map (right column) with the second-order a posteriori limited Lagrangian scheme at output times t = 0.375, t = 0.75, t = 1.125 and t = 1.5 (from top to bottom row) — Comparison of the deformed shape computed using the simpler cascade P1 → P0 in black line on the left panels only.

Frequently asked questions

Q1. What have the authors contributed in "A 3d cell-centered ader mood finite volume method for solving updated lagrangian hyperelasticity on unstructured grids" ?

In this paper, the authors 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. This strategy has been successfully tested in an hydrodynamics context and the present work aims at extending it to the case of hyperelasticity. 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. These test cases feature material bending, impact, compression, non-linear deformation and further bouncing/detaching motions. 

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.