Next Article in Journal
The Eulerian–Lagrangian Approach for the Numerical Investigation of an Acoustic Field Generated by a High-Speed Gas-Droplet Flow
Next Article in Special Issue
Radial Basis Functions Vector Fields Interpolation for Complex Fluid Structure Interaction Problems
Previous Article in Journal
Nanoparticle Delivery in Prostate Tumors Implanted in Mice Facilitated by Either Local or Whole-Body Heating
Previous Article in Special Issue
A Monolithic and a Partitioned, Reduced Basis Method for Fluid–Structure Interaction Problems
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

McCoy School of Engineering, Midwestern State University, 3410 Taft Blvd., Wichita Falls, TX 76308, USA
Submission received: 10 June 2021 / Revised: 16 July 2021 / Accepted: 22 July 2021 / Published: 3 August 2021
(This article belongs to the Special Issue Fluid Structure Interaction: Methods and Applications)

Abstract

:
With the development of mature Computational Fluid Dynamics (CFD) tools for fluids (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, fluid mesh adjustments according to the motion and finite 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 like to review the so-called monolithic approaches for the solution of coupled FSI systems as a whole in the context of the immersed boundary method. In particular, the focus is on the implicit monolithic algorithm for compressible solids immersed inside a compressible liquid. Notice here the main focus of this paper is on liquid or more precisely liquid phase of water as working fluid. Using the word liquid, the author would like to emphasize the consideration of the compressibility of the fluid and the assumption of constant density and temperature. It is a common practice to assume that the pressure variations are not strong enough to alter the liquid density in any significant fashion for acoustic fluid–solid interactions problems. Although the algorithm presented in this paper is not directly applicable to aerodynamics in which the density change is significant along with its relationship with the pressure and the temperature, the author did revisit his earlier work on merging immersed boundary method concepts with a fully-fledged compressible aerodynamic code based on high-order compact scheme and energy conservative form of governing equations. In the proposed algorithm, on top of a uniform background (ghost) mesh, a fully implicit immersed method is implemented with mixed finite element methods for compressible liquid as well as immersed compressible solids with a matrix-free Newton–Krylov iterative solution scheme. In this monolithic approach, with the simple modulo function, the immersed solid or structure points can be easily located and thus the displacement projections and force distributions stipulated in the immersed boundary method can be effectively and efficiently implemented. This feature coupled with the key concept of the immersed boundary method helps to avoid topologically challenging mesh adjustments and to incorporate parallel processing commands such as Message Passing Interface (MPI) and further vectorization of the numerical operation. Once these high-performance procedures are implemented coupled with the monolithic implicit matrix-free Newton–Krylov iterative scheme with immersed methods, effective and efficient reduced order modeling techniques can then be employed to explore phase and parametric spaces. The in-house developed programs are at the moment two-dimensional. Furthermore, based on the same approach implemented in one-dimensional test example with one continuum immersed in another continuum, such monolithic implicit matrix-free Newton–Krylov iterative approach can be extended for the study of composites with deformable aggregates and matrix.

1. 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].

2. Formulations and Theories

