A stabilized linear finite element method for anisotropic poroelastodynamics with application to cardiac perfusion

We propose a variational multiscale method stabilization of a linear finite element method for nonlinear poroelasticity. Our approach is suitable for the implicit time integration of poroelastic formulations in which the solid skeleton is anisotropic and incompressible. A detailed numerical methodology is presented for a monolithic formulation that includes both structural dynamics and Darcy flow. Our implementation of this methodology is verified using several hyperelastic and poroelastic benchmark cases, and excellent agreement is obtained with the literature. Grid convergence studies for both anisotropic hyperelastodynamics and poroelastodynamics demonstrate that the method is second-order accurate. The capabilities of our approach are demonstrated using a model of the left ventricle (LV) of the heart derived from human imaging data. Simulations using this model indicate that the anisotropicity of the myocardium has a substantial influence on the pore pressure. Furthermore, the temporal variations of the various components of the pore pressure (hydrostatic pressure and pressure resulting from changes in the volume of the pore fluid) are correlated with the variation of the added mass and dynamics of the LV, with maximum pore pressure being obtained at peak systole. The order of magnitude and the temporal variation of the pore pressure are in good agreement with the literature.


Introduction
Cardiac perfusion, a mechanism by which blood is removed and delivered to the myocardium, plays a significant role in heart function. Cardiovascular magnetic resonance imaging is the most commonly used experimental technique for acquiring perfusion data [1,2]. Computational methods such as the finite element method (FEM) provide an alternative to medical imaging and promise to provide insight into the mechanism of perfusion. Computational methods for modelling perfusion require a poroelastodynamic framework that describes the unsteady coupled interaction between the elastodynamics of a solid and fluid dynamics within the pores of the solid. Poroelasticity is also used in other areas of biomechanics, including in models of cartilage, liver, and cornea, and it also has been widely applied in geomechanics and hydrogeology [3][4][5]. The complete modelling of poroelastodynamics can be considered within a fluid-structure interaction framework that resolves each pore. However, it is difficult to obtain accurate information about the pores, such as the shape and connectivity between them, and simulations of models with this level of detail are generally computationally intractable. Thus, homogenized models are commonly employed [6][7][8][9][10]. In these reduced models, homogenization converts the small-scale phenomenon into macroscopic quantities by considering the poroelastic medium to be a homogenized mixture of solid and fluid. In such models, the solid part is often called the skeleton [11,12]. For cardiac perfusion, the myocardium is considered a homogeneous mixture of a solid skeleton compartment consisting of cardiac myocytes and collagen and a fluid compartment composed of coronary vessels to transport blood to and from the myocardium [11][12][13]. In homogenization, the length scale of the pores is characterized by capillary permeability and porosity [14]. Using the poroelastic approach, the Darcy flow assumption provides a continuous flow velocity and pressure in the pores [15]. Assuming Darcy flow also simplifies the model by considering the velocity as the pressure gradient multiplied by the permeability tensor [16].
Simulating poroelastodynamics in complex problems such as cardiac perfusion requires methods that can handle the incompressibility or near incompressibility of the anisotropic, nonlinear skeleton together with Darcy flow. In large deformation models, the need to satisfy the Ladyzhenskaya-Babuska-Brezzi (LBB) criterion [17] makes it challenging to treat the skeleton as either incompressible or nearly incompressible. Costa et al. [18] were among the first to perform a complete three-dimensional finite element analysis of the left ventricle (LV) of the heart with large deformation and a nonlinear and non-homogeneous incompressible model. They concluded that higher-order finite elements are required for a stable, accurate solution. Yang et al. [19] conducted an FEM analysis of poroviscoelastic soft tissue using biquadratic interpolation for the displacement of the skeleton and fluid pressure. They concluded that the same order of interpolation for displacement and pressure would not satisfy the LBB criterion. Chapelle et al. [11] studied the poroelastodynamics in an elliptical LV with near incompressibility of the skeleton using FEM with firstorder elements for the displacement and piecewise constant elements for the pressure. Richardson et al. [13] studied the poroelastodynamics in a realistic human LV using an Immersed Boundary/Finite Element method [20] that used first-order elements for both the displacement and velocity of the skeleton and the pore pressure together with stabilization using a volumetric energy term. Cookson et al. [21] used a mixed method with second-order elements for the displacement and first-order elements for the pressure for a comparably realistic LV. Lee et al. [22] used a similar approach for the analysis of poroelastodynamics in a porcine LV.
It is well known that in the absence of stabilization, using first-order elements for both the displacement and the pressure in incompressible and nearly incompressible material models leads to spurious pressure checkerboard modes and volumetric locking [23]. Prior work yielded a stable method for incompressible elasticity and poroelasticity through higherorder elements and reduced integration methods [24]. Because higher-order elements are computationally expensive, another choice is to use first-order elements for the displacement and piecewise constant elements for the pressure. However, such methods give only firstorder accuracy for the pressure and are limited to near incompressibility. Obtaining secondorder accuracy for both displacement and pressure requires at least linear elements for both displacement and pressure. Many stabilization techniques have been developed to enable the use of equal-order elements, such as formulations that use a modified deformation gradient [23] or mean dilation [25]. Developments in the finite element method for fluid mechanics led to a stabilized method for incompressible fluid with first-order elements for velocity and pressure [26,27]. Later, this method was extended to incompressible static elasticity [28]. Klaas et al. [29] presented the stabilized finite element method for incompressible and nearly incompressible static nonlinear hyperelasticity using linear elements for both displacement and pressure. Masud and Truster [30] used a variational multiscale framework for the stabilization using linear elements. Scovazzi et al. [31] and Zeng et al. [32] extended this method to unsteady formulations by adding a stabilization term corresponding to the inertia force using the Dynamic Variational Multiscale (D-VMS) method. Rossi et al. [33] presented a detailed formulation of the stabilized finite element method for unsteady hyperelasticity using an implicit time-stepping method.
This study extends the D-VMS method for large deformation hyperelasticity and isotropic materials presented in Rossi et al. [33] to a poroelastic anisotropic formulation suitable for complex biomechanical models, including the LV. The stabilization term in our formulation includes contributions from the anisotropic elastodynamics and poroelastodynamics. Grid convergence studies for anisotropic elastodynamics and poroelastodynamics test cases are presented. Further verification studies are reported for three-dimensional problems of shrinking and expanding cubes [11]. The method's capabilities are demonstrated by simulating LV dynamics and perfusion using an anatomical model derived from human image data across three consecutive cardiac cycles.

