An energy-stable mixed formulation for isogeometric analysis of incompressible hyper-elastodynamics

TL;DR

This paper introduces an energy-stable mixed formulation for isogeometric analysis of incompressible hyper-elastodynamics, linking fluid and solid dynamics, with demonstrated stability, convergence, and applicability to dynamic problems.

Contribution

It presents a novel energy-stable mixed formulation using NURBS-based elements for incompressible hyper-elastodynamics, improving stability and convergence over traditional methods.

Findings

Energy stability estimate established for the discretization.

Numerical assessment confirms inf-sup stability for various NURBS pairs.

The formulation demonstrates accurate results under different loading conditions.

Abstract

We develop a mixed formulation for incompressible hyper-elastodynamics based on a continuum modeling framework recently developed and smooth generalizations of the Taylor-Hood element based on non-uniform rational B-splines (NURBS). This continuum formulation draws a link between computational fluid dynamics and computational solid dynamics. This link inspires an energy stability estimate for the spatial discretization, which favorably distinguishes the formulation from the conventional mixed formulations for finite elasticity. The inf-sup condition is utilized to provide a bound for the pressure field. The generalized-$\alpha$ method is applied for temporal discretization, and a nested block preconditioner is invoked for the solution procedure. The inf-sup stability for different pairs of NURBS elements is elucidated through numerical assessment. The convergence rate of the proposed formulation with various combinations of mixed elements is examined by the manufactured solution method. The numerical scheme is also examined under compressive and tensile loads for isotropic and anisotropic hyperelastic materials. Finally, a suite of dynamic problems is numerically studied to corroborate the stability and conservation properties.

Tables

Table 1. Table 1: Three-dimensional tensile test: geometry setting and material properties.

Material properties: $G(\widetilde{𝑪},p)=G_{ich}^g\left(\widetilde{𝑪}\right)+{\sum }_{i=1,2}{G}_{ich}^{f_i}\left(\widetilde{𝑪}\right)+\frac{p}{{\rho }_0}$, $G_{ich}^g\left(\widetilde{𝑪}\right)=\frac{c_1}{2{\rho }_0}\left({\tilde{I}}_1-3\right)$, $G_{ich}^{f_i}\left(\widetilde{𝑪}\right)=\frac{k_1}{2k_2{\rho }_0}\left(e^{k_2{\tilde{E}}_i^2}-1\right)$, ${\tilde{E}}_i:={𝑯}_i:\widetilde{𝑪}-1$, ${𝑯}_i:=k_d𝑰+\left(1-3k_d\right)\left({𝒂}_i\otimes {𝒂}_i\right)$, ${\rho }_0=1.0\times {10}^3$ kg/m3, $c_1=7.64\times {10}^3$ Pa, $k_1=9.966\times {10}^5$ Pa, $k_2=524.6$, $k_d=0.226$.

Table 2. Table 2: The number of elements n e ​ n subscript 𝑛 𝑒 𝑛 n_{en} and the number of equations n e ​ q subscript 𝑛 𝑒 𝑞 n_{eq} in the system ( 30 ) for the five different meshes. Meshes 4 and 5 consist of linear tetrahedral elements.

Mesh	1	2	3	4	5
$n_{en}$	61440	120	120	5760	90000
$n_{eq}$	1785024	5091	7584	6396	75144

Table 3. Table 3: Three-dimensional beam bending: geometry setting, boundary conditions, and material properties.

Material properties: $G_{ich}\left(\widetilde{𝑪}\right)=\frac{c_1}{2{\rho }_0}\left({\tilde{I}}_1-3\right)+\frac{c_2}{2{\rho }_0}\left({\tilde{I}}_2-3\right)$, ${\rho }_0=1.1\times {10}^3$ kg/m3, $c_1=c_2=E/6$, $E=1.7\times {10}^7$ Pa. Reference scales: $L_0=1$ m, $M_0=1$ kg, $T_0=1$ s.

Equations

$${\mathbf{\phi }}_t(\cdot )=\mathbf{\phi }(\cdot ,t):{\Omega }_{\mathbf{X}}$$

$$\mathbf{X}$$

$$\mathbf{U}:=\mathbf{\phi }(\mathbf{X},t)-\mathbf{\phi }(\mathbf{X},0)=\mathbf{\phi }(\mathbf{X},t)-\mathbf{X},\mathbf{V}:={\frac{\partial \mathbf{\phi }}{\partial t}|}_{\mathbf{X}}={\frac{\partial \mathbf{U}}{\partial t}|}_{\mathbf{X}}=\frac{d\mathbf{U}}{dt}.$$

$$\mathbf{F}:=\frac{\partial \mathbf{\phi }}{\partial \mathbf{X}},J:=\text{det}\left(\mathbf{F}\right),\mathbf{C}:={\mathbf{F}}^T\mathbf{F}.$$

$$\tilde{\mathbf{F}}:=J^{-\frac{1}{3}}\mathbf{F},\tilde{\mathbf{C}}:=J^{-\frac{2}{3}}\mathbf{C},$$

$$G(\tilde{\mathbf{C}},p)=G_{ich}(\tilde{\mathbf{C}})+G_{vol}(p).$$

$$\rho (p):={\left(\frac{d{G}_{vol}}{dp}\right)}^{-1},{\beta }_{\theta }(p):=\frac{1}{\rho }\frac{d\rho }{dp}=-\frac{{\partial }^2{G}_{vol}}{\partial p^2}/\frac{\partial G_{vol}}{\partial p},{\mathbf{\sigma }}^{dev}:=J^{-1}\tilde{\mathbf{F}}\left(\mathbb{P}:\tilde{\mathbf{S}}\right){\tilde{\mathbf{F}}}^T,$$

$$\mathbb{P}:=\mathbb{I}-\frac{1}{3}{\mathbf{C}}^{-1}\otimes \mathbf{C},\tilde{\mathbf{S}}:=2\frac{\partial \left({\rho }_{0}G\right)}{\partial \tilde{\mathbf{C}}}=2\frac{\partial \left({\rho }_0{G}_{ich}\right)}{\partial \tilde{\mathbf{C}}},$$

$$0=\frac{d\mathbf{u}}{dt}-\mathbf{v},$$

$$0={\beta }_{\theta }(p)\frac{dp}{dt}+{\nabla }_{\mathbf{x}}\cdot \mathbf{v}$$

$$0=\rho (J)\frac{d\mathbf{v}}{dt}-{\nabla }_{\mathbf{x}}\cdot {\mathbf{\sigma }}^{dev}+{\nabla }_{\mathbf{x}}p-\rho (J)\mathbf{b},$$

$$\mathbf{u}=\mathbf{g},\text{mbox}on{\Gamma }_{\mathbf{x}}^{g,t},\mathbf{v}=\frac{d\mathbf{g}}{dt},\text{mbox}on{\Gamma }_{\mathbf{x}}^{g,t},({\mathbf{\sigma }}^{dev}-p\mathbf{I})\mathbf{n}=\mathbf{h},\text{mbox}on{\Gamma }_{\mathbf{x}}^{h,t}.$$

$$\mathbf{u}(\mathbf{x},0)={\mathbf{u}}_0(\mathbf{x}),p(\mathbf{x},0)=p_0(\mathbf{x}),\mathbf{v}(\mathbf{x},0)={\mathbf{v}}_0(\mathbf{x}).$$

$$I_1:=\text{tr}\mathbf{C},I_2:=\frac{1}{2}\left[{\left(\text{tr}\mathbf{C}\right)}^2-\text{tr}\left({\mathbf{C}}^2\right)\right].$$

$$G_{ich}(\tilde{\mathbf{C}})=\frac{c_1}{2{\rho }_0}\left({\tilde{I}}_1-3\right)+\frac{c_2}{2{\rho }_0}\left({\tilde{I}}_2-3\right),$$

$$G_{vol}(p)=\frac{\kappa }{{\rho }_0}\left(1-e^{-\frac{p}{\kappa }}\right),$$

$$\rho (p)={\rho }_0{e}^{\frac{p}{\kappa }},{\beta }_{\theta }(p)=1/\kappa .$$