2.1. 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 Figure 1 and Figure 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 [ H 0 , Γ v 1 ( Ω ) ] d
Ω [ w i ρ f ( v ˙ i g i ) + w i , j σ i j ] d Ω Γ f w i f i Γ f d Γ Ω s w i s f i s d Ω = 0 ,
with
Ω s w i s f i s d Ω = Ω s [ w i ( ρ s ρ f ) ( v ˙ i g i ) + w i , j ( σ i j s σ i j f ) ] d Ω .
In continuum model for Newtonian fluids, the stress components σ i j is a combination of the pressure p and the deviatoric stress components τ i j ,
σ i j = p δ i j + τ i j ,
with τ i j = μ ( v j , i + v i , j ) .
Furthermore, the continuity equation of the compressible viscous fluid can be expressed as
v i , i + p ˙ κ = 0 ,
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
d p d ρ = c 2 = κ ρ ,
we have
Δ ρ ρ = Δ p κ .
For water, the bulk modulus is around 2.1 × 10 9   N / m 2 and the density is around 1.0 × 10 3   k g / 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 σ i j s can be decomposed as the pressure p s and the deviatoric stress components τ i j s ,
σ i j s = p s δ i j + τ i j s .
In nonlinear mechanics, the solid deformation gradient F i j = x i s ( t ) / x j s ( 0 ) is introduced along with the Green-Lagrangian strain component ϵ i j , also defined in the tensor form as ( F T F I ) / 2 with an identity matrix I . Consequently, the energy conjugate stress S i j of the Green-Lagrangian strain, namely the second Piola–Kirchhoff stress can be derived from the elastic energy W ¯ , 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 energy W ¯ , 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 potential W ¯ at time t, is defined as
W ¯ = C 1 ( J 1 3 ) + C 2 ( J 2 3 ) + κ s ( J 3 1 ) 2 / 2 ,
where the invariants J i are functions of the invariants I i , namely J 1 = I 1 I 3 1 / 3 , J 2 = I 2 I 3 2 / 3 , and J 3 = I 3 1 / 2 , with I 1 = C k k , I 2 = [ ( I 1 ) 2 C i j C i j ] / 2 , and I 3 = d e t ( 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 term Q ¯ , or [ p s + κ s ( J 3 1 ) ] 2 / 2 κ s is added to the elastic energy W ¯ , 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,
J 3 1 + p s κ s = 0 ,
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 introduce
p ¯ = κ s ( J 3 1 ) .
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 as p ¯ . 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
S k l = W ϵ k l = W ¯ ϵ k l [ p s + κ s ( J 3 1 ) ] J 3 ϵ k l ,
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,
σ i j s = 1 det ( F ) F i , m S m n F j , n .
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 0 , Γ v 1 ( Ω ) ] d , which includes w s [ H 1 ( Ω s ) ] d , and find v and p Ω , p s Ω s , such that
Ω w i ρ ( v ˙ i g i ) d Ω + Ω ( w i , j τ i j p w i , i ) d Ω Γ f w i f i Γ f d Γ + Ω s [ w i s ( ρ s ρ ) ( v ˙ i g i ) + w i , j s ( τ i j s τ i j f ) ( p s p ) w i , i s ] d Ω + + Ω q ( v j , j + p , t κ ) d Ω + Ω s q s ( J 3 1 + p s κ s ) d Ω = 0 .
For the fluid domain the following interpolations are introduced for the entire domain Ω :
v h = N I v v I , w h = N I v w I , p h = N I p p I , q h = N I p q I ,
where N I v and N I p 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:
u s , h = N J u u J s , w s , h = N J u w J s , p s , h = N J p s p J s , q s , h = N J p s q J s ,
where N J u and N J p s are the interpolation functions at node J for the solid displacement vector dependent of the background fluid velocity vector and the independent solid pressure unknowns; and u J s , h , w J s , h , p J s , h , and q J s , h stand for the nodal values of the discretized solid displacement vector, admissible solid displacement variation, solid pressure, admissible solid pressure variation, respectively.
Now if we substitute both discretizations (13) and (14) into Equation (12), we can have the following discretization of the weak form: q h L 2 ( Ω h ) , q s , h L 2 ( Ω s h ) , w h [ H 0 , Γ v h 1 , h ( Ω h ) ] d , which includes w s , h [ H 1 , h ( Ω s h ) ] d ,
Ω h w i I N I v ρ v ˙ i h d Ω Γ f h w i I N I v f i Γ f h d Γ + Ω h ( w i I N I , j v τ i j p h w i I N I , i v ) d Ω + Ω s h [ w i J s N J u ( ρ s ρ ) ( v ˙ i h g i ) + w i J s N J , i u ( σ i j s σ i j f ) ] d Ω Ω h w i I N I v ρ g i d Ω + Ω h q I N I p ( v j , j h + p , t h κ ) d Ω + Ω s h q J s N J p s ( J 3 1 + p s , h κ s ) d Ω = 0 .
Therefore, in the entire domain Ω , the arbitrariness of the fluid velocity variations w i I , fluid pressure variations q I , and independent solid pressure variations q J s 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,
r i I v = 0 , r I p = 0 , r J p s = 0 ,
where the residuals are defined as
r i I v = Ω h N I v ρ v ˙ i h d Ω + Ω h [ N I , j v τ i j p h N I , i v ] d Ω Γ f h N I v , Γ f h f i Γ f h d Γ + Ω s h N ˜ [ N J u ( ρ s ρ ) ( v ˙ i h g i ) + N J , i u ( σ i j s σ i j f ) ] d Ω Ω h N I v ρ g i d Ω , r I p = Ω h N I p ( v j , j h + p , t h κ ) d Ω , r J p s = Ω s h N J p s ( J 3 1 + p s , h κ s ) d Ω .
Please note that in both Equations (15) and (17) an implicit function N ˜ is introduced to map from w I to w J s , 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 N ˜ 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.

2.2. 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:
( K u u * ) h ( K u s ) h ( K u s ) h T ( K s s ) h U ¯ h S ¯ h = R h * 0 ,
where R h * is an effective load vector, ( K u u * ) h and ( K u s ) 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:
V ¯ h T ( K u u ) h * V ¯ h c 1 V ¯ h V 2 f o r   a l l V ¯ h ker [ ( K u s ) h T ] ,
where c 1 > 0 , v V 2 = i , j v i x j L 2 ( V o l ) 2 and k e r [ ( K u s ) h T ] = V ¯ h | V ¯ h R n , ( K u s ) h T V ¯ h = 0 .
(b) Inf-Sup Condition:
inf S ¯ h sup U ^ h U ¯ h T ( K u s ) h S ¯ h | | U ¯ h | | | | S ¯ h | | c 2 > 0 ,
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 u s are often employed.

2.3. 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
u R ( x ) = + u ( y ) ϕ ( x y ) d y ,
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 t h moment depends on the number of support points. For instance, the discretized delta function with four support points is uniquely defined by the following:
  • ϕ is a continuous function, with ϕ ( r ) = 0 for | r | 2 ;
  • For all r, j   e v e n ϕ ( r j ) = j   o d d ϕ ( r j ) = 1 2 ;
  • For all r, j ( r j ) ϕ ( r j ) = 0 ; and
  • 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 ,
ϕ ( r ) = ( 2 r 3 ) + 4 r 4 r 2 + 1 8 , ϕ ( r 2 ) = 1 2 ϕ ( r ) , ϕ ( r 1 ) = 1 4 + r 2 + ϕ ( r ) , ϕ ( r + 1 ) = 3 4 r 2 ϕ ( r ) .
Naturally, in the three-dimensional case, one of the smoothed approximations to the delta function is given by
δ h ( x ) = 1 h 3 ϕ ( x 1 h ) ϕ ( x 2 h ) ϕ ( x 3 h ) .
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].