Mathematical formulation
For the poroelastic medium, the homogeneous mixture formulation represents both solid and fluid as a single continuum with porosity ϕ, which represents the fraction of fluid volume in the porous medium, such that the volume occupied by the fluid at time t is V f (t) = ∫ Ω t ϕ(x, t)dx, (1) in which Ω t is the physical region occupied by the poroelastic medium at time t. The volume of the skeleton is the difference between the total volume of Ω t and the fluid volume V f (t), The deformation gradient tensor F of the homogeneous mixture is computed with respect to the skeleton configuration and is defined by in which X are material coordinates of the skeleton in the reference configuration and x = χ (X, t) are the current coordinates of the material point X at time t. The Lagrangian displacement and velocity of the material are u(X, t) = x − X and v(X, t) = u(X, t). We consider the balance equations for the poroelastodynamics problem in Lagrangian form.
In the next section, we start by presenting the mathematical statement of the elastodynamics problem without the presence of the pore fluid (ϕ ≡ 0). Then, we extend the formulation to poroelastodynamics using the mass and momentum balance equations in the pore-fluid, coupled with the mass and momentum balance equations for the homogeneous mixture.

Elastodynamics
The dynamics of a hyperelastic solid are fully described by the equations of mass and momentum balance. Because there is no pore fluid, the homogeneous mixture is comprised solely of the skeleton and is referred to as a solid. The mass balance equation states that the mass of a solid M(t) does not change in time: M(t) = M(0). Introducing the mass per unit volume in the current configuration ρ s (x, t) and reference configuration M(t) = ∫ Ω t ρ s (x, t) dx = ∫ Ω 0 Jρ s (χ(X, t), t) dX = ∫ Ω 0 ρ 0 s (X) dX = M(0), (6) or, locally, Jρ s (χ(X, t), t) = ρ 0 s (X) .
For an incompressible solid, J ≡ 1 and ρ s (χ(X, t), t) ≡ ρ 0 s (X). With this, differentiating the above equation with respect to time yields the rate form, J˙= J ∇ ⋅ v = 0, in which ∇ · is the divergence in the current configuration. Because v = v(X, t) is the Lagrangian structural velocity, this is stated in fully Lagrangian form by H: ∇ X v = 0, (8) in which H = JF −T is the cofactor matrix of F. The cofactor matrix can be written more conveniently using the cross product between matrices [34], H = 1 2 F × F(In fact, because the cross-product is a linear operator, in this form, the linearization of the cofactor is straightforward.).
Momentum balance in the Lagrangian frame is in which b is the body force per unit mass, and S is the second Piola-Kirchhoff (PK2) stress.
To model the myocardium, we assume the PK2 stress can be additively decomposed into active and passive parts: S = S a + S p . The active stress S a is obtained from a contraction model, whereas the passive stress S p is obtained from a strain energy function Ψ = Ψ(C), such that S p = 2 ∂Ψ (C) ∂C , (10) in which C = F T F is the right Cauchy-Green strain tensor. Eq. (10) corresponds to a hyperelastic material description. To enforce the constraint (8), we introduce the Lagrange multiplier p(X, t), which is identified as the solid pressure [22]. Next, introducing the distortional component of the deformation gradient tensor F = J −1/3 F and the corresponding distortional Cauchy-Green strain C = F T F, the strain energy function is and the passive component of the PK2 stress becomes in which W (C) is the strain energy function corresponding to the isochoric deformations. Considering the Lagrange multiplier p(X, t) to include any spherical component of the total stress, the total PK2 stress is modified to include only the deviatoric part of the active stress, given as Here, DEV[S] is the PK2 stress tensor corresponding to the deviatoric Cauchy stress [35], given as where dev[σ] = σ − 1 3 tr(σ)I . (15) in which dev[σ] is the deviatoric component of the Cauchy stress σ. Because Using Eq. (16a), it is possible to write the momentum equation (Eq. (16b)) as a second-order differential equation in time. However, since a key part of the proposed stabilization method (explained later in Section 3.2) is that it acts on the velocity field, it is convenient to write the equation as a first order system. As shown in Rossi et al. [33], after discretization, this choice results in a two-step algorithm in which we first solve for the velocity and pressure fields and then update the displacement field. Note that, as demonstrated in Scovazzi et al. [31] for isotropic incompressible and nearly incompressible elasticity, we obtain a stable scheme for the dynamic problem only with velocity stabilization.

