A structure-preserving integrator for incompressible finite elastodynamics based on a grad-div stabilized mixed formulation with particular emphasis on stretch-based material models
Introduction
Preserving physical and geometrical properties in the fully discrete model has been an active research topic in computational science and engineering for decades. The foremost property is stability which guarantees numerical convergence. The pathological energy growth of the trapezoidal rule suggests stable schemes developed for linear problems are generally inapplicable to nonlinear analysis [1], and it thus stimulates the development of energy-conserving schemes. In continuum mechanics, the total linear and angular momenta are often demanded to be preserved in the solutions; from the perspective of mathematical physics, preserving the symplectic character in the phase space can be salubrious. Designing algorithms that preserve invariants has led to a research area known as the structure preserving integrators, and abundant numerical evidences have justified its advantage. Nevertheless, it is by no means a trivial task to preserve energy and symplecticity simultaneously [2], [3]. The superiority of either invariant over the other is controversial [4], [5], [6]. In this study, we restrict our discussion to schemes that preserve or dissipate energy and conserve momenta, and they are often referred to as the energy–momentum consistent schemes.
The development of energy–momentum consistent schemes heretofore is briefly summarized as follows. The initial attempt was made by enforcing the discrete conservation laws via the Lagrange multiplier method [7]. Its major drawback is that adding constraints engenders difficulty in achieving convergence for the nonlinear solver [8]. Inspired by the conserving algorithms developed for particle dynamics [9], an energy–momentum conserving method was introduced by modifying the mid-point rule with the stress evaluated in terms of a convex combination of the right Cauchy–Green tensors within a time interval [10]. The combination parameter needs to be determined iteratively through an algorithm for general materials based on the mean value theorem [11]. To alleviate the algorithmic complexity, the discrete gradient formula was proposed to provide an algorithmic stress definition, which conserves energy and momentum and conveniently works for general material models [12]. Following that, several alternate algorithmic stress formulas have been proposed [13], and they have been successfully applied to shell dynamics [14], multibody systems [15], contact and impact [16], [17], and incompressible elastodynamics [12]. A Petrov–Galerkin finite element formulation was utilized in time for Hamiltonian systems [18], and the aforementioned discrete gradient method can be realized within that framework by adopting a special quadrature formula in time. Recently, a framework based on polyconvexity and the notion of tensor cross product has been developed [19], [20], [21], and its combination with the discrete gradient method has further led to a multifield variational formulation [22], with applications in thermo-elasticity [23], viscoelasticity [24] and electro-mechanics [25], [26], to list a few.
The energy–momentum schemes have been further generalized with dissipation effects taken into account. Dissipation arises due to numerical and physical mechanisms. The former is often purposely introduced to damp the error caused by poorly resolved modes from spatial discretization. This goal was first achieved in linear dynamics by the HHT- [27] and generalized- schemes [28]. In nonlinear analysis, different approaches were made to introduce analogous damping effects on the high-frequency modes [29], [30], [31], [32]. On the other hand, using the energy–momentum method in non-reversible processes may also lead to the decay of the energy. Examples include elastoplasticity [33], [34], [35], viscoelasticity [36], [37], and thermo-elasticity [23], [38]. It is worth mentioning that the General Equations for Non-Equilibrium Reversible Irreversible Coupling (GENERIC) framework provides a systematic approach of extending the energy–momentum schemes to systems with physical dissipative mechanisms [24], [39].
Many engineering and biomedical materials exhibit volume-preserving behavior under large deformation. Devising a structure-preserving scheme that works well for incompressible elastodynamics is propitious yet non-trivial. It is non-trivial because the variational formulation typically needs to be modified to account for kinematic constraints. The invariants embedded in the variational formulation may take different forms accordingly. This further calls for new strategies in the design of the conserving integrator. In [12], Gonzalez resorted to the quasi-incompressible formulation based on the three-field variational principle [40] and identified a modified Hamiltonian1 as an invariant for the three-field variational formulation. Invoking the discrete gradient formula, the author designed a time-stepping scheme that preserves the modified Hamiltonian over time. In [17], the authors considered an alternate approach based on a two-field variational formulation [41], in which a pressure-like variable is introduced to be interpolated independently. The Hamiltonian of the form2 is conserved after applying the discrete gradient formula to design an algorithmic stress. However, that study was still restricted to the quasi-incompressible scenario. Recently, a seven-field variational formulation was proposed [22], in which the authors obtained the that the total energy, i.e., can be conserved by a careful design of the algorithmic stress. Later, the strategy was extended by focusing on quasi-incompressible stretch-based material models [42]. It is remarkable that the result is more physically relevant than the prior result of Gonzalez [12].
We mention that the above Hamiltonians suffer from a fundamental issue in their volumetric part. Under the setting of the Helmholtz free energy, the volumetric energy involves a bulk modulus and is written in terms of , the volume ratio. In the incompressible limit, the volumetric energy becomes an indeterminate form of the type as the bulk modulus approaches infinity. Thus, the strategies developed and conclusions made in prior studies cannot be straightforwardly applied to the fully incompressible regime.
There have been continued research efforts towards a structure-preserving integrator for nonlinear problems under constraints. Most of them exploited the Hu–Washizu type variational formulations [22], [42], [43], [44], while some leverage the Hellinger–Reissner formulation [45], [46]. To our knowledge, the finite elasticity in the fully incompressible regime has been rarely discussed. In this study, the variational formulation is based on a finite-strain extension of the Herrmann variational formulation [47], [48], [49], [50]. The derivation of the formulation was based on a complementary energy which was obtained from a Legendre transformation performed on the volumetric energy, and it was termed as the Gibbs free energy in the work of [48]. Based on the mixed formulation, we will show that the invariants of motion can be described in a detailed manner as the summation of the kinetic energy and the isochoric part of the potential energy, i.e., We mention that the isochoric potential energy is, in fact, identical to that of the Helmholtz strain energy [48]. Based on this, one primary goal of this work is to design a scheme that preserves the invariants in the fully discrete formulation. Inspired by the discrete gradient approach [12], an algorithmic stress formula is applied for the isochoric part of the stress to achieve this goal. The scaled mid-point formula [51], [52] is also considered as an alternate option for constructing the algorithmic stress. Although it is theoretically appealing [52], [53], its denominator is not a metric for the difference of the deformation tensors, and its numerical robustness becomes questionable. Our numerical evidence confirms this concern, and the scaled mid-point formula is not recommended for future investigation.
In this study, the material models are formulated in principal stretches to account for general isotropic hyperelastic materials. In particular, the stress and the elasticity tensor of the Ogden-type materials are carefully derived. A missing term in the documented formula for the elasticity tensor [54, Eqn. (6.197)] is identified, which has an impact on the convergence behavior for the Newton–Raphson iteration. We mention that, in practice, the quality of the discrete conservation laws is strongly influenced by the accuracy of the solution for the nonlinear systems [55]. Therefore, the missing term can be critical in certain applications. Moreover, for the eigenvalue-based models, the algorithm for the spectral decomposition impacts the solution quality substantially. To date, algorithms based on Cardano’s formula are still widely adopted. However, it necessitates a perturbation technique when two or three eigenvalues coincide [56], [57], which is known to have a detrimental impact on the eigenvector accuracy [58], [59]. In the context of conserving integrators, the perturbation technique further causes the loss of the discrete conservation properties and demands a specially-developed strategy to mitigate that issue [60]. The algorithm developed by Scherzinger and Dohrmann is a robust alternative for spectral decomposition [59]. Its accuracy has been demonstrated to be comparable with that of LAPACK [61]. In this study, we adopt this algorithm for the constitutive modeling of Ogden-type materials. Numerical tests indicate that it is as accurate as the invariant-based models in evaluating the algorithmic stresses, and there is thus no need to invoke the aforesaid strategy [60].
For incompressible elastodynamics, the volume ratio is an invariant for the motion as well. Its preservation in the integrator is worth pursuing and is known as an “unsolved issue” in computational elasticity [62], [63], [64]. Different from prior approaches that strive to satisfy the nonlinear kinematic constraint , we propose a plan of two essential steps within our framework. First, a spatial discretization strategy is needed to ensure the discrete satisfaction of the divergence-free condition. Second, with that velocity field, a time integration algorithm for the deformation state is demanded to ensure volume preservation. A time-stepping algorithm that satisfies this point is known as the volume-conserving algorithm [65], [66]. In this work, the Galerkin projection is made with a smooth generalization of the Taylor–Hood element using the spline technology [67], [68]. This type of mixed element is convenient for implementation, inf-sup stable, robust for large-strain analysis, and related to the concept of isogeometric analysis [69]. Yet, a drawback is that its constraint ratio is large [70, Sec. 4.3.7], especially when the polynomial degree is low. To enhance the solution quality, the grad-div stabilization technique is introduced in our proposed formulation. This stabilization technique has been found to be rather effective in improving the discrete mass conservation in fluid mechanics and transport problems [71], [72], [73] and can be interpreted as a sub-grid model [74]. Interestingly, it has been proved that, as the stabilization parameter approaches infinity, the solution of the Taylor–Hood element with the grad-div stabilization converges to that of the Scott–Vogelius element [75], which is discretely divergence-free [76]. In this work, we introduce the grad-div stabilization to the study of elastodynamics as the first step in our proposed roadmap. We also show that the grad-div stabilization term respects the momentum conservation and dissipates the energy, thereby rendering an energy-decaying and momentum-conserving scheme. It is anticipated that the proposed numerical framework may offer a robust and accurate approach for elastodynamic analysis.
The body of this work is organized as follows. In Section 2, the strong-form problem and the constitutive relations are presented. After that, we present the major contribution of this work in Section 3, including the energy–momentum consistent formulation, different algorithmic stress designs, grad-div stabilization, and a segregated predictor multi-corrector algorithm. In Section 4, the claimed numerical attributes are demonstrated by numerical examples. We draw conclusions in Section 5.
Section snippets
Kinematics
In this work, the colon-equal sign indicates the definition of quantities. We start by summarizing the notations to be used in the formulation of the problem. Let the bounded open sets and represent the initial and current configurations of a continuum body, respectively. Their boundaries are assumed to be Lipschitz and are denoted by and , with unit outward normals and , respectively. The boundary can be decomposed into two non-overlapping parts, that is, ,
Semi-discrete formulation
In this section, we present the spatially discrete formulation for the initial–boundary value problem. The spatial projection is made by adopting a smooth generalization of the Taylor–Hood element based on the spline technology [84] for the following reasons. It is convenient as the displacement and velocity are interpolated via isoparametric elements. The spline technology enables a straightforward way of integrating with the CAD system within the paradigm of isogeometric analysis [69]. In the
Numerical examples
In our numerical investigation, unless otherwise specified, the discrete pressure function space is generated by -refinement to achieve the highest possible continuity, and the degree elevation is adopted with and to construct the velocity function space; the Gonzalez discrete gradient is used to construct the algorithmic stress; we use Gaussian quadrature points in each direction; we use , , as the stopping criteria in the predictor multi-corrector
Conclusion
In this work, we first discussed the invariants for incompressible elastodynamics and designed a consistent scheme with respect to the energy and momenta. In particular, the Hamiltonian for fully incompressible materials is identified, in which the potential energy only involves the isochoric part of the energy. Considering that the previously existing studies are restricted to quasi-incompressible materials [12], [17], this newly identified Hamiltonian can be beneficial in analyzing dynamical
Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Acknowledgments
This work is supported by the National Natural Science Foundation of China [Grant Numbers 12172160, 12072143], Shenzhen Science and Technology Program [Grant Number JCYJ20220818100600002], Southern University of Science and Technology, China [Grant Number Y01326127], and the Department of Science and Technology of Guangdong Province [Grant Numbers 2020B1212030001, 2021QN020642]. Computational resources are provided by the Center for Computational Science and Engineering at the Southern
References (102)
- et al.
Constraint energy momentum algorithm and its application to non-linear dynamics of shells
Comput. Methods Appl. Mech. Engrg.
(1996) Conservative numerical methods for
J. Comput. Phys.
(1984)- et al.
A new solution procedure for application of energy-conserving algorithms to general constitutive models in nonlinear elastodynamics
Comput. Methods Appl. Mech. Engrg.
(2001) Exact energy and momentum conserving algorithms for general models in nonlinear elasticity
Comput. Methods Appl. Mech. Engrg.
(2000)- et al.
Dynamic analysis of rigid and deformable multibody systems with penalty methods and energy-momentum schemes
Comput. Methods Appl. Mech. Engrg.
(2000) - et al.
Energy-controlling time integration methods for nonlinear elastodynamics and low-velocity impact
Comput. Methods Appl. Mech. Engrg.
(2006) - et al.
A first order hyperbolic framework for large strain computational solid dynamics. Part I: Total Lagrangian isothermal elasticity
Comput. Methods Appl. Mech. Engrg.
(2015) - et al.
A computational framework for polyconvex large strain elasticity
Comput. Methods Appl. Mech. Engrg.
(2015) - et al.
A new mixed finite element based on different approximations of the minors of deformation tensors
Comput. Methods Appl. Mech. Engrg.
(2011) - et al.
A mixed variational framework for the design of energy-momentum schemes inspired by the structure of polyconvex stored energy functions
Comput. Methods Appl. Mech. Engrg.
(2018)
A new energy-momentum time integration scheme for non-linear thermo-mechanics
Comput. Methods Appl. Mech. Engrg.
An energy-momentum time integration scheme based on a convex multi-variable framework for non-linear electro-elastodynamics
Comput. Methods Appl. Mech. Engrg.
A thermodynamically consistent time integration scheme for non-linear thermo-electro-mechanics
Comput. Methods Appl. Mech. Engrg.
Formulation and analysis of conserving algorithms for frictionless dynamic contact/impact problems
Comput. Methods Appl. Mech. Engrg.
On the formulation of high-frequency dissipative time-stepping algorithms for nonlinear dynamics. Part I: low order methods for two model problems and nonlinear elastodynamics
Comput. Methods Appl. Mech. Engrg.
On the formulation of high-frequency dissipative time-stepping algorithms for nonlinear dynamics. Part II: second-order methods
Comput. Methods Appl. Mech. Engrg.
Generalized energy-momentum method for non-linear adaptive shell dynamics
Comput. Methods Appl. Mech. Engrg.
Energy-dissipative momentum-conserving time-stepping algorithms for finite strain multiplicative plasticity
Comput. Methods Appl. Mech. Engrg.
Volume-preserving energy-momentum schemes for isochoric multiplicative plasticity
Comput. Methods Appl. Mech. Engrg.
Energy consistent algorithms for dynamic finite deformation plasticity
Comput. Methods Appl. Mech. Engrg.
Variational and projection methods for the volume constraint in finite deformation elasto-plasticity
Comput. Methods Appl. Mech. Engrg.
A finite element formulation for nonlinear incompressible elastic and inelastic analysis
Comput. Struct.
Variational-based locking-free energy–momentum schemes ofhigher-order for thermo-viscoelastic fiber-reinforced continua
Comput. Methods Appl. Mech. Engrg.
Energy-decaying and momentum-conserving schemes for transient simulations with mixed finite elements
Comput. Methods Appl. Mech. Engrg.
A unified continuum and variational multiscale formulation for fluids, solids, and fluid-structure interaction
Comput. Methods Appl. Mech. Engrg.
On a variational theorem for incompressible and nearly-incompressible orthotropic elasticity
Int. J. Solids Struct.
On the stability of symplectic and energy-momentum algorithms for non-linear Hamiltonian systems with symmetry
Comput. Methods Appl. Mech. Engrg.
Algorithms for static and dynamic multiplicative plasticity that preserve the classical return mapping schemes of the infinitesimal theory
Comput. Methods Appl. Mech. Engrg.
A robust algorithm for finding the eigenvalues and eigenvectors of symmetric matrices
Comput. Methods Appl. Mech. Engrg.
Galerkin-based mechanical integrators for finite elastodynamics formulated in principal stretches–Pitfalls and remedies
Comput. Methods Appl. Mech. Engrg.
A stability study of some mixed finite elements for large deformation elasticity problems
Comput. Methods Appl. Mech. Engrg.
The importance of the exact satisfaction of the incompressibility constraint in nonlinear elasticity: mixed FEMs versus NURBS-based approximations
Comput. Methods Appl. Mech. Engrg.
A continuum and computational framework for viscoelastodynamics: I. finite deformation linear models
Comput. Methods Appl. Mech. Engrg.
Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement
Comput. Methods Appl. Mech. Engrg.
A low order Galerkin finite element method for the Navier–Stokes equations of steady incompressible flow: a stabilization issue and iterative methods
Comput. Methods Appl. Mech. Engrg.
Mixed finite element methods with convection stabilization for the large eddy simulation of incompressible turbulent flows
Comput. Methods Appl. Mech. Engrg.
Grad-div stabilization and subgrid pressure models for the incompressible Navier-Stokes equations
Comput. Methods Appl. Mech. Engrg.
An implicit thermomechanical theory based on a Gibbs potential formulation for describing the response of thermoviscoelastic solids
Internat. J. Engrg. Sci.
An extension of Herrmann’s principle to nonlinear elasticity
Appl. Math. Model.
Quasi-incompressible finite elasticity in principal stretches. continuum basis and numerical algorithms
Comput. Methods Appl. Mech. Engrg.
Robustness of isogeometric structural discretizations under severe mesh distortion
Comput. Methods Appl. Mech. Engrg.
Subdivision-stabilised immersed b-spline finite elements for moving boundary flows
Comput. Methods Appl. Mech. Engrg.
Multipatch isogeometric analysis for electrophysiology: Simulationin a human heart
Comput. Methods Appl. Mech. Engrg.
Two classes of mixed finite element methods
Comput. Methods Appl. Mech. Engrg.
Numerical studies of finite element variational multiscale methodsfor turbulent flow simulation
Comput. Methods Appl. Mech. Engrg.
Error analysis and iterative solvers for Navier-Stokes projection methods with standard and sparse grad-div stabilization
Comput. Methods Appl. Mech. Engrg.
Three-dimensional nonlinear dynamics of slender structures: Cosserat rod element approach
Int. J. Solids Struct.
Analysis of transient algorithms with particular reference to stability behavior
Lie–Poisson Hamilton–Jacobi theory and Lie-Poisson integrators
Phys. Lett. A
Symplectic-energy-momentum preserving variational integrators
J. Math. Phys.
Cited by (3)
A continuum and computational framework for viscoelastodynamics: II. Strain-driven and energy–momentum consistent schemes
2023, Computer Methods in Applied Mechanics and Engineering