scispace - formally typeset

Journal ArticleDOI

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

15 Jan 2022-Journal of Computational Physics (Academic Press)-Vol. 449, pp 110779

AbstractIn 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.

Topics: Finite volume method (54%), Discretization (54%), Order of accuracy (53%), Numerical analysis (52%), Solver (51%)

Summary (9 min read)

Jump to: [1 Introduction][2 Updated Lagrangian hyperelastic modeling for isotropic materials][2.1.1 Lagrange-Euler mapping.][2.1.2 Measures of deformation][2.1.3 Geometric conservation law (GCL)][2.2 Governing equations][2.2.1 Conservation laws][2.2.2 Constitutive law for isotropic materials][2.2.3 Volumetric shear strain decomposition][2.2.4 Examples of constitutive laws][2.3 Updated Lagrangian hyperelasticity for isotropic materials][3.1.1 Geometrical entities][3.1.2 Conservative and constitutive discrete variables][3.2 Discrete velocity gradient operators][3.3.1 Conservation laws - GCL, momentum and total energy][3.3.2 The subcell force expression ensuring thermodynamic consistency][3.3.3 The nodal solver ensuring the conservation of total energy and momentum][3.4.1 Physical conservation laws][3.4.2 Space-time discretization of the left Cauchy-Green tensor equation][3.5 Limiting: a posteriori MOOD loop][3.6 Internal consistency consideration][3.7 Time-step monitoring][3.8 Boundary condition treatments][4.1 Algorithm][4.2 Meshing and parallelization][5 2D and 3D test problems][5.1 Swinging plate][5.2 Elastic vibration of a Beryllium plate][5.3 Finite deformation of a cantilever thick beam][5.5 Twisting column][5.6 Rebound of a hollow circular bar][5.7 Impact of a jelly-like droplet][5.8 L-shaped block][5.9 Monitoring of the internal consistency] and [6 Conclusions and perspectives]

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

Content maybe subject to copyright    Report

HAL Id: hal-03215358
https://hal.archives-ouvertes.fr/hal-03215358
Preprint submitted on 3 May 2021
HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entic research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.
L’archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diusion de documents
scientiques de niveau recherche, publiés ou non,
émanant des établissements d’enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
A 3D cell-centered ADER MOOD Finite Volume
method for solving updated Lagrangian hyperelasticity
on unstructured grids
Walter Boscheri, Raphaël Loubère, Pierre-Henri Maire
To cite this version:
Walter Boscheri, Raphaël Loubère, Pierre-Henri Maire. A 3D cell-centered ADER MOOD Finite
Volume method for solving updated Lagrangian hyperelasticity on unstructured grids. 2021. �hal-
03215358�

A 3D cell-centered ADER MOOD Finite
Volume method for solving updated
Lagrangian hyperelasticity on unstructured
grids
Walter Boscheri
a
, Rapha¨el Loub`ere
b
, Pierre-Henri Maire
c
.
a
Dipartimento di Matematica e Informatica, Ferrara Italy
b
Institut de Math´ematiques de Bordeaux (IMB), Talence, France
c
CEA-CESTA, Le Barp, France
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 discretiza-
tion named EUCCLHYD and introduced in the context of Lagrangian hydrody-
namics. 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 or-
der schemes using DERivatives) approach is adopted to obtain second-order of ac-
curacy in time. This strategy has been successfully tested in an hydrodynamics
context and the present work aims at extending it to the case of hyperelasticity.
Here, the hyperelasticty equations are written in the updated Lagrangian frame-
work and the dedicated Lagrangian numerical scheme is derived in terms of nodal
solver, Geometrical Conservation Law (GCL) compliance, subcell forces and com-
patible discretization. The Lagrangian numerical method is implemented in 3D un-
der 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, com-
pression, non-linear deformation and further bouncing/detaching motions.
Key words: Cell-centered Lagrangian finite volume schemes, moving unstructured
Preprint submitted to Elsevier Science 30 April 2021

meshes, a posteriori MOOD limiting, ADER, hyper-elasticity
Contents
1 Introduction 2
2 Updated Lagrangian hyperelastic modeling for isotropic materials 4
2.1 Kinematics 5
2.2 Governing equations 8
2.3 Updated Lagrangian hyperelasticity for isotropic materials 13
3 Finite volume discretization 14
3.1 Mesh and notation 14
3.2 Discrete velocity gradient operators 16
3.3 Semi-discretization in space 17
3.4 Space-Time discretization ADER methodology 20
3.5 Limiting: a posteriori MOOD loop 27
3.6 Internal consistency consideration 28
3.7 Time-step monitoring 29
3.8 Boundary condition treatments 29
4 Implementation considerations 30
4.1 Algorithm 30
4.2 Meshing and parallelization 31
5 2D and 3D test problems 31
5.1 Swinging plate 32
5.2 Elastic vibration of a Beryllium plate 33
5.3 Finite deformation of a cantilever thick beam 35
5.4 Blake’s problem 37
5.5 Twisting column 38
5.6 Rebound of a hollow circular bar 39
5.7 Impact of a jelly-like droplet 42
5.8 L-shaped block 43
5.9 Monitoring of the internal consistency 45
6 Conclusions and perspectives 47
A Principal invariants of a tensor 48
B Boundary conditions (BCs) 49
References 50
1 Introduction
This work is concerned with the accurate multi-dimensional simulation of hyper-elasticity
models by updated Lagrangian Finite Volume (FV) numerical scheme.
Traditionally displacement-based Finite Element (FE) frameworks are employed when sim-
ulating engineering large strain transient situations, see for instance [34,27,6]. Tetrahedral
Email addresses: walter.boscheri@unife.it (Walter Boscheri),
raphael.loubere@u-bordeaux.fr (Rapha¨el Loub`ere), pierre-henri.maire@cea.fr
(Pierre-Henri Maire).
2