Poroelastodynamics
Although we consider both the skeleton and the pore fluid as incompressible, the homogeneous mixture may not be incompressible due to the fluid moving within the skeleton [11,36]. The total mass of the homogeneous mixture as a function of time is M(t) = M f (t) + M s (t) and is evaluated via The density of the mixture is ρ(x, t) = ϕ(x, t)ρ f (x, t) + (1 − ϕ(x, t))ρ s (x, t). For an unconfined poroelastic medium, even if the fluid and the solid are incompressible, the mixture is generally not incompressible (ρ ≠ ρ 0 ) due to the presence of an external source. For cardiac perfusion, the porous medium (myocardium) is connected to an external source (large blood vessels) and is unconfined.
in which m(X, t) is the (Lagrangian) added mass, i.e., the additional fluid mass per unit reference volume. For an incompressible solid, ρ s ≡ ρ 0 s . Taking the time derivative of Eq.
(18b), and using J˙= J ∇ ⋅ v, we find Jρ f ϕ( ∇ ⋅ v) + J ρ f ϕ = ṁ . This can be simplified as Using the definition of perfusion velocity w(x, t) = ϕ(x, t)(v f (χ −1 (x, t), t) − v(χ −1 (x, t), t)), and defining the Eulerian mass source ρ f s = ∂ ρ f ϕ ∂t + ∇ ⋅ ρ f ϕv f (20) in which s(x, t) is the volumetric source in the current configuration, we obtain, after algebraic manipulation, In the formulation used herein, the perfusion velocity is assumed to be governed by a Darcy flow model, such that w(x, t) = − K ∇p pore χ −1 (x, t), t , (22) in which p pore (X, t) is the pore pressure in the reference configuration and K is the permeability tensor. We assume the total pore pressure p pore (X, t) is the sum of three contributions: the Lagrange multiplier p(X, t) that impose skeleton incompressibility; the pressure derived from the pore pressure-volume relationship p PV (X, t); and the compaction pressure p c (X, t) that enforces a non-negative porosity ϕ at all times, so that p pore (X, t) = p(X, t) + p PV (X, t) + p c (X, t), (23) in which p PV and p c are functions (specified later) of the added mass m. Because the pore pressure and the added mass are defined in the reference configuration, Eq. (21) is pulled back to the reference configuration using the incompressibility assumption for the fluid, ρ f = ρ 0 f , along with Eqs. (22) and (23) to obtain the Lagrangian form of the added mass equation, (24) in which the Lagrangian permeability tensor is K 0 = JF −1 KF −T , the added mass dependent permeability tensor is K 0 m = K 0 ∂ p PV + p c / ∂m, and the Lagrangian volumetric source is S(X, t).
The constraint that the skeleton remain incompressible, J(1 − ϕ) = 1 − ϕ 0 , can be written in terms of the added mass m. In fact, by using Eq. (18) to eliminate ϕ, and using the fluid incompressibility assumption ρ f = ρ 0 f , we find the relationship J − 1 = m/ρ 0 f . Then, the incompressibility constraint in Lagrangian coordinates can be stated in rate form as Momentum balance of the mixture involves the inertia of both the skeleton and the pore fluid. Assuming the fluid and the skeleton have the same acceleration, we can simplify the equations of momentum balance, obtaining ρ 0 + m v = ∇ X ⋅ (Dev[P]) − ∇ X ⋅ (pH) + ρ 0 b, (26) in which ρ 0 = ρ(χ(X, 0), 0) = ρ 0 f ϕ 0 + ρ 0 s 1 − ϕ 0 The stress tensor S and the fluid pressures p PV and p c are derived from an energy function that depends on the strain tensor C, the Lagrange multiplier p and the added mass m, such that Ψ = Ψ(C, p, m).
in which W (C) characterizes the skeleton isochoric deformations, U(m) characterizes the changes in volume due to the fluid motion, and the last term enforces the incompressibility constraint of the solid skeleton. The above strain energy function allows for negative values of ϕ depending on the value of the source S. To enforce a non-negative value for the porosity ϕ, we split the energy U(m) into a pressure-volume relation U PV (m) and a compaction penalization U c (m), such that  m ρ 0 f . (28) With this definition of the energy, we obtain The final system of equations for the poroelastic mixture is Note that similar to Eq. (16), because the stabilization (Section 3.3) acts on the velocity field, it is convenient to write the balance equations as a first-order system.
Remark.-The density of the fluid ρ f and the pore pressure p pore are related by the free enthalpy of the fluid g m [11], given as Following an isothermal assumption, this relationship is simplified as [37] g m = p pore − p 0 pore ρ 0 f , (32) in which p 0 pore is the pore pressure in the initial reference configuration. The above choice holds if the fluid is incompressible, i.e.ρ f = ρ 0 f .

Constitutive equations
We verify that our formulation is stable and robust using several constitutive laws. For isotropic elastic solids, we use the neo-Hookean and Mooney-Rivlin models, Mooney−Rivlin: W (C) = C 1 I 1 − 3 + C 2 I 2 − 3 , (34) in which G, C 1 , and C 2 are material constants, and I 1 = tr(C) and I 2 = 1 2 (I 1 2 + tr(C 2 )) are invariants of C. For anisotropic elastic solids, we use the standard reinforced model, in which G f is a material constant and I 4f = f 0 ⋅ Cf 0 , with f 0 being the fibre direction in the initial reference configuration.
For the pore-fluid, the energy corresponding to the pressure-volume relationship U PV (m) and the compaction penalization U c (m) are given as in which κ s and c are constants, ϕ crit is the critical value of ϕ up to which compaction is permitted, and ϵ is a small value. Here, c is the maximum value of the pore pressure from the previous time step. When (m + ϕ 0 ) ≈ ϕ crit , the pore pressure component p c given as is very large to prevent further pore fluid extraction [11,22]. The functional form of p c is chosen this way to prevent m + ϕ 0 becoming smaller than ϕ crit . In fact, for very small value of ϵ ≈ 0, the contribution introduced by p c is negligible if m + ϕ 0 ≫ ϕ crit but becomes rapidly large as m + ϕ 0 approaches ϕ crit . However, very small values of ϵ will lead to divergence of the nonlinear solver. Preliminary numerical experiments not reported here suggest that the choice ϵ = 0.001 and ϕ crit = 0.001 ensures convergence of the nonlinear solver while enforcing the constraint m + ϕ0 > ϕcrit.
To model the biomechanics of the left ventricle, we use the Holzapfel-Ogden (HO) model [38], in which a, b, a i , b i , a fs , and b fs are material constants, I 4 s = s 0 ⋅ Cs 0 , and I 8fs = f 0 ⋅ Cs 0 , with s 0 the sheet direction in the reference configuration. For the pressure-volume relationship within the pores, we use the energy function proposed by Bruinsma et al. [39], U PV (m) = q 1 q 3 exp q 3 m + ϕ 0 + q 2 m + ϕ 0 ln q 3 m + ϕ 0 − 1 − q 1 exp q 3 ϕ 0 + q 2 ln ϕ 0 m + ϕ 0 , (39) in which q 1 , q 2 , and q 3 are material constants.
Finally, the active component of the stress tensor is S a = T t, I 4f f 0 ⊗ f 0 , (40) in which the scalar tension T = T(t, I 4f ) is a function of time and the current fibre elongation in which T a (t) is the active tension that is determined empirically.

