{'fraction_non_alphanumeric': 0.033939128530764175, 'fraction_numerical': 0.01126659633337978, 'mean_word_length': 4.750028613940712, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
STABILIZATION APPROACHES FOR THE HYPERELASTIC IMMERSED BOUNDARY METHOD FOR PROBLEMS OF LARGE-DEFORMATION INCOMPRESSIBLE ELASTICITY

The immersed boundary method is a model of fluid-structure interaction that describes a structure, or a collection of structures, immersed in fluid. This formulation uses Eulerian coordinates for the momentum, incompressibility, and viscosity of the fluidstructure system and Lagrangian coordinates for the structural deformations and resultant forces. Integral transforms with delta function kernels connect the two frames. In the continuum equations, both the fluid and the structure are typically modeled as incompressible. Upon discretization, however, the structure’s incompressibility is only approximately maintained. To obtain a robust method under large structural deformations, we introduce a volumetric energy in the solid region that stabilizes the formulation and improves the accuracy of the numerical scheme. This volumetric energy is added by decomposing the strain energy into isochoric and dilatational components, as in standard solid mechanics formulations. Using standard quasi-static solid mechanics benchmarks, we study the performance of the proposed stabilized method with various choices of the finite element basis employed for the structural discretization.


Introduction
The immersed boundary (IB) method is a framework for modeling fluid-structure interaction (FSI) introduced by Peskin [1,2] to model blood flow around the heart valve leaflets. Eulerian variables describe the momentum, velocity, and incompressibility of the fluid-solid system, and Lagrangian variables describe the deformation and resultant forces of the immersed structure. Integral transforms with delta function kernels mediate interaction between the two frames and maintain a continuous velocity field across the fluid-solid interface while avoiding the need for body-fitted descriptions of the fluid and structure. The method was originally formulated to describe thin structures occupying zero volume within the fluid [1,2] and was eventually extended to describe volumetric structures [3]. In the work of Boffi et al. [4], the IB equations are systematically derived in the framework of large-deformation continuum mechanics. This formulation is particularly useful for many biomedical applications, in which physiologically realistic and experimentally validated material models for soft tissue rely upon a continuum mechanics description (see, e.g., Holzapfel [5]). For this reason, in this work, we adopt the continuum formulation posed by Boffi et al. [4].
The numerical implementation used here follows the one described by Griffith and Luo [6]. This implementation uses a finite difference scheme to approximate the Eulerian equations and a finite element (FE) scheme for the Lagrangian equations, and it uses regularized delta 1 functions in approximations to the integral transforms. This stands in contrast to an earlier numerical method, the immersed finite element method (IFEM) [7], that uses the FE method to approximate both the Eulerian and Lagrangian equations. An advantage of our implementation is that the same regularized delta function is used over the entire domain, whereas IFEM constructs position specific delta functions using reproducing kernel particle methods (RKPM) [8]. IFEM, however, may use unstructured Eulerian grids, which is an advantage when working with complex computational domains.
Another hybrid IB method combining finite difference and finite element methods was proposed by Devendran and Peskin [9]. This method is only implemented for, and heavily relies upon, a representation of the structure using linear simplicial elements. Although different continuum material models may be selected, their implementation requires analytically calculating derivatives of the strain energy functional with respect to the coefficients of the finite element representation of the displacement field. Other implementations include an extension using radial basis functions to represent the structure [10], a particle based method to represent the structure [11], and the numerical implementation of Boffi et al. [4] and Roy et al. [12] that avoids the use of regularized delta functions, achieving regularization instead via the FE basis functions.
A key feature of the IB formulation is that one momentum equation is used for both the fluid and solid regions, and the Eulerian incompressibility constraint is imposed throughout the entire computational domain. In particular, because incompressibility is maintained in the Eulerian frame for all points in the computational domain, the structure's motion is automatically incompressible. The discretized IB equations, however, do not automatically maintain incompressibility within the solid region. This is because the operator that restricts the Eulerian velocity field to the Lagrangian mesh relies upon a regularized delta function with a non-zero gradient. This implies that the divergence of the Lagrangian velocity field will be non-zero as well. Further, once the velocity is restricted to the Lagrangian mesh, it is then projected onto the FE basis functions, and this projection is not guaranteed to preserve the incompressibility of the structural velocity field.
A systematic study of the loss of the structure's incompressibility was performed by Casquero et al. [13] using an IB-type method with divergence-conforming B-splines. Through the use of these basis functions, their method achieves negligible changes of volume in the Eulerian frame and reduced incompressibility errors in the Lagrangian frame. In the present work, we show that a simple stabilization method active only in the solid region can greatly reduce the loss of incompressibility. The proposed method can be used with Lagrangian elements without employing complicated function spaces. To assess the performance of the proposed method, we employ standard benchmark problems traditionally used for largedeformation incompressible elasticity.
The proposed stabilization method is rooted in the deviatoric-spherical decomposition of the Cauchy stress tensor. Specifically, the method penalizes changes in volume in the solid by adding an additional pressure contribution to the stress. This additional contribution depends on a stabilization parameter, which we call the numerical bulk modulus κ s . The method resembles standard displacement formulations of nearly incompressible structural mechanics, in which a volumetric penalty term is included in the structure's strain energy functional [14]. In nearly incompressible solid mechanics, the presence of this term is tuned via a physical parameter, the bulk modulus κ, representing the resistance of the solid to compression. As κ → ∞, the material can only experience incompressible deformations.
Unfortunately with many simple numerical methods, high values of the bulk modulus lead to volumetric locking, or sub-optimal convergence rates in the computed displacement [14].
In the proposed method, the penalization does not require κ s → ∞ since incompressibility is inherited from the Eulerian momentum equation, and we demonstrate that locking can be avoided even with simple linear finite elements. We modulate the numerical bulk modulus by using a numerical Poisson ratio ν s and standard linearized elasticity relations. We therefore limit the values of ν s to range between −1 and 1 2 . At ν s = −1 (equivalently κ s = 0), no volumetric stabilization is added and large changes in the solid's volume can be observed. We explore the effects of ν s as a volumetric stabilization parameter used to restore the incompressibility of the structure and show that a value of ν s = 0.4 is sufficient to achieve this goal. We emphasize that ν s is used to stabilize the exactly incompressible deformations of the structure, not to model a compressible material. Our tests clearly demonstrate that the numerical results obtained with the proposed stabilization are in good agreement with benchmark results from a fully incompressible solid mechanics formulation.
The regularized delta function in the IB method can be interpreted as a weighting function. In this sense, the IB method resembles particle based and mesh free methods, such as the element-free Galerkin method (EFG) [15] and RKPM [8], which use weighting functions to reconstruct data at given nodes. In fact, both the regularized delta function and the EFG and RKPM weighting functions are not interpolatory at the nodes. In the EFG method, it is known that if the support of the weighting functions is small, volumetric locking may result. A simple way to alleviate volumetric locking for hexahedral elements, known as selective reduced integration (SRI) [16], is to integrate the volumetric term associated with κ using a quadrature rule with reduced order of accuracy. A procedure equivalent to SRI has been shown to fix this issue [17]. We demonstrate herein that volumetric locking also occurs in our implementation of the IB method if we let ν s → 1 2 . As already explained, the proposed method can circumvent issues with locking and obtain accurate solutions by using values of ν s much smaller than 1 2 , even with low order elements. Although the proposed stabilization can improve simulations in which the stress tensor is not deviatoric, we also demonstrate that using a traceless Cauchy stress tensor leads to smaller errors in volume conservation. We analyze two different strategies to achieve a traceless stress: the first is based on the Flory decomposition of the deformation gradient tensor [18], and the second eliminates the volumetric contribution of the stress tensor using a deviatoric projection. Whereas the Flory decomposition is mainly used for hyperelastic materials, the deviatoric projection strategy is easily implemented for general elastic and hypoelastic materials. The Flory decomposition for hyperelastic materials is equivalent to a formulation that additively decouples the isochoric (volume-preserving) and dilatational (volume-changing) parts of the structure's energy. Such decompositions require making physical assumptions about the solid being studied, namely that uniform pressure only results in a change in size and does not result in changes in the structure's shape [19].
To test the performance of the proposed stabilization method to preserve incompressibility, we use standard quasi-static benchmarks with different material models. The first two tests, Cook's membrane [20] and the compressed block [21], are two-dimensional problems that invoke the plane-strain assumption. The final two tests, an anisotropic extension to Cook's membrane [22] and a torsion test [23], are fully three-dimensional. In each of these benchmarks, it will be shown that not using the proposed stabilization, corresponding to ν s = −1 (equivalently κ s = 0), leads to unphysical deformations. In addition, compared to other choices of the structure's strain energy functional, omitting the stabilization leads to the worst errors in the volume conservation. Furthermore, it is generally the case that using the Flory decomposition with a finite choice of κ s performs the best.

Continuous Equations of Motion.
We briefly outline the IB equations of fluidstructure interaction. In our formulation, Ω = Ω f t ∪ Ω s t is the computational domain, in which Ω f t and Ω s t are the disjoint regions occupied by the fluid and the structure at time t, respectively. We describe the reference configuration of the structure using Lagrangian reference coordinates X ∈ Ω s 0 , in which Ω s 0 is the solid domain at time t = 0. We describe the computational domain using Eulerian coordinates x ∈ Ω. We use the mapping χ(X, t) : Ω s 0 → Ω s t to connect the reference configuration of the structure to its configuration. The simplest version of the IB formulation defines the Cauchy stress on the full domain to be Because we use a Lagrangian description of the structure, it is convenient to use the first Piola-Kirchhoff stress P s to describe the elastic response of the structure. Let F = ∂χ ∂X be the deformation gradient, and let J = det(F). The first Piola-Kirchhoff stress is related to the corresponding Cauchy stress by σ s = 1 J P s F T . For the material models considered here, P s is determined from a strain energy functional Ψ(F) via P s = ∂Ψ ∂F . We consider a Newtonian fluid, which has a stress tensor given by σ f = µ ∇u + ∇u T + pI, in which I is the identity tensor. The IB form of the equations of motion, as derived by Boffi et al. [4], is: Here, u is the Eulerian velocity, U is the velocity of the structure, p is the pressure, ρ is the constant mass density, µ is the viscosity, f is the elastic body force of the immersed solid, and N is the outward unit normal along the solid boundary ∂Ω s 0 in the reference configuration. We use lower case variables (x and u) to denote spatial or Eulerian quantities and upper case variables (X and U ) to denote material or Lagrangian quantities. The operators ∇ 2 , ∇·, and ∇ are with respect to spatial coordinates, and D Dt = ∂ ∂t + u · ∇ is the material time derivative. The differential operator ∇ X · indicates the divergence with respect to the Lagrangian coordinates.
The equations (2) and (3) are the Navier-Stokes equations for an incompressible Newtonian fluid augmented by elastic forces in the solid region. In the IB formulation, interactions between the Eulerian and Lagrangian variables are communicated via integral transforms with Dirac delta function kernels, equations (4) and (5). The viscosity present in the system causes u to be continuous, and this fact, paired with (5), enforces the no-slip condition on the fluid-solid interface. Equations (3) and (5) imply that the solid motion is exactly incompressible. To demonstrate this, let ω s 0 ⊆ Ω s 0 be a subregion of the solid domain in the reference configuration, and let ω s t = χ(ω s 0 , t) be the current configuration of this subregion at time t. The volume of this subregion in the current configuration is ω s t dx, and its volume in the reference configuration is ω s 0 dX. Because ∂χ ∂t = U and ∇ · u = 0, the Reynolds transport theorem implies: which means the volume of material region ω s 0 does not change in time. Because ω s 0 ⊆ Ω s t is arbitrary, it also holds pointwise. It is also useful to recall that ω s t dx = ω s 0 J(X, t) dX, which requires that J(X, t) ≡ 1 for an incompressible solid. Deviations from J = 1 indicate a change in volume in compressible and nearly incompressible cases.
In practice, we solve a weak form of equation (4) that is amenable to discretization via standard nodal FE methods. Rather than prolonging the divergence of the stress as in equation (4), we instead prolong a force F that is weakly equivalent to ∇ X ·P s . To determine F , we introduce arbitrarily smooth test functions V (X) and require for all smooth V (X). From the divergence theorem, we obtain: We then define f (x, t) = Ω s 0 F (X, t)δ(x − χ(X, t)) dX. By construction, the deviatoric part will satisfy the property tr(dev[T]) = 0. In continuum mechanics, the Cauchy stress may be similarly decomposed: in which the pressure p is the scalar φ. When motions are incompressible, p is a Lagrange multiplier that enforces incompressibility. For compressible motions, pI encodes the material's volume change. In practice with the IB method, however, the following decomposition of the Cauchy stress is used: in which σ v = µ ∇u + ∇u T is the viscous stress and σ s is not necessarily deviatoric. We use dev[σ v ] for added clarity, even though σ v is already deviatoric due to equation (3). Although (10) and (11) are equivalent in the continuous case, using (11) in the discretized equations may lead to disastrous effects. We wish to introduce a stabilization to the structure's stress that corrects for the loss of incompressibility resulting from discretization. We first take the deviatoric component of σ s and then introduce a volumetric stabilization, such that the stress is in which p stab I is a stabilization term that acts like a pressure in the solid region. Such volumetric stabilization can also be included if the deviatoric solid stress is not considered: Similar to treatments of nearly incompressible elasticity, we define p stab as a volumetric penalization term. More specifically, p stab is derived from a volumetric energy U (J) that depends only on the structure's changes in volume: Thus, this formulation for stabilization parallels models of nearly incompressible elasticity. When defining U (J) in nearly incompressible elasticity, restrictions are placed on U (J) to achieve certain physically motivated properties. Specifically, it is necessary that this term satisfy U (1) = 0, so that the identity motion F = I introduces no extra energy. The contribution of the volumetric energy to the first Piola-Kirchhoff stress is JU (J)F −T , so it is beneficial for U (J) also to satisfy U (1) = 0. This implies no extra stress is introduced if J = 1. We also require that lim J→∞ U (J) = ∞ and lim J→0 U (J) = ∞ so that large dilatations and contractions are also energetically unfavorable. Finally, we want to control the effect of U (J) through the numerical bulk modulus κ s , such that κ s → ∞ represents exact incompressibility. A simple example of U (J) that satisfies the above conditions is: To modulate the κ s , we introduce the numerical Poisson ratio ν s . The two parameters are related via in which µ s refers to the linearized (F = I) shear modulus. Note that ν s = −1 yields p stab = 0, retrieving either (10) or (11). The previously mentioned models, (10) - (13), are equivalent in the continuous case. It is interesting to study all these formulations, though, because the discrete coupling operators are not guaranteed to preserve the discrete divergence of the Eulerian velocity. Therefore in the discretized equations, we may lose discrete incompressibility of the solid even if we maintain a discretely divergence-free Eulerian velocity field, and (10) - (13) are no longer equivalent.
2.2.1. Umodified Model. The models in which the solid's Cauchy stress is not necessarily deviatoric are derived from the following two energy functionals, which omit and include additional volumetric stabilization, respectively: The energy given by (17) yields the model given by (12), and (18) yields the model given by (13). Note that (17) is the formulation which may lead to the aforementioned disastrous effects.

Modified Model.
It is possible to obtain models (10) and (12) in different ways. One way is through the Flory decomposition F = J −1/3 F [18]. Note that det(F) = 1 by definition. We reformulate strain energy functionals (17) and (18), respectively, as: The energy given by (20) completely decouples energy associated to volume changing and volume preserving motions and achieves the desired split in the Cauchy stress. This decoupling is motivated by the physical assumption that a uniform pressure only produces changes in size but not changes in shape. Work by Sansour explores this physical assumption in depth [19], but we offer a brief explanation here. Specifically, with (20), we obtain an additive split in the Cauchy stress into purely deviatoric and dilatational stresses. This means that the only contributions to the stabilizing pressure will come from U (J).
We show that using the model with the Flory decomposition has a similar effect as using the deviator operator.
∂F is a fourth order tensor. Explicitly: with I denoting the fourth order identity tensor. By contracting (21) with ∂W ∂F , we obtain: Pushing forward (22), it is evident that using the Flory decomposition yields a traceless Cauchy stress:

Deviatoric Projection.
Another way to achieve a deviatoric Cauchy stress is through simply using the deviator operator (9). We refer to these stress models as deviatoric projections. In our tests, deviatoric projections for hyperelastic models will be constructed by using the deviator operator for the first Piola-Kirchhoff stress: Note that (24) resembles (22) with the exception of the J −1/3 pre-factor and that pushing forward (24) also yields a traceless tensor. The operator is applied to (17) and yields a Cauchy stress with the desired split. Herein, we study these models with and without the volumetric stabilization.

Constitutive Laws.
For incompressible isotropic solids, writing Ψ as a function of I 1 = tr(C) and I 2 = 1 2 (I 2 1 − tr(C 2 )) allows for material frame indifference. I 1 and I 2 are the first two tensor invariants of the right Cauchy-Green tensor C = F T F [24]. In incompressible cases, J ≡ 1, so there is no dependence on the third invariant, I 3 = det(C) = J 2 . In cases of compressibility, however, this is not the case; J ≡ 1, so the structure will undergo volumechanging deformations. Often, the energy for compressible materials is written as a function ofĪ 1 = J −2/3 I 1 andĪ 2 = J −4/3 I 2 , which are the invariants of C = F T F. By modifying the invariants in this way, we remove information about the volume change. Thus we refer to the invariants of C as the modified invariants.
Often, the Young's modulus is used when describing neo-Hookean models. To relate the Young's modulus E to µ s we use the formula µ s = E 2(1+ν) . Here we use ν = 1 2 because we are modeling a material whose motions are incompressible.
When using modified invariants the energy and stress are, respectively: Taking the deviatoric projection of (26), we have

2.3.2.
Mooney-Rivlin Models. If we add dependence on I 2 in a simple way, we acquire the Mooney-Rivlin material law. The unmodified invariant case is given by: in which c 1 and c 2 are material constants. By using modified invariants we have: Taking the deviatoric projection of (31), we have As with the noe-Hookean models, we study material models described in this work with and without volumetric stabilization. For the Mooney-Rivlin material law, this will require being able to relate material constants to κ s . For consistency between the small deformation (linear) and large deformation (nonlinear), we set µ s = 2(c 1 + c 2 ) when calculating κ s . This allows the use of the same formula, equation (16), that relates κ s and ν s to a material quantity.

Modified Standard Reinforcing Model.
To examine the effects of anisotropy, we use the modified standard reinforcing model [25]. This model is used to describe transversely isotropic materials with fibers given by a material vector A in the reference configuration and a = FA in the current configuration. The effect of the anisotropy appears through the pseudo-invariants I 4 and I 5 : in which B = FF T is the left Cauchy-Green strain. Because a is the stretched and rotated material vector, I 4 measures the stretch of the fiber, whereas I 5 encodes information related to the shear as well as the stretch [26]. The modified standard reinforcing model has the following strain energy functional and first Piola-Kirchhoff stress: Here, µ T is the shear modulus of the material in the plane transverse to the fibers, and µ L is the shear modulus along the length of the fibers. To determine κ s , µ T is used in equation (16) because this material model does not involve an isotropic shear modulus µ s . The parameter E L is similar to a Young's modulus but in the direction of the fiber. When we modify I 1 , we instead obtain: Likewise, for the modified standard reinforcing model, the deviatoric projection is: The anisotropic models considered herein take the following forms: with the exception of the deviatoric projection of the standard reinforcing model, which does not arise from an energy functional. If we modify both I 4 and I 5 by using C and use the volumetric stabilization, we arrive at a Cauchy stress with an additive deviatoricspherical split. However, in cases of uniform pressure, it is possible that a body will undergo a shape-changing deformation if the material is anisotropic. Thus, the volumetric split is not appropriate for the anisotropic part of the stress.
We remark that in biomechanics literature, the standard reinforcing model is often defined as Ψ = c 1 (I 1 − 3) + c 4 (I 4 − 1) 2 , without any dependence on I 5 [26]. It can be shown that omitting this pseudo-invariant implies that the linearized shear moduli in the direction of the fibers and perpendicular to the fibers must be the same. It can also be shown that the three modes of shear characteristic of transversely isotropic materials are not represented if I 5 is omitted [25]. The modified standard reinforcing model is arrived at by augmenting the standard reinforcing model in a way that allows for the consistency between the linear and finite regimes [26].

Numerical Methods
The numerical methods we use in this study consist of those for fluid-structure interaction and those for nearly incompressible structural mechanics, which we compare against our FSI results. The focus of this work is the volumetric stabilization to the numerical method of Griffith and Luo [6], and so the the descriptions of both methods are only outlined.
For FSI simulations we use IBAMR [27,28] which is an open-source adaptive and distributedmemory parallel implementation of the IB method. Specifically, we use the IBFE module in IBAMR, which allows the use volumetric structures. The quasi-static finite element benchmark solutions are computed using BeatIt [29]. Both IBAMR and BeatIt heavily rely on the parallel C++ finite element library libMesh [30] and on PETSc [31] linear solvers.
3.1. Fluid-Structure Interaction. The methods used for FSI are described in detail in two previous works [6,32]. Briefly, a staggered-grid finite difference method is used to discretize the Eulerian equations, and a nodal FE method is used to discretize the Lagrangian equations [33]. In this scheme, the Eulerian velocity u is approximated at the cell edges (faces in three spatial dimensions), and the pressure p is approximated at cell centers. We use standard second order accurate finite differences to discretize the Eulerian incompressible Navier-Stokes equations (2) and (3) [32,34]. In our computations, we use the unified formulation of the hyperelastic IB method [6].
We discretize the structure Ω s 0 via a triangulation T h = ∪ e K e , in which K e are the elements. On T h , we define Lagrangian basis functions {φ (X)} m =1 , in which m is the number of FE nodes in our triangulation. These functions belong to the common FE spaces of P1, P2, Q1, and Q2, which in two spatial dimensions denote the spaces of linear, quadratic, bilinear, and biquadratic basis functions, respectively; see [33]. In three dimensions, we only use P1 and Q1, which are the spaces of linear and trilinear basis functions. The mapping χ is approximated by χ h (X, t) = m =1 χ (t)φ (X). The Lagrangian force is also be approximated with the same finite element basis via F h (X, t) = m =1 F (t)φ (X). Gaussian quadrature rules are used to integrate the integral equations [14]. For each case, we use integration orders which exactly integrate the basis functions of our chosen FE space; we do not employ selective reduced integration.
Coupling between the Eulerian and Lagrangian variables is mediated by discretized integral equations with regularized delta function kernels. The regularization of the delta function is one place in the scheme where the solid's incompressibility is typically lost because the kernel functions used in practice generally have a non-zero gradient. Including a volumetric energy term in the solid's strain energy can help to overcome spurious volume changes that result from this construction. We use a regularized delta function of the form in which d is the spatial dimension. Specifically, we use the four point smoothed delta function of Peskin [3].
The prolongation and restriction operators of the coupling are constructed to be formal adjoints of one another. This means that there is conservation of power as we map data between the Eulerian grid and Lagrangian mesh. In two spatial dimensions, our discrete prolongation operator is: in which X e Q are quadrature points and w e Q are quadrature weights. We denote by x i− 1 2 ,j and x i,j− 1 2 the Eulerian grid point at (i∆x, (j + 1 2 )∆x) and ((i + 1 2 )∆x, j∆x), respectively. ∆x is the Eulerian grid spacing. (f 1 ) i− 1 2 ,j is the first component of the discrete Eulerian force evaluated at the corresponding Eulerian grid point, and (f 2 ) i,j− 1 2 is similarly defined. The discrete restriction operator is: in which (u 1 ) i− 1 2 ,j and (u 2 ) i,j− 1 2 are the components of the discrete Eulerian velocity field evaluated at their respective points. U IB is a Lagrangian velocity field that is in general not a sum of the basis functions {φ(X)} m =1 . To attain a velocity field U h that is a sum of the basis functions, U IB is projected onto {φ(X)} m =1 in an L 2 sense. For more details, we refer to the description of the method by Griffith and Luo [6].
All Dirichlet boundary conditions for the structure are imposed via a penalty method. Specifically, surface forces of the form T = κ D (χ − χ D ) are applied to the structure's boundary to approximate Dirichlet boundary conditions given by χ D . κ D denotes a stiffness used to penalize deviations from the desired value. We use the scaling κ D ∝ ∆x (∆t) 2 , in which ∆t is the time-step size, so that the the stiffness parameter increases as the mesh is refined. For all tests, we use a pre-factor of 2.5 for κ D .
3.2. Structural Mechanics. The results from the FSI calculations are compared to quasistatic fully incompressible elasticity finite element simulations. Specifically, we use a mixed displacement-pressure formulation to enforce the incompressibility constraint [35]. Piecewise linear polynomials are used to interpolate the displacement and pressure fields. This choice leads to an unstable numerical method, as the corresponding P1/P1 elements (indicating a piecewise linear approximation for both the displacement and pressure) do not satisfy the Ladyzhenskaya-Brezzi-Babuska (LBB) condition (also called the inf-sup condition) [33]. In 1986 Hughes, Franca, and Balestra circumvented the LBB condition for the Stokes equations using a simple stabilization method suitable for equal order approximations and recovering the optimal order of convergence [36]. Their method has been extended to linear and nonlinear elasticity [37,38], and it has been reinterpreted in a variational multiscale framework [39,40]. We use this stabilization method to guarantee the stability of the finite element simulations.

Benchmarks
We use standard benchmark problems for incompressible elasticity drawn from the solid mechanics literature, except that here the solid bodies are embedded in an incompressible Newtonian fluid. However, because the elastic part of the structural material response is hyperelastic, and thus path-independent, the steady states of the FSI problems are the same as those from pure solid mechanics formulations. In fact, the IB formulation used herein treats the solids as visco-elastic because the stress within the solid is σ = σ f + σ s . Recall that the fluid stress has the form σ f = µ ∇u + ∇u T + pI in the continuous case, so the equilibrium configuration defined by ∇ · σ = 0 and zero fluid velocity implies ∇ · σ s + ∇p = 0 as well.
We report the computed displacements of a point of interest and the total volume conservation for each benchmark. When listing the ranges for the volume conservation, we omit the coarsest discretizations of the structure. We also present the deformations of each benchmark for all relevant formulations of the energy functional (e.g. equations (18) - (19) for the isotropic models). These visualizations depict the average values of J for each element, and the extents of the color bar indicate cutoff values. The average of J for element e is calculated as follows: We use the same Gaussian quadrature rule here as for the approximation of the integrals in other parts of our method. We use material models for the structure with modified isotropic invariants, unmodified isotropic invariants, and the deviatoric projection of the isotropic part of the elastic stress. These models are studied with varying levels of volumetric stabilization that are tuned In all tests, the computational domain is Ω = [0, L] d , in which d = 2, 3 is the spatial dimension and L is the domain length. The spacing of the Eulerian grid is ∆x = L N , in which N is the number of cells in one dimension. We use zero velocity boundary conditions on the computational domain. This allows for the fluid velocity to settle down to zero as time goes on. The simulations are run until a final time T f , which is chosen such that the velocity is approximately zero. The loads on the structure are incrementally applied; at time t = 0 the load is zero and linearly increases in time until at time T l = αT f the load is fully applied. Here α ∈ (0, 1). Between times T l and T f , we let the structure relax to its resting configuration. Except where otherwise noted, the density is ρ = 1.0 kg m 3 , and the viscosity is µ = .16 N·s m 2 . Let ∆X be the Lagrangian mesh width. The mesh factor ratio M FAC = ∆X ∆x describes the relative grid spacing between the Eulerian and Lagrangian meshes. In our tests, the choice of M FAC = 1 is used for pressure driven cases and M FAC = 2 is used for shear driven cases. These choices of M FAC were made based on preliminary tests (data not presented). More specifically, we use M FAC = 1 for the compression block test and M FAC = 2 for the Cook's membrane, the anisotropic Cook's membrane, and the torsion tests. Results for this benchmark are summarized in figures (3) and (4). As shown in figure (3), most cases converge to the benchmark solution. With values of ν s close to 1 2 , the solution exhibits volumetric locking as in a typical displacement based FE formulation. Note the unmodified case in the final row of plots in figure (3) and that the displacement of the single point of interest converges under mesh refinement for these cases. However, it can be observed that qualitatively the deformation of the entire mesh is unphysical; see figure (2D). We emphasize that the case of unmodified invariants and zero volumetric energy are the only cases in which this unphysical behavior is observed. Figure (4) shows the percent change in the total area of the mesh after deformation. It is clear that the modified invariants and deviatoric projection yield improved results in terms of global area conservation in comparison to the unmodified invariants. This effect becomes more pronounced as the bulk modulus is decreased. It may appear as though the modified invariants and deviatoric projection cases have zero volume change, but this is not the case. The percent change in total volume for all elements considered ranged between .10% and .000021% for modified invariants and between .10% and 0% for the deviatoric projection. For the unmodified cases, this range was between 7.45% and .000065%. These ranges account for change in area in an absolute sense, whereas the plots display whether the change in area was a gain or loss. The quantity of interest is y-displacement at the center of the top face. Figure (5) reports these values at time T f = 100 s. The load time is T l = 40 s. Again, note the performance of the unmodified invariants in the final row of plots in figure (7). As before, the convergence behavior of the single recorded point is satisfactory although the overall deformations are unphysical in these cases; see figure (6). Particularly noticeable in this benchmark is the effect of using modified invariants versus unmodified invariants while using a nonzero bulk modulus. This is apparent in figures (6A) and (6B), where the deformations of the elements are smoothest in (6A) in the case in which modified invariants are used. As expected for values of ν s close to 1 2 , volumetric locking plagues the lower order elements, resulting in poor convergence. However, locking behavior is avoided for different values of ν s , corresponding to smaller bulk moduli, even for low order elements. Figure (8) reports the percent change in total volume. Modified invariants and deviatoric projection yield superior volume conservation. The percent change for all element types considered ranged between 2.1% and .0004% for the modified invariants, between 14% and .001% for the unmodified invariants, and between 2.4% and .0005% for the deviatoric projection.

Modified
Unmodified Avg J 0.95 1.01 Figure 2. Deformations of Cook's membrane using a neo-Hookean material model, equations (25) - (28). The first row is depicts cases with ν s = .4, and the second row depicts cases with ν s = −1 (equivalent to κ s = 0 and no volumetric-based stabilization). The first column depicts cases with modified invariants, and the second column depicts cases with unmodified invariants. Notice that the case with modified invariants with nonzero bulk modulus have the smoothest deformations whereas those of the case with unmodified invariants and zero bulk modulus behave unphysically.

4.3.
Anisotropic Cook's Membrane. This benchmark involves a fully three-dimensional and anisotropic Cook's membrane; see figure (9). It is similar to and based upon one studied by Wriggers et al. [22].  As in the other cases considered, the behavior for the case of zero volumetric penalization with unmodified invariants yields unphysical deformations. In this case, the poor behavior is located at one of the corners on the face where the traction is applied; see figures (10) and (11). Specifically, the element at this location collapses; two of the FE nodes are approximately in the same location. This was also the location where the y-displacement was measured for the plots in figure (12). Finally, as in the other cases considered, the case  of unmodified invariants is associated with worse volume conservation. Figure (13) depicts the percent change in total volume for these cases. The percent change for all element types considered ranged between 2.2% and .0087% for unmodified invariants. For the modified invariants and the deviatoric projection, both ranges were .11% and .00014%. Unlike the other cases considered, the computation with zero volumetric penalization seems to perform nearly as well as or better than the case with volumetric penalization; see figure (12). However, for the modified invariants, figures (9A) and (9C) show that using a nonzero bulk modulus produces a more uniform distribution of J and one that is closer to J = 1. Overall, the differences among all cases in the results presented for this test are fairly

Modified
Unmodified Avg J 0.90 1.05 Figure 6. Deformations of the compressed block test using a neo-Hookean material model, equations (25) - (28). The first row is depicts cases with ν s = .4, and the second row depicts cases with ν s = −1 (equivalent to κ s = 0 and no volumetric-based stabilization). The first column depicts cases with modified invariants, and the second column depicts cases with unmodified invariants. Again, notice that the case with modified invariants with nonzero bulk modulus has the smoothest deformations whereas those of the case with unmodified invariants and zero bulk modulus behave unphysically. Figure (15) show the computed displacements for modified invariants, unmodified invariants, and the deviatoric projection as well as for different values of the numerical Poisson ratio. The cases of unmodified invariants and zero bulk modulus lead to the most extremely unphysical deformations in all benchmarks studied; see figure (15 D). Figure (16) shows the displacement in the y-direction at the center point of the twisted face. Notice that in these plots, the cases with unmodified invariants and zero volumetric energy clearly delineate themselves from other cases. Unique to this test, the convergence of the computed displacement of this case is not deceptive; the convergence is poor and the deformations are also poor. Additionally, the effect of volumetric penalization is more drastic in this benchmark: the percent change in volume is generally much larger than the previous tests; see figure (17).  Figure 9. The anisotropic Cook's membrane benchmark. Traction in the y-direction is applied to the smallest face, and the opposite side is kept fixed. The quantity of interest is y-displacement as measured at the encircled point.
Specifically, the range of percent change for all element types considered is between 11% and .16% for the modified invariants, between 93% and 8.5% for the unmodified invariants, and between 61% and 1.5% for the deviatoric projection. The choice of numerical Poisson ratio also has a large effect on the displacement of the twisted face, which can be seen in figure (16). The differences between displacement curves for with and without volumetric penalization is more apparent, with the case of volumetric penalization performing much better here. We contrast that with the anisotropic Cook's membrane benchmark, which has only slight differences between the displacement curves for ν s = .4 and ν s = −1.
Finally, in this benchmark the deviatoric projection delineates itself from the modified invariants. However, in this test both the volume conservation and the displacement performed worse for the deviatoric projection than the modified invariants; see figures (16) and (17). For the other benchmarks, there was a negligible difference between the solution produced by the modified invariants and that produced by the deviatoric projection.

Discussion and Conclusion
In this work, we propose simple stabilization methods for the hyperelastic IB method that correct for changes in volume that the immersed structure may exhibit in the discretized equations. Mimicking nonlinear solid mechanics formulations, our strategy is to append a volumetric penalization term to the structure's stress. Additionally, we explore the effect of this term when added to a solid Cauchy stress which is deviatoric and find that use of a deviatoric Cauchy stress enhances the effect of the volumetric stabilization. Although these changes are standard in solid mechanics, the hyperelastic extension of the IB method does not incorporate them [4,6]. Prior benchmark studies in the solid mechanics literature have indicated that this decomposition is important for accurate simulations, especially in the nearly incompressible or fully incompressible limits. To date, however, the effect of the formulation of the elastic stress on the accuracy of the IB method does not appear to have been systematically studied.
In this work, we describe two ways to achieve a deviatoric Cauchy stress in the solid region: reformulation of the structure's strain energy functional in terms of modified invariants and use of the deviatoric operators (9) and (24) on the structure's stress. As presented here, the volumetric stabilization to the solid's stress comes from a term added to the solid's energy

Modified
Unmodified Avg J 0.85 1.05 Figure 10. Deformations of the compressed block test using the modified standard reinforcing model, equations (37) - (40). The first row is depicts cases with ν s = .4, and the second row depicts cases with ν s = −1 (equivalent to κ s = 0 and no volumetric-based stabilization). The first column depicts cases with modified invariants, and the second column depicts cases with unmodified invariants. We remark that the case with unmodified invariants and zero bulk modulus leads to collapsed element on the face where the traction is applied.

Modified Unmodified
Avg J 0.85 1.05 Figure 11. Deformations of the anisotropic Cook's membrane with zero and finite volumetric energy from a different view; see figure (10). The collapsed element in case (D) is clearly visible here.
functional. From this point of view, this work explores different forms of the energy functional of an incompressible elastic body on the accuracy of the IB method for benchmark problems of incompressible nonlinear elasticity. The effect of the volumetric penalization is explored through varying the numerical Poisson ratio. Setting ν s = −1 describes the situation of no volumetric penalization. As detailed herein, the IB method is exactly incompressible in the continuum limit, so the actual Poisson ratio of the immersed structure will automatically be ν = 1 2 . Thus, including a volumetric energy term makes no difference in the continuous equations. Upon discretization, however, incompressibility is exactly maintained for the Eulerian velocity field only. Nonetheless, even in the discrete case, it is assumed that the solid does inherit some incompressibility from the Eulerian variables. This explains why, if a deviatoric stress tensor is used, even for the case of ν s = −1 (κ s = 0) the structure's total volume remains mostly unchanged. If the solid stress is not fully deviatoric, we show that poor results can be expected if no volumetric  penalization is added. The numerical Poisson ratio thereby acts as a stabilization parameter in the discrete IB equations, rather than a material parameter. Further, as is desirable for a stabilization parameter, the effect of the numerical Poisson ratio vanishes in the case that the structural deformations are exactly incompressible. For this work we use four common benchmark problems found in structural mechanics literature. The two dimensional problems, Cook's membrane and the compressed block, are used to perform an in depth study of the proposed volumetric stabilization. For these problems we are able to use more element types and more choices of ν s . Additionally, we also use three-dimensional problems, an anisotropic extension to Cook's membrane and a torsion test. All benchmarks point broadly to the conclusions that some form of volumetric stabilization performs best. More specifically, the results suggest that a numerical Poisson ratio of ν s = .4 yields results as accurate as a fully incompressible FE approach.
Similar to low order displacement-based FE formulations, as ν → 1 2 , our method suffers  Reducing the Poisson ratio is clearly demonstrated to be a simple fix for these tests, in that the convergence behavior is improved for all elements, including the second order elements. Additionally, this change yields only a negligible sacrifice in total volume conservation. As indicated by figures (2), (6), (10), and (15), the reduction in pointwise volume conservation is also negligible. This change in numerical Poisson ratio is justified because ν s is a numerical parameter. Specifically, and unlike the typical situation in nearly incompressible elasticity, changing the numerical Poisson ratio does not affect the limiting material model obtained under grid refinement. The effect of anisotropy is explored through study of the anisotropic Cook's membrane benchmark. Broadly, the same conclusions which hold true for the other benchmarks also  Figure 14. Diagram depicting how the torsion is applied to the column. The face opposite the applied torsion is kept fixed. The quantity of interest is the y-displacement as measured at the encircled area.
hold for this benchmark. That is, unmodified isotropic invariants paired with zero volumetric penalization yield unphysical deformations. Introducing anisotropy opens the door for many more possible forms for the elastic energy functional such as modifying the pseudo-invariants or including the pseudo-invariants in the volumetric term (only for nearly incompressible materials). Though we argue that the use ofĪ 4 andĪ 5 is not appropriate for the cases considered here, the study of their use in an IB framework may still prove interesting in future work.
The effect of the proposed changes is most pronounced for the torsion test. Removing the volumetric term often has disastrous effects when unmodified invariants are used, and again the torsion test produces the most noticeable effects. The volume conservation is particularly poor, and the displacement of the point under study converges slowly. Even for relatively fine discretizations of this benchmark, our method demonstrates improvements in volume conservation of up to 59% when compared to the unstabilized case with unmodified invariants.
In closing, we emphasize that using elastic energy functionals with unmodified invariants and zero volumetric penalization, given by equations (17) and (45), consistently yield unphysical deformations in our IB simulations. Further, such choices of energy functional produces noticeably worse volume conservation. Going even further, modified invariants generally perform better than unmodified invariants in all our results. Including a volumetric penalization term also generally produces more accurate results, with the exception of the anisotropic Cook's membrane in which the differences in performance are minimal. Given its minimal computational cost and considerable improvements to the accuracy of the computed displacements and overall volume conservation, this study strongly indicates that volumetric energy-based stabilization with modified invariants should be the default structural formulation for IB-type methods similar to that used herein.

Modified
Unmodified Avg J 0.75 1.05 Figure 15. Deformations of the torsion test using the Mooney-Rivlin material law, equations (30) - (33). The first row is depicts cases with ν s = .4, and the second row depicts cases with ν s = −1 (equivalent to κ s = 0 and no volumetric-based stabilization). The first column depicts cases with modified invariants, and the second column depicts cases with unmodified invariants. Note the extremely poor deformations for unmodified invariants with no volumetric-based stabilization in case (D).
Carolina at Chapel Hill through the Research Computing division of UNC Information Technology Services.  Omitting the coarsest discretizations (m = 65), the largest deviations in total volume among all element types used are approximately 11% for the modified case, 93% for the unmodified case, and 61% for the deviatoric case.