mesh elements are often preferred for the discretization because the mesh generators are eas-
ily available. Nonetheless, such FE approach presents some well reported issues, for instance
a reduced order of convergence for strains and stresses [48], poor performance in nearly in-
compressible bending dominated scenarios [40] and spurious numerical instabilities (shear or
volumetric locking, parasitical hydrostatic pressure fluctuations, etc.) [5], reduction of the
accuracy when artificial stabilization terms are added. Notice that the introduction of high
order schemes usually helps but does not cure all of these issues.
The earliest attempt at solving solid dynamics equations in a FV framework is to be found
in the papers [56,55]. In this seminal approach, the solid dynamics equations are written un-
der Eulerian formalism as a first-order system of conservation laws. Their FV discretization
relies on Godunov-type methods.
In a series of papers [1,8,33,9] an original and promising FV computational framework has
been constructed for solving solid dynamics equations written under total Lagrangian formu-
lation. Impressive three-dimensional simulations of hyperelastic materials undergoing large
deformations are achieved within this framework. The underlying numerical fluxes are ob-
tained employing an acoustic approximate Riemann solver. 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. Here, the defor-
mation gradient is nothing but the Lagrange-Euler mapping which relates the initial and the
current configuration of the material under consideration. Let us point out that these three
supplementary geometrical conservation laws are redundant and nothing in this numerical
method guarantees their consistency at the discrete level.
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].
Others have explored the used of a mixed velocity-pressure approach using a Variational
Multi-Scale method [53], possibly with stabilized nodal stresses [19].
This work also considers the updated Lagrangian framework like for [39,14] contrarily to
[8,33,9]. Previously in [13] we have presented a second-order accurate cell-centered La-
grangian scheme on unstructured mesh restricted to the hydrodynamics system of conserva-
tion laws. This scheme is constructed upon a subcell discretization, popularized for staggered
Lagrangian schemes [15,16] and later extended to cell-centered ones [45], further associated
with a nodal solver relying on total energy conservation and Geometrical Conservation Law
(GCL) compliance. Second-order of accuracy is usually gained by a MUSCL-like approach
—piecewise linear reconstructions supplemented with limiters— and a predictor-corrector,
Runge-Kutta or a Generalized Riemann Problem (GRP) time discretization.
Contrarily, for the scheme in [13], we have relied on ADER (Arbitrary high order schemes
using DERivatives) methodology developed originally from an Eulerian perspective [54,10].
This approach is further supplemented with an a posteriori MOOD limiting strategy [21] to
stabilize and produce a fail-safe Lagrangian scheme. We have shown in [13] that such a cell-
centered numerical method is able to perform on classical and demanding hydrodynamics
test cases using unstructured mesh made of simplexes in 2D and 3D.
In this work we propose the extension of this numerical method to solve problems involving
2D and 3D elastic materials. We ought to solve an hyperelasticity model of PDEs (Partial
3

Differential Equations) [39,12,40,14]. Historically hypoelasticity models [57,7,58] have been
sometimes preferred, see for instance [59,30,47,52,20]. A parallel discussion about hypo- and
hyper-elastic models and their resolution can be found for instance in [49]. Here, we are
interested in simulating the large deformations undergone by isotropic hyperelastic materi-
als. The set of conservation laws under consideration is written under updated Lagrangian
representation and equipped with a constitutive law which satisfies not only the principle of
material frame indifference but is also thermodynamically consistent [35]. In this approach,
the Cauchy stress tensor is the derivative of the free energy with respect to the left Cauchy-
Green tensor which is a relevant measure of deformation. Moreover, one can show that the
Cauchy stress tensor might be expressed uniquely in terms of this deformation tensor. Tak-
ing advantage of the simplicial grid, we introduce a specific discretization of the time rate
of change of the left Cauchy Green tensor to ensure its compatibility with the Geometric
Conservation Law, refer to [28]. In this article are tackled several issues related to the 3D
extension of our ADER Lagrangian scheme, as well as the increase of complexity in the
modeling of hyper-elastic materials. First, the hyperelastic model demands the resolution
of a constitutive law which, in the framework of ADER methodology, requires some adap-
tation. Second, the a posteriori MOOD limiting strategy must consider new admissibility
criteria brought by the model in order to still ensure the robust and fail-safe characteristics
while maintaining an acceptable, if not optimal, accuracy. Third, the boundary conditions
(BCs) must be dealt with care as materials may balistically fly but also impact, bounce,
slide, spread, tear apart onto a wall or different materials. Fourth, in relation to the points
above, 3D Lagrangian numerical code requires extra-care as efficient simulations demand a
well designed parallelization methodology as well as appropriate BCs and robust limiting
strategy.
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.
For a broad and modern introductions on hypo- or hyper-elasticity we refer the readers in
particular to [39,12,14,49]. For 3D Lagrange updated FV numerical schemes among many
works we refer to [17,43,31,14]. In the forthcoming section the paper presents in details the
hyperelastic model and the governing equations to be solved. Then in the third section, the
Lagrangian numerical scheme is introduced with emphasis on the ADER approach, the nodal
solver, the a posteriori limiting strategy and some discussion on the internal consistency,
boundary conditions and implementation considerations. All numerical tests are gathered in
the fourth section. We propose the numerical results of our simulations for a large set of 2D
and 3D problems involving elastic materials impacting, detaching, compressing, swinging,
twisting, etc. Conclusions and perspectives are finally drawn in the last section.
2 Updated Lagrangian hyperelastic modeling for isotropic materials
In this section, we 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.
4