Numerical formulation
This section introduces a finite element method specialized to piecewise linear (P 1 ) finite elements for the incompressible poroelastic mixture. In particular, we extend a previously developed stabilization method for P 1 elements for isotropic elastodynamics [33] to anisotropic hyperelastic materials and anisotropic poroelastic materials. We use the second-order backward differentiation (BDF2) method to integrate the elastodynamic and poroelastodynamic equations in time. As in prior work [33], the BDF2 method with a time step size controlled by the CFL condition leads to accurate dynamics and dissipation of only the undesired high frequencies.
We denote by (v, w) Ω 0 = ∫ Ω 0 vw dX and (v, w) Ω 0 = ∫ Ω 0 v ⋅ wdX (42) the L 2 (Ω 0 ) and L 2 (Ω 0 ) d inner products on the interior of the domain Ω 0 and by v, w Γ 0 = ∫ Γ 0 vw dA (43) a functional on a regular subset Γ 0 of the boundary ∂Ω 0 . We also define the element-wise inner products in which K is an element of the triangulation T ℎ , such that Ω = ∪ K ∈ T ℎ K. Here we consider K to represent either a triangle (in two spatial dimensions) or a tetrahedron (in three dimensions). To simplify the notation, from this point on, we omit the subscript X if the

Author Manuscript
Author Manuscript Author Manuscript Author Manuscript inner products are taken over Ω 0 or Ω 0 ′ , because all quantities are Lagrangian quantities described using reference coordinates.

Abstract variational multiscale framework for nonlinear dynamics
We present here the abstract framework for nonlinear problems [33] using the D-VMS method. Denote by V an infinite dimensional space, V* its dual, and N a nonlinear operator associated with a generic nonlinear problem whose variational statement can read as Find y ∈ V such that, ∀ φ ∈ V, in which y is the solution vector and N(y) ∈ V *. The variational multiscale framework introduces a decomposition of the space V into a finite dimensional subspace V and its compliment V′ such that V = V + V ′. With this, we can now decompose into a resolvable coarse scale and unresolvable fine scale vectors y = y + y′ and φ = φ + φ′ such that (45) can be recast as Find y ∈ V and y′ ∈ V′ such that, ∀φ ∈ V and ∀φ′ ∈ V′, Assuming the existence of the Fréchet derivative ℒ[y]( ⋅ ) of N( ⋅ ), as in Rossi et al. [33], defined as we can approximate the nonlinear problem as follows Find y ∈ V and y′ ∈ V′ such that, ∀φ ∈ V and ∀φ′ ∈ V′, As we will demonstrate in the numerical test, the linearized problem above results in a second-order accurate method when using linear finite elements.
Introducing the formal adjoint of ℒ*[y] of ℒ[y], the above problem is recast as Find y ∈ V and y′ ∈ V′ such that, ∀φ ∈ V and ∀φ′ ∈ V′, The solution to (49b) can be formally expressed in terms of the fine-scale Green's function integral operator ℳ′ in which g′(·, ·) is the fine-scale Green's function. Approximating the Green's function using localization arguments as g′(X, Y) ≈ τ y δ(X − Y), in which δ is the multidimensional Dirac delta function, we obtain the multiscale variational formulation, Find y ∈ V such that, The matrix τ y contains the problem dependent parameters that we use to stabilize the variational formulation.

Variational multiscale method for hyperelastodynamics
We now apply the variational multiscale framework (51) to the hyperelastodynamic problem (16). Using the vectors y = [u, v, p] and φ = [w, w, q], we have Introducing the multiscale decomposition and imposing u = v, we find that u′ = v′, and in which F, H and P are, respectively, the deformation gradient tensor, its cofactor, and the first Piola-Kirchhoff stress tensors evaluated using the coarse scale displacements u. Integration by parts leads to Using a minimalist approach to avoid spurious pressure oscillations [33], we approximate This approximation is equivalent to choosing Introducing the finite dimensional test function spaces and the corresponding trial function spaces the stabilized weak form of the hyperelastodynamic equations reads Find u h , v ℎ ∈ S v ℎ , and p ℎ ∈ S q ℎ , such that, ∀w ℎ , w ℎ ∈ V v ℎ , and q ℎ ∈ V q ℎ , in which the terms corresponding to the deviatoric component of the stress tensor disappear because of the choice of linear shape functions. We use the BDF2 scheme for the time discretization and Newton's method for the nonlinear equations. After discretization, as shown in Rossi et al. [33], the system can be solved as a two-step algorithm in which we first solve for the velocity and pressure fields, followed by an explicit update of the displacement field. The numerical implementation follows that given by Rossi et al. [33].

Variational multiscale method for poroelastodynamics
We now apply the variational multiscale framework (51) to the poroelastodynamic problem (30). Using the vectors y = [u, v, p, m] and φ = [w, w, q, r], we have Introducing the multiscale decomposition as before, we find in which K 0 = JF −1 KF −T and K 0 m = K 0 ∂ p PV + p c / ∂m. We remark that in (62), we consider the permeability tensor K 0 to be a model parameter defined in the reference configuration. If K 0 is defined in current configuration, then Eq. (62) should contain its linearization in terms of the displacements. Integration by parts leads to in which we assume K 0 = K 0 T . Using again a minimalist approach to avoid spurious pressure oscillations, we choose  Using the function spaces defined in (58) and (59), the stabilized weak form of the poroelastodynamic equations reads, Find u h , v ℎ ∈ S v ℎ , and p h , m ℎ ∈ S q ℎ , such that, ∀w ℎ , in which the terms corresponding to the deviatoric component of the stress tensor disappear because we have specialized the formulation to piecewise linear shape functions, and the vector W corresponds to the Lagrangian perfusion velocity. Similar to Eq. (60), we discretize Eq. (65) in time using the BDF2 scheme and solve the resulting system of equations for the velocity, pressure, and added mass using Newton's method. We then update the displacement explicitly.