3. 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 d h as c p d T 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 d u as c v d T 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 = ρ R T 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
d u = p d ( 1 / ρ ) + d q ,
where d q is the specific heat transferred to the system from the environment.
For isentropic (quasistatic, adiabatic, and reversible) gas, with d q = 0 , we have
c v d T = p d ( 1 / ρ ) ,
Thus, combine with the ideal gas law, p v = R T , we have
c v / R d ( p / ρ ) = p d ( 1 / ρ ) .
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 = d p / d ρ or γ p / ρ = γ R T . 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
( ρ e ) t + x i [ ( ρ e + p ) v i ] x i ( k T x i ) = x i ( τ i j v j ) + f i v i + q s ,
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.

4. Matrix-Free Newton–Krylov Iteration

4.1. 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
A ϕ j = A 1 , A 2 , , A ( n 1 ) , A n ϕ 1 j ϕ 2 j ϕ 3 j ϕ ( n 1 ) j ϕ n j = i = 1 n ϕ i j A i ,
where A i represents the ith column of the matrix A .
Denote B as a matrix with m rows and n columns
ψ 1 T ψ 2 T ψ 3 T ψ ( k 1 ) T ψ k T B = ψ 1 T B ψ 2 T B ψ ( k 1 ) T B ψ k T B
the multiplication ψ i T 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,
ψ i T B = ψ i 1 , ψ i 2 , , ψ i ( m 1 ) , ψ i m B 1 B 2 B ( m 1 ) B m = j = 1 m ψ i j B j ,
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
E A = I 1 F 0 1 0 2 .
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.