$$0=\frac{1}{\kappa }\frac{dp}{dt}+\frac{d}{dt}{\nabla }_{\mathbf{x}}\cdot \mathbf{u}=\frac{d}{dt}\left(\frac{p}{\kappa }+{\nabla }_{\mathbf{x}}\cdot \mathbf{u}\right).$$

$$0=\frac{p}{\kappa }+{\nabla }_{\mathbf{x}}\cdot \mathbf{u},$$

$${\nabla }_{\mathbf{x}}\cdot {\mathbf{\sigma }}^{dev}-{\nabla }_{\mathbf{x}}p={\rho }_0\mathbf{b}.$$

$$0={\mathbf{B}}^k\left(\dot{\mathbf{y}},\mathbf{y}\right):=\frac{d\mathbf{u}}{dt}-\mathbf{v},$$

$$0={\mathbf{B}}^p\left(w_p;\dot{\mathbf{y}},\mathbf{y}\right):={\int }_{{\Omega }_{\mathbf{x}}^t}{w}_p{\nabla }_{\mathbf{x}}\cdot \mathbf{v}d{\Omega }_{\mathbf{x}},$$

$$0={\mathbf{B}}^m\left({\mathbf{w}}_{\mathbf{v}};\dot{\mathbf{y}},\mathbf{y}\right):={\int }_{{\Omega }_{\mathbf{x}}^t}{\mathbf{w}}_{\mathbf{v}}\cdot \rho (J)\frac{d\mathbf{v}}{dt}+{\nabla }_{\mathbf{x}}{\mathbf{w}}_{\mathbf{v}}:{\mathbf{\sigma }}_{dev}-{\nabla }_{\mathbf{x}}\cdot {\mathbf{w}}_{\mathbf{v}}p-{\mathbf{w}}_{\mathbf{v}}\cdot \rho (J)\mathbf{b}d{\Omega }_{\mathbf{x}}$$

$$-{\int }_{{\Gamma }_{\mathbf{x}}^{h,t}}{\mathbf{w}}_{\mathbf{v}}\cdot \mathbf{h}d{\Gamma }_{\mathbf{x}},$$

$$W_P(\mathbf{X},t):=w_p({\phi }_t(\mathbf{X}),t),$$

$$\mathbf{B}(\mathbf{X},t):=\mathbf{b}({\phi }_t(\mathbf{X}),t),$$

$$0={\mathbf{B}}^k\left(\dot{\mathbf{Y}},\mathbf{Y}\right):=\frac{d\mathbf{U}}{dt}-\mathbf{V},$$

$$0={\mathbf{B}}^p\left(W_P;\dot{\mathbf{Y}},\mathbf{Y}\right):={\int }_{{\Omega }_{\mathbf{X}}}{W}_P{\nabla }_{\mathbf{X}}\mathbf{V}:\left(J{\mathbf{F}}^{-T}\right)d{\Omega }_{\mathbf{X}},$$

$$0={\mathbf{B}}^m\left({\mathbf{W}}_{\mathbf{V}};\dot{\mathbf{Y}},\mathbf{Y}\right):={\int }_{{\Omega }_{\mathbf{X}}}{\mathbf{W}}_{\mathbf{V}}\cdot {\rho }_0\frac{d\mathbf{V}}{dt}+{\nabla }_{\mathbf{X}}{\mathbf{W}}_{\mathbf{V}}:\left(J{\mathbf{\sigma }}_{dev}{\mathbf{F}}^{-T}\right)-{\nabla }_{\mathbf{X}}{\mathbf{W}}_{\mathbf{V}}:\left(J{\mathbf{F}}^{-T}\right)P$$

$$-{\mathbf{W}}_{\mathbf{V}}\cdot {\rho }_0\mathbf{B}d{\Omega }_{\mathbf{X}}-{\int }_{{\Gamma }_{\mathbf{X}}^H}{\mathbf{W}}_{\mathbf{V}}\cdot \mathbf{H}d{\Gamma }_{\mathbf{X}}.$$