Choice of the stabilization parameter
For simplicity, the stabilization matrix τ is chosen to be isotropic [33], i.e., τ = τI. The stabilization parameter τ is an intrinsic time scale that can be obtained in various ways [31]. Here, h e is the characteristic length of element e, c μ e is the shear wave velocity, and T ℎ is the general triangulation. Although the stabilization term adds some numerical dissipation, its contribution is small, thanks to the choice of the stabilization parameter in the D-VMS method. This was demonstrated in [31,33] and shown in the numerical results in the following sections. Further, for the D-VMS method, numerical dissipation is necessary to satisfy the LBB criterion. The expression for τ (Eq. (66)) was obtained in previous studies [31,33] after intensive numerical experiments to minimize the numerical dissipation while maintaining stabilization. In the expression (Eq. (66)), the term min(Δt μ ,Δt) is chosen so that the numerical dissipation is not too large to cause catastrophic failures in computation in the limit Δt → ∞. The term Δt μ /100 is introduced so that the numerical dissipation is stable in the limit Δt → 0. Further, to reduce excessive numerical dissipation for incompressible case, Rossi et al. [33] introduced the value of c τ in the range [0.01, 0.03] after numerical experimentation to maintain the stabilization, and we follow this approach here.

Choice of time step size
Even though the implicit time integration is unconditionally stable, for strongly nonlinear problems, the time step size choice is closely related to the convergence of nonlinear iterations. For the nonlinear convergence, an empirical relationship is obtained for the time step size as proportional to the velocity of the shear wave, in which α CFL is the CFL number. In our numerical experiments, when we use 0.1 ≤ α CFL ≤ 1, the nonlinear solver generally requires at most five Newton iterations. Note that stabilization method and the assumptions Eqs.

Results
Numerical tests are presented for benchmark test cases in two and three dimensions for both hyperelastic and poroelastic problems. In addition, the method's capability is demonstrated by simulating LV dynamics and perfusion in an anatomical model derived from human image data.

Elastodynamics
This section details numerical tests of the anisotropic hyperelastic formulation, including a two-dimensional anisotropic compressed block, three-dimensional shear deformation, threedimensional bending column, and three-dimensional twisting column. The results obtained from the D-VMS method are first compared with the standard Q 1 finite element with reduced selective integration (Q 1 −B-bar) method [23], P 1 − P 1 (piecewise linear interpolation for displacement, velocity and pressure) method, and the P 1 − P 0 (piecewise linear interpolation for displacement and velocity and piecewise constant interpolation for pressure) method. Fig. 2 compares the grid convergence for the displacement, velocity, and pressure for the D-VMS and Q 1 −B-bar methods. It is clear from the figure that the displacement, velocity and pressure converge with grid refinement. Further, Fig. 3 shows that the pressure field obtained using the D-VMS method is smooth and in agreement with the Q 1 −B-bar method. Thus, it is evident that D-VMS produces a solution that is free of locking with a non-oscillatory pressure field. Note that the P 1 − P 0 locks, and the P 1 − P 1 formulation appears to avoid locking but generates an oscillatory pressure field.
where u X * , u Y * , and u Z * . are the displacements in X, Y, and Z directions, and α and β are constants.
The exact solution for the pressure is   6 shows the pressure obtained from the numerical simulation using the D-VMS and P 1 − P 1 methods and the error in the pressure at t = 0.25s. The simulations are carried out for meshes with mesh spacing h = 0.167m and time step size Δt = 0.05s. The pressure field obtained using the D-VMS method is less oscillatory than the P 1 − P 1 method. Fig.  6 also shows the fibre stretches obtained using the D-VMS and P 1 − P 1 methods and the error in fibre stretch obtained using the D-VMS method. The maximum value of the error in pressure and the fibre stretch are found to be 0.0033Pa and 0.005, respectively. Similar to the pressure, the fibre stretch is also oscillatory for the P 1 − P 1 method compared to the D-VMS method.
A convergence study was performed using a set of five grids with the coarsest mesh spacing of h = 0.25m (with time step size Δt = 0.1s) and the finest mesh spacing of h = 0.015625m (with time step size Δt = 0.00625s). Fig. 7 shows the convergence of displacement, velocity, pressure, and fibre stretch. Similar to an earlier isotropic convergence study [33], a secondorder convergence is observed for all variables. This illustrates that the D-VMS method can be extended to anisotropic cases without loss in accuracy.

Hyperelastic anisotropic bending column-We
next consider a threedimensional anisotropic bending column. Fig. 8 shows the computational domain for the bending column. The bending axis is slightly inclined such that the solution is not symmetric. A fixed displacement boundary condition is used for the bottom plane Y = 0, and a traction-free boundary condition is used for all the other faces. The initial velocity is set to v(X, 0) = 5 3 Y , 0, 0 m/s .
All other physical parameters are the same as that of the bending column. The mesh for the grid convergence is also the same as the bending problem. Fig. 10 compares the pressure fields obtained from the various methods. Similar to the bending column, the pressures obtained from the P 1 − P 1 method are highly oscillatory. Fig. 11 shows that the result for the finer mesh obtained using the D-VMS method is in good agreement with that obtained from the Q 1 -B-bar method. Fig. 11 also shows the deformed shape of the twisting column without the fibres (G f = 0). It is clear from the figure that the anisotropicity introduced asymmetric twisting due to the fibre reinforcement. Fig. 12 shows clear convergence of the pressure and the fibre stretch with refined grids.

Poroelastodynamics
This section presents numerical tests of the poroelastic formulation, including an anisotropic two-dimensional compressing block, isotropic three-dimensional shear deformation, and isotropic three-dimensional swelling and shrinking cubes.