4.2. 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 N ˜ 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
p = r m + 1 , k 1 r , v f m + 1 , k 1 Δ V f , 0 r , u c m + 1 , k 1 Δ U c , 0 r , u s m + 1 , k 1 Δ U s , 0 r , p f m + 1 , k 1 Δ P f , 0 r , p c m + 1 , k 1 Δ P c , 0 r , p s m + 1 , k 1 Δ P s , 0 .
with
J = ( r , v f m + 1 , k 1 , r , u c m + 1 , k 1 , r , u s m + 1 , k 1 , r , p f m + 1 , k 1 , r , p c m + 1 , k 1 r , p s m + 1 , k 1 ) .
A set of orthonormal vectors v i , with 1 i n in the Krylov space K n and an ( n + 1 ) × n upper-Hessenberg matrix H ¯ n . If we define V n = ( v 1 , v 2 , , v n ) and V n + 1 = ( v 1 , v 2 , , v n + 1 ) ,
J V n = V n + 1 H ¯ n ,
and the remaining process in the GMRES method is to solve the least square problem:
min z K n p J z or min y R n p J V n y .
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,
min y R n γ e 1 H ¯ n y ,
γ = 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 = J q i , and for j = 1 to i , we have h j i = w T v j and w is updated with w h j i v j . Consequently, we obtain h ( i + 1 ) i = w 2 and v i + 1 = w / h ( i + 1 ) i . w = J q i with
J q i [ r ( Θ m + 1 , k 1 + e q i ) r ( Θ m + 1 , k 1 ) ] / e .
After we establish the elements of an upper n × n Hessenberg matrix H n as well as an upper ( n + 1 ) × n Hessenberg matrix H ¯ n , for j = 1 to n , and i = 1   to j 1 , a factorization of H n is carried out
h i j = c i h i j + s i h ( i + 1 ) j   and   h ( i + 1 ) j = s i h i j + c i h ( i + 1 ) j ,
where the entities of the rotation processes are calculated as
r = h j j 2 + h ( j + 1 ) j 2 , c j = h j j / r ,   and   s j = h ( j + 1 ) j / r .
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 j j = r , p j = c j b j ,   and   p j + 1 = s j b j .
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
Δ Θ k , n = Δ Θ k , 0 + i = 1 n y i q i ,
or
< Δ V f , Δ U c , Δ U s , Δ P f , Δ P c , Δ P s > = Δ Θ k , 0 + i = 1 n y i q i .
Notice here for time increments, we use the trapezoidal temporal discretization is as the following,
b ( t + Δ t ) = b ( t ) + Δ t 2 b ˙ ( t ) + b ˙ ( t + Δ t ) b ˙ ( t + Δ t ) = b ˙ ( t ) + Δ t 2 b ¨ ( t ) + b ¨ ( t + Δ t )
and consequently,
b ( t + Δ t ) = b ( t ) + b ˙ ( t ) Δ t + Δ t 2 4 b ¨ ( t ) + b ¨ ( t + Δ t )
and more importantly,
δ b ( t + Δ t ) = Δ t 2 δ b ˙ ( t + Δ t ) δ b ˙ ( t + Δ t ) = Δ t 2 δ b ¨ ( t + Δ t )
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
A = 1 1 0 0 0 0 0 1 1 0 0 0 0 1 0 1 0 0 0 0 1 1 0 0 0 0 1 0 1 0 0 0 0 0 1 1 0 0 0 1 0 1 ;
H = h 11 h 12 h 13 h 14 h 1 k h 21 h 22 h 23 h 24 h 2 k 0 h 32 h 33 h 34 h 3 k 0 0 h 43 h 44 h 4 k 0 0 0 0 0 0 0 h ( k + 1 ) k .
therefore, for k = 1
c 1 s 1 s 1 c 1 h 11 ( 1 ) h 21 ( 1 ) y 1 = c 1 s 1 s 1 c 1 p 1 ( 1 ) 0 .
h 11 ( 2 ) 0 y 1 = p 1 ( 2 ) p 2 ( 2 ) .
with p 1 ( 1 ) = β , p 2 ( 1 ) = 0 , p 1 ( 2 ) = c 1 β , p 2 ( 2 ) = s 1 β . Naturally, if p 2 ( 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
c 1 s 1 0 s 1 c 1 0 0 0 1 h 11 ( 1 ) h 12 ( 1 ) h 21 ( 1 ) h 22 ( 1 ) 0 h 32 ( 1 ) y 1 y 2 = c 1 s 1 0 s 1 c 1 0 0 0 1 p 1 ( 1 ) 0 0 .
1 0 0 0 c 2 s 2 0 s 2 c 2 h 11 ( 2 ) h 12 ( 2 ) 0 h 22 ( 2 ) 0 h 32 ( 2 ) y 1 y 2 = 1 0 0 0 c 2 s 2 0 s 2 c 2 p 1 ( 2 ) p 2 ( 2 ) 0 .
h 11 ( 2 ) h 12 ( 2 ) 0 h 22 ( 3 ) 0 0 y 1 y 2 = p 1 ( 2 ) p 2 ( 3 ) p 3 ( 3 ) .
with p 1 ( 1 ) = β , p 2 ( 1 ) = 0 , p 3 ( 1 ) = 0 ; p 1 ( 2 ) = c 1 β , p 2 ( 2 ) = s 1 β , p 3 ( 2 ) = 0 ; p 1 ( 3 ) = p 1 ( 2 ) , p 2 ( 3 ) = c 2 p 2 ( 2 ) , p 3 ( 3 ) = s 2 p 2 ( 2 ) . Naturally, if p 3 ( 3 ) or rather s 2 p 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
c 1 s 1 0 0 s 1 c 1 0 0 0 0 1 0 0 0 0 1 h 11 ( 1 ) h 12 ( 1 ) h 13 ( 1 ) h 21 ( 1 ) h 22 ( 1 ) h 23 ( 1 ) 0 h 32 ( 1 ) h 33 ( 1 ) 0 0 h 43 ( 1 ) y 1 y 2 y 3 = c 1 s 1 0 0 s 1 c 1 0 0 0 0 1 0 0 0 0 1 p 1 ( 1 ) 0 0 0 .
1 0 0 0 0 c 2 s 2 0 0 s 2 c 2 0 0 0 0 1 h 11 ( 2 ) h 12 ( 2 ) h 13 ( 2 ) 0 h 22 ( 2 ) h 23 ( 2 ) 0 h 32 ( 1 ) h 33 ( 1 ) 0 0 h 43 ( 1 ) y 1 y 2 y 3 = 1 0 0 0 0 c 2 s 2 0 0 s 2 c 2 0 0 0 0 1 p 1 ( 2 ) p 2 ( 2 ) 0 0 .
1 0 0 0 0 1 0 0 0 0 c 3 s 3 0 0 s 3 c 3 h 11 ( 2 ) h 12 ( 2 ) h 13 ( 2 ) 0 h 22 ( 3 ) h 23 ( 3 ) 0 0 h 33 ( 3 ) 0 0 h 43 ( 1 ) y 1 y 2 y 3 = 1 0 0 0 0 1 0 0 0 0 c 3 s 3 0 0 s 3 c 3 p 1 ( 2 ) p 2 ( 3 ) p 3 ( 3 ) 0 .
h 11 ( 2 ) h 12 ( 2 ) h 13 ( 2 ) 0 h 22 ( 2 ) h 23 ( 2 ) 0 0 h 33 ( 4 ) 0 0 0 y 1 y 2 y 3 = p 1 ( 2 ) p 2 ( 3 ) p 3 ( 4 ) p 4 ( 4 ) .
with p 1 ( 1 ) = β , p 2 ( 1 ) = 0 , p 3 ( 1 ) = 0 , p 4 ( 1 ) = 0 ; p 1 ( 2 ) = c 1 β , p 2 ( 2 ) = s 1 β , p 3 ( 2 ) = 0 , p 4 ( 2 ) = 0 ; p 1 ( 3 ) = p 1 ( 2 ) , p 2 ( 3 ) = c 2 p 2 ( 2 ) , p 3 ( 3 ) = s 2 p 2 ( 2 ) , p 4 ( 3 ) = 0 ; p 1 ( 4 ) = p 1 ( 2 ) , p 2 ( 4 ) = p 2 ( 3 ) , p 3 ( 4 ) = c 3 p 3 ( 3 ) , p 4 ( 4 ) = s 3 p 3 ( 3 ) . Naturally, if p 4 ( 4 ) or rather s 2 p 3 ( 3 ) is sufficiently small, the third equation will be approximately satisfied, and the iteration will stop at the level when k = 3 .

5. 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
x ( E f u x ) = ρ f 2 u t 2 ρ f g , x ( E s u x ) = ρ s 2 u t 2 ρ s g .
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
ρ f E f g x 2 2 + c 1 E f x + c 4 ; ρ f E f g x 2 2 + c 2 E f x + c 5 ; and ρ s E s g x 2 2 + c 3 E s x + c 6 .
Therefore, six constants c i will be derived based on the following six boundary conditions
u ( 0 ) = 0 , u ( L ) = 0 , u ( L 1 ) = u ( L 1 + ) , u ( L 2 + ) = u ( L 2 ) , E f u x ( L 1 ) = E s u x ( L 1 + ) , E f u x ( L 2 + ) = E s u x ( L 2 ) .
In this model study, for simplicity, we assign L 1 = L / 3 = l , L 2 = 2 L / 3 = 2 l , E s = a E f = a E , 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 β = ρ g l , the solution of the six unknown constants c can be determined from the boundary conditions (41) as
c = β 2 b + 2 , b + 4 , 3 b , 0 , 3 b 3 , 1 + b 2 b / a
The comparisons with analytical solutions as well as existing computational fluid and solid software package such as ADINA as illustrated in Figure 6 and Figure 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 π / 40 t )   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   k g / 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   P a · s , the bulk modulus κ f = 2.1 × 10 9   N / m 2 , and the density ρ f = 1000   k g / m 3 .
It is evident that even with the coarse mesh used, the developed formulation with high-order 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 / 4 c mixed finite elements as illustrated in detail in Refs. [4,64]. A typical immersed solid is represented by 4 × 4 9 / 4 c 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 Figure 9 and Figure 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 Figure 12 and Figure 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.

