A structure-preserving integrator for incompressible finite elastodynamics based on a grad-div stabilized mixed formulation with particular emphasis on stretch-based material models

Abstract

We present a structure-preserving scheme based on a recently-proposed mixed formulation for incompressible hyperelasticity formulated in principal stretches. Although there exist several different Hamiltonians introduced for quasi-incompressible elastodynamics based on different multifield variational formulations, there is not much study on the fully incompressible materials in the literature. The adopted mixed formulation can be viewed as a finite-strain generalization of Herrmann variational formulation, and it naturally provides a new Hamiltonian for fully incompressible elastodynamics. Invoking the discrete gradient and scaled mid-point formulas, we are able to design fully-discrete schemes that preserve the Hamiltonian and momenta. Our analysis and numerical evidence also reveal that the scaled mid-point formula is non-robust numerically. The generalized Taylor–Hood element based on the spline technology conveniently provides a higher-order, robust, and inf-sup stable spatial discretization option for finite strain analysis. To enhance the element performance in volume conservation, the grad-div stabilization, a technique initially developed in computational fluid dynamics, is introduced here for elastodynamics. It is shown that the stabilization term does not impose additional restrictions for the algorithmic stress to respect the invariants, leading to an energy-decaying and momentum-conserving fully discrete scheme. A set of numerical examples is provided to justify the claimed properties. The grad-div stabilization is found to enhance the discrete mass conservation effectively. Furthermore, in contrast to conventional algorithms based on Cardano’s formula and perturbation techniques, the spectral decomposition algorithm developed by Scherzinger and Dohrmann is robust and accurate to ensure the discrete conservation laws and is thus recommended for stretch-based material modeling.

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-$\alpha$ [27] and generalized-$\alpha$ 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 Hamiltonian¹

$$\tilde{H}≔{\int }_{{\Omega }_{\mathit{X}}}\frac{1}{2}{\rho }_0{\Vert \mathit{V}\Vert }^2+W\left({\Theta }^{2/3}\tilde{\mathit{C}}\right)+p\left(J-\Theta \right)+\lambda \left(\Theta -1\right)d{\Omega }_{\mathit{X}}$$

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 $\tilde{H}$ 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 form²

$$\hat{H}≔{\int }_{{\Omega }_{\mathit{X}}}\frac{1}{2}{\rho }_0{\Vert \mathit{V}\Vert }^2+W\left(\mathit{C}\right)+\frac{\varepsilon }{2}{p}^{2}d{\Omega }_{\mathit{X}}$$

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

$$\check{H}≔{\int }_{{\Omega }_{\mathit{X}}}\frac{1}{2}{\rho }_0{\Vert \mathit{V}\Vert }^2+W\left(\mathit{C}\right)d{\Omega }_{\mathit{X}}$$

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 $J$, the volume ratio. In the incompressible limit, the volumetric energy becomes an indeterminate form of the type $0\times \infty$ 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., 

$$H≔{\int }_{{\Omega }_{\mathit{X}}}\frac{1}{2}{\rho }_0{\Vert \mathit{V}\Vert }^2+G_{\mathrm{ich}}\left(\tilde{\mathit{C}}\right)d{\Omega }_{\mathit{X}}.$$

We mention that the isochoric potential energy $G_{\mathrm{ich}}$ 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 $J$ 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 $J=1$, 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.