{'fraction_non_alphanumeric': 0.03800256380375248, 'fraction_numerical': 0.01683952919240182, 'mean_word_length': 4.655133781468301, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 22, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A unified continuum and variational multiscale formulation for fluids, solids, and fluid-structure interaction

We develop a unified continuum modeling framework using the Gibbs free energy as the thermodynamic potential. This framework naturally leads to a pressure primitive variable formulation for the continuum body, which is well-behaved in both compressible and incompressible regimes. Our derivation also provides a rational justification of the isochoricvolumetric additive split of free energies in nonlinear elasticity. The variational multiscale analysis is performed for the continuum model to construct a foundation for numerical discretization. We first consider the continuum body instantiated as a hyperelastic material and develop a variational multiscale formulation for the hyper-elastodynamic problem. The generalized-α method is applied for temporal discretization. A segregated algorithm for the nonlinear solver, based on the original idea introduced in [107], is carefully analyzed. Second, we apply the new formulation to construct a novel unified formulation for fluid-solid coupled problems. The variational multiscale formulation is utilized for spatial discretization in both fluid and solid subdomains. The generalized-α method is applied for the whole continuum body, and optimal high-frequency dissipation is achieved in both fluid and solid subproblems. A new predictor multi-corrector algorithm is developed based on the segregated algorithm. The efficacy of the new formulations is examined in several benchmark problems. The results indicate that the proposed modeling and numerical methodologies constitute a promising technology for biomedical and engineering applications, particularly those necessitating incompressible models.


Introduction
Continued advancement in the variational multiscale (VMS) method for computational fluid dynamics (CFD) [5,59,131], multiscale boundary conditions to model the distal vasculature [126], numerical optimization methods [84,129], novel coupling procedures for fluid-structure integration (FSI) [7,30,92], and new solver technologies [90,91] have led to increasingly sophisticated simulation technologies of three-dimensional patient-specific cardiovascular problems [83,85,117]. These simulation methods have been applied to investigate a wide range of cardiovascular problems with increasing clinical utility, such as coronary artery disease [100], aneurysms [72], and congenital heart disease [130]. These advances have also led to new open problems and a call for new computational technologies in biomedical modeling [118]. In particular, there is a pressing need to accurately predict transmural stresses for nonlinear, anisotropic, nearly incompressible, viscoelastic materials in the setting of FSI. This technology will benefit the investigation of long-term vascular growth and remodeling driven by mechanical forces and mechanobiological response. This work represents our first step towards developing a robust, stable, accurate, and efficient finite element technology to address the aforementioned need. We focus on (1) constructing a unified modeling framework for solid and fluid dynamics, (2) developing a new numerical formulation to handle both compressible and incompressible isotropic hyperelastic materials in a consistent manner, and (3) formulating a new coupling procedure for FSI problems. It is worth mentioning that the technology proposed in this work, although motivated by biomedical problems, is also applicable to a wide range of more general engineering problems. In this section, we will first review traditional methods for nonlinear incompressible solid dynamics and then introduce our new formulation based on a new continuum basis. After that, we will discuss a new coupling approach for FSI problems. Lastly, we will provide an outline of the body of this work.

Continuum mechanics and numerical methods for solid dynamics
Incompressible solids refer to materials that undergo very small volume changes under mechanical loading. This property is ubiquitous in elastic-plastic materials, rubber-like materials, and biomaterials. From a numerical point of view, incompressibility poses a strict restriction on the choice of suitable methods. The classical treatment of incompressibility boils down to the use of a displacement-pressure element pair that satisfies the celebrated Ladyzhenskaya-Babuška-Brezzi (LBB) condition [12]. An improper choice of elements may lead to mesh locking and spurious pressure modes [48,Chapter 4]. In practical calculations, the auxiliary pressure variable in the mixed formulation is not favored due to the additional computational cost. The cost-effectiveness of an algorithm was emphasized during the early days of finite element method development. Consequently, there has been a strong motivation for developing displacement elements that circumvent the volumetric locking phenomenon in structural analysis. One possible solution is invoking higher-order finite element methods [43]. Nevertheless, low-order elements still enjoy significant popularity in nonlinear problems for at least two reasons. First, they are simple to implement and to use in conjunction with adaptive mesh generation; second, they are robust in nonlinear dynamic analysis. One popular approach for the low-order displacement elements is theB andF projection methods [23,28,47]. The idea of these projection methods goes back to the mean-dilatation approach [93] and the reduced and selective integration method [80]. In theB projection method, the dilatational part of the strain is projected onto an appropriate function space to alleviate the volumetric constraint. The equivalence of this approach to the mixed finite element is now well-understood [46,80]. TheF projection method is a generalization of theB method for large deformation problems [61,115]. Without a doubt, the Q 1 /Q 0 element, as the low-order quadrilateral/hexahedral element, is the most widely used element in structural analysis, with Q 1 for the displacement discretization space and Q 0 for the projected dilatational strain space. Interested readers are referred to [48,Chapter 4] for a comprehensive review on this topic. For problems with complex geometries, mesh generation poses an additional constraint on the choice of numerical methods. Automated hexahedral mesh generation still poses significant challenges for practical problems and remains a labor-intensive and time-consuming process [111]. On the other hand, the algorithms for tetrahedral mesh generation are mature [35,112] and are routinely used for complex problems [125]. However, the aforementioned projection methods cannot be directly applied to low-order simplicial elements, since the P 1 /P 0 element does not satisfy the LBB condition and will suffer from locking [50]. To overcome the above issues, particularly for problems with incompressible materials and complex geometries, it is imperative to enable low-order simplicial elements for nearly-incompressible materials. There are several approaches that have claimed success in this regard. TheF projection method is generalized by projecting the dilatational part of the deformation gradient onto a patch-wise constant space with predefined non-overlapping element patches [24]. This method requires one to group elements into patches, and hence introduces additional complication in mesh generation as well as in the matrix assembly. In the mixed-enhanced approach, the compatible strain fields are augmented by an additional strain field. In doing so, the volumetric locking can be avoided with linear pressure interpolation [119]. This method bears some similarity with the MINI element and the stabilized method based on residual-free bubbles. However, those methods seem to be not completely satisfactory. It is still common to see that practitioners in the area of soft tissue modeling make compromise between the incompressibility material property and the geometrical complexity.
Over the years, the stabilized methods have gained tremendous success in handling numerical instabilities in the finite element modeling of fluid flows [13,52]. Later, the VMS method was introduced to provide a rationale for the stabilized methods and a framework for designing new computational methods [5,49,51,59,60,106]. Among the stabilization techniques, the instability associated with the LBB condition was resolved by invoking a Petrov-Galerkin formulation and adding a perturbation term to the test function [52]. It was proven that many implementationally convenient elements, such as the equal-order simplicial elements, are convergent for the Stokes problem by employing this stabilized formulation. Since the Stokes equations are formally identical to the displacement-pressure formulation for the incompressible isotropic elasticity equations, this stabilization technique can be directly applied to the incompressible elasticity equations. It is worth noting that the Petrov-Galerkin method has also been developed for the stress-displacement formulation of nearly-incompressible elasticity based on the classical Hellinger-Reissner principle [17,33]. In [68,81,82,22], the stabilized method was applied to finite elasticity and plas-ticity by introducing a pressure stabilization term in an ad hoc way. In [87], the same static finite elasticity model was analyzed in the VMS framework. It should be pointed out that those formulations are restricted to static calculations and are not well-suited for dynamic fully-incompressible problems [107].
In recent years, there is a growing interest in extending the stabilized and VMS methods to solid dynamics with particular interests in nearly incompressible and bending dominated scenarios. In [70], the linear momentum equation was augmented by a dynamic equation for the deformation gradient to improve the stress accuracy. The authors designed a Petrov-Galerkin finite element formulation for such system to provide stabilization for nearly incompressible materials. In [36], J, the determinant of the deformation gradient, was added as one additional independent variable to enhance the performance in the incompressible limit. In [11], the authors introduced a computational framework for polyconvex hyperelasticity by treating the deformation gradient, its cofactor, and its determinant as independent kinematic variables. This framework was numerically investigated in [10,37] by invoking the entropy variables (conjugate stresses) of these kinematic variables. The polyconvex hypothesis guarantees the existence of entropy variables and leads to a new Hellinger-Reissner type variational principle. In general, this computational framework enjoys second-order accuracy for the stress and strain and shows robust performance in compressible, nearly incompressible, and fully incompressible regimes. A drawback associated with this formulation is that there are nineteen independent variables for describing kinematics in three-dimensional problems, which makes the algorithm quite noneconomical. In parallel to the aforementioned development, a pressure-rate equation was proposed in [107]. It was shown that this rate-type equation for the pressure field significantly improves the calculations of nearly incompressible elastodynamic problems. In [102], a VMS analysis was performed on this pressure-rate type formulation, and the authors coined it as "dynamic variational multiscale method". This formulation shows robust performance in nearly incompressible and bending dominated scenarios and has been recently generalized to viscoelasticity [132]. However, due to the usage of J as a variable in the model derivation, the authors had to treat the nearly incompressible and the fully incompressible cases separately [11,37,102,107]. More importantly, considering several different mixed formulations have been proposed in the aforementioned works, we feel the role of those equations needs to be elucidated. In our opinion, the pressure rate equation [102,107,132] or the equation for J [10, 36] is a differential mass equation. Consequently, one should judiciously use the classical relation ρJ = ρ 0 , and, in our opinion, it is redundant to solve ρJ = ρ 0 together with the differential mass equation in the discrete problem [107,132]. This point will be clarified in Section 2.
In traditional approaches, the constitutive relations of finite elasticity are derived based on the Helmholtz free energy or the strain energy [45]. In non-equilibrium thermodynamics, one can alternatively derive equivalent formulations based on other thermodynamic potentials. These thermodynamic potentials are related through Legendre transformations. Among all the thermodynamic potentials, the Helmholtz free energy is favored for the discussion of constitutive modeling, likely because it is a function of temperature, specific volume, and the number of molecules. These variables are easy to measure and control in laboratories. For example, it is easy to measure the temperature while it is very hard, if not impossible, to measure its conjugate variable, the entropy. However, it should be pointed out that there are certain drawbacks of formulations based on the Helmholtz free energy. For example, in the incompressible limit, the Helmholtz free energy ceases to serve as a valid thermodynamic potential, since one of the state variables, the density, is constrained to be a constant [79]. Therefore, one has to perform a Legendre transformation with respect to the specific volume to transform the independent variable to the thermodynamic pressure. The resulting thermodynamic potential serves as a proper starting point for the discussion of incompressible materials. Deviating from traditional approaches, we propose that one should flexibly choose an appropriate thermodynamic potential for the discussion based on the particular problem.
In this work, we consider the Gibbs free energy as the thermodynamic potential and derive the constitutive relations based on the Coleman-Noll approach [40,86,124], with both elastic and viscous material responses taken into account. A by-product of this derivation is that the isochoric-volumetric additive split of the free energy, which was regarded as a postulate based on physical intuition, can be justified with solid arguments. The resulting system provides a unified framework for the viscous incompressible Navier-Stokes equations, the compressible Navier-Stokes equations, the compressible hyperelasticity, and the fully incompressible hyperelasticity. In particular, the resulting compressible Navier-Stokes equations derived from the Gibbs free energy are written in the pressure primitive variables. It is well-known that the Navier-Stokes equations are well-behaved in the pressure primitive variables, but not for variable sets where the density is an independent variable [42]. This fact partly corroborates our argument that the Helmholtz free energy degenerates in the incompressible limit. One attractive feature of a unified theory for fluids and solids is that this theory may serve as a bridge for generalizing the VMS method from computational fluid dynamics to computational solid dynamics. Indeed, the VMS method was originally proposed to provide "a paradigm for computational mechanics" [51]. Yet, over the years, progress was made primarily in the area of CFD [6,59,60]. We hope that this work will spark further research on developing VMS formulations for computational mechanics and enable CFD experts to readily bring their expertise to bear on solid mechanics.

Fluid-structure interaction
Fluid-structure interaction problems refer to the coupling of the fluid and structural equations defined on non-overlapping domains with appropriate interface conditions. Over the past several decades, tremendous advancement has been made in computational FSI problems [6,14,66,78]. In terms of how the fluid-structure interface is treated, the numerical methods for FSI problems may be categorized into two major groups: the boundary-fitted approach and the immersed boundary approach. In the boundary-fitted approach, the fluid problem is typically posed in an arbitrary Lagrangian-Eulerian (ALE) coordinate system [55,108], and the fluid domain is represented by a mesh that is deformable with the Lagrangian solid mesh. The boundary-fitted approach enjoys the exact coupling conditions on the interface and accurate stress calculations near the interface, with the expense of moving the mesh. For problems involving very large deformations, it may be necessary to regenerate the mesh to maintain the mesh quality and the solution accuracy, which can be quite expensive for three-dimensional calculations [63].
The immersed boundary method was introduced as an alternative approach for FSI simulations [97]. This approach releases the requirement of the boundary-fitting condition and hence can be quite attractive for problems with very large deformations of the solid boundary, such as the heart valve dynamics [39,66]. One major shortcoming of the immersed boundary approach is the loss of accuracy near the interface. Hence, its applicability is limited when the interfacial physics is important. In this work, we adopt the boundaryfitted approach for our FSI formulation.
In this work, one critical feature that differs from traditional FSI formulations is that the solid problem is written in pressure primitive variables, and its governing equations constitute a set of first-order equations, including a differential mass balance equation, a set of linear momentum balance equations, and a set of kinematic equations. This formulation for the solid problem is formally similar to the ALE formulation of the incompressible Navier-Stokes equations, where the kinematic equations are replaced by the equations governing the mesh motion. Indeed, one can view our FSI formulation as a unified continuum body governed by the mass and linear momentum balance equations, where, in the solid subdomain, the deviatoric part of the Cauchy stress is elastic and the problem is written in the Lagrangian reference frame; in the fluid subdomain, the deviatoric part of the Cauchy stress is viscous and the problem is written in the ALE coordinate system. This unified formulation naturally allows one to construct a uniform VMS formulation for spatial discretization, and the resulting numerical scheme provides a residual-based turbulence model for the fluid subproblem and a VMS scheme for the solid subproblem with the purpose of stabilizing the pressure instability arising from equal-order interpolations.
Writing the FSI problem in a unified formulation also allows us to perform time integration in a uniform way. The generalized-α method has shown to be an accurate and robust temporal scheme for structural dynamics [19], fluid dynamics [62], and FSI problems [6,26,69]. One issue associated with the traditional FSI formulation is that there is a mismatch for the choice of the parameters in the generalized-α scheme [9, p. 120]. The structural dynamics problem is typically written in the pure displacement formulation and hence involves a second-order time derivative. In contrast, the fluid dynamics problem involves only a first-order time derivative. To achieve controllable high-frequency dissipation, the parameters in the generalized-α method are parametrized by the spectral radius of the amplification matrix at an infinitely large time step. This parametrization is different for first-order [19] and second-order systems [62]. In traditional FSI formulations, this leads to a dilemma for the choice of parameters in the temporal scheme. In [6], the parametrization was chosen based on the first-order system to provide optimal dissipation for the fluid problem. That choice sacrifices the dissipation for the solid problem. In [69], the optimal parametrization was chosen for the fluid and solid subdomains separately. Yet, that leads to a kinematic inconsistency on the fluid-solid interface. It has been observed that this kinematic inconsistency may lead to failure in FSI simulations, and an interpolation procedure along the fluid-solid interface was proposed to address this issue [64]. In this work, since the solid dynamics is written as a first-order system, the aforementioned numerical challenge is conveniently resolved. The generalized-α method is applied for the unified continuum problem using the parametrization for the first-order problem [19]. One can get optimal numerical dissipation for both the fluid and solid subproblems.
Methods for solving the discretized FSI problem can be categorized into two families: staggered and monolithic methods. The latter can be further categorized into block-iterative, quasi-direct, and direct strategies [123]. A segregated algorithm for the solid dynamics, based on the original idea introduced in [107], naturally leads to a new solution procedure for FSI problems. First, one solves the pressure and the velocity for the continuum body in the matrix problem. Second, the segregated algorithm for the solid is invoked to obtain the displacement field in the solid subdomain. Third, the solid displacement is applied as a boundary condition for the mesh motion problem in the fluid subdomain. This coupling procedure can be categorized as the quasi-direct solution strategy [9,Chapter 6]. In comparison with traditional FSI formulations, the additional cost in this FSI formulation is mainly due to the introduction of the pressure variable in the solid problem. Considering the number of degrees of freedom in the solid subdomain is typically much smaller than that of the fluid subdomain, the new formulation in fact does not significantly increase the cost of the solution of the system.

Structure and content of the paper
The body of this work is organized as follows. In Section 2, a unified continuum model is derived by choosing the Gibbs free energy as the thermodynamic potential. In Section 3, we perform VMS analysis for the resulting continuum model. In Section 4, we derive a fully discrete scheme for hyper-elastodynamics. In Section 5, the algorithm developed for solid dynamics is coupled with fluid dynamics and constitutes a novel FSI formulation. The coupling procedure and the implementation details for FSI problems are discussed. In Section 6, benchmark problems are studied to examine the numerical formulations. We draw conclusions in Section 7.

Continuum Mechanics
In this section, we begin by presenting the ALE formulation of the balance laws. Following that, we derive constitutive relations based on the Gibbs free energy using the Coleman-Noll type analysis. An interesting result is that the additive split of free energies can be justified in this derivation. We recover several familiar models within our modeling framework. In the last part, for various well-known volumetric energies, their Legendre transformations are derived and discussed.

Continuum mechanics on moving domains
In this section, we discuss the kinematics of a deformable body and present the balance equations defined in an ALE frame of reference. Let Ω X , Ω x , and Ω χ be bounded open sets in R n d , where n d represents the number of space dimensions. They represent the domain occupied by the continuum body in the material (Lagrangian), the current (Eulerian), and the referential frames, respectively. The Lagrangian-to-Eulerian map at time t is a diffeomorphism defined as This map satisfies that ϕ(X, 0) = X, which implies that x is the current location of a material particle whose initial location is X. The displacement and the velocity of the material particle are given by In this work, d (·) /dt designates a total time derivative. The deformation gradient and the Jacobian determinant are defined as The referential-to-Eulerian map at time t is a diffeomorphism defined aŝ and satisfiesφ(χ, 0) = χ. The mesh displacement and the mesh velocity are defined aŝ The mesh deformation gradient and the mesh Jacobian determinant are defined aŝ It also proves convenient to introduce the Lagrangian-to-referential mapping at time t as ϕ(·, t) : Ω X → Ω t χ =φ(Ω X , t) ∀t ≥ 0, X → χ =φ(X, t), ∀X ∈ Ω X .
One can analogously define the displacement, velocity, and deformation gradient for the motion of the referential frame relative to the Lagrangian reference frame [105]. The three mappings introduced above are related through an operator composition This operator composition is illustrated in Figure 1. Let ρ = ρ(x, t) denote the density. The conservation of mass represented on the referential domain Ω χ is Let σ denote the Cauchy stress, and let b denote the body force per unit mass. The linear momentum balance equation on the referential domain Ω χ is The balance of angular momentum is given by (2.6) It can be shown that (2.6) is equivalent to the symmetry of the Cauchy stress, Let E := ι + v · v/2 denote the total energy per unit mass, where ι is the internal energy per unit mass; let q denote the heat flux; let r denote the heat source per unit mass. The balance of total energy is The second law of thermodynamics can be stated as Using the Piola identity and the localization argument, one can derive the advective form of the balance equations and the second law of thermodynamics as follows, (2.14) Moreover, the displacement-velocity kinematic relation can be written in the advective form as Details of the derivation of (2.10)-(2.14) can be found in [108].
Remark 1. ρJ = ρ 0 is a widely used algebraic mass balance equation. However, we posit this is really a special case due to the Lagrangian description. In the general scenario, such as the ALE description, one has to adopt the differential equation (2.10) to describe mass conservation [77].

Constitutive relations
There are multiple concepts of pressure in the literature. The mechanical pressure is defined as the dilatational (i.e., hydrostatic) part of the Cauchy stress. The thermodynamic pressure is an intensive state variable. Its negative value conjugates to the specific volume [15]. Sometimes, mathematicians introduce pressure as a Lagrangian multiplier to enforce the incompressibility constraint. In this work, we adopt the thermodynamic definition for the pressure, and the derivation will reveal how these three concepts are related. In continuum mechanics, the Helmholtz free energy (or the strain energy) is often used as the thermodynamic potential to derive and describe the constitutive relations. Such a derivation often follows the classical Coleman-Noll approach [21,45,76]. Among these relations, the thermodynamic pressure is often expressed as a function of the density and the temperature. However, this relation only remains valid for compressible materials. In the incompressible limit, the pressure-density curve becomes a vertical line of infinite slope since the density is constrained as a constant. We believe this pathological behavior of the pressure-density curve near the incompressible limit is the bane of incompressible solid solvers. A similar argument was made in [40, p. 319] based on a thought experiment, with the conclusion that "it would seem unreasonable to allow a constitutive relation for an incompressible elastic body to involve the pressure." To remedy this degeneracy, one can derive the constitutive relations based on the Gibbs free energy [79]. The Gibbs free energy and the Helmholtz free energy are related by a Legendre transformation [15]. For compressible materials, both the free energies are valid; for incompressible materials, the Helmholtz free energy degenerates, and the Gibbs free energy remains valid. In this section, we will derive the system of equations for a continuum mechanics model based on the Coleman-Noll approach. For simplicity, we take the referential frame of reference to be identical to the material frame of reference within this section. Then we directly have Consequently, the advective form of the balance equations and the second law of thermodynamics can be written as In the above, D represents the dissipation. The specific Gibbs free energy per unit mass is defined as Taking material time derivatives on both sides of (2.21) results in Substituting the internal energy balance equation (2.19), the second law of thermodyanmics (2.20), and the mass balance equation (2.16) into the above relation, one obtains To facilitate our discussion, we introduce the following notations.
1. The right Cauchy-Green tensor C is defined as C := F T F .
2. The deformation gradient F can be multiplicatively decomposed into dilatational and distortional parts [31] as In the above relation, I represents the second-order identity tensor, J 1 3 I represents the volume-changing (dilatational) part of the deformation, andF = J − 1 3 F is the volume-preserving (distortional) part of the deformation. Correspondingly, the right Cauchy-Green tensor can be decomposed as 3. We can obtain the following differentiation relation wherein I is the fourth-order identity tensor and δ IJ is the Kronecker delta. It is straightforward to show that PP = P, which implies that P is a projection.
4. The Cauchy stress σ can be additively split into deviatoric and hydrostatic parts, The second Piola-Kirchhoff stress S is defined by S := JF −1 σF −T , and one can show that dev[σ] = J −1 F (P : S) F T .

5.
The spatial velocity gradient can be additively split into the rate of deformation tensor d and the spin tensor w as Furthermore, we split d into deviatoric and hydrostatic parts as The inner product of σ and ∇ x v can be written as (2.24) 6. Algebraic manipulations can show that the time rate of C can be expressed in terms of d in the following relation The time rate ofC can be derived as Hence, we have (2.25) Using (2.24) and (2.25), the relation (2.22) can be rewritten as (2.26) In the above relation, the time rate of G is on the left-hand side, and the time rates ofC, θ, and p appear on the right-hand side. Invoking Truesdell's principle of equipresence [124], we demand that the Gibbs free energy is a function ofC, p, and θ: Taking material time derivatives on both sides, we obtain the relation

Substituting (2.27) into (2.26) leads to
We make the following choices for the Cauchy stress, the heat flux, the entropy, and the density.
In (2.29),μ is the dynamic shear viscosity; in (2.30),λ is the second viscosity coefficient, and 2 3μ +λ is the bulk viscosity; in (2.31),κ is the termal conductivity. Combining the constitutive relations (2.29) and (2.30), one obtains the Cauchy stress as The first term in (2.34) represents the isochoric elastic stress [45], the second term is the pressure, the third term gives the viscous shear stress, and the last term gives the bulk viscous stress. Observing that σ is symmetric, the angular momentum balance law (2.12) is automatically satisfied. The constitutive relation for the heat flux (2.31) is the Fourier's law. The constitutive relation for the entropy density s (2.32) and the density ρ (2.33) coincides with the classical thermodynamic definitions [15]. Invoking (2.21), (2.32), and (2.33), one can get the constitutive relation for the internal energy per unit mass as Proposition 1. Given the constitutive relations (2.29)-(2.33), the dissipation D defined in (2.20) takes the form Proof. The constitutive relations (2.32) and (2.33) make the last two terms in (2.28) vanish. Therefore, one has The constitutive relation for the heat flux (2.31) leads to The constitutive relation (2.30) leads to Using the constitutive relation (2.29), the first term in (2.36) can be simplified as In summary, one has which completes the proof.
The relation (2.35) suggests that the dissipation D is guaranteed to be non-negative if the dynamic shear viscosity, the bulk viscosity, and the thermal conductivity are non-negative. Proof. One only needs to verify that (2.29) is material frame indifferent. Considering a proper orthogonal tensor Q, one has det(QF ) = det(Q)det(F ) = det(F ) = J. Therefore, QF = QF . It is known that under the rigid-body motion described by Q, dev[d] transforms to Qdev[d]Q T . Consequently, the right-hand side of (2.29) transforms to Proposition 3. The Gibbs free energy takes the following additive decoupled form, Proof. The density ρ is independent ofC sinceC is volume-preserving (det(C) = 1). Consequently, the constitutive relation (2.33) leads to Integrating the above partial derivative gives (2.37), where G vol (p, θ) = ρ −1 dp.
It is worth pointing out that the free energy adopted in the above discussion is the specific free energy, which means that it is the energy per unit mass. It is sometimes useful to introduce the free energy per unit volume in the material configuration [40,45], which is denoted as G R . The two energies are linked by the relation G R = ρ 0 G. Correspondingly, we denote G R iso := ρ 0 G iso and G R vol := ρ 0 G vol . Let H R denote the Helmholtz free energy per unit volume in the material configuration. It can be obtained through a Legendre transformation of −G R with respect to p, namely, Invoking the additive split of the Gibbs free energy (2.37), one has This relation can be summarized as the following proposition.
Proposition 4. The Helmholtz free energy H R admits the additive decoupled form Remark 2. The additive split of the Helmholtz free energy was introduced as a postulate based on the multiplicative decomposition of the deformation gradient [31]. Above, we have given a justification of the additive split (2.38) for hyperelastic materials. Extension of the arguments to inelastic materials is under assessment.
With the energy split (2.37), the constitutive relations (2.33) and (2.34) can be rewritten as Here, we introduce the isobaric thermal expansion coefficient α p and the isothermal compressibility coefficient β θ [15]. They are defined as Making use of the constitutive relation (2.39), they can be expressed explicitly as

Examples of closed systems of equations
Before proceeding further, we notice that the time derivative of ρ can be expanded as Hence, the governing equations in the ALE frame of reference become Ideal gas model The Gibbs free energy for the ideal gas is wherein R is the specific gas constant, C v is the specific heat at constant volume, p ref and θ ref are the reference values of the pressure and temperature. With this choice, the constitutive relations are The balance equations (2.43)-(2.45) together with the above constitutive equations constitute the pressure primitive variable formulation for the compressible Navier-Stokes equations [42]. It was shown that the set of pressure primitive variables is among the sets of variables that are well-behaved for both compressible and incompressible flows [41,42,53,75]; on the contrary, any set of variables involving density (e.g. conservation variables or density primitive variables) becomes ill-defined in the incompressible limit [41]. That observation, in part, justifies our motivation of deriving a continuum mechanics model based on the Gibbs free energy.
Incompressible viscous flow We introduce the Gibbs free energy for incompressible viscous fluid flow as With this choice, the constitutive relations are Compressible hyperelastic model The Gibbs free energy for compressible hyperelastic materials takes the following general form.
If we ignore the viscous effect, the constitutive relations are Incompressible hyperelastic model Noticing that incompressibility impliesC = C, the Gibbs free energy for incompressible hyperelastic materials takes the following form.
Ignoring the viscous effect, the constitutive relations can be written as Remark 3. It is worth pointing out that the displacement-velocity relation (2.15) is not the unique choice for describing kinematics. The deformation gradient transport relation [108, p. 24], can be utilized to calculate the strain and stress as well [27,109]. In fact, this strategy is expected to give second-order spatial accuracy in the calculation of the strain and the stress with linear elements. A trade-off is that n 2 d additional differential equations for F need to be solved. Based on the polyconvexity hypothesis, it seems natural to introduce kinematic relations for F , JF −T , and J [10,11,36]. In that approach, the kinematic equations involve 2n 2 d + 1 degrees of freedom. In this work, we choose to solve the simple displacement-velocity equation (2.15). Its simple structure leads to an additional benefit in the design of the algorithm.

Legendre transformation of volumetric energies
Here, we want to point out that, in the literature, there exist derivations of constitutive relations based on a different Gibbs free energy [99,116]. In those works, the Gibbs free energy is taken as a function of the stress, which implicitly requires the Helmholtz free energy to be convex with respect to the strain. Convexity is a strong requirement and rules out many important nonlinear materials. This is a key difference with the present work, which does not impose the requirement. Most hyperelastic materials are postulated to be polyconvex [4,86]. A direct consequence of polyconvexity is that, with the isochoric-volumetric split shown in Proposition 4, the volumetric energy is convex. Thus, it is indeed legitimate to perform a Legendre transformation on the volumetric part of the energy. The derivations in [99,116] correspond to the Hellinger-Reissner variational principle; while our derivation corresponds to the Herrmann variational principle [44,67]. In the following, we present a few examples of the Lengendre transformations for several classical volumetric energies, and we use κ to denote the bulk modulus.
Quadratic volumetric energy The first example is the quadratic energy The corresponding specific energy is The conjugate function of the above H vol (v) is According to (2.39) and (2.41), we have

ST91 volumetric energy
The quadratic volumetric energy (2.46) has been under criticism because it approaches a finite value as the volume goes to zero. This may lead to numerical instabilities [113]. To circumvent this issue, the following volumetric energy was proposed in [114] and is now widely used.
The specific energy is The Taylor expansion of the above G vol is Clearly, the G vol (p) for the ST91 volumetric energy (2.47) can be viewed as a high-order modification of the one associated with the quadratic energy. From the formula of G vol (p), one obtains M94 volumetric energy The volumetric free energy proposed in [89] and its corresponding specific free energy are

Its Taylor expansion is
Based on the above formula, one arrives at Here, the constitutive relation for the density is the linear barotropic relation, which is usually used to describe small density variations near a reference value [38].

L94 volumetric energy
The volumetric free energy proposed in [73] and the corresponding specific free energy are

Its Taylor expansion is
Interestingly, this volumetric free energy gives a constant isothermal compressibility coefficient.

Remark 4.
For compressible materials, one can derive pressure by taking derivative of the volumetric energy with respect to J as It can be verified that (2.51) is compatible with the constitutive equation (2.39). At this stage, one can clearly see that the equation (2.51) is an equation of state, just like p = ρRθ for the ideal gas model. In CFD, it is very rare to see one put p = ρRθ into a weak form and solve it by finite element or finite volume methods. By analogy, we feel it is worth raising a question on the validity of using (2.51) in the mixed formulation for finite elasticity, although such an approach is regarded as well-established [22,45,68,82].
Remark 5. In [107], it has been observed that, even for small-strain elasticity, it is not advisable to use 0 = p/κ + ∇ x · u in the VMS formulation for transient analysis. This fact confirms our doubts about using (2.51) in the VMS formulation for finite elastodynamics.
Remark 6. In the incompressible limit, the bulk modulus approaches infinity, and we have G vol (p) → p/ρ 0 , ρ → ρ 0 , β θ → 0. Therefore, the above constitutive relations are well-defined in both compressible and incompressible regimes. In contrast, the constitutive relation (2.51) based on the Helmholtz free energy H R vol will blow up in the incompressible limit.
Remark 7. It is interesting to notice that, although the volumetric energies H R vol are different, their conjugate counterparts G vol are very similar. The first two terms of their Taylor expansions are identical. Notice that, besides the convexity condition, physical intuition also suggests that the volumetric energy H R vol achieve its minimum value at J = 1 and blows up to infinity as J → 0 and J → ∞. We feel that these conditions may imply some mathematical properties for G vol . These properties may help design constitutive relations directly based on G vol . We believe future work in this area will be useful both theoretically and practically.

Variational Multiscale Analysis
The variational multiscale method was introduced as a general framework for subgrid-scale modeling in computational mechanics [49,51]. As a generalization of the stabilized methods, it improves the stability bound for singularly perturbed problems and overcomes the infsup condition for saddle-point problems. Interested readers are referred to [1,20,59] for comprehensive reviews. In this section, we invoke the residual-based VMS method [5,95] to construct a formulation for the continuum problem derived in Section 2. Here, and in what follows, we restrict our discussion to the isothermal condition 1 . The system of equations (2.42)-(2.45) 2 are simplified as The constitutive relations for the density ρ and the isothermal compressibility β θ can be simplified as univariate functions of the pressure: In this section, we consider the strong-form problem endowed with proper initial conditions and periodic boundary conditions. To simplify the notation, we choose V to denote both the trial and the test function spaces, which are assumed to be identical in this section. Let (·, ·) Ωx denote the L 2 inner product over the domain Ω x . The variational formulation for the equations (3.1)-(3.3) can be stated as follows.
Now we introduce a projection operator P : With the aid of the projection operator, we have a well-defined direct-sum decomposition of the function space V as where I is the identity operator. Here, V represents the unresolved fine scales. With this space decomposition, we can decompose the trial solution y and the test function w as y =ȳ + y ,ȳ = Py, y = (I − P) y, w =w + w ,w = Pw, w = (I − P) w.
With the decomposition of w and by virtue of the linear dependency of the variational formulation B(w, y) in w, we can decompose the original variational formulation into a coupled system as (3.7) The above two equations are usually referred to as the coarse-scale and the fine-scale equations [5]. We assume that B is Fréchet differentiable with respect to y ∈ V up to the n-th derivative. Using the Taylor's Formula in the Banach space [2], the left-hand side of (3.7) can be expanded as In the above, D k y B represents the k-th derivative of B in terms of the second argument y. It is a k-linear functional on V × · · · × V . Moving B (w ,ȳ) to the right-hand side, we have One can represent the right-hand side of (3. Inserting (3.9) into (3.6), a closed, finite-dimensional system can be obtained forȳ, Given the analytic form of the mapping F , one may obtainȳ from (3.10) and y from (3.9). The resulting y =ȳ + y is the exact solution of the original problem (3.5). However, obtaining an analytic form for F is as hard as solving the original problem analytically, if not harder. A practical approach is to systematically design an approximated mapping F . Replacing F in (3.9)-(3.10) by the approximated mapping, one may obtain a suite of computable formulations for the fine-and coarse-scale components. Similar to the residual-based VMS modeling approach for turbulence [5], we introduce a perturbation series to represent y and derive a detailed pathway to constructF . The difference between our approach and the one adopted in [5] is that, in addition to the approximation of the fine-scale Green's operator and the truncation of the perturbation series, we introduce one additional approximation procedure, i.e. the truncation of the Taylor expansion formula in (3.8). This additional step is due to the general nonlinear term that may appear in finite elasticity. Since the derivation of the fine-scale approximation goes deeper into functional analysis, we give the detailed derivation in A. In our model, the fine-scale component is approximated as In the above, the choice τ = diag(τ K , τ C , τ M ) implies that the fine-scales are postulated to be decoupled. The precise formulas for the stabilization parameters τ K , τ C , and τ M depend on the specific problem considered. For simple linear problems, τ can be computed as a local mean-value of the fine-scale Green's function [49,58]; sometimes, error estimates provide a guidance for the design of τ [32,54]; for complex problems, scaling arguments are usually made for the design of τ [121]. The detailed formula of τ for solid and fluid dynamics will be given in the subsequent sections. With (3.11), we can complete our VMS formulation as follows, This formulation provides a basis for our finite element formulations in the subsequent sections.

Formulation for solid dynamics
In this section, we restrict our discussion to hyper-elastodynamics. Within the general VMS framework developed in Section 3, the problem is spatially discretized using the VMS formulation. The generalized-α method is utilized for temporal discretization. As is shown in [107], a block decomposition of the tangent matrix reveals that the problem can be solved in a segregated manner without loosing consistency in the nonlinear solver. Lastly, we discuss the choice of the stabilization parameters.

Initial-boundary value problem
We consider the hyper-elastodynamic problem written in the Lagrangian reference frame, The time interval of interest is denoted as (0, T ), with T > 0. The boundary Γ x = ∂Ω x can be partitioned into two non-overlapping subdivisions: Γ x = Γ g x ∪ Γ h x , wherein, Γ g x represents the Dirichlet part of the boundary, and Γ h x represents the Neumann part of the boundary. Boundary conditions for this problem are imposed as Given the initial data u 0 , p 0 , and v 0 , the initial conditions for the strong-form problem (4.1)-(4.3) can be stated as Here we only consider hyperelastic materials, and the deviatoric part of the Cauchy stress takes the following specific form The constitutive relations ρ = ρ(p) and β θ = β θ (p) are given by the relations (3.4). To simplify our subsequent discussion, we introduce the first Piola-Kirchhoff stress P := JσF −T , andP := Jσ dev F −T .
Remark 8. Although upon first glance, the mass equation (4.2) may look similar to the pressure rate equation with artificial compressibility proposed by A.J. Chorin in [18], the concept is different. The isothermal compressibility coefficient β θ here is a real physical quantity, not a numerical artifact. It is zero for fully incompressible materials.

Variational multiscale formulation
Based on the VMS formulation derived in Section 3, the formulation for the strong-form problem (4.1)-(4.3) can be constructed conveniently. Consider a discretization of the current domain into finite elements. The union of element interiors is denoted by Ω x . Let us denote the finite dimensional trial solution spaces for the solid displacement, pressure, and velocity in the current domain as S u h , S p h , and S v h , respectively. We assume that functions in the trial spaces satisfy the Dirichlet boundary conditions (4.4)-(4.5) on Γ g x . Let V u h , V p h , and V v h denote the corresponding test function spaces. The VMS formulation can be stated as follows. Find y h (t) : Here u h0 , p h0 , and v h0 are the L 2 projections of the initial data onto the finite dimensional spaces S u h , S p h , and S v h , respectively. The subscript h denotes a mesh parameter. In the above and henceforth, the formulations for the kinematic equations, the mass equation, and the linear momentum equations are indicated by the subscripts k, p and m, respectively.
Remark 9. The VMS formulation (4.8)-(4.12) is derived directly from (3.12) by adopting the finite dimensional trial solution spaces as V and taking the following assumptions, (3) u = v = 0 on the boundary ; (4) (w p h , β θ (p h + p )dp/dt) Ωx = (w p h , β θ (p h )dp h /dt) Ωx ; The assumptions are adopted to simplify the numerical model. Similar assumptions have been made in the residual-based VMS modeling for turbulence [5].
The VMS formulation (4.8)-(4.12) can be pulled back to the material frame of reference. To obtain that, we define the test functions defined in the material frame of ref- H(X, t) := h(ϕ(X, t), t), G(X, t) := g(ϕ(X, t), t).
The corresponding trial solution spaces on Ω X are denoted as S U h , S P h , and S V h ; the corresponding test function spaces are denoted as V U h , V P h , and V V h . The VMS formulation in the material frame of reference can be stated as follows.
Here U h0 , P h0 , and V h0 are the L 2 projections of the initial data onto the finite dimensional spaces S U h , S P h , and S V h , respectively.
Remark 10. If piecewise linear elements are used for the spatial discretization, the term, in (4.9) and the term in (4.14) vanish.

Temporal discretization
A desirable time integration algorithm for structural dynamics should satisfy at least three requirements: unconditional stability, second-order accuracy, and dissipation on high-frequency modes. The generalized-α method is one type of time integration schemes that satisfy these requirements [19,62,65]. For this reason, we adopt this method for time integration in this work. The time interval [0, T ) is divided into a set of n ts subintervals of size ∆t n := t n+1 − t n delimited by a discrete time vector {t n } nts n=0 . The solution vector and its first-order time derivative evaluated at the time step t n are denoted as Y n andẎ n ; the basis function for the discrete function spaces is denoted asN A . With those, the residual vectors can be represented as The fully discrete problem can be stated as follows. At time step t n , givenẎ n , Y n , the time step ∆t n , and the parameters α m , α f , and γ, findẎ n+1 and Y n+1 such that The parameters α m , α f , and γ define the time integration scheme. It has been shown in [62] that, for linear problems, second-order accuracy in time can be achieved, provided and unconditional stability can be attained if The generalized-α method allows one to control the damping effect on the high frequency modes while maintain the accuracy of the temporal scheme. This desired property is achieved for first-order linear equations by choosing the parameters as wherein ∞ denotes the spectral radius of the amplification matrix at the highest mode [19,62]. In long-time nonlinear calculations, it is advisable to choose ∞ strictly less than one to avoid detrimental effects from the high frequency modes [6,62,74]. Unless otherwise specified, we choose ∞ = 0.5 in this work.
Remark 11. Using the generalized-α method for the first-order structural dynamic problem has recently been shown to enjoy improved numerical properties in comparison with the generalized-α method applied for the second-order structural dynamic systems [65].

A segregated algorithm
It can be shown that the fully discrete kinematic equation (4.16) is equivalent tō Invoking the relations (4.19)-(4.21), the left-hand side of above equation can be written explicitly as Equations (4.17), (4.18), and (4.23) constitute a system of nonlinear algebraic equations and can be solved by the Newton-Raphson method. The solution vector Y n+1 at the time step t n+1 with n = 0, · · · , n ts − 1 is solved iteratively by means of a predictor multi-corrector scheme. We define as the solution vector at the Newton-Raphson iteration step i = 0, · · · , i max . We denote that the residual vectors evaluated at the iteration stage i as The tangent matrix of the above nonlinear system is denoted as K k,(i),Ṗ := 0, (4.29) Notice that (4.28) and (4.30) are diagonal matrices. This special structure of the two submatrices can be exploited to perform a block decomposition of the tangent matrix K: With the above decomposition, the original linear system of equations for the Newton-Raphson iteration, can be solved in a two-stage segregated algorithm. In stage one, the and (4.38), we have the following observations, With the above relations, the linear system (4.37) can be reduced to a smaller problem, . (4.39) The linear system (4.38) can be reduced to the relation Therefore, the solution procedure can be consistently rewritten into two smaller problems. One first solves the equation (4.39) to obtain ∆Ṗ n+1,(i) , ∆V n+1,(i) T . Then the displacement increment ∆U n+1,(i) can be obtained through the relation (4.40). We can summarize the above discussion as the following segregated predictor multi-corrector algorithm.
Predictor stage: Set: Multi-corrector stage: Repeat the following steps i = 1, . . . , i max : 1. Evaluate the solution vectors at the intermediate stages: 2. Assemble the residual vectors R (i) based on (4.24)-(4.27) using the solution evaluated at the intermediate stages.
3. Let R (i) l 2 denote the l 2 -norm of the residual vector. If either one of the following stopping criteria is satisfied for prescribed tolerances tol R , tol A , set the solution vector at time step t n+1 asẎ n+1 =Ẏ n+1,(i−1) and Y n+1 = Y n+1,(i−1) , and exit the multi-corrector stage; otherwise, continue to step 4.
and return to step 1.
Since the kinematic equation is linear, it is reasonable to expect that this equation can be solved in one Newton-Raphson iteration. Indeed, we have the following proposition.
The proof of this proposition is given in B. Due to this fact, the right-hand side of (4.39) becomes − R p,(i) , R m,(i) T , and the relation (4.40) reduces to Remark 12. In general,R k (1) = 0. In the algorithm we presented, the same-Y predictor [62] is adopted, and it is straightforward to show that, with this predictor, If the time-stepping parameters are α m = α f = γ = 1 and the predictors are chosen as Setting α m = α f = γ = 1 corresponds to the backward Euler method. In [102], a non-traditional predictor is chosen to enforceR k (1) = 0. Remark 13. This segregated algorithm based on block decomposition was first proposed by G. Scovazzi and co-authors in [107]. This solution procedure (4.39)-(4.40) significantly reduces the size of the linear system for the implicit solver. In comparison with the second-order pure-displacement formulation, writing the system into a first-order system does not lead to a significant increase in computational cost. As will be revealed in Section 5.4, the first-order formulation is very appealing for FSI problems.

Stabilization parameters
The design of the stabilization parameter is of critical importance for the behavior of the stabilized finite element formulation. In this work, a practical setting of the stabilization parameters is wherein ∆x is the diameter of the circumscribing sphere of the tetrahedral element, c m and c c are two non-dimensional parameters, c is the maximum wave speed in the solid body. In this work, the formula for c is chosen based on a small-strain isotropic linear elastic material [54,107]. For compressible materials, c is given by the bulk wave speed and for incompressible materials, c is given by the shear wave speed Remark 14. The stabilization parameter is the crux of the design of the VMS formulation. In [54], the authors proposed two choices of the stabilization parameter for the smallstrain elastodynamic problem: 0.5∆x/c and 0.5∆t. In [37,70], the stabilization parameter is designed based on ∆t. In [107], the authors proposed a set of stabilization parameters as c s15 ∆t/2c CF L or c s15 ∆x/2c, wherein c s15 is a non-dimensional parameter, and c CF L is the global CFL number. In [102], the stabilization parameter for compressible materials is designed as Remark 15. In the current design, the elastic wave speed c is based on the small-strain elastic theory. By exploiting the eigenvalue structure of the large deformation problem, one can obtain the elastic wave speed for nonlinear materials [10,36]. This will surely lead to a better design of the stabilization parameters.

Formulation for fluid dynamics and fluid-structure interaction
In this section, we present a new framework for FSI problems based on the unified continuum model derived in Section 2. The formulation for solid dynamics is directly adopted from the one developed in Section 4; the formulation for fluid dynamics is constructed based on the general VMS formulation given in Section 3. The finite element formulation, time integration, and a new solution procedure are discussed in this section. In all discussions related to FSI problems, we use a superscript f to indicate quantities related to fluids, a superscript m to indicate quantities related to the mesh motion in the fluid subdomain, and a superscript s to indicate quantities related to solids.

Initial-boundary value problem
In this section, we present the initial-boundary value problem for the FSI problem. The time interval of interest is denoted as (0, T ), with T > 0. The domain occupied by the continuum Figure 2: Illustration of the fluid-structural interaction problem setting.
body in the referential frame Ω χ admits the decomposition Ω f χ and Ω s χ are the two subdomains occupied by the fluid and solid respectively. The fluidsolid interface is a R n d −1 -dimensional manifold and is denoted as Γ I χ . The current domain Ω x can be decomposed correspondingly as The boundary Γ x = ∂Ω x can be partitioned into four non-overlapping subdivisions: In the above decomposition, Γ gs x represents the Dirichlet part of the solid boundary; Γ g f x represents the Dirichlet part of the fluid boundary; Γ hs x represents the Neumann part of the solid boundary; Γ h f x represents the Neumann part of the fluid boundary. The unit outward normal vector to Γ x is denoted as n. The fluid-solid interface in the current configuration is denoted as Γ I x ; n s and n f represent the unit outward normal vector to Ω s x and Ω f x on the interface Γ I x , respectively. The configurations and the boundary decomposition are illustrated in Figure 2. The referential configuration can be defined by the Lagrangian-to-referential map χ =φ (X, t) . We define the mapφ in the solid subdomain Ω s X as the identity map χ =φ (X, t) = X, for X ∈ Ω s X and t ≥ 0. (5.1) Due to the composition relation (2.3), one has ϕ =φ • id Ω s X . This implies that the material description is adopted in the solid subdomain. In the fluid subdomain, the referential configuration can be determined by the referential-to-Eulerian mapφ. Typically, one determineŝ u m in Ω f χ first, and the referential-to-Eulerian map in the fluid subdomain is given by the definition (2.1). The motion of the referential domain needs to maintain the regularity for functions defined on Ω χ [94]. Oftentimes, the construction of this mapping is tailored to specific problems. In this work, we consider two options: the harmonic extension algorithm and the pseudo-linear-elasticity algorithm.
The harmonic extension algorithm One simple and effective way of constructing the mapφ is by solving forû m as a harmonic extension of the trace of U s on Γ I χ [131], viz.
Onceû m is obtained, the mesh velocityv m in the fluid subdomain can be calculated aŝ In all, we define the mesh displacementû and the mesh velocityv aŝ The pseudo-linear-elasticity algorithm One may also model the motion of the fluid subdomain by solving a succession of pesudo-linear-elastostatic equations [63,131]. We consider a time instantt < t such thatt is close to t. In numerical computations,t is often conveniently chosen as the previous time step. The mesh displacement at timet < t is given byû χ,t :=x − χ =φ χ,t − χ. Consequently, one has the following identity, Freezing the time instantst and t, we defineũ m (x) asũ m •φ(χ,t) :=û (χ, t) −û χ,t . Given the fictitious Lamé parameters λ m and µ m ,û m is determined by solving the following linear elastostatic equations: In the above, u s t and u s t are the solid displacements in the current configurations at time t andt, respectively. In our implementation, the fictitious Lamé parameters are divided by the Jacobian determinant of the element mapping [63]. This effectively increases the robustness of the mesh moving algorithm.
The fluid-solid body can be viewed as a single continuum body governed by the following balance equations.
in Ω x , (5.7) In the fluid subdomain, we consider viscous incompressible fluid flow, with the constitutive relations given in Section 2.3. Due to the fact that β f θ (p) = 0, the strong-form problem (5.7)-(5.8) in the fluid subdomain can be further simplified as In the solid subdomain, we havev = v due to the Lagrangian description (5.1) adopted in the solid subproblem. Consequently, the strong-form problem (5.7)-(5.8) for the solid can be explicitly written as Given the Dirichlet data g s and g f and the boundary tractions h s and h f , the boundary conditions can be stated as In the fluid subdomain, the initial conditions are given by a divergence-free velocity field In the solid subdomain, the initial conditions are specified as

Variational multiscale formulation
In this section, we apply the VMS formulation developed for the general continuum problem in Section 3 to the FSI problem. The VMS formulation for the solid problem has been derived in Section 4.2. For completeness of the FSI formulation, it is stated as follows. Let S s u h , S s p h , and S s v h denote the finite dimensional trial solution spaces for the solid displacement, pressure, and velocity in the current domain, and let V s u h , V s p h , and V s v h denote the corresponding test function spaces. We assume that the Dirichlet boundary conditions are built into the definitions of the trial solution spaces. The variational multiscale formulation is stated as follows. Find y s In the above formulation, Ω s x represents the union of solid element interiors. The above formulation can be conveniently pulled back to the material frame of reference, as was done in Section 4.2.
Next, we give the formulation for the motion of the fluid subdomain. We only present the variational formulation for the harmonic extension algorithm. The formulation based on the pseudo-linear-elastic algorithm can be similarly constructed. Let S m u h denote the trial solution space of the mesh displacementû m h , and let V m u h denote the corresponding test function space. The kinematic boundary conditionû m Lastly, we present the VMS formulation for the fluid problem. Let S f p h and S f v h denote the trial solution space of the fluid pressure and velocity, and let V f p h and V f v h denote the test function spaces. We assume that functions in the space S f v h satisfy the Dirichlet boundary condition on Γ g f x . The VMS formulation is stated as follows. Find In the above, ξ = {ξ i } n d i=1 are the coordinates of an element in the parent domain, and C I is a positive constant derived from an element-wise inverse estimate [32]. Unless otherwise specified, the value of C T is taken to be 4. The value of C I is independent of the mesh size but relies on the polynomial order of the interpolation basis functions. For linear interpolations, C I is suggested to be 36 [29, p. 65].

Formulation for the coupled problem
Based on the individual subproblems given in the previous section, the semi-discrete FSI problem can be stated as follows. Find y s wherein the last term is an additional stabilization term defined as This stabilization term (5.38) provides additional stabilization to control oscillations [56,57,110]. Due to the choice ofτ , the term (5.38) is proportional to the residual, and therefore the consistency of the semi-discrete formulation (5.37) is maintained. In this work, the whole fluid-solid domain is treated as a single continuum body and is discretized by a single set of mesh with continuous basis functions across the fluid-solid interface. This mesh choice and equal-order interpolations directly imply the satisfaction of the following kinematic constraint relations on the interface, It also implies the following relations for the test functions, Invoking standard variational arguments [6], the relations (5.40) lead to the weak continuity of the traction across the fluid-solid interface, viz.
The above relation, together with (5.39), ensures the correct coupling between the fluid and the solid sub-systems.

Temporal discretization
In this section, we apply the generalized-α method to the semi-discrete formulation (5.37) to construct the fully discrete scheme. The time interval [0, T ) is divided into a set of n ts subintervals of size ∆t n := t n+1 − t n delimited by a discrete time vector {t n } nts n=0 . Let y s n := {u s n , p s n , v s n } andẏ s n := {u s n ,ṗ s n ,v s n } denote the solution vector and its first-order time derivative of the solid displacement, pressure, and velocity evaluated at time t n ; let y f n := û m n , p f n , v f n andẏ f n := v m n ,ṗ f n ,v f n denote the solution vector and its first-order time derivative of the mesh displacement in the fluid subdomain, fluid pressure, and fluid velocity evaluated at time t n . Let N A denote the basis function on the current configuration and letN A denote the corresponding basis function defined on the reference configuration. They are related by the mappingφ as N A =N A •φ −1 . We denote the residual vectors as The fully discrete scheme can be stated as follows. At time step t n , givenẏ n and y n , finḋ y n+1 and y n+1 such that The choice of the parameters α m , α f and γ is given by (4.22).

Remark 16.
It is shown in [19] that for second-order ordinary differential equations, the parametrization for α m in terms of ∞ takes the form Only when ∞ = 0, the parametrizations of α m for the first-order and the second-order system coincide. In FSI simulations, oftentimes one solves the solid problem using the pure displacement formulation and solves the fluid problem using the Navier-Stokes equations [9]. In staggered FSI solvers, one can choose separate, optimal parameters for the generalizedα methods [26]; in monolithic FSI solves, the mismatch of the α m parameter engenders a dilemma [9, pp. 119-120]. In [6], the parametrization (4.22) is adopted for the whole FSI system since the Navier-Stokes equations are more challenging than the structural problems. However, a drawback with that approach is the non-optimal high-frequency dissipation in the structural equations. In [69], the optimal choices of α m are applied to the fluid and the solid equations separately. That leads to incompatibility of the acceleration at the fluid-solid interface and has been observed to engender failure in FSI simulations [64]. In our proposed approach, this issue is properly addressed since the solid mechanics problem is written as a first-order system by introducing the displacement-velocity kinematic equations. This allows one to achieve optimal numerical dissipation in the monolithic FSI solver without violating the kinematic compatibility.

A predictor multi-corrector algorithm based on the segregated algorithm
The coupled nonlinear system of equations (5.41)-(5.46) can be solved monolithically using the Newton-Raphson method [98]. For FSI problems, a monolithic solver can be implemented using the block-iterative [120], quasi-direct [122], or the direct coupling techniques [6]. In the quasi-direct coupling approach, the fluid and solid equations are treated as a block and the mesh is treated as a block. One solves a block of equations using the most recent unknowns from the other block. The quasi-direct coupling approach enjoys a good balance of robustness and computational cost and has been recommended for general FSI problems [8]. In this work, we invoke the quasi-direct coupling methodology and the segregated algorithm for the solid equations developed in Section 4.4. This leads to a novel algorithm for the nonlinear solver in FSI problems. This algorithm is stated as follows.

Evaluate the solutions at the intermediate stage aṡ
There are 2n d + 1 degrees of freedom associated with each computational node in both subdomains; in step 3, n d + 1 degrees of freedom need to be solved by the implicit solver. In the implicit solver, the mass and linear momentum equations are discretized by VMS in space and the generalized-α method in time. In the assembly routine, one only needs to call the proper constitutive relations for the material. In doing so, one can implement the solution vector and the tangent matrix problem over the whole continuum body in a unified approach, which may simplify the data structure management in the numerical implementation. The linear system of equations (5.50)-(5.53) can be solved by means of the Generalized Minimal Residual (GMRES) algorithm [103] with an additive Schwarz preconditioner; the linear system of equations for the mesh motion (5.54) is symmetric positive definite and hence can be solved by means of the conjugate gradient method. The PETSc package [3] is adopted to provide an interface for a wide range of solver options.

Benchmark computations
In this section, we first use manufactured solutions to assess the convergence behavior of the algorithm developed in Section 4 for both compressible and incompressible materials. Following that, a classical benchmark problem for incompressible solids is studied, and the results for the VMS formulation with P 1 /P 1 and Q 1 /Q 1 elements are compared with the results of theF -projection method. In the last, two FSI benchmark problems are studied to demonstrate the effectiveness of the new FSI formulation proposed in Section 5.

Manufactured solution for compressible hyperelasticity
In the first example, the displacement field is given in a closed-form formulation for a unit cube (1 cm × 1 cm × 1 cm). In the prescribed displacement field, T 0 = 1.0 × 10 −3 s, β 1 = 10 −3 π rad/cm. This manufactured solution describes rotation of the cube with the bottom surface Z = 0 clamped. A static version of this manufactured solution has been used in [71] for code verification purposes. The material model is chosen as the Neo-Hookean model with the volumetric energy given by (2.47). Its Gibbs free energy is The  Figure 4. Note immediately that the L 2 -norms of the errors in displacement and velocity converge quadratically; the asymptotic convergence rate for the pressure is 1.8; the deformation gradient and the deviatoric part of the Cauchy stress converge linearly. In [88], a similar sub-optimal convergence rate for the pressure field has been observed, using a different stabilized formulation.

Manufactured solution for fully incompressible hyperelasticity
In the second example, we consider a fully incompressible Neo-Hookean material model. Its Gibbs free energy is The geometrical domain is again a unit cube (1 cm × 1 cm × 1 cm). The shear modulus is chosen as 100 KPa; the density ρ 0 is chosen to be 1.0 × 10 3 kg/m 3 . The analytic forms of the displacement and the pressure fields are given as follows.    Figure 5. For the fully incompressible material, the convergence rates for the displacement and velocity are of order two; the convergence rates for the deformation gradient and the deviatoric part of the Cauchy stress are firstorder; the convergence rate for the pressure becomes first-order for the fully incompressible material. Although there is no convergence proof for this nonlinear problem, the convergence rate of the pressure in the L 2 -norm for linear elasticity is proved to be first-order [34,52]. Therefore, the first-order convergence rate for the pressure field is indeed expected.

Nearly incompressible block under compression
This example was initially proposed in [101] as a benchmark problem for nearly incompressible solids. The original problem was proposed as a quasi-static problem. In this work, we pose the problem in the dynamic setting and adopt the Neo-Hookean model described by Notice that the material model we used here is slightly different from the one used in the original paper [101] because we demand that the Gibbs free energy takes the decoupled form (2.37) and the isochoric part of the energy is a function ofC. Following [101], the material parameters are chosen as µ s = 80.194 MPa, κ = 400889.806 MPa, and ρ 0 = 1.0 × 10 3 kg/m 3 . The corresponding Poisson's ratio is 0.4999. The problem setting is illustrated in Figure 6. Symmetry boundary conditions are applied on the X = Y = Z = 0 planes. A 'dead' surface load H is applied on a quarter portion of the top of the block, and the load assumes the negative z-direction in the reference configuration. The magnitude of H is measured by |H| = ΛH 0 , and the reference value H 0 = 4 MPa. We calculate the compression level of the upper center point with Λ = 20, 40, 60, and 80. The load force is gradually increased as a linear function of time. The problem is integrated in time with 200 time steps. We performed simulations using theF -projection method with Q 1 /Q 0 element and the VMS formulation with Q 1 /Q 1 element and P 1 /P 1 element. The stabilization parameters are chosen as c m = 0.1 and c c = 0.1. In Figure 7, the compression levels of the upper center point with mesh refinement for different loading ratios and different numerical methods are illustrated. The VMS formulation with P 1 /P 1 element tends to give very stiff response with coarse meshes. In contrast, theF -projection method with Q 1 /Q 0 element and the VMS formulation with Q 1 /Q 1 element tend to give very soft response with coarse meshes. With mesh refinement, convergent results are obtained for all three different methods.

Flow over an elastic cantilever
The two-dimensional flow-induced oscillation of an elastic cantilever attached to a fixed square block was initially designed in [127] as a benchmark problem for FSI algorithms. It has been used extensively in the literature to assess the quality of FSI algorithms [6,26,128]. In this work, the original two-dimensional problem is extended to a three-dimensional problem by extruding the original problem in the third direction [128]. On the inflow surface, a uniform flow in the x-direction is imposed with magnitude 51.3 cm/s; on the outflow surface, a zero traction boundary condition is applied; on the four lateral boundary surfaces, slip boundary conditions (zero normal velocity and zero tangential traction) are applied. The geometry and boundary conditions are illustrated in Figure 8. The fluid density and dynamic shear viscosity are set to be ρ f 0 = 1.18 × 10 −3 g/cm 3 andμ = 1.82 × 10 −4 poise, respectively. The solid is modeled as the Neo-Hookean material with the volumetric free energy given by (2.48), The solid referential density ρ s 0 is 0.1 g/cm 3 , the shear modulus µ s and the Poisson's ratio are 9.2593 × 10 5 dyn/cm 2 and 0.35, respectively.      The fluid subdomain is discretized with 74696 tetrahedral elements, and the solid subdomain is discretized with 5369 tetrahedral elements (Figure 9). The mesh motion is governed by the pseudo-linear-elasticity algorithm, wherein the mesh Young's modulus is set to be unity and the mesh Poisson's ratio is set to be 0.3. The mesh is fixed at the inflow and outflow boundary surfaces as well as the square block surface. On the lateral surfaces, the mesh is only allowed to move in the tangential direction. In the algorithm, the mesh Young's modulus is multiplied with the inverse of the Jacobian determinant of the element mapping to provide the Jacobian stiffening. This procedure is found to be necessary for this simulation. The problem is simulated with ∆t = 1.0 × 10 −3 s up to T = 10.0 s. We set c m = c c = 0 in the study of this example. In our numerical experiences, for compressible materials, setting c m = c c = 0 gives more accurate results. The vertical displacement of the cantilever over time is plotted in Figure 10, and the comparison of the obtained results with those in the literature is listed in Table 1. The results of the new FSI computational framework are in good agreement with the reported results. Snapshots of the simulation results are depicted in Figure 11.

The Greenshields-Weller numerical benchmark
The FSI benchmark example we consider next describes wave propagation in an elastic tube [38]. This benchmark example has been adopted to verify several existing cardiovascular FSI solvers [7,96]. The computational domain consists of a right circular hollow cylinder representing the elastic tube and an inner right circular cylinder representing the fluid domain. The length of the tube is 10 cm, the inner radius for the fluid domain is 1 cm, and the outer radius is 1.2 cm. For the elastic material, the reference density ρ s 0 = 1 g/cm 3 , the Young's modulus E s = 1.0 × 10 7 dyn/cm 2 , and the Poisson's ratio is 0.3. In the original benchmark problem [38], a small-strain linear elastic material is used. Apparently, the smallstrain model is not incorporated in our theory. We adopt the Neo-Hookean model with the volumetric free energy given by (2.48), In doing so, the density is given by ρ = ρ 0 (1 + p/κ), which is identical to the density-pressure relation used in [38]. The fluid is described by the incompressible Navier-Stokes equations with fluid density ρ f 0 = 1 g/cm 3 and dynamic shear viscosityμ = 0.04 poise. The initial conditions for all fields are set to be zero. A step change in pressure is imposed on the inlet surface and the pressure value p in is set to be 5 kPa. At the fluid outlet surface and the lateral surface of the solid, stress-free boundary conditions are applied. The solid velocity in the axial direction on the inlet and outlet surfaces are set to be zero. The geometry and the boundary conditions are illustrated in Figure 12 (a). The solid domain is discretized by structured tetrahedral elements. In the fluid domain, a boundary layer mesh is created with thickness 0.2 cm, and the rest fluid domain is unstructured. The mesh we used is depicted in Figure 12 (b). There are 2.19 × 10 6 elements and 3.76 × 10 5 nodal points in the finite element mesh. The problem is simulated with a fixed time step ∆t = 2.0 × 10 −7 s up to T = 8.0 × 10 −3 s. Since the material is compressible, we set c m = c c = 0 in the study of this example. The mesh motion is given by the harmonic extension algorithm. The analytic value of the wave speed is 8.77 m/s. To obtain the numerically predicted wave speed, we define the wave front as the location where the pressure is 2.5 kPa. Using a linear function to fit the locations of the wave front with respect to time, we obtain that the numerical wave speed is 8.49 m/s, which is 3% lower than the analytic value and 1% lower than the numerical prediction given in [38] (8.58 m/s). With the pressure wave at time 8.0 ms, we may obtain the distance between two peaks is 2.84 cm, and consequently the frequency of the wave is 298.9 Hz. The analytic values for the wave frequency using three different formulas are 269 Hz, 308 Hz, and 336 Hz, respectively [38], and the numerical prediction given in [38] is 318 Hz. Our numerical predicted wave frequency is 10 % higher, 3% lower, and 12% lower than the analytic values, and 6 % lower than the numerical prediction given by [38]. The same problem is also simulated with the Young's modulus E s = 1.0 × 10 6 dyn/cm 2 and p in = 0.5 kPa. In this case, the analytic value of the wave speed 2.77 m/s. Our numerical prediction gives 2.70 m/s, which is 3% lower than the theoretical value. These results demonstrate satisfactory agreement with published results.

Conclusions
We present a comprehensive suite of theoretical and numerical methodologies for solid dynamics, fluid dynamics, and FSI problems. We summarize the contributions in detail as follows.
A novel continuum modeling framework is developed. Our derivation differs from traditional approaches in three aspects. (1) We start with a dynamic model for the continuum body. In the thermodynamic configuration space, a quasi-static process is constrained on an equilibrium manifold in a lower dimension, representing a succession of equilibrium states. This is really an idealized scenario, since it ignores several critical physical effects, and hence does not represent real physics [15,Chapter 4]. (2) We choose the thermodynamic potential as the Gibbs free energy instead of the Helmholtz free energy. The two potentials both give equivalent descriptions of material behavior for any unconstrained process. However, for constrained processes, one should judiciously choose the potential. In our discussion, we are particularly interested in the isochoric process, and the Helmholtz free energy degenerates in this case. A Legendre transformation is performed to transform the independent variable from the specific volume to the pressure. This leads to a theory based on the Gibbs free energy and consequently a pressure primitive variable formulation. This formulation is known to be well-behaved in both compressible and incompressible regimes. (3) Our pressure equation is derived from the mass balance equation instead of the equation of state. In the traditional two-field variational principle, the equation for the pressure field is often introduced as an algebraic equation of state [25,45,81]. We feel that it is unnecessary and improper to put an algebraic equation in a weak form and discretize it by the finite element method (see Remark 5).
Our numerical formulation is designed based on VMS. We provide a formal derivation of the fine-scale model based on the general continuum model. The VMS framework is applied to construct a numerical formulation for finite hyper-elastodynamics; the generalizedα method is applied for time integration; the segregated algorithm introduced in [ The new formulation for hyper-elastodynamics is naturally extended to FSI problems. Due to the general continuum framework we developed, the coupling between viscous fluids and hyper-elastodynamics can be viewed as a uniform continuum problem, where the body is subdivided into subdomains with different material behavior. Over the whole continuum domain, VMS is applied for spatial discretization, and the generalized-α method is applied for time integration. The uniform treatment of FSI problems enjoys several desirable attributes.
(1) The formulation is well-behaved in both compressible and incompressible regimes for both fluids and solids. (2) The pressure instability arising from the equal-order interpolation is handled by the fine-scale modeling in both subdomains. (3) The generalized-α method is applied for the first-order FSI system and naturally achieves optimal high-frequency dissipation in both subdomains. (4) The resulting matrix problems in both subdomains take similar structure, which simplifies the data structure management in the code implementation. Two benchmark problems are simulated to examine the effectiveness of the new FSI formulation, and the results are in good agreement with published results.
There are several promising directions for future work. (1) This new computational framework will be extended to handle nonlinear anisotropic viscoelastic materials to account for more realistic tissue behavior. (2) The FSI computational framework is well-suited to patient-specific cardiovascular biomedical modeling with more physically realistic material models and complex geometries. (3) The concept of the immersed finite element method [16,133] can be utilized with the new FSI formulation developed in this work, which may open a door to novel approaches for complex FSI problems. In summary, a new computational methodology has been established, and it may provide an effective approach to handle problems that cannot be addressed with the previously existing methods.

A Fine-scale approximation
Physically sound assumptions [5,104] suggest that one may represent y by a perturbation series: Here we have to introduce our first assumptions to handle with the residual term o ( y n V ) in the above formulation. We assume that y V < 1 such that the residual term o ( y n V ) is negligible. Grouping the coefficients of the powers of ε leads to a recurrence formula , for 3 ≤ k < n.
Although D y B (w ,ȳ) [·] is a linear operator, it still a challenging task to derive an analytic form of the Green's operator. Therefore, one more approximation needs to be introduced. The most simple but effective approach is to approximate G by wherein τ ∈ R 7×7 [5]. With this formula, we can solve the equations (A.3)-(A.5) and obtain the perturbation series y m = G (G m ) ≈G (G m ) = −τ G m . Therefore, we have y = ∞ m=1 ε m y m ≈ −τ ∞ m=1 ε m G m . The infinite series y needs to be truncated to numerically calculate y . It was shown that y 1 carries most of the subgrid energy [104]. This suggests that one may truncate the perturbation series as y ≈ εy 1 . Then the approximation of the fine-scale component can be explicitly written as y = F (ȳ, Res (ȳ)) ≈F (ȳ, Res (ȳ)) := εy 1 = −τ Res (ȳ) .
In the above, the fine-scale component y is systematically approximated. This approximation formula is utilized to construct the VMS formulations in this work.