6. 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 N ˜ 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].

Funding

Not Applicable.

Institutional Review Board Statement

Not Applicable.

Informed Consent Statement

Not Applicable.

Data Availability Statement

Not Applicable.

Acknowledgments

The author is very grateful for ASEE Air Force Summer Faculty Fellowship in 2008 and 2009 and the opportunity of collaboration with Raymond Gordnier, Air Vehicle Directory, Air Force Research Laboratory, WPAFB, OH 45433.

Conflicts of Interest

The author declares no conflict of interest.

References

  1. Belytschko, T. Fluid-structure interaction. Comput. Struct. 1980, 12, 459–469. [Google Scholar] [CrossRef]
  2. Liu, W.; Belytschko, T.; Chang, H. An arbitrary Lagrangian-Eulerian finite element method for path-dependent materials. Comput. Methods Appl. Mech. Eng. 1986, 58, 227–245. [Google Scholar] [CrossRef]
  3. Shyy, W.; Udaykumar, H.S.; Rao, M.M.; Smith, R.W. Computational Fluid Dynamics with Moving Boundaries; Taylor & Francis: London, UK, 1996. [Google Scholar]
  4. Wang, X.S. Fundamentals of Fluid-Solid Interactions-Analytical and Computational Approaches; Elsevier Science: Amsterdam, The Netherlands, 2008. [Google Scholar]
  5. Bathe, K.J.; Zhang, H.; Wang, M.H. Finite element analysis of incompressible and compressible fluid flows with free surfaces and structural interactions. Comput. Struct. 1995, 56, 193–213. [Google Scholar] [CrossRef] [Green Version]
  6. Bathe, K.J.; Zhang, H. Finite element developments for general fluid flows with structural interactions. Int. J. Numer. Methods Eng. 2004, 60, 213–232. [Google Scholar] [CrossRef] [Green Version]
  7. Tornberg, A.; Shelley, M. Simulating the dynamics and interactions of flexible fibers in stokes flows. J. Comput. Phys. 2004, 196, 8–40. [Google Scholar] [CrossRef]
  8. Zhang, J.; Childress, S.; Libchaber, A.; Shelley, M. Flexible filaments in a flowing soap film as a model for one-dimensional flags in a two-dimensional wind. Nature 2000, 408, 835–838. [Google Scholar] [CrossRef]
  9. Huerta, A.; Liu, W.K. Viscous flow with large free surface motion. Comput. Methods Appl. Mech. Eng. 1988, 69, 277–324. [Google Scholar] [CrossRef] [Green Version]
  10. Tezduyar, T. Stabilized finite element formulations for incompressible-flow computations. Adv. Appl. Mech. 1992, 28, 1–44. [Google Scholar]
  11. Tezduyar, T. Finite element methods for flow problems with moving boundaries and interfaces. Arch. Comput. Methods Eng. 2001, 8, 83–130. [Google Scholar] [CrossRef]
  12. Wall, W.A.; Schott, B.; Ager, C. A monolithic approach to fluid-structure interaction based on a hybrid eulerian-ALE fluid domain decomposition involving cut elements. Int. J. Numer. Methods Eng. 2019, 119, 208–237. [Google Scholar]
  13. Zhang, L.; Wagner, G.J.; Liu, W.K. Modeling and simulation of fluid structure interaction by meshfree and FEM. Commun. Numer. Methods Eng. 2003, 19, 615–621. [Google Scholar] [CrossRef]
  14. Oshe, S.J.; Fedki, R.P. Level Set Methods and Dynamic Implicit Surfaces; Springer: Berlin/Heidelberg, Germany, 2002. [Google Scholar]
  15. Sethian, J.A. Level Set Methods and Fast Marching Methods; Cambridge University Press: Cambridge, UK, 1996. [Google Scholar]
  16. LeVeque, R.; Li, Z. Immersed interface methods for Stokes flow with elastic boundaries or surface tension. SIAM J. Sci. Comput. 1997, 18, 709–735. [Google Scholar] [CrossRef] [Green Version]
  17. Li, Z.; LeVeque, R. The immersed interface methods for elliptic equations with discontinuous coefficients and singular sources. SIAM J. Sci. Comput. 1994, 31, 1019–1044. [Google Scholar]
  18. Glimm, J.; Grove, J.W.; Li, X.L.; Tan, D.C. Robust computational algorithms for dynamic interface tracking in three dimensions. SIAM J. Sci. Comput. 2000, 21, 2240–2256. [Google Scholar] [CrossRef]
  19. Li, X.L.; Jin, B.X.; Glimm, J. Numerical study for the three-dimensional Rayleigh-Taylor instability through the TVD/AC scheme and parallel computation. J. Comput. Phys. 1996, 126, 343–355. [Google Scholar] [CrossRef] [Green Version]
  20. Tryggvason, G. Numerical simulations of the Rayleigh-Taylor instability. J. Comput. Phys. 1988, 75, 253–282. [Google Scholar] [CrossRef] [Green Version]
  21. Farhat, C.; Geuzaine, P.; Grandmont, C. The discrete geometric conservation law and the nonlinear stability of ALE schemes for the solution of flow problems on moving grids. J. Comput. Phys. 2001, 174, 669–694. [Google Scholar] [CrossRef]
  22. Hughes, T.; Liu, W.K.; Zimmerman, T.K. Lagrangian-Eulerian finite element formulations for incompressible viscous flows. Comput. Methods Appl. Mech. Eng. 1981, 29, 329–349. [Google Scholar] [CrossRef]
  23. Kamakoti, R.; Shyy, W. Fluid-structure interaction for aeroelastic applications. Prog. Aerosp. Sci. 2004, 40, 535–558. [Google Scholar] [CrossRef]
  24. Lian, Y.; Shyy, W. Numerical simulations of membrane wing aerodynamics for micro air vehicle applications. J. Aircr. 2005, 42, 865–873. [Google Scholar] [CrossRef]
  25. Shyy, W.; Berg, M.; Ljungqvist, D. Flapping and flexible wings for biological and micro air vehicles. Prog. Aerosp. Sci. 1999, 35, 455–505. [Google Scholar] [CrossRef]
  26. Chessa, J.; Belytschko, T. The extended finite element method for two-phase fluids. J. Appl. Mech. 2003, 70, 10–17. [Google Scholar] [CrossRef]
  27. Wagner, G.; Moës, N.; Liu, W.K.; Belytschko, T. The extended finite element method for rigid particles in stokes flow. Int. J. Numer. Methods Eng. 2001, 51, 293–313. [Google Scholar] [CrossRef]
  28. Glowinski, R.; Pan, T.; Hesla, T.I.; Joseph, D.D.; Périaux, J. A fictitious domain approach to the direct numerical simulation of incompressible viscous flow past moving rigid bodies: Application to particulate flow. J. Comput. Phys. 2001, 169, 363–426. [Google Scholar] [CrossRef] [Green Version]
  29. Wang, X.; Yang, Y.; Wu, T. Model studies of fluid-structure interaction problems. Comput. Model. Eng. Sci. 2019, 119, 5–34. [Google Scholar]
  30. McQueen, D.; Peskin, C. Computer-assisted design of butterfly bileaflet valves for the mitral position. Scand. J. Thor. Cardiovasc. Surg. 1985, 19, 139–148. [Google Scholar] [CrossRef] [PubMed]
  31. Peskin, C. Numerical analysis of blood flow in the heart. J. Comput. Phys. 1977, 25, 220–252. [Google Scholar] [CrossRef]
  32. Peskin, C. The immersed boundary method. Acta Numer. 2002, 11, 479–517. [Google Scholar] [CrossRef] [Green Version]
  33. Beyer, R. A computational model of the cochlea using the immersed boundary method. J. Comput. Phys. 1992, 98, 145–162. [Google Scholar] [CrossRef]
  34. Fauci, L.; Peskin, C. A computational model of aquatic animal locomotion. J. Comput. Phys. 1988, 77, 85–108. [Google Scholar] [CrossRef]
  35. Stockie, J.M.; Wetton, B.R. Analysis of stiffness in the immersed boundary method and implications for time-stepping schemes. J. Comput. Phys. 1999, 154, 41–64. [Google Scholar] [CrossRef] [Green Version]
  36. Dillon, R.; Fauci, L.; Fogelson, A.; Gaver, D., III. Modeling biofilm processes using the immersed boundary method. J. Comput. Phys. 1996, 129, 57–73. [Google Scholar] [CrossRef]
  37. van Loon, R.; Anderson, P.D.; van de Vosse, F.N.; Sherwin, S.J. Comparison of various fluid-structure interaction methods for deformable bodies. Comput. Struct. 2007, 85, 833–843. [Google Scholar] [CrossRef] [Green Version]
  38. Yu, Z. A DLM/FD method for fluid/flexible-body interactions. J. Comput. Phys. 2005, 207, 1–27. [Google Scholar] [CrossRef]
  39. Wang, X.; Zhang, L.T.; Liu, W.K. On computational issues of immersed finite element methods. J. Comput. Phys. 2009, 228, 2535–2551. [Google Scholar] [CrossRef]
  40. Sulsky, D.; Kaul, A. Implicit dynamics in the material-point method. Comput. Methods Appl. Mech. Eng. 2004, 193, 1137–1170. [Google Scholar] [CrossRef]
  41. Lim, S.; Ferent, A.; Wang, X.; Peskin, C. Dynamics of a closed rod with twist and bend in fluid. SIAM J. Sci. Comput. 2008, 31, 273–302. [Google Scholar] [CrossRef]
  42. Gay, M.; Zhang, L.T.; Liu, W.K. Stent modeling using immersed finite element method. Comput. Methods Appl. Mech. Eng. 2005, 195, 4358–4370. [Google Scholar] [CrossRef]
  43. Liu, W.; Liu, Y.; Farrell, D.; Zhang, L.; Wang, X.; Fukui, Y.; Patankar, N.; Zhang, Y.; Bajaj, C.; Chen, X.; et al. Immersed finite element method and its applications to biological systems. Comput. Methods Appl. Mech. Eng. 2006, 195, 1722–1749. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  44. Liu, W.; Jun, S.; Zhang, Y.F. Reproducing kernel particle methods. Int. J. Numer. Methods Fluids 1995, 20, 1081–1106. [Google Scholar] [CrossRef]
  45. Boffi, D.; Gastaldi, L. A finite element approach for the immersed boundary method. Comput. Struct. 2003, 81, 491–501. [Google Scholar] [CrossRef]
  46. Mori, Y.; Peskin, C. Implicit second-order immersed boundary methods with boundary mass. Comput. Methods Appl. Mech. Eng. 2008, 197, 2049–2067. [Google Scholar] [CrossRef]
  47. Newren, E.; Fogelson, A.L.; Guy, R.D.; Kirby, R.M. Unconditionally stable discretizations of the immersed boundary equations. J. Comput. Phys. 2007, 222, 702–719. [Google Scholar] [CrossRef]
  48. Liu, Y.; Liu, W.K.; Belytschko, T.; Patankar, N.; To, A.C.; Kopacz, A.; Chung, J.H. Immersed electrokinetic finite element method. Int. J. Numer. Methods Eng. 2007, 71, 379–405. [Google Scholar] [CrossRef]
  49. Liu, Y.; Zhang, L.; Wang, X.; Liu, W.K. Coupling of Navier-Stokes equations with protein molecular dynamics and its application to hemodynamics. Int. J. Numer. Methods Fluids 2004, 46, 1237–1252. [Google Scholar] [CrossRef]
  50. Fogelson, A.; Wang, X.; Liu, W.K. Special issue: Immersed boundary method and its extensions—Preface. Comput. Methods Appl. Mech. Eng. 2008, 197, 25–28. [Google Scholar]
  51. Liu, W.; Kim, D.W.; Tang, S. Mathematical foundations of the immersed finite element method. Comput. Mech. 2007, 39, 211–222. [Google Scholar] [CrossRef]
  52. Mittal, R.; Iaccarino, G. Immersed boundary methods. Annu. Rev. Fluid Mech. 2005, 37, 239–261. [Google Scholar] [CrossRef] [Green Version]
  53. Wang, X. From immersed boundary method to immersed continuum method. Int. J. Multiscale Comput. Eng. 2006, 4, 127–145. [Google Scholar] [CrossRef]
  54. Wang, X. An iterative matrix-free method in implicit immersed boundary/continuum methods. Comput. Struct. 2007, 85, 739–748. [Google Scholar] [CrossRef]
  55. Sulsky, D.; Brackbill, J.U. A numerical method for suspension flow. J. Comput. Phys. 1991, 96, 339–368. [Google Scholar] [CrossRef]
  56. Boffi, D.; Gastaldi, L.; Heltai, L. On the CFL condition for the finite element immersed boundary method. Comput. Struct. 2007, 85, 775–783. [Google Scholar] [CrossRef]
  57. Boffi, D.; Gastaldi, L.; Heltai, L.; Peskin, C. On the hyperelastic formulation of the immersed boundary method. Comput. Methods Appl. Mech. Eng. 2008, 197, 2210–2231. [Google Scholar] [CrossRef]
  58. Wang, X.; Liu, W.K. Extended immersed boundary method using FEM and RKPM. Comput. Methods Appl. Mech. Eng. 2004, 193, 1305–1321. [Google Scholar] [CrossRef]
  59. Zhang, L.; Gerstenberger, A.; Wang, X.; Liu, W.K. Immersed finite element method. Comput. Methods Appl. Mech. Eng. 2004, 193, 2051–2067. [Google Scholar] [CrossRef]
  60. Shu, C.W. Essential non-oscillatory and weighted essentially non-oscillatory schemes. Acta Numer. 2020, 1–60. [Google Scholar]
  61. Wang, X.; Gordnier, R. A combination of immersed methods and compact scheme based computational fluid dynamics for fluid-structure interactions. In Air Force ASEE Summer Faculty Fellow Report; Midwestern State University: Wichita Falls, TX, USA, 2009. [Google Scholar]
  62. Wang, X.; Bathe, K.J. Displacement/pressure based finite element formulations for acoustic fluid-structure interaction problems. Int. J. Numer. Methods Eng. 1997, 40, 2001–2017. [Google Scholar] [CrossRef]
  63. Wang, X.; Bathe, K.J. On mixed elements for acoustic fluid-structure interactions. Math. Model. Methods Appl. Sci. 1997, 7, 329–343. [Google Scholar] [CrossRef]
  64. Bao, W.; Wang, X.; Bathe, K.J. On the Inf-Sup condition of mixed finite element formulations for acoustic fluids. Math. Model. Methods Appl. Sci. 2001, 11, 883–901. [Google Scholar] [CrossRef] [Green Version]
  65. Bathe, K. Finite Element Procedures; Prentice Hall: Englewood Cliffs, NJ, USA, 1996. [Google Scholar]
  66. Li, S.; Liu, W.K. Reproducing kernel hierarchical partition of unity part I: Formulation and theory. Int. J. Numer. Methods Eng. 1999, 45, 251–288. [Google Scholar] [CrossRef]
  67. Liu, W.; Chen, Y.; Chang, C.T.; Belytschko, T. Advances in multiple scale kernel particle methods. Comput. Mech. 1996, 18, 73–111. [Google Scholar] [CrossRef]
  68. Liu, W.; Chen, Y.; Uras, R.A.; Chang, C.T. Generalized multiple scale reproducing kernel particle methods. Comput. Methods Appl. Mech. Eng. 1996, 139, 91–157. [Google Scholar] [CrossRef]
  69. Liu, W.; Chen, Y.J. Wavelet and multiple scale reproducing kernel methods. Int. J. Numer. Methods Fluids 1995, 21, 901–932. [Google Scholar] [CrossRef]
  70. Wang, X. On the discretized delta function and force calculation in the immersed boundary method. In Computational Fluid and Solid Mechanics; Bathe, K., Ed.; Elsevier: Amsterdam, The Netherlands, 2003; pp. 2164–2169. [Google Scholar]
  71. Griffith, B.; Patankar, N.A. Immersed methods for fluid-structure interaction. Annu. Rev. Fluid Mech. 2020, 52, 421–448. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  72. Wang, X.; Zhang, L. Immersed Methods for Coupled Systems of Continua; World Scientific Publishing Co.: Singapore, 2021. [Google Scholar]
  73. Pan, T. Private Communications. 2021. Available online: http://www.math.uh.edu/ (accessed on 4 May 2021).
  74. Vincent, S.; Caltagirone, J.P. A monolithic approach of fluid–structure interaction by discrete mechanics. Fluids 2021, 6, 95. [Google Scholar] [CrossRef]