$${\mathsf{N}}_i^0(\xi )=\{\begin{array}{ll}1& \text{mbox}if{\xi }_i\le \xi <{\xi }_{i+1},\\ 0& \text{mbox}otherwise.\end{array}$$

$${\mathsf{N}}_i^{\mathsf{p}}(\xi )=\frac{\xi -{\xi }_i}{{\xi }_{i+\mathsf{p}}-{\xi }_i}{\mathsf{N}}_i^{\mathsf{p}-1}(\xi )+\frac{{\xi }_{i+\mathsf{p}+1}-\xi }{{\xi }_{i+\mathsf{p}+1}-{\xi }_{i+1}}{\mathsf{N}}_{i+1}^{\mathsf{p}-1}(\xi ).$$

$${\mathsf{R}}_i^{\mathsf{p}}(\xi ):=\frac{{\mathsf{w}}_i{\mathsf{N}}_i^{\mathsf{p}}(\xi )}{{\sum }_{j=1}^n{\mathsf{w}}_j{\mathsf{N}}_j^{\mathsf{p}}(\xi )}.$$

$$\mathbf{\alpha }:=\{{\alpha }_1,{\alpha }_2,\cdots ,{\alpha }_{\mathsf{m}-1},{\alpha }_m\}=\{-1,{\alpha }_2,\cdots ,{\alpha }_{\mathsf{m}-1},-1\}$$

$$\mathbf{\alpha }-1:=\{{\alpha }_1,{\alpha }_2-1,\cdots ,{\alpha }_{\mathsf{m}-1}-1,{\alpha }_{\mathsf{m}}\}=\{-1,{\alpha }_2-1,\cdots ,{\alpha }_{\mathsf{m}-1}-1,-1\}.$$

Full text

An energy-stable mixed formulation for isogeometric analysis of incompressible hyper-elastodynamics

Ju Liu${}^{\textup{ a}}$, Alison L. Marsden${}^{\textup{ a}}$, Zhen Tao${}^{\textup{ b}}$

a *Department of Pediatrics (Cardiology), Bioengineering,

and Institute for Computational and Mathematical Engineering, Stanford University,

Clark Center E1.3, 318 Campus Drive, Stanford, CA 94305, USA

b* *Institute for Computational Engineering and Sciences, The University of Texas at Austin,

201 East 24th Street, 1 University Station C0200, Austin, TX 78712, USA

E-mail address: liuju@stanford.edu, amarsden@stanford.edu, taozhen.cn@gmail.com *

Abstract

We develop a mixed formulation for incompressible hyper-elastodynamics based on a continuum modeling framework recently developed in [36] and smooth generalizations of the Taylor-Hood element based on non-uniform rational B-splines (NURBS). This continuum formulation draws a link between computational fluid dynamics and computational solid dynamics. This link inspires an energy stability estimate for the spatial discretization, which favorably distinguishes the formulation from the conventional mixed formulations for finite elasticity. The inf-sup condition is utilized to provide a bound for the pressure field. The generalized-$\alpha$ method is applied for temporal discretization, and a nested block preconditioner is invoked for the solution procedure. The inf-sup stability for different pairs of NURBS elements is elucidated through numerical assessment. The convergence rate of the proposed formulation with various combinations of mixed elements is examined by the manufactured solution method. The numerical scheme is also examined under compressive and tensile loads for isotropic and anisotropic hyperelastic materials. Finally, a suite of dynamic problems is numerically studied to corroborate the stability and conservation properties.

Keywords: Incompressible elasticity, Mixed formulation, Inf-sup condition, Energy stability, Generalized-$\alpha$ method, Anisotropic arterial wall model

1 Introduction

1.1 Motivation and literature survey

Over the past few decades, significant progress has been achieved in the finite element modeling of solid mechanics problems. A central topic is to devise a numerical scheme that works well in the incompressible limit. Under the small-strain assumption, this issue is well-understood, and it boils down to interpolating the displacement and pressure with elements that satisfy the Ladyzhenskaya-Babuška-Brezzi (LBB) or the inf-sup condition [4]. Under large strains, most materials exhibit volume-preserving behavior, which makes it imperative to appropriately handle the incompressibility constraint. This issue is particularly relevant for modeling biological tissues, which are largely incompressible due to their high water content. In fact, the nonlinear nature of large strain analysis, together with the kinematic constraint, makes the numerical analysis of incompressible materials quite challenging. Classical treatments of this class of problems include the $\bar{\textup{F}}$-projection method [9, 11, 24], the enhanced assumed strain method (EAS) [46, 44], and the mixed $u/p$ formulation [48].

The $\bar{\textup{F}}$-projection and EAS methods share some similarities. Both methods are developed based on heuristic splits of the deformation gradient; the geometrically linear versions of the two methods are linked with the mixed finite element method [4, 19]. Nevertheless, there are drawbacks of both. For the $\bar{\textup{F}}$-projection method, its implementation requires a nonlocal matrix inversion if the projection is onto a continuous finite element space. The EAS method relies crucially on a static condensation procedure to maintain the pure displacement code structure. The penalty nature of the pure displacement formulation inevitably induces an ill-conditioned stiffness matrix, which imposes a severe constraint on the choice of linear solvers. It has long been known that both methods suffer from mesh instability or the hourglass mode [49] and hence necessitate further refinements to numerical technologies for the hourglass control.

The mixed $u/p$ formulation introduces a pressure-like variable as the Lagrange multiplier for the incompressibility constraint in the strain energy [48]. The resulting scheme necessitates interpolating the displacement and pressure fields independently. Performing a linearization of this formulation provides a justification for the use of inf-sup stable elements [3]. Yet, for nonlinear problems, linearized stability is often insufficient to guarantee nonlinear stability [16]. It remains unclear whether there is any a priori nonlinear stability estimate for the mixed $u/p$ formulation.

In the meantime, the stabilized finite element method, as a technique initially developed for computational fluid dynamics, has been extended to solid mechanics based on various variational formulations [1, 6, 30, 36, 40, 42, 51]. Using the stabilized formulation allows one to interpolate physical quantities with equal-order interpolations. This feature gives practitioners maximum flexibility in mesh generation and numerical implementation, and allows for low-order elements which are more robust than their higher-order counterparts. Equal-order interpolations always give an optimal constraint ratio [20, Chapter 4], which may be regarded as another appealing feature for incompressible elasticity. The stabilization term can be interpreted as a subgrid scale model within the variational multiscale framework [21, 23, 36, 40]. It has been observed that for inelastic models, the subgrid scale model requires careful design [51]. This observation partly motivates this work, in which we aim to design a stable numerical formulation for incompressible hyperelasticity that does not rely on subgrid scale numerical models with tunable parameters.

1.2 Overview of the proposed method

It is well-known that a finite element scheme is based on the formulation (i.e., the variational principle) and the discrete function spaces (i.e., the elements). Both components need to be properly accounted for in the design of numerical schemes. In this work, we introduce a mixed variational formulation different from the existing mixed $u/p$ formulation [48]. In that formulation, the momentum balance equations are coupled with an algebraic equation of state, which relates the pressure with $J$, the determinant of the deformation gradient [17, Chapter 8]. In the incompressible limit, this relation reduces to $J=1$. In the new mixed formulation, the momentum equations are coupled with the differential mass equation written in terms of the pressure primitive variable set. The volumetric behavior is reflected through the so-called isothermal compressibility factor [36]. In the incompressible limit, this term approaches zero, and the mass equation degenerates to the divergence-free constraint for the velocity field. Although $J=1$ is equivalent to the divergence-free constraint for the velocity field at the continuum level, they lead to different schemes at the discrete level. Based on the new mixed formulation, an a priori energy stability estimate can be obtained, and the inf-sup condition leads to a bound for the pressure solution. We regard these estimates as critical numerical properties embedded in the formulation that guarantee reliable results.

It should be pointed out that there are some existing formulations [15, 26, 38] that bear some similarity to ours, the key difference being that the Cauchy stress was expressed in a rate form in prior formulations. It is known that the rate constitutive equations are not built from free energies and cannot account for reversible elastic behavior [45]. Therefore, prior formulations cannot have an a priori energy stability estimate. Additionally, the rate constitutive equation requires special numerical considerations [25]. We aim to address these issues through the proposed formulation.

The choice of elements plays an equally critical role in numerical design for large-strain elasticity problems. Here, we attempt to provide a numerical technique that can be conveniently and robustly extended to the higher-order regime. The NURBS elements have been shown to enjoy superior robustness for large strain analysis [8, 33]. We adopt the same set of NURBS basis functions for the description of the geometry and approximation of the displacement field, aligning the proposed numerical formulation with the paradigm of isogeometric analysis [22]. The unique concept of $k$-refinement in isogeometric analysis allows one to generate higher-continuity basis functions without proliferation of degrees of freedom. However, it should be pointed out that in the setting of mixed finite elements, although the $k$-refinement leads to a pair of velocity-pressure elements that enjoy nearly the optimal constraint ratio [20, Chapter 4], it has been observed that such element types are not always inf-sup stable [41]. To remedy this issue, it has been proposed to use subdivision technology to generate a NURBS analogue for the Q1-iso-Q2 element [10, 28, 41]. In this work, we adopt an alternative approach, the inf-sup stable smooth generalizations of the Taylor-Hood element. In our opinion, the Taylor-Hood element is more convenient for implementation, especially in the parallel setting. We numerically assess the inf-sup stability for different combinations of the $p$- and $k$-refinements for generating the Taylor-Hood elements. It will be observed that the elements pass the numerical test if the polynomial degree is elevated at least once by the $p$-refinement to generate the discrete velocity space. Using the above new mixed formulation and the stable smooth generalizations of the Taylor-Hood element offer a new approach for incompressible large strain elastodynamics with several appealing features: it is well-behaved in the incompressible regime, the semi-discrete formulation respects energy stability, it does not involve tunable parameters or subgrid scale numerical models, it can achieve improved accuracy, especially for stress calculations, by employing higher-order smooth basis functions.

The remainder of the work is organized as follows. In Section 2, we state the governing equations and weak formulation for hyper-elastodynamics. In Section 3, the numerical scheme is presented and its numerical properties are analyzed. Following that, we numerically assess the inf-sup stability of different pairs of mixed NURBS elements. The elements that pass the test are used in the simulations for benchmark problems in Section 4. We draw conclusions in Section 5.

2 Hyper-elastodynamics

2.1 The initial boundary-value problem

Let $\Omega_{\bm{X}}$ and $\Omega_{\bm{x}}$ be bounded open sets in $\mathbb{R}^{d}$ with Lipschitz boundaries, wherein $d$ represents the number of spatial dimensions. The motion of the body is described by a family of diffeomorphisms parameterized by the time coordinate $t$,

[TABLE]

In the above, $\bm{x}$ represents the current position of a material particle originally located at $\bm{X}$, which implies $\bm{\varphi}(\bm{X},0)=\bm{X}$. The displacement and velocity of the material particle are defined as

[TABLE]

In this work, we use $d\left(\cdot\right)/dt$ to denote a total time derivative. The spatial velocity is defined as $\bm{v}:=\bm{V}\circ\bm{\varphi}_{t}^{-1}$. Analogously, we define $\bm{u}:=\bm{U}\circ\varphi_{t}^{-1}$. The deformation gradient, the Jacobian determinant, and the right Cauchy-Green tensor are defined as

[TABLE]

We define $\tilde{\bm{F}}$ and $\tilde{\bm{C}}$ as

[TABLE]

which represent the distortional parts of $\bm{F}$ and $\bm{C}$. We denote the thermodynamic pressure of the continuum body as $p$ and the density as $\rho$. The mechanical behavior of an elastic material can be described by a Gibbs free energy $G(\tilde{\bm{C}},p)$. It is shown in [36] that the Gibbs free energy can be additively split into an isochoric part and a volumetric part,

[TABLE]

The constitutive relations for the density $\rho$, the isothermal compressibility factor $\beta_{\theta}$, and the deviatoric part of the Cauchy stress can be described in terms of the Gibbs free energy as follows,

[TABLE]

wherein the projector $\mathbb{P}$ and the fictitious second Piola-Kirchhoff stress $\tilde{\bm{S}}$ are defined as

[TABLE]

$\mathbb{I}$ is the fourth-order identity tensor, and $\rho_{0}$ is the density in the referential configuration. Interested readers are referred to [36] for a detailed discussion of the governing equations and the constitutive relations. It is known that $\rho J=\rho_{0}$ due to mass conservation in the Lagrangian description. We can therefore introduce $\rho(J)=\rho_{0}/J$ as an alternative way of defining the density in the Lagrangian framework. In fact, we will adopt this choice in the following discussion. Under the isothermal condition, the energy equation is decoupled, and it suffices to consider the following equations for the motion of the continuum body,

[TABLE]

In the above system, the equations (1) describe the kinematic relation, and the equations (2) and (3) describe the conservation of mass and the balance of linear momentum. The boundary $\Gamma_{\bm{x}}^{t}=\partial\Omega^{t}_{\bm{x}}$ can be partitioned into two non-overlapping subdivisions: $\Gamma_{\bm{x}}^{t}=\Gamma_{\bm{x}}^{g,t}\cup\Gamma_{\bm{x}}^{h,t},$ wherein $\Gamma^{g,t}_{\bm{x}}$ is the Dirichlet part of the boundary, and $\Gamma^{h,t}_{\bm{x}}$ is the Neumann part of the boundary. Boundary conditions can be stated as

[TABLE]

Given the initial data $\bm{u}_{0}$, $p_{0}$, and $\bm{v}_{0}$, the initial conditions can be stated as

[TABLE]

The equations (1)-(5) constitute an initial-boundary value problem for elastodynamics.

Remark 1**.**

It is known that $J=1$ is equivalent to $\nabla_{\bm{x}}\cdot\bm{v}=0$ due to the identity $dJ/dt=J\nabla_{\bm{x}}\cdot\bm{v}$. However, the usage of $\nabla_{\bm{x}}\cdot\bm{v}=0$, or more generally (2), is uncommon in the literature. A reason is that the constraint $J=1$ is fitted into the elastostatic model, and the usage of $\bm{v}$ inevitably necessitates an elastodynamic model, which needs additional considerations in the numerical formulation. Another reason could be the missing link between $\beta_{\theta}$ and the strain energy. The constitutive relation for $\beta_{\theta}$ allows compressible materials and is recently derived in [36].

Since the above system looks different from the existing theory for hyperelasticity, we give an example of the constitutive model here. Let $I_{1}$ and $I_{2}$ designate the first and second invariants of the right Cauchy-Green tensor, that is,

[TABLE]

For isotropic materials, the isochoric part of the free energy can be conveniently expressed in terms of $\tilde{I}_{1}:=J^{-2/3}I_{1}$ and $\tilde{I}_{2}:=J^{-4/3}I_{2}$. The Mooney-Rivlin model can be expressed as

[TABLE]

where $c_{1}$ and $c_{2}$ are parameters that have the same dimension as pressure. The volumetric part of the Gibbs free energy can be built as a Legendre transformation of the Helmholtz volumetric free energy [36]. Here, we give an example

[TABLE]

which is transformed from the energy proposed in [34]. In (6), $\kappa$ designates the bulk modulus. This free energy leads to the relation

[TABLE]

As the bulk modulus $\kappa$ approaches infinity, the material becomes incompressible, and we have $G_{vol}(p)=p/\rho_{0}$ in the limit. This volumetric energy leads to $\rho(p)=\rho_{0}$ and $\beta_{\theta}(p)=0$.

2.2 Reduction to the small-strain theory

Assuming the strain is infinitesimally small, we have $\nabla_{\bm{x}}=\nabla_{\bm{X}}$ and $\rho(J)=\rho_{0}$. We also assume that $G_{vol}$ adopts the form given in (6). Then the mass equation (2) can be written as

[TABLE]

Integrating the above relation in time results in

[TABLE]

with a proper choice of the reference value for the pressure. Assuming further that the we are seeking a static equilibrium solution, the momentum equation (3) becomes

[TABLE]

The equations (8)-(9) constitute the classical mixed formulation for the small strain elastostatics [20, Chapter 4].

Remark 2**.**

For elastodynamics, one may instinctively add an inertial term to (9) and couple it with (8). However, numerical simulations indicate that this system is probably ill-posed. It is suggested to couple (9) with (7) rather than (8) for dynamic calculations [42]. A potential mathematical explanation is that (8) does not provide the proper coercive structure in the dynamic setting. This point will be further clarified in Proposition 1.

2.3 Weak formulation

Henceforth, we restrict our discussion to fully incompressible materials. Let us denote the trial solution spaces for the displacement, velocity, and pressure in the current domain as $\mathcal{S}_{\bm{u}}$, $\mathcal{S}_{\bm{v}}$, and $\mathcal{S}_{p}$, respectively. The Dirichlet boundary condition defined on $\Gamma^{g}_{\bm{x}}$ is properly built into the definitions of the $\mathcal{S}_{\bm{u}}$ and $\mathcal{S}_{\bm{v}}$. Let $\mathcal{V}_{\bm{v}}$ and $\mathcal{V}_{p}$ denote the corresponding test function spaces. The mixed formulation on the current configuration can be stated as follows. Find $\bm{y}(t):=\left\{\bm{u}(t),p(t),\bm{v}(t)\right\}^{T}\in\mathcal{S}_{\bm{u}}\times\mathcal{S}_{p}\times\mathcal{S}_{\bm{v}}$ such that for $t\in[0,T]$,

[TABLE]

for $\forall\left\{w_{p},\bm{w}_{\bm{v}}\right\}\in\mathcal{V}_{p}\times\mathcal{V}_{\bm{v}}$, with $\bm{y}(0)=\left\{\bm{u}_{0},p_{0},\bm{v}_{0}\right\}^{T}$. Here $\bm{u}_{0}$, $p_{0}$, and $\bm{v}_{0}$ are the $\mathcal{L}^{2}$ projections of the initial data onto the trial solution spaces. It is worth pointing out that although the material is fully incompressible, we still use $\rho(J)=\rho_{0}/J$ in (2.3), since the resulting discrete scheme cannot guarantee pointwise satisfaction of $J=1$. In the above and henceforth, the formulations for the kinematic equations, the mass equation, and the linear momentum equations are indicated by the superscripts $k$, $p$ and $m$, respectively. The equations (10)-(2.3) constitute the weak form of the problem. Performing integration by parts and using the localization argument, one can show the equivalence between the weak-form problem and the initial-boundary value problem. Let us define the following quantities on the material frame of reference via a pull-back operator:

[TABLE]

Correspondingly, the trial solution spaces are denoted as $\mathcal{S}_{\bm{U}}$, $\mathcal{S}_{P}$, and $\mathcal{S}_{\bm{V}}$; the test function spaces are denoted as $\mathcal{V}_{P}$ and $\mathcal{V}_{\bm{V}}$. The weak formulation can be alternatively stated as follows. Find $\bm{Y}(t):=\left\{\bm{U}(t),P(t),\bm{V}(t)\right\}^{T}\in\mathcal{S}_{\bm{U}}\times\mathcal{S}_{P}\times\mathcal{S}_{\bm{V}}$ such that for $t\in[0,T]$,

[TABLE]

for $\forall\left\{W_{P},\bm{W}_{\bm{V}}\right\}\in\mathcal{V}_{P}\times\mathcal{V}_{\bm{V}}$, with $\bm{Y}(0)=\left\{\bm{U}_{0},P_{0},\bm{V}_{0}\right\}^{T}$. Here $\bm{U}_{0}$, $P_{0}$, and $\bm{V}_{0}$ are the $\mathcal{L}^{2}$ projections of the initial data onto the spaces $\mathcal{S}_{\bm{U}}$, $\mathcal{S}_{P}$, and $\mathcal{S}_{\bm{V}}$ respectively.

3 Numerical formulation

In this section, we discuss the numerical procedures for the solution of the incompressible hyper-elastodynamics based on the weak formulation given in Section 2.3.

3.1 Spline spaces on the parametric domain

We start by reviewing the construction of B-splines and NURBS basis functions. Given the polynomial degree $\mathsf{p}$ and the dimensionality of the B-spline space $\mathsf{n}$, the knot vector can be represented by $\Xi:=\left\{\xi_{1},\cdots,\xi_{\mathsf{n}+\mathsf{p}+1}\right\}$, wherein $0=\xi_{1}\leq\xi_{2}\leq\cdots\leq\xi_{\mathsf{n}+\mathsf{p}+1}=1$. With the knot vector, the B-spline basis functions of degree $\mathsf{p}$, denoted as $\mathsf{N}_{i}^{\mathsf{p}}$ for $i=1,\cdots,\mathsf{n}$, can be defined recursively. The definition starts with the case of $\mathsf{p}=0$, in which the basis functions are defined as piecewise constants,

[TABLE]

For $\mathsf{p}\geq 1$, the basis functions are defined through the Cox-de Boor recursion formula,

[TABLE]

The NURBS basis functions of degree $\mathsf{p}$ are defined by the B-spine basis functions and a weight vector $\{\mathsf{w}_{1},\cdots,\mathsf{w}_{\mathsf{n}}\}$ as

[TABLE]

If we ignore the repetitive knots, the knot vector can be defined by a vector $\{\zeta_{1},\cdots,\zeta_{\mathsf{m}}\}$ representing the distinctive knots and a vector $\{r_{1},\cdots,r_{\mathsf{m}}\}$ recording the corresponding knot multiplicities. In this work, we consider open knot vectors, meaning $r_{1}=r_{\mathsf{m}}=\mathsf{p}+1$. We further assume that $r_{i}\leq\mathsf{p}$ for $i=2,\cdots,\mathsf{m}-1$. At the point $\zeta_{i}$, the B-spline basis functions have $\alpha_{i}:=\mathsf{p}-r_{i}$ continuous derivatives. The vector

[TABLE]

is referred to as the regularity vector. We adopt the notation

[TABLE]

When $\alpha_{i}$ takes the value $-1$, the basis functions are discontinuous at $\zeta_{i}$. The spaces $\mathcal{N}^{\mathsf{p}}_{\bm{\alpha}}$ and $\mathcal{R}^{\mathsf{p}}_{\bm{\alpha}}$ are defined as

[TABLE]

The notations $\mathcal{N}^{\mathsf{p}}_{\alpha}$ and $\mathcal{R}^{\mathsf{p}}_{\alpha}$ are used to indicate that $\alpha_{i}=\alpha$ for $i=2,\cdots,\mathsf{m}-1$, meaning the spline function spaces have continuity $C^{\alpha}$. The construction of multivariate B-spline and NURBS basis functions follows a tensor-product manner. Consider a unit cube $\hat{\Omega}:=(0,1)^{d}$, which is referred to as the parametric domain. Given $\mathsf{p}_{l}$, $\mathsf{n}_{l}$ for $l=1,\cdots,d$, we denote the knot vectors as $\Xi_{\mathsf{l}}=\{\xi_{1,l},\cdots,\xi_{\mathsf{n}_{l}+\mathsf{p}_{l}+1,l}\}$. Associated with each knot vector, the univariate B-spline basis functions $\mathsf{N}^{\mathsf{p}_{l}}_{i_{l},l}$ for $i_{l}=1,\cdots,\mathsf{n}_{l}$ are defined. Consequently, the tensor-product B-spline basis functions can be defined as

[TABLE]

Given the weight vectors $\{\mathsf{w}_{1,l},\cdots,\mathsf{w}_{\mathsf{n},l}\}$ for $l=1,\cdots,d$, the univariate NURBS basis functions $\mathsf{R}^{\mathsf{p}_{l}}_{i_{l},l}$ are defined. Correspondingly, the multivariate NURBS basis functions are defined as

[TABLE]

The tensor product NURBS space is denoted as

[TABLE]

3.2 Semi-discrete formulation and a priori estimates

In this work, we always consider three-dimensional problems (i.e. $d=3$). Two discrete function spaces $\hat{\mathcal{S}}_{h}$ and $\hat{\mathcal{P}}_{h}$ can be defined on $\hat{\Omega}=(0,1)^{3}$ as

[TABLE]

where $1\leq\mathsf{a}$ and $0\leq\mathsf{b}\leq\mathsf{a}$ are integers. We assume that the referential configuration of the body can be exactly parametrized by a geometrical mapping $\bm{\psi}:\hat{\Omega}\rightarrow\Omega_{\bm{X}}$. The discrete functions on the referential domain are defined through the pull-back operators,

[TABLE]

This pair of elements can be viewed as a generalization of the Taylor-Hood element [18], where the polynomial degree and the continuity can achieve arbitrarily high order. With the discrete function spaces $\mathcal{S}_{h}$ and $\mathcal{P}_{h}$ defined, we define the trial solution spaces for the displacement, pressure, and velocity on the referential configuration as

[TABLE]

Given the displacement $\bm{U}_{h}\in\mathcal{S}_{\bm{U}_{h}}$, one may obtain $\bm{\varphi}_{h}=\bm{U}_{h}(\bm{X},t)+\bm{X}$. Consequently, the trial solution spaces for the displacement, pressure, and velocity on the current configuration can be defined as

[TABLE]

and the test function spaces are defined as

[TABLE]

The semi-discrete formulation can be stated as follows. Find $\bm{y}_{h}(t):=\left\{\bm{u}_{h}(t),p_{h}(t),\bm{v}_{h}(t)\right\}^{T}\in\mathcal{S}_{\bm{u}_{h}}\times\mathcal{S}_{p_{h}}\times\mathcal{S}_{\bm{v}_{h}}$ such that for $t\in[0,T]$,

[TABLE]

for $\forall\left\{w_{p_{h}},\bm{w}_{\bm{v}_{h}}\right\}\in\mathcal{V}_{p_{h}}\times\mathcal{V}_{\bm{v}_{h}}$, with $\bm{y}_{h}(0):=\left\{\bm{u}_{h0},p_{h0},\bm{v}_{h0}\right\}^{T}$. Here $\bm{u}_{h0}$, $p_{h0}$, and $\bm{v}_{h0}$ are the $\mathcal{L}^{2}$ projections of the initial data onto the finite dimensional trial solution spaces. In the following, we demonstrate that the above semi-discrete formulation is embedded with energy stability and momentum conservation properties. The properties guarantee that the numerical solutions preserve critical structures of the original system. In contrast, to the best of the authors’ knowledge, there is no such stability estimate for the conventional mixed $u/p$ formulation [48] or the formulations based on rate constitutive equations [15, 26, 38].

Proposition 1** (A priori energy stability estimate).**

For fully incompressible materials, assuming the boundary data $\bm{g}$ is time independent, we have

[TABLE]

Proof.

Since the Dirichlet boundary data $\bm{g}$ is independent of time, one is allowed to choose $w_{p_{h}}=p_{h}$ in (17) and $\bm{w}_{\bm{u}_{h}}=\bm{v}_{h}$ in (3.2), and this leads to the following,

[TABLE]

Rearranging terms in the above equality leads to

[TABLE]

∎

Remark 3**.**

For compressible materials, one may analogously obtain a stability bound where a pressure-squared term enters into the integral on the left-hand side of (19). This gives a mathematical reason for the success of equal-order interpolations when the material is compressible. However, we do not favor this type of ‘energy’ estimates because the pressure-squared term does not carry physical meanings. To remedy this issue, an entropy variable can be introduced by leveraging the convexity of the volumetric energy, and a physically relevant entropy stability is expected [35, 43]. This is beyond the scope of this work and remains an area of future research.

Proposition 2** (Semi-discrete momentum conservation).**

Considering the pure Neumann boundary condition, we have the following conservation properties of the semi-discrete formulation,

[TABLE]

Proof.

The above conservation properties are direct consequences of choosing $\bm{w}_{\bm{v}_{h}}=\bm{e}_{i}$ and $\bm{w}_{\bm{v}_{h}}=\bm{e}_{i}\times\bm{\varphi}_{h}$ respectively in (3.2), where $\bm{e}_{i}$ is a unit vector in the $i$-th direction. ∎

Due to the incompressibility, the pressure force does not contribute to the energy. Therefore, the energy stability estimate (19) does not involve the pressure field. The inf-sup condition needs to be utilized to provide a bound for the pressure field. We assume that there exists a positive constant $\beta$ such that

[TABLE]

wherein $\|\cdot\|_{0}$ and $\|\cdot\|_{1}$ denote the $\mathcal{L}^{2}$ and $\mathcal{H}^{1}$ norm over $\Omega_{\bm{x}}$. Using the semi-discrete equation (3.2), the above inequality implies

[TABLE]

If we further assume that $\rho(J_{h})$ is uniformly bounded, using the Cauchy-Schwarz inequality, we may get

[TABLE]

with $\tilde{C}$ being a constant. Therefore, given the velocity, the deformation state, and the external forces, the pressure field is bounded. We note that the assumption on the boundedness of the density cannot be rigorously justified based on the current numerical formulation. It is anticipated that this issue can be resolved by invoking the structure-preserving discretization technique [12], which results in discrete solutions with pointwise divergence-free velocity field. With the exact satisfaction of the incompressibility constraint, the density remains as a constant.

Remark 4**.**

The linearization of $J-1=0$ results in a divergence operator acting on the virtual displacement field. This fact has been frequently used to justify the usage of inf-sup stable elements in the two-field variational principle [3]. However, we feel this may not be a good interpretation. First, the linearization argument cannot recover the compressible case (8). Second, the solvability of the Newton-Raphson procedure does not provide a bound for the solution.

3.3 Temporal discretization

We invoke the generalized-$\alpha$ method [27] for the temporal discretization of the weak form problem (10)-(2.3). The time interval $[0,T]$ is divided into a set of $n_{ts}$ subintervals of size $\Delta t_{n}:=t_{n+1}-t_{n}$ delimited by a discrete time vector $\left\{t_{n}\right\}_{n=0}^{n_{ts}}$. The solution vector and its first-order time derivative evaluated at the time step $t_{n}$ are denoted as $\bm{y}_{n}$ and $\dot{\bm{y}}_{n}$. The fully discrete scheme can be stated as follows. At time step $t_{n}$, given $\dot{\bm{y}}_{n}$, $\bm{y}_{n}$, the time step size $\Delta t_{n}$, and the parameters $\alpha_{m}$, $\alpha_{f}$, and $\gamma$, find $\dot{\bm{y}}_{n+1}$ and $\bm{y}_{n+1}$ such that for $\forall\left\{w_{p},\bm{w}_{\bm{v}}\right\}\in\mathcal{V}_{p}\times\mathcal{V}_{\bm{v}}$,

[TABLE]

The choice of parameters $\alpha_{m}$, $\alpha_{f}$ and $\gamma$ determines the accuracy and stability of the temporal scheme. Importantly, the high-frequency dissipation can be controlled via a proper parametrization of these parameters, while maintaining second-order accuracy and unconditional stability (for linear problems). For the above first-order dynamic problems, the parametrization is

[TABLE]

wherein $\varrho_{\infty}\in[0,1]$ denotes the spectral radius of the amplification matrix at the highest mode [27]. Setting $\varrho_{\infty}=1$ recovers the mid-point rule. For nonlinear structural dynamics, the mid-point rule is observed to have a pile-up effect for the energy error and often leads to diverged results for long-time simulations. In this study, the value of $\varrho_{\infty}$ is fixed to be $0.5$.

Remark 5**.**

Interested readers are referred to [7] for the parametrization of $\alpha_{m}$, $\alpha_{f}$, and $\gamma$ for second-order structural dynamics. A recent study shows that using the generalized-$\alpha$ method for first-order structural dynamics enjoys improved dissipation and dispersion properties and does not suffer from overshoot [29]. Moreover, using a first-order structural dynamic model is quite propitious for the design of an FSI scheme [36].

Remark 6**.**

It is tempting to apply the discrete energy-momentum methods [47] to the semi-discrete system. Those algorithms yield fully discrete systems that inherit the energy stability and momentum conservation properties and are thence particularly well-suited for transient analysis. For problems we are interested in, the solution may be driven to a static equilibrium by external forces, and the stress formula in the energy-momentum methods will become ill-defined. Because of this, we retain the generalized-$\alpha$ method in this work.

3.4 A Segregated predictor multi-corrector algorithm

The equations (21)-(26) constitute a system of nonlinear algebraic equations to be solved in each time step, and we invoke the Newton-Raphson method with consistent linearization. At time step $t_{n+1}$, the solution vector $\bm{y}_{n+1}$ is solved by means of a predictor multi-corrector algorithm. We denote $\bm{y}_{n+1,(l)}:=\left\{\bm{u}_{n+1,(l)},p_{n+1,(l)},\bm{v}_{n+1,(l)}\right\}^{T}$ as the solution vector at the Newton-Raphson iteration step $l=0,\cdots,l_{max}$. The residual vectors evaluated at the iteration stage $l$ are denoted as

[TABLE]

The consistent tangent matrix associated with the above residual vectors is

[TABLE]

wherein

[TABLE]

$\bm{\mathrm{I}}$ is the identity matrix, and $\bm{\mathrm{O}}$ is the zero matrix. The above diagonal structure of the two blocks can be utilized to construct a block factorization of $\bm{\mathrm{K}}_{(l)}$, with which the solution procedure of the linear system of equations in the Newton-Raphson method can be consistently reduced to a two-stage algorithm [36, 42]. In the first stage, one obtains the increments of the pressure and velocity at the iteration step $l$ by solving the following linear system,

[TABLE]

In the second stage, one obtains the increments for the displacement by

[TABLE]

To simplify notations in the following discussion, we denote

[TABLE]

Readers are referred to the Appendix of [37] for the explicit formulas of the block matrices in (29).

Remark 7**.**

In [36], it was shown that $\bm{\mathrm{R}}^{k}_{(l)}=\bm{0}$ for $l\geq 2$ for general predictor multi-corrector algorithms; in [40], a special predictor is chosen so that $\bm{\mathrm{R}}^{k}_{(l)}=\bm{0}$ for $l\geq 1$. In our experience, setting $\bm{\mathrm{R}}^{k}_{(l)}=\bm{0}$ for $l\geq 1$, regardless of the predictor chosen, simplifies the implementation and does not deteriorate the convergence rate of the Newton-Raphson solution procedure.

Based on the above discussion, a predictor multi-corrector algorithm for solving the nonlinear algebraic equations in each time step can be summarized as follows.

Predictor stage: Set:

[TABLE]

Multi-corrector stage: Repeat the following steps for $l=1,\dots,l_{max}$:

Evaluate the solution vectors at the intermediate stages:

[TABLE] 2. 2.

Assemble the residual vectors $\bm{\mathrm{R}}^{m}_{(l)}$ and $\bm{\mathrm{R}}^{p}_{(l)}$ using $\bm{y}_{n+\alpha_{f},(l)}$ and $\dot{\bm{y}}_{n+\alpha_{m},(l)}$. 3. 3.

Let $\|\bm{\mathrm{R}}_{(l)}\|_{\mathfrak{l}^{2}}$ denote the $\mathfrak{l}^{2}$-norm of the residual vector. If either one of the following stopping criteria

[TABLE]

is satisfied for two prescribed tolerances $\textup{tol}_{\textup{R}}$, $\textup{tol}_{\textup{A}}$, set the solution vector at the time step $t_{n+1}$ as $\bm{y}_{n+1}=\bm{y}_{n+1,(l-1)}$ and $\dot{\bm{y}}_{n+1}=\dot{\bm{y}}_{n+1,(l-1)}$, and exit the multi-corrector stage; otherwise, continue to step 4. 4. 4.

Assemble the tangent matrices (29). 5. 5.

Solve the following linear system of equations for $\Delta\dot{p}_{n+1,(l)}$ and $\Delta\dot{\bm{v}}_{n+1,(l)}$,

[TABLE] 6. 6.

Obtain $\Delta\dot{\bm{u}}_{n+1,(l)}$ from the relation (28). 7. 7.

Update the solution vector as

[TABLE]

and return to step 1.

In our experience, the choice of the linear solver for (30) critically impacts the overall numerical efficiency and robustness, especially for three-dimensional problems. Linear solvers based on algebraic factorizations (such as incomplete LU) are prone to fail due to the appearance of a zero sub-matrix $\bm{\mathrm{O}}$ in (30), which may lead to zero-pivoting. Hence, an iterative solution procedure for (30) is specifically designed based on a nested block preconditioning technique. Readers are referred to [37] for more details.

4 Numerical results

In this section, we perform numerical investigations using the proposed scheme. Unless otherwise specified, we use $\mathsf{p}+\mathsf{a}+1$ Gauss quadrature points in each direction; the pressure function space is generated by the $k$-refinement to achieve the highest possible continuity.

4.1 Numerical Inf-Sup test

The inf-sup condition for the discrete problem states that there exists a constant $\beta$ independent of the mesh size such that

[TABLE]

We examine the inf-sup condition for the proposed discrete spaces $\mathcal{S}_{\bm{v}_{h}}$ and $\mathcal{S}_{p_{h}}$ using the numerical inf-sup test [4]. Let $N_{A}$ and $M_{\tilde{A}}$ denote the velocity and pressure basis functions on the current configuration where $A$ and $\tilde{A}$ are the node number. The following matrices are defined.

[TABLE]

We consider the following eigenvalue problem: Find $\gamma^{h}_{i}$ and $\bm{\psi}_{i}$ such that

[TABLE]

The value of $\beta^{h}$ is determined as the square root of the smallest non-zero eigenvalue. The regularity vector $\bm{\alpha}=\{-1,\alpha,\cdots,\alpha,-1\}$ is the same in all three directions. The numerical integration is performed by the Gauss quadrature rule with $\mathsf{p}+\mathsf{a}+2$ quadrature points in each direction to ensure accuracy. The eigenvalues are calculated by the SLEPc package [14]. The trend of $\beta^{h}$ is examined as we progressively refine the mesh for $0\leq\mathsf{b}\leq\mathsf{a}$. We consider a curved geometry for the domain, which is exactly represented by NURBS and illustrated in Figure 1. The computed values of $\beta^{h}$ for $\mathsf{p}=2$, $3$, and $4$ with $0\leq\mathsf{b}\leq\mathsf{a}\leq 2$ are presented in Figure 2. It can be observed that $\beta^{h}$ approaches zero with mesh refinement when $\mathsf{a}=\mathsf{b}$. To confirm this observation, we investigate the cases of $\mathsf{a}=3$ and $\mathsf{a}=4$ with $\mathsf{p}$ fixed to be $2$, with results reported in Figure 3. Again, we observe that $\beta^{h}$ shows a clear trend of approaching zero with mesh refinement only when $\mathsf{a}=\mathsf{b}$. To further validate this finding, we also study a unit cube for the domain, which allows us to start the test with $\mathsf{p}=1$. Again, the same trend of $\beta^{h}$ is observed. Based on the collected results, we make the following salient observations. For the smooth generalizations of the Taylor-Hood element, if the velocity space is generated by pure $k$-refinement (i.e., $\mathsf{a}=\mathsf{b}$), the resulting element pair is not inf-sup stable. If the velocity space is generated by pure $p$-refinement from the pressure space (i.e., $\mathsf{b}=0$), the smallest eigenvalues are bounded below from zero. Also, if $\mathsf{a}\geq 2$, the velocity spaces generated with $1\leq\mathsf{a}-\mathsf{b}$ also pass the numerical inf-sup test. This suggests that one may still perform $k$-refinement to increase the regularity of the velocity space if it is followed by a $p$-refinement of order at least one.

4.2 Convergence studies

In this example, we investigate the convergence behavior of the proposed numerical scheme. We consider an incompressible Neo-Hookean material model

[TABLE]

The geometrical domain is a unit cube with dimension 1m $\times$ 1m $\times$ 1m. The modulus $c_{1}$ is chosen as $1$ Pa, and the density $\rho_{0}$ is 1 kg/m3. The analytic forms of the displacement and pressure fields on the referential configuration adopt the following forms,

[TABLE]

In this example, the reference values are chosen as $L_{0}=1$ m, $M_{0}=1$ kg, $T_{0}=1$ s; both $\beta$ and $\gamma$ are chosen to be $2\pi$ rad; $c$ and $d$ are non-dimensional parameters that take the value 0.2. On the faces $Y=Z=0$ m and $Y=Z=1$ m, the body is fully clamped, and traction boundary conditions are applied on the rest faces. For the simulations, we use $\textup{tol}_{\textup{R}}=10^{-10}$ and $\textup{tol}_{\textup{A}}=10^{-12}$ as the stopping criteria in the predictor multi-corrector algorithm. Two different time step sizes are used to ensure that the temporal error does not pollute the spatial convergence rate. The relative errors of the displacement and pressure fields are reported in Figure 4 for varying values of $\mathsf{p}$ with $\mathsf{a}=1$, $\mathsf{b}=0$. We notice immediately that all the errors decrease with the optimal rates. In Figure 5, we report the convergence rates for $\mathsf{a}=2$, which resembles a smooth generalization of the spectral element [50]. In Figure 5 (a), we note that the increase of the value of $\mathsf{a}$ does not improve the convergence rate, regardless of the value of $\mathsf{b}$. Yet, the velocity error is smaller than that of the $\mathsf{a}=1$ case. From Figure 5 (b), we can see that the pressure errors are almost indistinguishable for $\mathsf{a}=1$ and $\mathsf{a}=2$.

4.3 Three-dimensional compression of a block

In this example, we examine the performance of the new formulation using the benchmark problem initially designed in [39]. On the boundary faces $X=Y=Z=0$, we apply symmetry boundary conditions, and we disallow horizontal displacement on the top surface. A ‘dead’ load with magnitude $3.2\times 10^{8}$ Pa is applied on a quarter portion of the top surface, which assumes the negative Z-direction in the referential configuration. The block is initially stress free with zero displacement. The surface traction load is applied as a linear function of time and reaches the prescribed magnitude at time $T=1$ s. We adopt an incompressible Neo-Hookean model given by the following energy function,

[TABLE]

The material properties are chosen as $\rho_{0}=1.0\times 10^{3}$ kg/m3 and $c_{1}=8.0194\times 10^{7}$ Pa. We simulate the problem with a fixed time step size $\Delta t=5.0\times 10^{-3}$ s. We fix the values of $\mathsf{a}$ and $\mathsf{b}$ to $1$ and [math] in this example and choose $\textup{tol}_{\textup{R}}=10^{-3}$ and $\textup{tol}_{\textup{A}}=10^{-6}$ as the stopping criteria. For comparison purposes, we also simulate the problem with the variational multiscale (VMS) formulation [36] using equal-order interpolations. As a classical benchmark problem, the primary quantity of interest is the displacement at the upper center point (i.e. the point at $X=Y=0$, $Z=1$ in the reference configuration). In Figure 6 (b), the compression levels at this point calculated by different methods are illustrated. For the coarsest mesh (two elements per side), the stable element with $\mathsf{p}=2$ gives a very good prediction of the compression level. It is interesting to note that the equal-order interpolation using the $Q_{1}/Q_{1}$ element with the VMS formulation gives a fairly good result for a finer mesh with four elements per side. Using the same mesh, the stable elements with $\mathsf{p}=1$ and $\mathsf{p}=2$ gives slightly softer predictions, which is due to oscillations of the higher-order methods at the tip. Using a mesh with eight elements per side, both stable elements give indistinguishable results in comparison with the reference value. In Figure 7, we further compare the pressure profiles at the current configuration calculated by a coarse mesh (two elements per side) with the value of $\mathsf{p}$ varying from $2$ to $6$. The pressure profile calculated by a fine mesh (16 elements per side and $\mathsf{p}=2$) is depicted to serve as a reference solution profile. It can be observed that the increase of the polynomial degree $\mathsf{p}$ improves of the solution quality. For the case of $\mathsf{p}=6$, the calculated result essentially captures the major feature of the pressure field.

4.4 Tensile test of an anisotropic fiber-reinforced hyperelastic soft tissue specimen

In this example, we examine the performance of the proposed formulation for an anisotropic hyperelastic material, which has been designed to describe arterial tissue layers with distributed collagen fibers [13]. The geometry set-up and the material model are summarized in Table 1. The groundmatrix is modeled as an isotropic Neo-Hookean material, with $c_{1}$ being the shear modulus. The $i$th family of collagen fibers is modeled by an exponential function $G^{f_{i}}_{ich}$. The unit vector $\bm{a}_{i}$ characterizes the mean orientation of the fiber, and $\kappa_{d}$ is a dispersion parameter that characterizes the distribution of the collagen fibers. In this study, we assume the mean orientation of the two families of fibers has no component in the radial direction and is completely determined by $\varphi$, the angle between the fiber orientation and the loading direction. For a circumferential specimen, the tensile load is along the circumferential direction and $\varphi=49.98^{\circ}$; correspondingly, for an axial specimen, the value of $\varphi$ is $40.02^{\circ}$. We consider only one-eighth of the specimen by applying symmetry boundary conditions. On the loading surface, a master-slave relation is enforced for the nodes to ensure that the surface moves only in the loading direction. The loading traction is applied gradually and reaches $2$ N in 200 seconds. We simulate the problem with a fixed time step size $\Delta t=2.0\times 10^{-2}$ s. Again, we fix the value of $\mathsf{a}$ and $\mathsf{b}$ to be $1$ and [math] and use $\textup{tol}_{\textup{R}}=10^{-3}$ and $\textup{tol}_{\textup{A}}=10^{-6}$ as the stopping criteria. Three different meshes are used for the proposed formulation: mesh 1 consists of $61440$ elements with $\mathsf{p}=2$, mesh 2 consists of $120$ elements with $\mathsf{p}=1$, and mesh 3 consists of $120$ elements with $\mathsf{p}=2$. In Figure 8, the load-displacement curves calculated by the three different meshes for the circumferential and axial specimen are plotted. It is hard to distinguish the results in Figure 8 (a). In Figure 8 (b), we provide a detailed comparison near the tensile load $0.35$ N. The curve obtained from mesh 3 is still very close to the reference solid line, indicating improved accuracy with increasing polynomial degree. For comparison purposes, we present the stress results calculated by the VMS formulation [36] with linear tetrahedral elements using two different spatial resolutions (see Table 2). From Figures 9 and 10, we observe that the essential feature of the Cauchy stress is captured in mesh 2, although there are slight oscillations near the corners. The results calculated from the mesh 1 and mesh 3 are almost indistinguishable, indicating that increasing the polynomial degree improves the accuracy of the stress results. In contrast, the stress is poorly resolved in mesh 4 due to the low-order elements. The results of mesh 5 illustrate that mesh refinement helps improve the quality of the stress results. Yet, one can still observe a discontinuous pattern and oscillations of the stress profile.

4.5 Three-dimensional beam bending

In this example, we present a three-dimensional beam vibration problem to evaluate the performance of the elastodynamics formulation in a bending dominated scenario [5]. The problem configuration as well as the material properties are illustrated in Table 3. The beam is fully clamped at the base, and the other faces are specified by zero tractions. The body is initially stress free with zero displacement. The vibration is initiated through an initial velocity

[TABLE]

This initial condition leads to an oscillatory motion of the beam. For the simulations, we choose $\mathsf{p}=1$, $\mathsf{a}=1$, and $\mathsf{b}=0$ for the discrete function spaces. We use $\textup{tol}_{\textup{R}}=10^{-8}$ and $\textup{tol}_{\textup{A}}=10^{-8}$ as the stopping criteria. The numerical results show the deformation state of the beam calculated from the two different meshes are indistinguishable, suggesting a coarse mesh with $\Delta x=L_{0}/2$ is capable of accurately describing the beam dynamics (Figure 11). Since the boundary data is time independent and the body force and surface tractions are zero, the total energy of the beam is conserved according to Proposition 1. We observe that the total energy is well-preserved up to $T=10$ s (Figure 12 (a)). From the periodic pattern of the kinetic and potential energies, we obtain an average period of the oscillation is 0.9018 s. To better illustrate the energy conservation, we plot the relative errors of the energy in Figure 12 (b), using three different spatial meshes. Interestingly, the error of the total energy achieves its maximum value when the beam reaches its largest deformation. For the coarsest mesh ($\Delta x=L_{0}/2$), the error accumulates slightly over time, and we can see that the relative error reaches about one percent at around $9.5$ s. We also observe that the spatial mesh refinement helps reduce the error of the total energy. For the meshes with $\Delta x=L_{0}/4$ and $\Delta x=L_{0}/6$, we do not observe a pile-up effect of the energy error. Also, the magnitude of the relative error is reduced with mesh refinement. In comparison with the previously published results [2, 32], the new formulation enjoys a better discrete energy conservation property.

4.6 A spinning annular disk

In this example, we study a spinning annular disk with zero traction boundary condition imposed on all boundary faces. In Figure 13 (a), the geometrical setting is illustrated. The inner radius of the disk is $0.5$ m, the outer radius is $1.5$ m, and the thickness is $1$ m. The material of the disk is Neo-Hookean with density $\rho_{0}=10$ kg/m3 and shear modulus $c_{1}=7.5$ Pa. Both the geometrical and material settings follow the benchmark example in [31]. The initial displacement is zero, and the spinning motion is initiated by an initial angular velocity of $1$ rad/s in the x-y plane, that is

[TABLE]

We choose the reference scales as $L_{0}=1$ m, $M_{0}=1$ kg, and $T_{0}=1$ s. The geometry of the domain can be exactly parametrized by quadratic NURBS, and hence we choose $\mathsf{p}=2$ for the discrete function spaces with $\mathsf{a}=1$ and $\mathsf{b}=0$. A coarse mesh is generated with 32 elements in the circumferential direction, 4 elements in the radial direction, and 4 elements in the axial direction; a fine mesh is generated with 64 elements in the circumferential direction, 8 elements in the radial direction, and 8 elements in the axial direction. The time step size is $\Delta t=2\times 10^{-4}$ s, and the problem in integrated up to $T=10.0$ s. For the simulations, we use $\textup{tol}_{\textup{R}}=10^{-8}$ and $\textup{tol}_{\textup{A}}=10^{-8}$ as the stopping criteria. In Figure 13 (b), a snapshot of the simulated velocity in the annular disk is depicted. Due to the zero traction boundary condition and the zero body force, this problem serves as a benchmark for examining the energy stability as well as the momentum conservation properties. In Figure 14 (a), we can see that the kinetic energy and the total energy are nicely conserved. In Figure 14 (b), the relative errors of the total energy over time are plotted, which are uniformly smaller than $3\times 10^{-6}$. The exact value of the linear momentum is zero, and we see that the absolute errors are less than $1.5\times 10^{-13}$ in Figure 14 (c). The x- and y-components of the angular momentum are zero, with numerical values having absolute errors less than $10^{-13}$ (Figure 15). The analytic value of the z-component of the angular momentum is $78.5$ kg$\cdot$m2/s, and we depict its relative error from the simulation with the coarse mesh. Note that the error of the z-component of the angular momentum is highly oscillatory and is bounded by $8\times 10^{-9}$. The numerical results corroborate the estimates given in Section 3.2.

5 Conclusions and future work

In this work, we presented a new numerical formulation for incompressible hyper-elastodynamics. We have revealed that the proposed formulation possesses a physically compatible notion of numerical stability, and the inf-sup condition can be utilized to give a bound for the pressure. These properties favorably distinguish the proposed formulation from previously existing ones [15, 26, 38, 48]. We use smooth generalizations of the Taylor-Hood element based on NURBS for the spatial discretization, aiming to provide a higher-order method that is stable, robust, and implementationally convenient. The inf-sup stability for the elements is elucidated through numerical assessment. A variety of benchmark examples are simulated to investigate the effectiveness of the method in different loading conditions and for different material models. In particular, two dynamic problems are studied to verify the numerical stability and conservation properties.

In addition to the superior accuracy in stress calculations, the adoption of NURBS elements makes the description of material anisotropy convenient because the mesh naturally aligns along the axial, circumferential, and radial directions. These attributes make the proposed formulation a promising candidate for biomedical problems. Based on the proposed formulation, the anisotropic arterial wall model will be further refined with detailed stress-driven mass production and removal for individual constituents that comprise the tissue. This will lead to a three-dimensional patient-specific predictive tool for vascular growth and remodeling.

Acknowledgements

This work is supported by the National Institutes of Health (NIH) under the award numbers 1R01HL121754 and 1R01HL123689, the National Science Foundation (NSF) CAREER award OCI-1150184, and computational resources from the Extreme Science and Engineering Discovery Environment supported by the NSF grant ACI-1053575. The authors acknowledge TACC at the University of Texas at Austin for providing computing resources that have contributed to the research results reported within this paper.