A Revisit of Implicit Monolithic Algorithms for Compressible Solids Immersed Inside a Compressible Liquid

: With the development of mature Computational Fluid Dynamics (CFD) tools for ﬂuids (air and liquid) and Finite Element Methods (FEM) for solids and structures, many approaches have been proposed to tackle the so-called Fluid–Structure Interaction or Fluid–Solid Interaction (FSI) problems. Traditional partitioned iterations are often used to link available FEM codes with CFD codes in the study of FSI systems. Although these procedures are convenient, ﬂuid mesh adjustments according to the motion and ﬁnite deformation of immersed solids or structures can be challenging or even prohibitive. Moreover, complex dynamic behaviors of coupled FSI systems are often lost in these iterative processes. In this paper, the author would of composites with deformable aggregates and matrix.


Introduction
Fluid-Structure interactions (FSI) are ubiquitous in many engineering problems. Many research efforts have been invested in the development of modeling methods over the past few decades [1][2][3][4]. In practice, iterations are used between existing finite element methods (FEM) for solids and computational fluid dynamics (CFD) for fluids [5,6]. In engineering practices, such iterative approaches are effective and convenient, yet complex dynamic behaviors of coupled systems can also be filtered out [7,8]. Furthermore, mesh adaption processes such as the arbitrary Lagrangian Eulerian (ALE) and other mesh updating methods can be topologically challenging and computationally prohibitive [9][10][11][12][13]. In this paper, the focus is on the immersed boundary method and its concepts which are initially established in applied mathematics communities. Interface tracking procedures also include level set methods [14,15], immersed interface methods [16,17], and front tracking methods [18][19][20] which have been proposed to eliminate or bypass numerical problems associated with large motions of immersed objects [21][22][23][24][25]. Immersed methods do provide a very direct and simple coupling between immersed deformable solids with the surrounding fluid without specific interface tracking, which enable effective and efficient exploration of phase and parametric spaces in the design of complex FSI systems [26][27][28][29]. Since its inception [30][31][32], the immersed boundary method has been employed for the solutions of wave propagation in cochlea [33], swimming motions of marine worms [34], wood pulp fiber dynamics [35], and biofilm processes [36]. In fact, it has been discovered that the modeling philosophy of immersed methods is similar to the fictitious domain method [28,[37][38][39]. In fact, such a principle can also be implemented to handle deformable solids immersed in another deformable solid medium as discussed in Ref. [40]. In early developments of immersed boundary methods, immersed elastic fibers are employed to construct various objects. As proposed in Ref. [41], torsion, bending, and transverse shear effects have also been added to immersed structures, namely beams instead of fibers. In this paper, the author would like to revisit immersed methods for compressible FSI systems ranging from acoustic fluid-solid interactions to compressible (or rather nearly incompressible) liquid interacting with deformable structures. In particular, mixed finite element formulations are used for both immersed solids and background fluids [42,43]. The communication between immersed solids and background fluids is accomplished with various kernels or radial-based interpolation functions widely implemented for meshless methods [44]. Such delta functions provide a higher order smoothness and can be extended to non-uniform fluid grids [45][46][47]. The focus of this paper is on some earlier results presented by the author and his collaborators [41,43,48,49], in particular those of acoustic fluid-solid interaction systems for which mixed finite element formulations and mixed finite elements which either satisfy the Inf-Sup conditions or pass the numerical Inf-Sup tests are essential in lieu of immersed methods. Similar immersed methods have already been implemented to electrokinetics and mass and heat transfers [48,50,51]. A fairly comprehensive review from the perspective of a merger between traditional aerodynamics and the immersed boundary method is available in Ref. [52].

Mixed Finite Element Formulations
The key feature of immersed methods is to allow independent Lagrangian mesh of immersed solids to move on top of the background fluid which consists of a fixed Eulerian mesh. However, in order to be consistent with the actual physics, it is important to combine the immersed domain (Ω s ) and the physical exterior domain (Ω f ) into a complete background domain (Ω) as illustrated in Figures 1 and 2. It is this complete domain that stays the same for any instantaneous position of the immersed continuum. Thus, there will be no need for frequent mesh adjustments. In fact, a background uniform mesh, namely ghost mesh is often used to effectively locate the vicinity of the immersed solids with modulo functions.   As illustrated in Refs. [53,54], the underlining principle of immersed methods is the energy conservation. In summary, the energy and power inputs to the fluid domain introduced by the immersed solids or immersed structures are identical to those from the equivalent body forces. To secure this key feature, the same delta function or kernel function must be used in both the distribution of the resultant nodal forces and the interpolation of the solid velocities based on the surrounding or rather background fluid velocities. In fact, this treatment in immersed methods can be viewed as the synchronization of the fluid motion with the solid motion within the immersed solid domain Ω s , namely v s = v f . This important feature is also very similar to what is employed in the fictitious domain methods [37,38] along with the so-called distributed Lagrange multipliers as the equivalent body forces. Although in theory, the equivalent body forces can be directly calculated along with independent fluid and solid velocity vectors [55], the coupling of equivalent body forces and solid motions is implicit and nonlinear. Therefore, a matrix-free Newton-Krylov iterative procedures must be established to link both fluid and solid domains. Moreover, if the surrounding fluid is viscous and incompressible, the immersed solid or immersed structure must be incompressible in immersed methods and vice versa [56][57][58][59].
In this work, we review the implemented mixed finite element formulations for both compressible immersed solids or rather nearly incompressible solids and compressible liquid or rather nearly incompressible fluid. First, the weak form of the fluid-solid system can be written as ∀w ∈ [ with In continuum model for Newtonian fluids, the stress components σ ij is a combination of the pressure p and the deviatoric stress components τ ij , Furthermore, the continuity equation of the compressible viscous fluid can be expressed as where κ is the bulk modulus of the fluid. Notice here we focus on liquid or more precise liquid phase of water as working fluid in the context of acoustic fluid-solid interactions namely the pressure variations are not strong enough to alter the water density in any significant fashion as in the case of underwater detonation where strong shock will trigger both phase transition of water and density changes must be included. Moreover, from dp dρ we have ∆ρ For water, the bulk modulus is around 2.1 × 10 9 N/m 2 and the density is around 1.0 × 10 3 kg/m 3 . It is clear that as long as the pressure variation is within the ambient pressure which is 1.0 × 10 5 N/m 2 . The change of the water density is very minimum and should be ignored, so is the density of the immersed solid. Notice that when we deal with the air as the working fluid, fully-fledged energy conservation form of governing equations must be used along with non-oscillatory time incremental schemes [60,61].
Likewise, for immersed solids, the solid stress components σ s ij can be decomposed as the pressure p s and the deviatoric stress components τ s ij , In nonlinear mechanics, the solid deformation gradient F ij = ∂x s i (t)/∂x s j (0) is introduced along with the Green-Lagrangian strain component ij , also defined in the tensor form as (F T F − I)/2 with an identity matrix I. Consequently, the energy conjugate stress S ij of the Green-Lagrangian strain, namely the second Piola-Kirchhoff stress can be derived from the elastic energyW, which is often related to the three invariants of the Cauchy-Green deformation tensor C also defined in the tensor form as F T F.
We must also point out that the solid domain Ω s and the material point position x s all refer to the current solid configurations, and for clarity they could be denoted as Ω s (t) and x s (t), respectively. To properly introduce the material constitutive laws, we must first introduce elastic energyW, which is often related to the invariants of the Cauchy-Green deformation tensor C. A typical hyperelastic material can be described by the Mooney-Rivlin material model with large displacements and large strains. Hence, the strain energy potentialW at time t, is defined as where the invariants J i are functions of the invariants I i , namely J 1 = I 1 I −1/3 3 , J 2 = I 2 I −2/3 3 , and J 3 = I 1/2 3 , with I 1 = C kk , I 2 = [(I 1 ) 2 − C ij C ij ]/2, and I 3 = det(C). In the formulations for acoustoelastic FSI systems, in order to avoid locking issues due to high bulk modulus for almost compressible solids, an additional elastic energy termQ, or −[p s + κ s (J 3 − 1)] 2 /2κ s is added to the elastic energyW, along with solid unknown pressure p s . In the mixed finite element formulation, p s is an independent pressure unknown completely different from the background fluid pressure p, where κ s is the solid bulk modulus and J 3 stands for the determinant of the deformation gradient.
To delineate this pressure unknown from the pressure calculated based on the volumetric strain and the bulk modulus, we introducē Naturally, the continuity equation becomes a constraint between the pressure unknown expressed as p s and the pressure as the function of the volumetric strain expressed asp. In the displacement-based finite element formulation, such a constraint is strictly enforced, with p s =p. However, for almost incompressible materials with large bulk modulus, which covers many soft materials such as biological materials, the displacement-based finite element formulation will introduce so-called checkerboard pressure distributions and need to be replaced with mixed finite element formulations or combined with reduced integration techniques. Therefore, the second Piola-Kirchhoff stress can be expressed as with W =W +Q.
Of course, to relate to the expression in Equation (7), the solid Cauchy stress can often be recovered from the second Piola-Kirchhoff stress, which is an energy conservative dual of the Green-Lagrangian strain, In immersed methods, within the immersed domain, the solid displacement is dependent on the velocity of the fluid occupying the same immersed domain. Hence, the primary unknowns for the coupled fluid-solid system are the fluid velocity v, the fluid pressure p, and the solid pressure p s . Define the Sobolev spaces, the weak form of governing equations can be expressed as: ∀q ∈ L 2 (Ω), q s ∈ L 2 (Ω s ), w ∈ [H 1 0,Γ v (Ω)] d , which includes ∀w s ∈ [H 1 (Ω s )] d , and find v and p ∈ Ω, p s ∈ Ω s , such that For the fluid domain the following interpolations are introduced for the entire domain Ω : where N v I and N p I are the interpolation functions at node I for the fluid velocity vector and the fluid pressure; and v I , w I , p I , and q I represent the nodal values of the discretized fluid velocity vector, admissible fluid velocity variation, fluid pressure, and admissible fluid pressure variation, respectively.
Please note that the interpolation functions for the velocity vector unknowns in general are denoted by the superscript v, the independent pressure unknowns for immersed solids are denoted by p s and the independent pressure unknowns for background fluids are denoted by p.
Likewise for the solid domain Ω s , the discretization is based on the following: Therefore, in the entire domain Ω, the arbitrariness of the fluid velocity variations w iI , fluid pressure variations q I , and independent solid pressure variations q s J yields four equations for three velocity components and one pressure at each fluid node denoted as I and one equation for the independent solid pressure at each solid node denoted as J, where the residuals are defined as Please note that in both Equations (15) and (17) an implicit functionÑ is introduced to map from w I to w s J , namely from the fluid mesh to the solid mesh, with the kernel functions. It is the implicit and nonlinear mapping between the immersed solid and the background fluid that consolidates the solid displacement unknowns with the fluid velocity unknowns. This nonlinear mapping can be better illustrated in specific nestings of subroutines although it is difficult to assign explicit analytical expression. Furthermore, a mapping functionÑ represents the mapping from the fluid domain implemented with background ghost mesh to the immersed sold domain the position of which is identified by the modulo function and projection will make that mapping efficient and straightforward. Furthermore, both Equations (12) and (15) are fully-fledged nonlinear form of the governing equation before any linearization. In fact, in the monolithic approach, the linearization will be coupled with the geometrical increment of the immersed solids in the program. The geometrical nonlinearity of the moving structures/solids are hidden in the so-called solid domain Ω s which is also occupied by the fluid. Therefore, the immersed solid in computation is the original solid subtract the background fluid which is physically not present as depicted in Figure 1 and in our monolithic implementation will also be constrained to follow the same kinematic motions as the immersed solid though with distinct fluid material constitutive laws identical to the surrounding fluid.
The isentropic compressible fluid formulation adopted in this paper with a constant density is in fact the same mixed finite element formulations developed for acoustoelastic FSI systems [62,63]. With sufficient mesh resolution and time step, pressure wave propagation and related scattering and radiation phenomena can be modeled. In general, small strain and small deformation assumptions are introduced for the acoustoelastic FSI systems. Therefore, immersed methods which have the advantage of automatically tracking the FSI interfaces is not necessary for this type of linear FSI systems [53,54]. However, with immersed solids, in nonlinear acoustic FSI model with finite interfacial motions, the monolithic implicit Newton-Krylov iterative procedure with immersed methods becomes essential. Furthermore, in the acoustoelastic FSI with free surface as illustrated in Figure 3, we have also demonstrated that with the same mixed finite element formulation for both acoustic fluid and immersed two-dimensional solid mimicking the immersed flexible structure, free surface modes No. 1 and 2 representing low frequency sloshing modes, solid modes No. 3 and 4 with intermediate frequency coupled or wet structural modes, along with clearly high and low pressure on top and bottom of the immersed structure, and acoustoelastic modes No. 5 and 6 with constant pressure on free surface can be evaluated in one monolithic coupled system eigenvalue problem. Moreover, with numerically stable mixed finite element formulation and mixed elements satisfying the Inf-Sup conditions, even a very coarse mesh can yield insightful results for fully coupled FSI systems. Naturally, the mesh density must be significantly higher if our aim in the simulation is to capture the high frequency range such as the acoustoelastic FSI models instead of the lower frequency ranges such as the wet structural modes or slowing modes as illustrated in Figure 3. The simple fact that a fairly coarse mesh as depicted in Figure 3 could capture coupled frequencies and corresponding modes from low to high ranges demonstrates the promise of this monolithic coupled code with mixed finite element formulations for immersed compressible solids and compressible (or rather nearly incompressible) liquid with constant density.

Stability Analysis
For nonlinear FSI analysis discussed in this paper, the incremental iteration within each time step involves a linearized transient analysis, hence, after transformation, the following equation is considered: where R * h is an effective load vector, (K * uu ) h and (K us ) h represent the effective stiffness matrix block corresponding to displacement unknowns and pressure unknowns, respectively [4].
The stability of such mixed finite element formulations with corresponding mixed elements is based on ellipticity and Inf-Sup conditions, which are necessary and sufficient conditions for well-posedness.
(a) Ellipticity Condition: where where the constant c 2 is independent of the mesh size h and the bulk modulus.
The Inf-Sup condition as depicted in (20) for mixed finite element formulations is elegant and useful, although the analytical proof for specific mixed elements and systems can be very difficult [64]. In practice, numerical Inf-Sup tests as discussed in [4,65] with eigenvalues or singular values associated with the coupling matrix K us are often employed.

Kernel Functions
The discretized delta functions implemented in the original Immersed Finite Element Methods proposed in Refs. [39,58,59] are the same as those in the so-called reproducing kernel particle methods (RKPM) or meshless finite element methods [44,[66][67][68][69]. In immersed methods, as pointed out in Refs. [53,54], to ensure the energy conservation, we must employ the same discretized delta function for the distribution of solid nodal forces and the interpolation from the background fluid velocities to the immersed solid velocities. Moreover, the modified window function can also be used in the discretized delta function for non-uniform meshing for the background fluid domain. Both wavelet and smooth particle hydrodynamics (SPH) methods belong to the same class of reproducing kernel methods where the "reproduced" function u R (x) is derived as with a projection operator or a window function φ(x).
In essence, the window function is required to be flatter at ω = 0 in the Fourier domain as the order of reproducing gets higher. In Figure 4, we provide a comparison of various kernel functions in function and spectrum domains. With the increase of the reproducing order, for the same finite support region, the kernel functions approach to the ideal low pass filter.  As presented in detail in Ref. [58], consider first a one-dimensional case, in accordance with the translation invariance [31], for all r, where r is the parameter representing the position of the submerged boundary point and is scaled with respect to the grid size h, the discretized delta function satisfies ∑ j φ 2 (r − j) = C, where C is a numerical constant.
In addition, to uniquely define the discretized delta function for all r, we also have ∑ j (r − j) m φ(r − j) = 0, where the selection of the m th moment depends on the number of support points. For instance, the discretized delta function with four support points is uniquely defined by the following:
For all r, ∑ j even 3.
For all r, ∑ j φ 2 (r − j) = C, where C is a numerical constant.
In general, for 0 < r < 1, the discretized delta function φ(r − j) covers four non-zero support points. However, for the degenerate case of the 4-point discretized delta function centered at r = 0, we have five support points, namely r − j = −2, −1, 0, 1, 2. From the degenerate case, we can easily derive the constant C. Hence, we obtain the following four admissible branches of solutions for 0 < r < 1, Naturally, in the three-dimensional case, one of the smoothed approximations to the delta function is given by It is interesting to point out that the discretized delta function in Equation (22) is very close to (1 + cos πr/2)/4, with r ∈ [−2, 2]. Moreover, it is easy to confirm that the discretized delta function φ(r), with r ∈ [−2, 2], defined in Equation (22), has C 1 continuity [70].
It is noted that in the current version of the implicit code the background ghost mesh is uniform which enables a very speedy estimate of the immersed solid location, and the order of accuracy is no more than two. Further work on the improvement of accuracy with adaptive meshing and much sharper FSI interface capturing techniques is discussed in Ref. [71].

Compressible Fluid Model for General Grids
In this section, unlike the early presentation of compressible fluid with working fluid as liquid, more specifically liquid form of water, in which the density is kept as constant due to the nearly compressible condition or rather large bulk modulus in the range of GPa. In this section, we revisit some early attempt to expand the immersed boundary method to truly compressible fluid with working fluid as compressible air. Naturally, it is a common practice to adopt the conservative forms of governing equations for background compressible fluid domain. Notice that we adopt the same thermodynamic concepts with an increment of the specific enthalpy dh as c p dT where c p is the specific heat for constant pressure and T is the temperature in absolute scale, namely Rankine or Kelvin along with an increment of the specific internal energy du as c v dT where c v is the specific heat for constant volume. Define the total specific mechanical energy as e, we have e = u + v i v i /2. Notice for ideal gas, we have p = ρRT with c p − c v = R and the gas constant R is defined as R u /M where R u is the universal gas constant and M is the molar mass of the gas. Hence, the pressure p can also be written as p = ρ(γ − 1)e = ρ(γ − 1)(e − v i v i /2) with γ = c p /c v . Define the reciprocal of the density ρ as the specific volume v, the first law of thermodynamics yields where dq is the specific heat transferred to the system from the environment. For isentropic (quasistatic, adiabatic, and reversible) gas, with dq = 0, we have Thus, combine with the ideal gas law, pv = RT, we have Consequently, for isentropic gas, the pressure p is governed by p = ρ γ C, where C is a constant. Furthermore, the speed of sound c is defined as c 2 = dp/dρ or γp/ρ = γRT. Finally, use total specific enthalpy form, h = e + p/ρ = c 2 /(γ − 1) + v i v i /2, and the Fourier law for thermal conduction, the conservation of energy can be written as where q s stands for the rate of heat transfer from the surrounding environment to the system or control volume per unit volume, a so-called heat source.

Matrix Operation Specifics
To depict clearly the monolithic implicit immersed method with Newton-Krylov iterations, we start with the basic linear algebra concepts for the solution of general linear system of equations and the solution of eigenvalue problems. Assume A is a matrix with m rows and n columns A{φ 1 , φ 2 , . . . , φ (k−1) , φ k } = {Aφ 1 , Aφ 2 , . . . , Aφ (k−1) , Aφ k } the multiplication Aφ j represents the linear combination of n columns of the matrix A using n entities of the vector φ j with φ j ∈ R n and 1 ≤ j ≤ k; namely where A i represents the ith column of the matrix A. Denote B as a matrix with m rows and n columns the multiplication ψ T i B represents the linear combination of m rows of the matrix B using m entities of the vector ψ i with ψ i ∈ R m and 1 ≤ i ≤ k; namely, where B j represents the jth row of the matrix B.
Consequently, reduced row echelon form (rref) of the matrix A after the left multiplication of the elementary matrix E, namely the row combination and manipulation operation. Assume the rank of the matrix A is r, the number of pivot variables is r, the number of free variables is n − r, we derive In general, if we have r < m, n, the identity matrix corresponding to the pivot variables I 1 has the size r × r, the free variable matrix F has the size r × (n − r), the zero matrix 0 1 has the size (m − r) × r, and the zero matrix 0 2 has the size (m − r) × (n − r). Introduce the new identity matrix I 2 with the size (n − r) × (n − r), the null space vectors can be expressed as −F I 2 ; where the last m − r rows of the transformation matrix E, often called elementary matrix, are the left null vectors. If r = m ≤ n, we have the full column rank, zero matrices 0 1 and 0 2 do not exist, namely there is no left null vector, and F has the size r × (n − r), or m × (m − r), whereas if r = n ≤ m, we have the full row rank and the zero matrix 0 1 has the size (m − r) × r, or (m − n) × n, namely there is no null vector.

Matrix-Free Newton-Krylov Iteration
In the matrix-free Newton-Krylov iteration, the target is the solution vector ∆Θ, or rather ∆V for fluid velocity unknowns, ∆P for fluid pressure unknowns, and ∆P s for independent solid pressure unknowns. In this monolithic implicit immersed method with iterative solutions based on Newton-Krylov and generalized minimal residual method (GMRES) along with the numerical construction of the Jacobian matrix in the Newton-Raphson iteration, no specific analytical equation can be established other than an abstract operatorÑ which is implemented in recursive, iterative, and nested subroutines, the details of which are also presented and elaborated in the research manuscript by the author and his collaborator Professor Lucy Zhang in the Department of Mechanical Engineering of Rensselaer Polytechnic Institute [72]. We start our matrix-free Newton-Krylov with A set of orthonormal vectors v i , with 1 ≤ i ≤ n in the Krylov space K n and an (n + 1) × n upper-Hessenberg matrixH n . If we define V n = (v 1 , v 2 , . . . , v n ) and V n+1 = (v 1 , v 2 , . . . , v n+1 ), and the remaining process in the GMRES method is to solve the least square problem: min z∈K n p − Jz or min y∈R n p − JV n y . (29) Assuming γ is the length of the initial residual vector p and e 1 is the unit vector, the first column of (n + 1) × (n + 1) identity matrix, γ = p 2 , b = γe 1 , preconditioning matrix Λ, for i = 1 to n, using the modified Gram-Schmidt orthogonalization process; we then have q i = Λ −1 v i and w = Jq i , and for j = 1 to i, we have h ji = w T v j and w is updated with w − h ji v j . Consequently, we obtain After we establish the elements of an upper n × n Hessenberg matrix H n as well as an upper (n + 1) × n Hessenberg matrixH n , for j = 1 to n, and i = 1 to j − 1, a factorization of H n is carried out where the entities of the rotation processes are calculated as Through this rotation process, the upper Hessenberg matrix is converted to a diagonal matrix with the coefficients defined as: for j = 1 to n h jj = r, p j = c j b j , and p j+1 = −s j b j . (34) Finally, the termination criteria: if |b n+1 | < , the solution vector ∆Θ, or rather ∆V f , ∆U c , ∆U s , ∆P f , ∆P c , and ∆P s is expressed as or Notice here for time increments, we use the trapezoidal temporal discretization is as the following, and consequently, and more importantly, where b can stand for the fluid interior nodal unknowns V f , the fluid-solid interface displacement unknowns U c , the solid interior nodal unknowns U s , the fluid interior pressure unknowns P f , the FSI coupled interface pressure unknowns P c , and the solid interior pressure unknowns P s . Furthermore, the solid velocities are based on the projections of the background fluid velocities within the finite support region. Thus, the solid displacements at the FSI interface and the interior are integrated with the time integration scheme as stipulated in Equations (36) and (38). To be more specific, the update processes occur in very incremental iteration within one time step, which means, the positions of the immersed structures/solids are updated in all iteration within the matrix-free Newton-Krylov iteration.
For specific elaboration of the combination of Newton-Raphson iteration with Newton-Krylov iteration, we introduce therefore, for k = 1 2 or rather −s 1 β is sufficiently small, the second equation will be approximately satisfied, and the iteration will stop at the level when k = 1.
Moreover, for k = 2 2 is sufficiently small, the third equation will be approximately satisfied, and the iteration will stop at the level when k = 2.
Finally, for k = 3 3 is sufficiently small, the third equation will be approximately satisfied, and the iteration will stop at the level when k = 3.

Results and Verifications
In this section, we present several numerical examples to illustrate the efficacy and accuracy of the monolithic implicit immersed method with matrix-free Newton-Krylov iterations. The software developed through this work could be available upon request.
A model problem is introduced with a one-dimensional continuum within [L 1 , L 2 ] immersed in another one-dimensional continuum within [0, L] as shown in Figure 5. The physical materials occupying region [0, L 1 ] ∪ [L 2 , L] for one continuum and region [L 1 , L 2 ] for another continuum have a constant cross-section area A and are subjected to the gravitational acceleration g in x direction. Thus, the governing equations can be written as (39) Figure 5. One-dimensional model problem of one continuum immersed in another continuum.
Equation (39) represent a set of two second-order partial differential equations (PDE) for both continuum with three connected regions, namely [0, L 1 ], [L 2 , L], and [L 1 , L 2 ]. If we consider the boundaries marked by L 1 and L 2 dependent of time t and the displacement u, this problem will be non-trivial. For extreme conditions, let's look at the corresponding static sets of equations, namely the right hand sides of Equation (39) which represent a set of two second-order ordinary differential equations (ODE). Finally, six constants c i will be introduced for all three regions, namely (0, L 1 ), (L 2 , L), and (L 1 , L 2 ), and the solution u(x) can be simply expressed as Therefore, six constants c i will be derived based on the following six boundary conditions In this model study, for simplicity, we assign L 1 = L/3 = l, L 2 = 2L/3 = 2l, E s = aE f = aE, and ρ s = bρ f = bρ. Moreover, we scale the three constants c i , i = 1, 2, 3 with the same parameter l/E. Finally, with β = ρgl, the solution of the six unknown constants c can be determined from the boundary conditions (41) as The comparisons with analytical solutions as well as existing computational fluid and solid software package such as ADINA as illustrated in Figures 6 and 7 confirm the validity of the monolithic implicit immersed method with matrix-free Newton-Krylov iterations [39]. In addition, this simple one-dimensional model problem also confirms the subtle shift in math and engineering fields in which the background compressible fluids are often replaced with the background compressible solids [55].
A clear validation and verification example is the driven cavity for viscous fluid with a dimension of 1 × 1 m as illustrated in Ref. [39]. To compare the dynamical behaviors, the top shear velocity of the cavity is set as 0.1 sin(2π/40t) m/s. The deformable solid is situated initially in the center of the cavity with zero velocity. To avoid the overlap of the mesh for the immersed solid and the background fluid mesh, we choose the typical immersed solid with the dimension of 0.13 × 0.13 m. The submerged solid is made of a compressible rubber material with the material constant κ s = 1 × 10 7 N/m 2 and the density ρ s = 1000 kg/m 3 . For hard solids we use C 1 = 200 and C 2 = 100 N/m 2 . The case with hard solids resembles the same driven cavity problem with immersed rigid bodies [73]. In addition, the viscous fluid is represented with compressible model with a constant wave speed. Notice that in this acoustoelastic FSI example, the Mach number is much smaller than 0.01 just as the case presented in Ref. [74]. Moreover, we have the dynamic viscosity µ = 1 Pa · s, the bulk modulus κ f = 2.1 × 10 9 N/m 2 , and the density ρ f = 1000 kg/m 3 .  It is evident that even with the coarse mesh used, the developed formulation with highorder elements provides reasonable results comparable to a reference solution. In addition, no spatial oscillation and checkerboard pressure bands are observed, as demonstrated in Figure 8. Of course, the proposed description can be further improved with refined meshes. As one of the main messages, we must point out that in engineering practice, before a large number of finite elements are used, it is always beneficial to employ coarse meshes with high-order elements to obtain a reasonable estimation of complicated problems.
In the implementation, the background fluid mesh is constructed for the entire cavity, which includes the space occupied by the immersed solids, consists of 20 × 20 9/4c mixed finite elements as illustrated in detail in Refs. [4,64]. A typical immersed solid is represented by 4 × 4 9/4c mixed finite elements. In the post computation presentation, we split each 9-node element into four 4-node elements. As a result, the velocity mesh will be two times denser than the pressure mesh. As shown in Figure 8, the mesh construction clearly demonstrates the philosophical difference between the monolithic implicit immersed method with matrix-free Newton-Krylov iterations and the traditional ALE approaches. In our method, the solid mesh is floating right on top of the background fluid mesh, whereas in the ALE formulation, the solid mesh is surrounded by the fluid mesh with a different mesh density. For the case with one immersed square solid, as long as the immersed solid does not move and deform significantly, it is still possible to solve this fluid-solid system with traditional modeling techniques such as the arbitrary Lagrangian-Eulerian (ALE) formulation. For the comparison of fixed immersed solid, we compare the velocity components of the midpoint (0.5, 0.5) of the immersed solid. With the same time step size, as shown in Figure 9, even with the coarse mesh of high-order mixed elements, the results from the monolithic implicit immersed method with matrix-free Newton-Krylov iterations and the ADINA FSI solver are very close. Considering the fact that these results are derived from two completely different approaches with entirely different meshes, these comparisons are very assuring. Moreover, for the case with the center of the immersed solid tethered with a soft spring (4 N/m), despite the difference of the fluid mesh shown in Figure 8, the trajectories of the block center derived from the immersed continuum method and the ADINA FSI solver do match well with each other as depicted in Figures 9 and 10. Further studies with different periods also demonstrate in Figure 10 that a slightly larger numerical viscous exists in comparison with the traditional computational approaches.
It is very important to point that the discretized delta function provides the possibility of linking the immersed solid node with the surrounding fluid nodes as done in meshless finite element methods. However, as pointed out in Ref. [56], the immersed solid node can also be restricted to communicate with the very finite element for the background fluid domain. As long as the virtual power input to the surrounding fluid is preserved, the concept of immersed methods will stay the same. This procedure is the same as the standard finite element procedure to distribute concentrated forces and to interpolate displacement/velocity within the element. Whether or not this procedure is as efficient as the meshless type of communication depends very much on the search algorithm for the purpose of locating the immersed solid nodes.  The benefit of the immersed methods is clear for the case with five immersed solids as illustrated in Figure 11. For this case, also depicted in Ref. [39], it is no longer feasible to use the ALE formulation, whereas it is relatively straightforward to add a few more deformable solids in the monolithic implicit immersed method with matrix-free Newton-Krylov iterations. The only possible comparison came from Professor Tsorng-Whay Pan in the Department of Mathematics of the University of Houston. The simulation animation can be accessed in Ref. [73].
It is very important to emphasize that in this monolithic implicit matrix-free Newton-Krylov iterative solution method, the fluid velocity, the immersed solid velocity, the fluid pressure, and the solid pressure must converge simultaneously. In fact, since we obtain the immersed solid velocity based on the projection of the velocity of fluid occupying the same position, namely their simultaneously convergence is guaranteed. However, the fluid pressure and the solid pressure are calculated with two completely different procedures and the concept of addition and subtraction as illustrated in Figure 1, which ultimately marks the success of the entire monolithic implicit immersed method with matrix-free Newton-Krylov iterative solution method. As shown in Figures 12 and 13, a typical converged data sheet along with two sets of pressure and velocity contour and band plots are presented for a simple sedimentation and rotation of a compressible solid in a viscous and slightly compressible fluid.  In Figure 14, we present two time snapshots taken from a case study reported in Ref. [61] in which a submerged flexible structure is interacting with the surrounding compressible fluid, more specifically, air. In this case, fully-fledged energy conservation form of governing equations in aerodynamics based on conformal mapping and compact schemes have been coupled with non-oscillatory time incremental schemes [60,61]. The specific details are imbedded in the proprietary code in Wright-Patterson Air Force Base. In order to illustrate the efficacy of immersed boundary methods for compact schemes, we intentionally set the boundaries to reflect the shock waves. It is clear that shocks are captured nicely as well as their interactions with the immersed structure.
In these implementations, we have a uniform background (ghost) mesh which allow us to employ the simple modulo function to locate the immersed solid or structure points and to implement the displacement projections and force distributions. This feature along with the concept of immersed methods enable us to avoid the topologically challenging mesh adjustments which are often prohibitive. The uniform background (ghost) mesh adopted for the compact scheme is also the important attributes for high accuracy and high computational efficiency. These features also naturally fit into parallel processing commands such as MPI and further vectorization of the operation.

Conclusions and Discussion
In this study, we confirm that immersed boundary methods can be applicable to compressible liquid interacting with immersed compressible solid in acoustic FSI problems. It seems equally valid for shock wave propagating in compressible fluid and acoustic waves interacting with both compressible fluid and immersed solid with the help of compact schemes.
In the monolithic implicit immersed method with Newton-Krylov iterations, in order to satisfy energy conservation, namely the energy input to the fluid domain from the immersed solid is the same as that from the equivalent body force, the same delta function must be used in both the distribution of the resultant nodal force and the interpolation of the solid velocity based on the surrounding fluid velocities. In fact, the key treatment in this fully implicit and monolithic immersed method, the synchronization of the fluid motion occupying the same geometric locations marked by the domain Ω s with the immersed solid, namely v s = v f . Although no explicit mathematical expression is possible in the implicit monolithic algorithms for compressible solid immersed inside a compressible fluid. However, an operatorÑ is introduced in Equation (17) to summarize such a nested fully implicit Newton-Krylov iterative procedure. In this procedure, the iterative solution based on Newton-Krylov and generalized minimal residual method (GMRES) along with the numerical construction of the Jacobian matrix in the Newton-Raphson iteration are combined as a single iterative procedure which is elaborated in detail in this paper. Furthermore, preliminary results derived with two-and three-dimensional compressible solver based on compact schemes are also included.
The constraint of Equation (42) introduces the (distributed) Lagrange multiplier as the equivalent body force. In this case, the equivalent body forces can be directly calculated along with independent fluid and solid velocity vectors. Of course, such a procedure will introduce a set of new unknowns equal to the number of velocity unknowns for solids. The early success of immersed finite element formulations for the type of immersed methods is only the beginning of this type of modeling treatment for complex FSI systems. Recent mathematical studies and extensions to electrokinetics and mass and heat transfers can also be found in Refs. [48,51].