Figures (23)
References
More filters

Journal ArticleDOI
Abstract: Gmsh is an open-source 3-D finite element grid generator with a build-in CAD engine and post-processor. Its design goal is to provide a fast, light and user-friendly meshing tool with parametric input and advanced visualization capabilities. This paper presents the overall philosophy, the main design choices and some of the original algorithms implemented in Gmsh. Copyright (C) 2009 John Wiley & Sons, Ltd.

4,113 citations


Book
12 Sep 2000
Abstract: Preface. List of Boxes. Introduction. Lagrangian and Eulerian Finite Elements in One Dimension. Continuum Mechanics. Lagrangian Meshes. Constitutive Models Solution Methods and Stability. Arbitrary Lagrangian Eulerian Formulations. Element Technology. Beams and Shells. Contact--Impact. Appendix 1: Voigt Notation. Appendix 2: Norms. Appendix 3: Element Shape Functions. Glossary. References. Index.

3,818 citations


Proceedings ArticleDOI
01 Jan 1989
TL;DR: Cell-centered and mesh-vertex upwind finite-volume schemes are developed which utilize multi-dimensional monotone linear reconstruction procedures which differ from existing algorithms (even on structured meshes).
Abstract: Solution and mesh generation algorithms for solving the Euler equations on unstructured meshes consisting of triangle and quadrilateral control volumes are presented. Cell-centered and mesh-vertex upwind finite-volume schemes are developed which utilize multi-dimensional monotone linear reconstruction procedures. These algorithms differ from existing algorithms (even on structured meshes). Numerical results in two dimensions are presented.

2,043 citations


Journal ArticleDOI
TL;DR: This paper presents and study a class of graph partitioning algorithms that reduces the size of the graph by collapsing vertices and edges, they find ak-way partitioning of the smaller graph, and then they uncoarsen and refine it to construct ak- way partitioning for the original graph.
Abstract: In this paper, we present and study a class of graph partitioning algorithms that reduces the size of the graph by collapsing vertices and edges, we find ak-way partitioning of the smaller graph, and then we uncoarsen and refine it to construct ak-way partitioning for the original graph. These algorithms compute ak-way partitioning of a graphG= (V,E) inO(|E|) time, which is faster by a factor ofO(logk) than previously proposed multilevel recursive bisection algorithms. A key contribution of our work is in finding a high-quality and computationally inexpensive refinement algorithm that can improve upon an initialk-way partitioning. We also study the effectiveness of the overall scheme for a variety of coarsening schemes. We present experimental results on a large number of graphs arising in various domains including finite element methods, linear programming, VLSI, and transportation. Our experiments show that this new scheme produces partitions that are of comparable or better quality than those produced by the multilevel bisection algorithm and requires substantially smaller time. Graphs containing up to 450,000 vertices and 3,300,000 edges can be partitioned in 256 domains in less than 40 s on a workstation such as SGI's Challenge. Compared with the widely used multilevel spectral bisection algorithm, our new algorithm is usually two orders of magnitude faster and produces partitions with substantially smaller edge-cut.

1,619 citations


Journal ArticleDOI
TL;DR: The basic explicit finite element and finite difference methods that are currently used to solve transient, large deformation problems in solid mechanics are reviewed.
Abstract: Explicit finite element and finite difference methods are used to solve a wide variety of transient problems in industry and academia. Unfortunately, explicit methods are rarely discussed in detail in finite element text books. This paper reviews the basic explicit finite element and finite difference methods that are currently used to solve transient, large deformation problems in solid mechanics. A special emphasis has been placed on documenting methods that have not been previously published in journals.

1,109 citations


Frequently Asked Questions (2)
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. 

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.