Anisotropic porous compressed block-
We now consider a porous version of the anisotropic compressed block discussed previously with an initial porosity ϕ 0 = 0.1 and a permeability K = 1.0 × 10 −5 I cm 4 dyne −1 s −1 . The source term is proportional to the pore pressure (S = β so p pore ). A source constant of β so = 1.0 × 10 −5 cm 2 dyne −1 s −1 is used. Fig. 13 compares the fields of pore pressure and added mass obtained from the D-VMS method and the P 1 − P 1 method. The P 1 − P 1 results are clearly oscillatory, although they are less oscillatory as compared to the hyperelastodynamics case. This may be because the skeleton is slightly compressible for the poroelastic case due to the presence of the fluid as compared to the full incompressibility in the hyperelastic case. Fig. 14 shows that the pore pressure field p pore and the fields of added mass m converge under mesh refinement. Fig. 15 shows the convergence of the pressure p, the pore pressure p pore , and the added mass m at the midpoint on the top surface of the compressed block.  (1 + αZsin (t)) 8 3 + κρ 0 f αZsin (t) .

Isotropic porous shear-Next
(77) The detailed derivation of this expression is given in Appendix B.
For the numerical simulations, we use the initial condition v (X, 0) = v* (X, 0) and m (X, 0) = 0. The boundary conditions are the same as that of the hyperelastic shear deformation (Section 4.1.2). For the pore-fluid, c = 0 is used in the compaction penalization (Eq. (36)). Fig. 16 compares the pore pressure fields and the added mass obtained from the P 1 − P 1 method and the D-VMS method. Similar to the hyperelastic case, the D-VMS method results in a smooth pressure field compared to the oscillatory pressure field obtained using the P 1 − P 1 method. Fig. 16 also shows a negligible difference between the results obtained from the simulation and the exact solution. Fig. 17 shows that the displacement, velocity, pore pressure, and added mass obtained using the D-VMS method are second-order accurate.

Swelling and shrinking-
Here we show the performance of the D-VMS method on the two poroelastic test cases proposed by Chapelle et al. [11]. Fig. 18 shows the computational domains for the two cases: a swelling cube and a shrinking cube. For the swelling and shrinking cases, the permeabilities of the fluid are taken as K = 10 −4 I m 2 kPa −1 s −1 and K = 2.5 × 10 −3 I m 2 kPa −1 s −1 , respectively. The simulations are performed using a mesh with mesh spacing h = 0.05cm (and time step size Δt = 0.001s) for the swelling cube and h = 0.05 mm (and time step size Δt = 0.0001s) for the shrinking cube. For the first test case, the swelling of the cube is due to the pore pressure difference across two opposite faces of the cube (Fig. 18a). The pore pressure p pore is gradually increased on the left surface (X = 0), and p pore = 0 is applied to the right surface. The other four faces are set to be impermeable to constrain the flow to be in one direction, and no source/sink is used. Note that no external forces are acting on the cube for this problem. The deformation is entirely driven by the Darcy flow in the pore fluid. Chapelle et al. [11] used a pressure-based equation for Darcy flow, whereas we used the added mass-based formulation. Thus, for the cube swelling case, the pore pressure boundary conditions used by Chapelle et al. [11] are converted to boundary conditions for m. The added mass m on the left surface (X = 0cm) is set to increase gradually m(0, Y , Z, t) = 0.5 1 − exp −t 2 /0.25 ρ 0 f , and m on the right surface (X = 1.0cm) is kept fixed at zero (Fig. 18a).
For the second case, the shrinking of the cube is driven by external forces applied on all six faces of the cube. In this case, normal displacement is constrained on three faces (X = 0,Y = 0, Z = 0), and normal forces are applied on the other three faces (X = 1.0 mm,Y = 1.0 mm, Z = 1.0 mm). In addition, a pressure-dependent sink term, S = −β si (p pore − p si ), is imposed, and a zero normal gradient of m is applied on all the faces of the cube. For the numerical simulation, β si = 0.01 kPa −1 s −1 and p si = 0 kPa are used. Note that for this case, the fluid added mass and the porosity are coupled with the normal solid pressure and the solid pressure and its deformation depends on the porosity (Eq. (30)). Thus, the problem serves to verify the two-way coupling between the solid deformation and the Darcy flow. Fig. 19 shows the deformed position and fields for the total pore pressure p pore (for the swelling cube) and added mass m (for the shrinking cube). For the swelling case, Fig. 19c shows that the gradient in m and pore pressure p pore resulted in swelling of the cube near the left surface (X = 0) and fluid flow from the left to the right surface. The swelling gradually decreases from the left to the right surface. Fig. 19c shows an excellent agreement between our results and the literature for the transient variation of pore pressure p pore and m at the middle point (0 cm,0 cm,0 cm) and the right corner point (0.5 cm,0.5 cm,0.5 cm).
For the shrinking case, Fig. 19d shows that the fluid in the cube is completely drained at a steady-state, resulting in a reduction in the volume of the cube by almost 10%. Fig. 19d also shows the transient variation of the pore pressure p pore and m at the mid-point. The pressure increases gradually as per the applied pressure on the three surfaces. When the fluid is almost completely drained out (m ≈ −ϕ 0 ), the penalty pressure p c takes a very large negative value, and the total pore pressure p pore suddenly drops to zero, restricting further drainage of fluid (Fig. 19d). The results are in very good agreement with the results from Chapelle et al. [11].