Figure 1. The concept of fictitious fluid domain vs. immersed domain.
Figure 1. The concept of fictitious fluid domain vs. immersed domain.
Fluids 06 00273 g001
Figure 2. An illustration of a Lagrangian mesh on top of a Eulerian mesh for immersed methods.
Figure 2. An illustration of a Lagrangian mesh on top of a Eulerian mesh for immersed methods.
Fluids 06 00273 g002
Figure 3. Pressure distribution of an acoustic fluid-solid interaction system (from the author’s earlier work [29,62,63]).
Figure 3. Pressure distribution of an acoustic fluid-solid interaction system (from the author’s earlier work [29,62,63]).
Fluids 06 00273 g003
Figure 4. A comparison of various discretized delta functions and their spectra (from the author’s earlier work [58,70]).
Figure 4. A comparison of various discretized delta functions and their spectra (from the author’s earlier work [58,70]).
Fluids 06 00273 g004
Figure 5. One-dimensional model problem of one continuum immersed in another continuum.
Figure 5. One-dimensional model problem of one continuum immersed in another continuum.
Fluids 06 00273 g005
Figure 6. Analytical solutions.
Figure 6. Analytical solutions.
Fluids 06 00273 g006
Figure 7. Comparisons with the monolithic implicit immersed method with matrix-free Newton–Krylov iterations and ADINA.
Figure 7. Comparisons with the monolithic implicit immersed method with matrix-free Newton–Krylov iterations and ADINA.
Fluids 06 00273 g007
Figure 8. Snapshots of results from the immersed continuum method and the ADINA FSI solver.
Figure 8. Snapshots of results from the immersed continuum method and the ADINA FSI solver.
Fluids 06 00273 g008
Figure 9. A comparison of horizontal and vertical velocities derived from the traditional ALE formulation and the immersed methods.
Figure 9. A comparison of horizontal and vertical velocities derived from the traditional ALE formulation and the immersed methods.
Fluids 06 00273 g009
Figure 10. Extra numerical dissipation due to immersed methods.
Figure 10. Extra numerical dissipation due to immersed methods.
Fluids 06 00273 g010
Figure 11. Driven cavity with five compressible immersed solids.
Figure 11. Driven cavity with five compressible immersed solids.
Fluids 06 00273 g011
Figure 12. Simultaneous convergence of the monolithic implicit immersed method with matrix-free Newton– Krylov iteration.
Figure 12. Simultaneous convergence of the monolithic implicit immersed method with matrix-free Newton– Krylov iteration.
Fluids 06 00273 g012
Figure 13. Fully converged velocity and pressure contours and bandplots.
Figure 13. Fully converged velocity and pressure contours and bandplots.
Fluids 06 00273 g013
Figure 14. Shockwave propagating within compressible air with compact schemes (from the author’s earlier work Ref. [61]).
Figure 14. Shockwave propagating within compressible air with compact schemes (from the author’s earlier work Ref. [61]).
Fluids 06 00273 g014
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Wang, S. A Revisit of Implicit Monolithic Algorithms for Compressible Solids Immersed Inside a Compressible Liquid. Fluids 2021, 6, 273. https://0-doi-org.brum.beds.ac.uk/10.3390/fluids6080273

AMA Style

Wang S. A Revisit of Implicit Monolithic Algorithms for Compressible Solids Immersed Inside a Compressible Liquid. Fluids. 2021; 6(8):273. https://0-doi-org.brum.beds.ac.uk/10.3390/fluids6080273

Chicago/Turabian Style

Wang, Sheldon. 2021. "A Revisit of Implicit Monolithic Algorithms for Compressible Solids Immersed Inside a Compressible Liquid" Fluids 6, no. 8: 273. https://0-doi-org.brum.beds.ac.uk/10.3390/fluids6080273

Article Metrics

Back to TopTop