Perfusion in left ventricle
This section presents numerical simulations of perfusion in the left ventricle using a previously developed image-based geometry [40] shown in Fig. 20. Like many previous studies [11,13,22], we only consider a portion of the LV. As illustrated in Fig. 20a, we apply zero normal and circumferential displacements at the base, allowing only radial expansion along the plane of the base of the LV.
The blood pressure inside the LV is modelled by a uniform pressure throughout the endocardial surface with a prescribed waveform p endo (t). The HO model (Eq. (38) in the myocardial skeleton are modelled using myofibres and collagen sheets as specified in Richardson et al. [13]. The myocardium is considered porous, with source/sink terms dependent on the pore pressure [11], S(X, t) = β so p so − p pore (X, t) − β si p pore (X, t) − p si , (78) in which β so and β si are constants and p so and p si are the source and sink pressures. For the numerical simulation, β so = β si = 3 × 10 −6 cm 2 dyne −1 s −1 , p so = 2.7 × 10 4 dyne/cm 2 , and p si = 1.3 × 10 4 dyne/cm 2 are used.
We simulate three consecutive cardiac cycles, including diastolic and systolic phases. The endocardial pressure is applied with a maximum values of 1.067 kPa and 14.53 kPa during the diastole and systole, respectively, given by p endo (t) =    Fig. 21a shows the temporal variation of the volume of the myocardium for both hyperelastic and poroelastic cases for three complete cycles. The hyperelastic simulations shown here are carried out using the same methodology but with zero porosity and external source. For both cases, a periodic state is obtained by the third cycle. The volume remains almost constant throughout the cardiac cycle for the hyperelastic case because the myocardium is treated as incompressible. For the poroelastic case, the volume changes depending on the added mass. Note that the volume of the skeleton J − m/ρ 0 f remains almost constant throughout the cycle.
In our poroelastic formulation, the change in the added mass m is governed by the Darcy flow pore fluid model, and the fluid flow into or out of the pores depends on the pore pressure-dependent source/sink. Fig. 21b shows the temporal variation of the space-averaged m for the endocardial and epicardial surfaces. During the diastolic phase, m increases on both surfaces because of the smaller pore pressure, which results in fluid flow from the source to the pores. Notice that at the end of diastole, the average m is larger at the endocardium than the epicardium. Furthermore, m falls at both the endocardium and epicardium during the systole. However, the rate of decrease in m is greater near the endocardium. The variation of m at the epicardium and endocardium in the diastole and systole could be attributed to the larger pore pressure at the epicardium (endocardium) during the diastolic (systolic) phase (Fig. 22). Fig. 22 shows the fields of p, p PV , and p pore at various cross-sections of the LV at the periodic state. Fig. 22a, b, and q show that, the pressure is negative everywhere at the end of the diastole. This is because the inflation elongates the fibres that are tangential to the surfaces at the epicardium and endocardium. Thus, the tensile forces along the fibres make the volumetric component of the total fibre stress and hence the pressure negative [8,42,43]. Further, the magnitude of the negative pressure is larger near the endocardium than the epicardium. This results from larger fibre elongation and the associated fibre stress near the endocardium as compared to the epicardium (Fig. 22q). In addition, the pressure is larger near the base than the middle and apex portions. This can be attributed to the constraints imposed on the displacements of the base and associated reductions in the fibre stretch and fibre stress. The fibre stretch and the fibre stress are also smaller near the apex than the middle portion because of the restricted motion near the apex, resulting in a slightly larger pressure near the apex. Fig.  22b shows that the pressure variation in the transmural direction is different at different angular locations because of the asymmetric deformation. During the systolic phase, Fig.  22c, d, and r show that p achieves a maximum near the endocardium and a minimum near the epicardium. Further, the pressure is positive everywhere because the fibres are compressed during contraction, and the fibre stress is positive. The fibre stress is more near the endocardium than the epicardium (Fig. 22e). Thus, the compressive force and associated pressure are greater near the endocardium than the epicardium. Fig. 23 shows temporal variations of the pore pressure p pore , its components p and p PV , the added mass m, and the source S. Note that when m = 0, p PV is zero. Thus, p is the pressure driving the source S and the added mass m in the beginning of the cycle. Because of the negative pressure during the diastole (Fig. 22b), the associated smaller pore pressure results in fluid flow from the source to the pores and an increase in the added mass m. With increasing m, p PV increases exponentially, leading to a positive total pore pressure p pore at both surfaces (Fig. 22j). However, the total pore pressure p pore is less than the source pressure of 2.7 kPa, so that there is net flow into the myocardium. Thus, the added mass m increases at the epicardium and the endocardium throughout the diastole (Fig. 23). Note that, the space averaged p and the pore pressure p pore are slightly larger near the epicardium than the endocardium at the end of the diastole (Fig. 22q). Thus, the added mass m is slightly larger near the endocardium than the epicardium. By the end of the diastole, the total pore pressure p pore nearly balances with the source and sink pressures, and the source S is nearly zero (Fig. 23), leading to a constant added mass m. Fig. 23 shows that p increases rapidly at both the endocardium and the epicardium during systole. The pressure p is very large near the endocardium as compared to the epicardium (Fig. 22r). At both surfaces, p pore increases and the associated source S rapidly decreases ( Fig. 23), leading to the shrinking of the pores, similar to that reported in the literature [11,22,42]. This leads to a rapid decrease in the added mass m. Note that the rate of decrease of m is very large near the endocardium as compared to the epicardium. This is because the pore pressure p pore is larger near the endocardium than the epicardium (Fig. 22r). Since p PV also increases exponentially with increasing m (Fig. 23), the rate of increase of p pore is less than that of p. Furthermore, the source S reaches its minimum value for both the endocardium and the epicardium when the pore-fluid pressure p pore reaches its maximum (Fig. 23). However, the added mass m continues to drop until the source S becomes positive (Fig. 23). We remark that the order of magnitude and the type of variation of the pore pressure are in good agreement with the literature [11,22].

Conclusion
This paper develops a stabilized equal-order mixed finite element method for anisotropic incompressible hyperelastodynamics and poroelastodynamics using linear finite elements. Our approach extends a dynamic variational multiscale method, presented by Rossi et al. [33] for incompressible isotropic hyperelastodynamics, to anisotropic hyperelastodynamics and poroelastodynamics. We verify the method's convergence, stability, and robustness for various anisotropic hyperelastic problems and isotropic/anisotropic poroelastic problems. We find that our method results in second-order accuracy in space and time, and the results are in excellent agreement with the standard selective reduced integration method (Q 1 -Bbar). For poroelastodynamics, our results are in excellent agreement with the benchmark results of Chapelle et al. [11]. In addition, poroelastic simulations of ventricular perfusion are performed using a realistic left ventricular geometry. Good qualitative agreement is obtained with prior results reported in the literature. The temporal variation of the various components of the pore pressure are correlated with the variation of the added mass m and the dynamics of the myocardium. The anisotropic behaviour greatly influences the magnitude of the pore pressure and its distribution across the myocardium. The skeleton pressure during the diastole (systole) was found to be negative (positive) in most of the places due to tension (compression) on the myofibres.  The right Cauchy-Green strain C is C(X, t) = 1 + αZ sin 2 (t) 0 (1 + αZsin (t))αXsin (t) 0 1 βsin(t) (1 + αZsin (t))αXsin (t) βsin (t) −α 2 cos 2 (t)X 2 + α 2 X 2 − β 2 cos 2 (t) + β 2 + 1 , (B.3) and the first invariant of C is I 1 = − α 2 cos 2 (t)X 2 − α 2 cos 2 (t)Z 2 + α 2 X 2 + α 2 Z 2 − β 2 cos 2 (t) + 2αZsin (t) + β 2 + 3.

(B.4)
Using the neo-Hookean model, the divergence of the deviatoric part of PK1 stress is given as 3(1 + αZsin (t)) 5 3 e X + − 2Gβsin 2 (t)α 3(1 + αZsin (t)) 5 3 e Y + (1 + αZsin (t)) −5 3 9 −sin (t)Gα 11α 2 cos 2 (t)X 2 − α 2 cos 2 (t)Z 2 + 5β 2 cos 2 (t) −11α 2 X 2 + α 2 Z 2 + 2αZsin (t) − 5β 2 − 3 e Z . Using the following expression for the body force, b = (0) e X + (0) e Y + sin (t) 18(1 + αZsin (t)) −9 Z 2 α α 2 X 2 + 2β 2 sin 2 (t) + 2sin (t)β 2 Z − αX 2 − 2βY ρ 0 f + 22sin (t) sin 2 (t)G Zα α 2 X 2 + (1 + αZsin (t)) 5 3 , (1 + αZsin (t)) 5 3 11α 2 cos 2 (t)X 2 − α 2 cos 2 (t)Z 2 + 5β 2 cos 2 (t) − 11α 2 X 2 + α 2 Z 2 + 2αZsin (t) − 5β 2 − 3 + b Z . (B.8c) The above system of equations can be solved to obtain  (1 + αZsin (t)) 8 3 . (B.9) The pore pressure p pore is  Convergence for the displacement u A , velocity v A , and pressure p A at the midpoint on the top surface of the anisotropic compressed block obtained using the D-VMS method and Q 1 -B-bar method at t = 1.0s.  Pressure fields for the two-dimensional anisotropic compressed block obtained using the D-VMS method, Q 1 -B-bar method, P 1 − P 1 method, and P 1 − P 0 method at t = 0.1s for the finest mesh spacing h = 0.2 cm.      Pressure fields (top) and fibre stretch fields (bottom) obtained using the D-VMS method and P 1 − P 1 method along with the absolute errors in the pressure and fibre stretch at t = 0.25s for the shear deformation of the anisotropic cube (Section 4.1.2).    is fixed and the axis of bending is asymmetric so that the solution is not axisymmetric. Right: Pressure fields for the three-dimensional anisotropic bending column problem obtained using the D-VMS method, Q 1 -B-bar method, and P 1 − P 1 method at t = 0.2s.  Deformed shapes of the three-dimensional anisotropic bending column obtained using the D-VMS method (shown with mesh) with various mesh spacing superimposed with the shape obtained using the Q 1 -B-bar method (with h = 0.167m) at t = 0.5s.  Pressure fields for the three-dimensional anisotropic twisting column obtained using the D-VMS method, Q 1 -B-bar method, and P 1 − P 1 method at t = 0.02s (Section 4.1.4).  Deformed shapes of the three-dimensional anisotropic twisting column problem obtained using the D-VMS method (shown with mesh) with various mesh spacing superimposed with the shape obtained using the Q 1 -B-bar (with h = 0.167m) method and deformed shape of three-dimensional twisting column without fibres (G f = 0) obtained using the Q 1 -B-bar method at t = 0.1s.  Pressure fields (left) and fibre stretch fields (right) for the three-dimensional anisotropic twisting column obtained using the D-VMS method with various meshes and the Q 1 -B-bar method at t = 0.1s.  Pore pressure fields (left) and added mass fields (right) for the anisotropic porous compressed block obtained using the D-VMS method and P 1 − P 1 method at t = 0.3s with mesh spacing h = 0.2cm (Section 4.2.1).  Pore pressure fields (top) and added mass fields (bottom) for the anisotropic porous compressed block obtained using different meshes at t = 1.0s.  Pore pressure fields (top) and added mass fields (bottom) obtained using the D-VMS method and P 1 − P 1 method and absolute errors in the pore pressure and added mass at t = 0.25s for the shear deformation of the poroelastic isotropic cube (Section 4.2.2).    Computational setup along with boundary conditions for the verification of poroelastic problems of (a) a swelling poroelastic cube resulting from an applied pore pressure gradient, and (b) a shrinking poroelastic cube resulting from applied compressive forces (Section 4.2.3).       Temporal variation of the volume of the myocardium LV vol for the poroelastic and hyperelastic cases and added mass m at the epicardial and endocardial surfaces, for three complete cardiac cycles.  Fields of (a) -(d) p, (e) -(h) p PV , (i) -(l) p pore , and (m) -(p) m at different cross sections of the LV (X = 0cm, Y = 0cm, and Z = 3.0cm) at the end of the diastole t = 2.1s and the peak of systole t = 2.25s, and the variation of p, p PV , p pore , and m across the line marked in the plane Z = 3.0cm, at the (q) end of diastole and (r) the peak of systole.  Temporal variation of p, p PV , p pore , m, and S at the endocardial and epicardial surfaces during the third cardiac cycle.