A sharp interface Lagrangian-Eulerian method for rigid-body fluid-structure interaction
Ebrahim M. Kolahdouz Amneet P. S. Bhalla Lawrence N. Scotten Brent A. Craven Boyce E. Griffith
math.NA cs.NA


Corella: A Private Multi Server Learning Approach
 based on Correlated Queries
    Hamidreza Ehteram,
	Mohammad Ali Maddah-Ali,
	and Mahtab Mirmohseni
	Hamidreza Ehteram and Mahtab Mirmohseni are with the Department of Electrical Engineering, Sharif University of Technology, Tehran, Iran (emails: ehteram.hamidreza@ee.sharif.edu, mirmohseni@sharif.edu).
	Mohammad Ali Maddah-Ali is with Nokia Bell Labs, Holmdel, NJ, USA (email: mohammad.maddahali@nokia-bell-labs.com).

    December 30, 2023
========================================================================================================================================================================================================================================================================================================================================================================================================




This paper introduces a sharp interface method to simulate fluid-structure interaction (FSI)
involving rigid bodies immersed in viscous incompressible fluids. The capabilities of this methodology 
are benchmarked using a range of test cases and demonstrated using large-scale models of biomedical FSI.
The numerical approach developed herein, which we refer to as an immersed Lagrangian-Eulerian (ILE) method, integrates aspects of
partitioned and immersed FSI formulations by solving separate momentum equations for the fluid and solid subdomains, 
as in a partitioned formulation, while also 
using non-conforming discretizations of the dynamic fluid and structure regions, as in an immersed formulation.
A simple Dirichlet-Neumann coupling scheme is used, in which the motion of the immersed solid is driven by fluid traction forces evaluated along 
the fluid-structure interface, and the motion of the fluid along that interface is constrained to match the solid velocity and thereby satisfy the no-slip condition.
 To develop a practical numerical method, we adopt a penalty approach that approximately imposes the no-slip condition along the fluid-structure interface.
 In the coupling strategy, a separate discretization of the fluid-structure interface is tethered to the volumetric 
 solid mesh via stiff spring-like penalty forces. 
 Our fluid-structure coupling scheme relies on an immersed interface method (IIM) for discrete geometries, which enables the accurate determination of both velocities and stresses 
 along complex internal interfaces.
Numerical methods for FSI can suffer from instabilities related to the added mass effect, 
but the computational tests indicate that the methodology introduced here remains stable for selected test cases 
across a range of solid-fluid density ratios, including extremely small, nearly equal, equal, and large density ratios.
Biomedical FSI demonstration cases include results obtained using this method
to simulate the dynamics of a bileaflet mechanical heart valve in a pulse duplicator, and to model transport of
blood clots in a patient-averaged anatomical model of the inferior vena cava.



Keywords:	Fluid-structure interaction, immersed methods, immersed interface method, low density ratios, rigid body dynamics, mechanical heart valve, inferior vena cava, clot transport








§ INTRODUCTION




Predictive numerical models of fluid-structure interaction (FSI) have long been of major interest in the scientific computing community. 
Numerical simulations of FSI problems can be characterized by the solution approach taken for the coupled system of momentum equations associated with the fluid and structure.
Partitioned formulations of FSI 
describe a fluid-structure system using distinct, non-overlapping fluid and structure regions.
Commonly used numerical methods of this type include well-known arbitrary Lagrangian-Eulerian (ALE) schemes <cit.>.
These formulations can yield outstanding resolution of flows and stresses up to the fluid-structure interface. 
Despite their high accuracy, however, ALE methods for FSI can be constrained by the difficulties of body conforming grid regeneration and mesh morphing, 
which can make it challenging to use these approaches to simulate systems involving very large structural displacements or deformations, and 
to handle transient contact between moving structures.
Overset methods using overlapping Chimera grid systems are another class of partitioned methods that have been used
to simulate moving rigid and flexible bodies <cit.>.
In these approaches, a complex domain is decomposed into multiple geometrically simple overlapping grids, and boundary information is exchanged between
these grids through interpolation.
Immersed formulations of FSI <cit.> 
are alternatives to body-fitted methods. 
Many immersed approaches to FSI have been developed, including Peskin's immersed boundary (IB) method <cit.> and various sharp-interface 
IB methods <cit.>. These
methods commonly describe the fluid in Eulerian form (i.e. using fixed physical coordinates) 
and the structure in Lagrangian form (i.e. using material coordinates attached to the structure), 
and they use non-conforming
discretizations along the fluid-structure interface.  
Because these methods avoid using body-conforming discretizations of the interface, they are readily able to treat models with very large structural deformations or displacements, 
and they facilitate simulations that fundamentally involve contact or near-contact between structures <cit.>.
The method presented herein, which we call an immersed Lagrangian-Eulerian (ILE) method, combines a partitioned approach to FSI
with an immersed coupling strategy.

The key challenge in developing immersed methods for FSI is linking the Eulerian and Lagrangian variables. 
Peskin's IB methods, for example, regularize singular forces and stress discontinuities along the 
fluid-structure interface, which enables straightforward discretization approaches but can yield low accuracy 
in the flows and stresses near those interfaces. Efforts have been made to improve the accuracy of the method, including 
the development of formally second-order IB methods
that realize second-order accuracy when applied to specific problems <cit.>
and IB methods that use Cartesian grid adaptive mesh refinement to enhance spatial resolution near fluid-structure interfaces <cit.>.
For general FSI problems, however, formally second-order accurate IB methods 
 still only realize first-order convergence rates <cit.>. 
Motivated by the goal of improving the accuracy of the original IB method, 
the immersed interface method (IIM) was introduced by LeVeque and Li <cit.> for  
elliptic PDEs with discontinuous coefficients and singular forces. 
The IIM subsequently was extended to the 
incompressible Stokes <cit.> and Navier-Stokes 
<cit.> equations, 
and it was combined with level set methods to represent 
the interface <cit.>. 
When applied to the incompressible Navier-Stokes equations, the IIM sharply imposes interfacial stress discontinuities 
through an extended finite difference discretization that accounts for jump conditions induced 
by singular forces at the interface.
Modern IIMs use generalized Taylor series expansions to extend the physical jump conditions from the interface to the finite difference
discretization in the Eulerian domain while permitting the use of efficient linear solvers based on the unmodified finite difference
discretizations <cit.>.
The IIM has been used to simulate various phenomena, including acoustics and elastodynamic wave propagation <cit.>,
fluid interfaces with insoluble surfactants <cit.>, the osmotic swelling of a deforming 
capsule <cit.>, and vesicle electrohydrodynamics <cit.>.
Other sharp interface FSI methods have been developed 
<cit.>,
and most of these methods achieve higher-order accuracy by adopting approaches that are similar
 to body-fitted discretization methods, such as local modifications to the finite difference stencils,
 to allow for the accurate reconstruction of boundary conditions in the vicinity of the immersed interface.


Many different IB methods have been developed to treat FSI problems involving rigid bodies. 
Differences between these approaches are 
mainly related to the way that the rigidity constraint is enforced to account for the effect of the structure in the fluid region.
Previous IB approaches to rigid-body FSI 
that use regularized coupling operators include Lagrange-multiplier-based fictitious-domain methods <cit.>, 
direct forcing IB methods <cit.>,
projection-based methods <cit.>,
immersed finite element methods <cit.>,
methods based on computing exact Lagrange multipliers for the rigidity constraint <cit.>,
penalty immersed boundary methods <cit.>,
immersed boundary lattice Boltzmann methods <cit.>,
and level set based approaches <cit.>.
Sharp-interface approaches designed for rigid body FSI include embedded boundary methods <cit.>, cut-cell methods <cit.>, and
the curvilinear immersed boundary method <cit.>.
Unlike the IB approach introduced by Peskin, these methods all solve the fluid equations in the exterior regions
surrounding the immersed object.
Of the sharp-interface IB methods considered here, immersed interface methods are the most similar to Peskin's IB 
method because they treat thin interfaces that are fully immersed in fluid, which enables the use of fast Cartesian 
grid fluid solvers. IIMs have been developed for bodies with prescribed motion 
<cit.>,
but they are more commonly used for thin flexible interfaces <cit.>.
To our knowledge, the few IIM models that treat volumetric rigid structures 
do not simulate fluid-structure interaction per se, but instead prescribe the motion of the immersed body <cit.>.
Xu and Wang used a feedback control to construct the singular force density <cit.>.
Recent work by some of us introduced an IIM for discrete surfaces described by a general finite element mesh 
to sharply resolve fluid dynamics for problems with prescribed motion <cit.>.
The present study uses this method to develop a new sharp interface approach to rigid-body FSI.


The present ILE method introduces a partitioned approach to FSI with an immersed coupling strategy that sharply resolves flow features up to the fluid-structure interface. 
Like partitioned formulations, the present approach uses distinct momentum equations for the fluid and solid regions. However, 
like immersed methods, and unlike typical partitioned methods, our ILE approach uses a non-conforming discretization of the dynamic fluid-structure interface that
does not require any grid regeneration or mesh morphing to treat large structural motions.
The fluid and solid subproblems are solved in a partitioned manner using independent, non-conforming discretizations and are coupled only through interface conditions. 
The ILE equations of motion are first introduced using an exactly constrained formulation that exactly imposes kinematic interface conditions through 
a Lagrange multiplier force distribution applied along the fluid-structure interface.
Solving the exactly constrained equations would require 
the solution of an extended saddle-point system involving an exact Lagrange multiplier force along with
the Eulerian velocity and pressure fields.
Developing efficient linear solvers for such systems is challenging even for conventional 
IB formulations with regularized delta functions <cit.>. Consequently, to obtain a practical numerical method, 
we next reformulate this scheme using a penalty approach that relaxes the kinematic constraint, and we use this penalty ILE method in all of our numerical tests. Specifically, our penalty formulation uses two representations of the fluid-structure interface, including 
a surface mesh and the boundary of a bulk volumetric structural mesh, that are connected by forces that impose kinematic and dynamic interface conditions. 
The dynamics of the 
volumetric structural mesh are driven by the exterior fluid traction obtained from solving the equations of fluid dynamics.
The surface mesh moves according to the local fluid velocity and locally exerts an approximate Lagrange multiplier force distribution back 
to the fluid generated from stiff spring-like penalty forces that link the surface mesh to the boundary of its volumetric counterpart.
At least formally, in the limit of infinite spring stiffness, the two interface representations become exactly conformal in their motion.
Results demonstrate that this approach is able to control these discrepancies effectively, 
and for sufficiently large penalty spring stiffnesses, the penalty formulation has little impact on the computed dynamics.
To discretize the jump conditions, we leverage our recently developed IIM
for discrete surfaces <cit.>, which allows us to impose stress jump conditions along 
complex interfaces within a Cartesian grid framework
and to use fast structured-grid solvers for the incompressible Navier-Stokes equations.
This IIM formulation describes fluid dynamics on both sides of the fluid-structure interface. However, only the fluid forces exerted by 
the “exterior” fluid have a physical effect on the structural dynamics, and the motion of the structure determines the fluid velocity at the 
fluid-structure interface. 

Our approach also can be viewed as a sharp implementation of the distributed Lagrange multiplier (DLM) technique first introduced by Glowinski, Patankar, and coworkers for immersed rigid structures <cit.>. 
In the DLM approach to FSI, a Lagrange multiplier force field is introduced to 
impose the kinematic condition at the fluid-structure interface. In the DLM literature, the Lagrange multiplier field (either exact or approximate) 
has typically been smoothed using, e.g., via regularized delta functions like those used in Peskin's IB method. 
From this standpoint, our approach is different from other sharp interface immersed methods, 
in which the velocity matching condition is imposed directly, 
for example, through velocity reconstruction <cit.> or through cut-cell approaches <cit.>. 
These approaches forgo Lagrange multipliers entirely and instead solve the fluid momentum equations only within the fluid subdomain. 
In contrast, immersed approaches, including DLM methods, typically extend the fluid domain into the solid domain, so that the fluid momentum equation 
is solved on the entire computational domain. 
The fluid velocity 
determined within the solid region can be different from the actual solid velocity, although volumetric DLM formulations 
can impose the constraint that the fictitious fluid velocity matches
the solid velocity within the overlapping region occupied by the structure.
Ultimately, however, it is necessary only for the fluid and solid velocities to match along the fluid-structure interface. 
To our knowledge, the present method is the first DLM-type formulation to sharply impose this constraint.
In addition, the current scheme uses a constant-coefficient flow solver and yet can readily treat immersed bodies that are 
lighter than the ambient fluid because it only imposes the kinematic condition along the fluid-solid interface. This is in contrast 
to volumetric DLM schemes, such as the approach of Nangia et al. <cit.>, which approximately constrains 
fluid and solid velocities to match in the extended fluid domain and allows for solid-fluid mass density ratios less than one only 
through the use of a variable-density fluid solver <cit.>. 
 
 


Instabilities due to artificial added mass effect have been 
observed in weakly coupled FSI schemes in which the fluid and solid equations are linked via explicit time stepping schemes.
Such instabilities can occur if the mass density of the solid  is comparable to or less than the mass density of the fluid .
Added mass effect instabilities have been discussed in both
sharp-interface IB-type methods for FSI  <cit.> 
as well as body-fitted methods, including ALE methods for FSI <cit.>.
Various approaches, including subcycling or using modified coupling conditions, have been developed to
maintain stability <cit.>.
Strong coupling schemes, in which the governing equations for the fluid and solid subdomains are simultaneously 
integrated in time <cit.>, have been shown to improve the stability of these FSI formulations.
In one common strong coupling approach, solutions are transferred between the fluid and structure 
multiple times within each time step (i.e. through subiterations) until convergence is achieved in the forces and displacements 
<cit.>.
Using subiterations substantially increases the computational cost per time step, however. Further, 
instabilities at low density ratios have still been reported in some situations even if using strong coupling <cit.>.
Substantial work has been devoted to understanding the sources of these instabilities and to developing methods to overcome these instabilities 
<cit.>. 
We computationally examine the performance of the ILE method across a broad range of density ratios and, at least 
for the examples considered herein, we do not detect added mass-related instabilities.
Specifically, as described in the following, we are able to use a simple Dirichlet-Neumann
coupling scheme  <cit.> while achieving reasonable time step sizes and 
avoiding the use of subcycling <cit.> or other iterative techniques 
for models involving extremely small, nearly equal, equal, and large mass density ratios.


To assess the robustness and accuracy of the proposed
algorithm, results obtained using our ILE method are compared to experimental and computational data 
from benchmarks widely used to test numerical methods for FSI <cit.>. 
We consider test cases in two and three spatial dimensions involving both smooth and sharp geometries, in various fluid conditions with Reynolds numbering 
reaching up to =1147, 
and with different numbers of translational and rotational
degrees of freedom.
Finally, cases demonstrating the application of this methodology to biomedical models are presented, including
to simulating the dynamics of a bileaflet mechanical heart valve in a pulse duplicator 
with a peak Reynolds number of 22600, and to modeling the transport of blood clots in a patient-averaged anatomical model of the inferior vena cava
with a peak Reynolds number of 1500.








§ CONTINUOUS EQUATIONS OF MOTION



This section outlines our ILE approach to fluid-structure interaction. 
Our methodology builds on a conventional partitioned formulation for FSI, detailed in Sec. <ref>, while 
leveraging a coupling scheme based on the immersed interface method.




 §.§ Partitioned formulation of FSI


Typical partitioned formulations describe the fluid-structure system occupying a computational domain Ω
via moving fluid and structure subdomains, respectively Ω^f_t and Ω^s_t
and indexed by time t, so that Ω=_t ∪_t; see Fig. <ref>.
The regions meet along the fluid-structure interface, Γ^_t = Ω^f_t∩Ω^s_t.
Fixed physical coordinates are ∈Ω.
We describe the structural kinematics in Lagrangian form via reference coordinates ∈Ω^s_0 attached to the solid, and we use the motion
map : (Ω^s_0 , t)↦Ω^s_t to determine the physical position of 
a solid material point  at time t.
In the absence of additional loading terms, the equations of motion are
   
    D/Dt(,t)     = ·(,t),     ∈Ω^f_t,
    ·(̆,t)    = 0,     ∈Ω^f_t,
    ∂/∂ t(,t)    = (̆(,t),t),     ∈Γ^_0,
    d/dt∫_Ω^s_0 ∂/∂ t(,t)   d   =∫_Γ^_t(,t)  d a,
    d/dt∫_Ω^s_0( - _) ×( ∂/∂ t(,t)) d   = ∫_Γ^_t( - _(t)) ×(,t)   d a,
 
in which (̆,t) is the fluid velocity,  is the mass density of the fluid,  is the mass density 
of the structure, 
_(t) and _ are, respectively, the center of mass of the structure in the current and reference configurations,
(,t) is the fluid stress tensor,

    (,t) = -p(,t)   + ( (̆,t) + ^̆T(,t) ),   ∈Ω^f_t,

p(,t) is the fluid pressure,  is the dynamic viscosity of the fluid,  (,t) = (,t) ·n⃗(,t) is the fluid traction, and
(,t) is the unit normal vector pointing into Ω^f_t along Γ^_t.
Eq. (<ref>) describes the fluid momentum in Eulerian form, Eq. (<ref>) is the incompressibility constraint, 
Eq. (<ref>) is the kinematic condition along the fluid-structure interface, which implies the no-slip and no-penetration conditions,
and Eqs. (<ref>) and (<ref>)
describe the dynamics of the linear and angular momentum of the immersed rigid body in Lagrangian form.
Eqs. (<ref>) and (<ref>) also account for
the dynamic conditions at the fluid-structure interface because the rigid body forces are balanced by
the fluid traction.
 




 §.§ Immersed Lagrangian-Eulerian (ILE) formulation



We now introduce an immersed formulation of FSI that describes the same physical model as the partitioned 
formulation detailed in Sec. <ref>. As in that formulation, the
computational domain
is Ω, with ∈Ω indicating fixed physical coordinates.
As in Sec. <ref>, the structural kinematics are described in Lagrangian form via reference coordinates ∈Ω^s_0 attached to the solid, and we use the 
same motion map : (Ω^s_0 , t)↦Ω^s_t to determine the physical position of solid material point  at time t.
In the immersed formulation, however, we solve the incompressible Navier-Stokes equations on the full computational domain Ω, 
including both the fluid and solid subdomains. 
We split the computational domain Ω  into an exterior 
fluid region Ω^f,+_t and an interior fluid region Ω^f,-_t,
each parameterized by time t, with superscripts ‘+’ (‘-’) indicating values obtained from the ‘exterior’ (‘interior’) side of the fluid-structure interface.
Using this notation, we have Ω^f,-_t ≡_t and Γ^_t = Ω^f,+_t∩Ω^f,-_t.
We extend the definition of the fluid velocity $̆, pressurep, viscosity, and stress tensorto hold in the entire computational domainΩ, so that the extended fluid stress tensoris

    (,t) = -p(,t)   + ( (̆,t) + ^̆T(,t) ),  ∈Ω =Ω^f_t∪Ω^s_t.

Our approach applies a singular surface force density along the fluid-structure interface to impose the 
kinematic constraint, which implies a discontinuity in the traction associated with the extended fluid 
stress,, alongΓ^_t.
A jump in a scalar fieldψ(,t)at position= (,t)along the interface is denoted by

    ψ(x⃗,t)  = lim_ϵ↓ 0ψ(x⃗ + ϵn⃗(,t),t) - lim_ϵ↓ 0ψ(x⃗ - ϵn⃗(,t),t)
    	= ψ^+(x⃗,t) - ψ^-(x⃗,t),

in which·indicates the jump in the value across the interface,n⃗(,t)is the outward unit normal vector (into the exterior fluid region) 
alongΓ^_t,
andψ^+(x⃗,t)andψ^-(x⃗,t)are the limiting values approaching the interface positionx⃗from the exterior 
fluid regionΩ^f,+_tand interior fictitious fluid regionΩ^f,-_t, respectively.
By considering the jump in the extended fluid stress, the governing equations are 
   
    D/Dt(,t)    = - p(,t) +  ^2 (̆,t),    ∈Ω,
    ·(̆,t)    = 0,    ∈Ω,
    ((,t),t) ·n⃗((,t),t)    = -J^-1(,t)  (,t),     ∈Γ^_0,
    ∂/∂ t(,t)    = u⃗((,t),t),     ∈Γ^_0,
    d/dt∫_Ω^s_0 ∂/∂ t(,t)   d   =∫_Γ^_t^+(,t)  d a,
    d/dt∫_Ω^s_0( - _) ×( ∂/∂ t(,t)) d   = ∫_Γ^_t( - _(t)) ×^+(,t)   d a,

in whichJ(,t)is the surface Jacobian determinant <cit.>,(,t)is an interfacial surface force density that 
is the Lagrange multiplier to maintain the kinematic condition for the constraint in Eq. (<ref>) applied along the fluid-solid interfaceΓ^_t,
and^+(,t)is the exterior fluid traction.
Importantly, Eqs. (<ref>) and (<ref>), which account for the dynamic interface condition,  imply that only 
the fluid momentum and stresses from the exterior fluid subregion have any physical effect
in driving the dynamics of the structure.


The jump discontinuity in Eq. (<ref>) can be decomposed into discontinuities in the pressure and viscous stress, which in current coordinates are
   
    p(x⃗,t)    =  J^-1(^-1(,t),t)  (^-1(,t),t) ·n⃗(x⃗,t),     ∈Γ^_t,
    ∂u⃗/∂n⃗(,t)    = - ( - n⃗(,t)⊗n⃗(,t))   J^-1(^-1(,t),t)  (^-1(,t),t),     ∈Γ^_t.

Higher order jump conditions, including those associated with the first normal derivative of the pressure and the second normal derivative of the velocity, 
can be also derived <cit.>.





 §.§ A penalty approach to the ILE formulation

The formulation introduced in Sec. <ref> requires the solution of a saddle-point system that couples the Eulerian and Lagrangian variables <cit.>.
To develop a practical numerical method, we relax the kinematic constraint, Eq. (<ref>), by introducing two representations of 
the fluid-structure interface 
and applying penalty forces to penalize displacements between the two representations. 
This penalty method determines an approximate Lagrange multiplier force instead of solving for the Lagrange multiplier to exactly impose the condition <cit.>.
Specifically, along with the mapping(,t)that determines the kinematics of the
structure, we introduce an explicit representation of the fluid-structure interface that is parameterized via(,t)and that moves 
with the fluid, so that∂(,t)/∂t = u⃗((,t),t); see Fig. <ref>.
In this study, we use a penalty formulation similar to one proposed
 by Goldstein et al. <cit.>,
in which the rigidity constraint is approximately imposed through a linear spring-like force via

    (,t) = κ((,t)-(,t)),  ∈Γ^_0.

Here,κis the spring stiffness penalty parameter. This force penalizes deviations from
the constraint (,t) =  (,t)and, in the discretized equations, acts to ensure that the volumetric and surface meshes are at least approximately conformal in their motion.
Note that
it is possible to control the discrepancy between the two configurations because
 asκ→∞, the formulation exactly imposes the constraint that the two interfaces representations move together.
 



 §.§ Rigid body dynamics

The general formulation for description of kinematics of a rigid body includes both 
translational and rotational motions of material points in the body.
The kinematics in the current configuration coordinates of the solid can be written 
in terms of reference coordinates as

    (,t)=(t)  ( -_)   +  _  +  _(t),

in which_is the position of the center of mass in the reference coordinates,_(t)is the displacement of the center of mass,
and(t)is the rotation matrix.
The rotation matrix(t)can be expressed using the Euler angles,ϕ= ϕ(t),θ= θ(t), andψ= ψ(t), as

    (t)=
    [                     cos(ψ)cos(θ) -sin(ϕ)cos(ϕ)+cos(ψ)sin(θ)sin(ϕ)  sin(ϕ)sin(ϕ)+cos(ψ)sin(θ)cos(ϕ);                     sin(ψ)cos(θ)  cos(ϕ)cos(ϕ)+sin(ψ)sin(θ)sin(ϕ) -cos(ϕ)sin(ϕ)+sin(ψ)sin(θ)cos(ϕ);                          -sin(θ)                     cos(θ)sin(ϕ)                     cos(θ)cos(ϕ) ],

here using thex-y-zconvention in the rotation order. Following the approach of Akkerman et al. <cit.>, we avoid
explicitly forming the Euler angles
except for cases where there is only one rotational degree of freedom.
As will be presented in Sec. <ref>, 
in the time integration of the rigid body dynamics, the rotation matrix is an additional problem unknown
that is integrated in time along with the displacement and rotational degrees of freedom.
This significantly reduces the complexity of the calculation of the rotation angle for three degree of freedom (3-DOF) problems in two spatial dimensions and six degree of freedom (6-DOF) problems in three spatial dimensions.
Using Eq. (<ref>), the rigid body displacement and velocity are respectively defined as

    (,t) = ((t) - )(-_)   +  _(t)

and

    (,t)=(t)  (-_)   +  _(t).

Denoting the rigid body center of mass in the current configuration as_(t) = (_,t)=_   +   _(t), 
Eq. (<ref>) becomes

    (,t) =(t)^-1(t)((,t)-_(t))   +  _(t) =(t)  ((,t)-_(t))   +  _(t),

in which

    (t)=(t)^-1(t)=
    [       0 -ω_3(t)  ω_2(t);  ω_3(t)       0 -ω_1(t); -ω_2(t)  ω_1(t)       0 ]

is the skew-symmetric angular velocity matrix. The pseudovectorω⃗(t)=(ω_1(t),ω_2(t),ω_3(t))can be extracted from this tensor matrix.
Using the angular velocity vector, the rigid body velocity field in Eq. (<ref>) can be re-written in terms ofω⃗(t)as,

    (,t) = ω⃗(t)  × ((,t)-_(t))   +  _(t).

For a three-dimensional problem, the three components of the translational velocity of the center of mass_(t), together
with the three components of the angular velocityω⃗(t), form the six degrees-of-freedom that completely determine
the kinematics of the rigid body.

Linear and angular momentum conservation in the rigid body are described by a system of six ordinary differential equations,

    d/dt(m  _(t))    = ^net(t) ,
    d/dt((t)  ω⃗(t))    =^net(t) ,

in whichmis the mass of the solid object,^net(t)is the global force including the sum of all the forces exerted on the rigid body, and^net(t)is the net torque vector.(t)is the inertia tensor of the solid body in the current configuration, which is defined in terms of the inertia tensor in the reference configuration,_0, via

    (t) = (t)  _0  ^T(t),

in which_0is

    _0 = ∫__0 ( - _)·(-_)  d -∫__0 (-_)⊗(-_)  d .

In our fluid-structure interaction framework, the net force and torque vectors in Eqs. (<ref>) and (<ref>) 
are computed as

    ^net   = ( - )V  +  ∫_Γ^_t^+(,t)  d a,
    ^net   =∫_Γ^_t( - _(t))  × ^+(,t )  d a,

in whichVis the volume of the solid object,is the gravity vector, and^+(,t)is 
the exterior fluid traction vector exerted on the solid object by the fluid.
In the penalty formulation of the ILE method, quantities on the left hand sides of Eqs. (<ref>) and (<ref>) 
are evaluated onΓ^_t, which moves with the local fluid velocity.





§ DISCRETE EQUATIONS OF MOTION


This section introduces numerical methods for the penalty formulation of the immersed Lagrangian-Eulerian method detailed in Sec. <ref>. 
This approach leverages our immersed interface method (IIM) for discrete surface representations <cit.>. 
We include only the key aspects of this method; for additional details and benchmarking studies of problems with prescribed motion, see Kolahdouz et al. <cit.>. 
Standard methods are used for the rigid body dynamics. 
A second order accurate Strang time step splitting approach <cit.>  is used to obtain systems of equations that can be treated via efficient linear solvers.



 §.§ Eulerian discretization

The incompressible Navier-Stokes equations are discretized on an adaptively refined Cartesian grid using
a staggered-grid finite difference discretization <cit.> in which the pressure is approximated at the centers of the Cartesian grid cells and the components of the velocity
are approximated at the centers of the edges (in two spatial dimensions) or faces (in three spatial dimensions) of the grid cells.
Standard compact second-order accurate differencing schemes are used for the divergence, gradient, and Laplace operators. 
The discrete divergence of the velocityD⃗ ·$̆ is evaluated at the cell centers, whereas the discrete pressure gradient G⃗ p and 
the components of the discrete Laplacian of the velocity L u⃗ are 
evaluated at the cell edges (in two spatial dimensions) or faces (in three spatial dimensions).
For the nonlinear advection terms,
a staggered-grid version <cit.> of the xsPPM7 variant <cit.> of the piecewise parabolic method (PPM) <cit.> is used.
Physical boundary conditions are prescribed along the boundaries of the computational domain Ω as described previously <cit.>. 
Adaptive computations use a discretization approach described by Griffith <cit.> that employs Cartesian grid adaptive mesh refinement (AMR).

To account for the jump conditions along the fluid-solid interface that occur in the ILE formulation, we modify the definitions of
G⃗ p  and L u⃗ for those stencils
that cross the immersed interface.
Using generalized Taylor series expansions <cit.>, it can be shown that if the interface
 cuts between two Cartesian grid points at location _∘=(x_∘,y_∘,z_∘), such that x_i,j,k≤ x_∘ < x_i+1,j,k, with _i,j,k∈Ω^f,- and _i+1,j,k∈Ω^f,+, then for a piecewise differentiable quantity ψ, we have

    ∂ψ/∂ x(_i+,j,k)    =ψ_i+1,j,k-ψ_i,j,k/Δ x+sgn{n^x}/Δ x∑_m=0^2(d^+)^m/m!∂^mψ/∂ x^m__∘ + O(Δ x^2), 
    ∂^2ψ/∂ x^2(_i,j,k)    =ψ_i+1,j,k-2ψ_i,j,k+ψ_i-1,j,k/Δ x^2+sgn{n^x}/Δ x^2∑_m=0^3(d^+)^m/m!∂^mψ/∂ x^m__∘ + O(Δ x^2),

in which Δ x is the grid spacing in the x direction, ψ_i,j,k=ψ(_i,j,k), d^+= x_i+1,j,k-x_∘ > 0, and n^x is the x-component of the normal vector n⃗=(n^x,n^y,n^z) at the intersection point _∘.
The full implementation of this approach to the three-dimensional incompressible Navier-Stokes equations,
including the application of the jump  corrections to the stencils of the pressure and the viscous terms and algorithms for identifying intersections between the finite difference stencils and 
the discrete interface representation, is detailed in our earlier work <cit.>.






 §.§ Lagrangian discretization



Let _h be a triangulation of _0, the reference configuration of the volumetric rigid body, 
composed of elements U^e such that _h = ∪_e U^e, in which e indexes the mesh elements. We take
{_l}_l=1^M to be the positions of the M nodes of the mesh in the reference 
configuration, {_l(t)}_l=1^M 
to be the current positions of the nodes, and {ϕ_l()}_l=1^M 
to be the corresponding interpolatory nodal (Lagrangian) basis functions.
A continuous description of the configuration of the structure is defined by

    _h(,t) =∑_l=1^M_l(t)ϕ_l(),  ∈_0.


The configuration of the fluid-structure interface representation that moves with the fluid is described by the 
mapping :(Γ_0^,t)↦Γ_t^. 
To obtain a discrete representation of this interface, we use a surface mesh that corresponds to the restriction of the volumetric 
solid mesh to ∂Ω^s_0 ≡Γ_0^. 
For the discrete representation of the fluid-structure interface we have

    _h(,t) =∑_l=1^M_l(t)ϕ_l(),   ∈Γ_0^,

except that in practice, we only need to evaluate this sum over the lower-dimensional subset of nodes that are located on surface 
mesh, and not over all of the nodes in the volumetric representation. The reason is that the interpolatory (nodal) basis functions associated with interior nodes 
vanish on the surface. Similarly, the surface force density is determined by
 
    _h (,t) =∑_l=1^M_l(t)ϕ_l(),   ∈Γ_0^.

Again, this sum only needs to be evaluated using the M^fs surface nodes. In an implementation, it is straightforward to use separate data structures 
for the volumetric and surface structural meshes. For the remainder of the paper, we adopt the convention that all computations involving the 
surface representation are performed using only a surface mesh with 
appropriate nodal degrees of freedom and surface-restricted basis functions.


Stress jump conditions are imposed by evaluating the correction terms from the interfacial forces and interface configuration 
(i.e. generalizations of Eqs. (<ref>) and (<ref>)).
Geometrical quantities, including the surface normals and surface Jacobian determinant, that are needed by 
the IIM discretization, are obtained by directly differentiating Eq. (<ref>).
Note, however, that the standard nodal basis functions are C^0 but not C^1 at element boundaries, and so quantities that are obtained in terms of ∂_h/∂ are discontinuous in both the reference and current configurations.
In particular, the pointwise jump conditions determined from the mesh geometry and the surface Jacobian J are generally discontinuous between the elements.
Following the approach introduced in our prior work <cit.>, we obtain a continuous approximation to the jump conditions through the L^2 projection.
Briefly, given a function ψ∈ L^2(Γ_0^), its L^2 projection P_hψ onto the subspace 
S_h = span{ϕ_l()}_l=1^M^fs is defined by requiring P_hψ to satisfy

    ∫_Γ_0^(ψ() - P_hψ())  ϕ_l()  d A = 0,   ∀ l=1,…,M^fs.

The L^2 projection of a vector-valued quantity is determined component-wise.
Because the L^2 projection is defined via integration, the function ψ does not need to be continuous or even to have well-defined nodal values.
By construction, however, P_h ψ will inherit any smoothness provided by the subspace S_h.
In particular, for C^0 Lagrangian basis functions, P_h ψ will be at least continuous.
In our numerical scheme, we separately compute the projection of the normal and tangential components of the surface force per unit current area, J^-1_h(,t), onto S_h,
as needed to specify the conditions for the pressure and the viscous stress.
We drop the subscript “h" in the remainder of the paper to simplify notation.
To solve for the projected jump conditions, linear systems of equations involving the mass matrix ℳ need to be solved, in which ℳ has
components ℳ_kl = ∫ϕ_k ϕ_l.
Eq. (<ref>) is evaluated using seventh-order Gaussian quadrature. 
Notice that these projections are computed only along the fluid-solid interface and involve only surface degrees of freedom. Consequently, the computational cost of evaluating these projections is 
asymptotically smaller than the solution of the fluid equations.
Note also that similar to the conventional IB method, a force-spreading operator = [] can be defined to evaluate and apply 
the correction terms [] to the Eulerian discretization via a discrete Eulerian 
force density = [] <cit.>.

The velocity of the fluid-structure interface representation that moves with the fluid is determined by evaluating the Eulerian velocity (̆,t) 
on the interface. 
As detailed previously <cit.>, 
it is possible to interpolate the discretized Eulerian velocity field $̆ to the Lagrangian 
interface mesh using a corrected  bilinear (or, in three spatial dimensions, trilinear) interpolation scheme that accounts for the known 
discontinuities in∂u⃗/∂n⃗. In general, however, the basic interpolation scheme will produce 
an interface velocity field that is not in the space of the nodal basis functions, which implies that it cannot be used directly 
to update the configuration of the interface. 
To obtain a suitable surface mesh velocity field, we project the interpolated velocity field 
onto the space spanned by the nodal basis functions using Gaussian quadrature. This procedure implicitly 
defines a velocity-restriction operator= [,], so that= [,]   $̆.



The FSI coupling approach used herein crucially relies on the accurate evaluation of the exterior fluid traction.
This requires evaluating the exterior fluid pressure and exterior viscous shear stress.
To evaluate the exterior pressure at a position ∈Γ^_t, we use

    p^+_h(x⃗,t) =   p(x⃗,t)  + [p](x⃗^-,t),

in which p^-=[p](x⃗^-,t) is the interior pressure interpolated to a position ^- away from the interface in the opposite direction 
of the normal vector  and at a distance equal to 1.3 times the diagonal size of one grid cell. 
This factor has been chosen on an empirical basis <cit.>.
Here,  is the unmodified bilinear (or trilinear) interpolation operator involving quantities on one side of the interface.
To evaluate the exterior wall shear stress, a one-sided approximation to the normal derivative of the velocity is calculated using the same interfacial velocity reconstruction 
procedure that is used to determine the interface velocity along with the velocity value at a neighboring location in the direction of the normal vector x⃗^+.
As with the pressure, unmodified bilinear (or trilinear) interpolation is used to obtain the velocity away from the interface.
A one-sided finite difference formula is used to calculate the normal derivative,

    (∂/∂n⃗)^+_h(x⃗,t) = []̆(x⃗^+,t)-(̆x⃗,t)/ĥ,

in which the distance ĥ is chosen to be slightly larger than the diagonal size of the Cartesian mesh (1.05 times the diagonal size), 
so that regular bilinear (or trilinear) interpolation can be used to evaluate []̆(x⃗^+,t, ĥ) ensuring that the interpolation only uses values on one side of the interface.
It is possible to use a second-order formula with a three-point stencil which requires interpolating 
an additional point in the normal direction, but preliminary numerical experiments (data not shown) suggest the computation using only two points 
is more stable.
Moreover, as shown previously <cit.>, this simple scheme is adequate to achieve a point-wise first-order accurate approximation to the wall shear stress.
Note that as with velocity interpolation, the pressure and wall shear stress can be evaluated at arbitrary locations along the 
interface. We also use a surface-restricted L^2 projection to obtain nodal values of these interfacial quantities.





 §.§ Time integration scheme




In advancing from time step n at time t to time step n+1 at time t +,
we define a vector of variables Υ⃗ that includes all of the Eulerian and Lagrangian quantities that need to be updated
as Υ⃗=[,̆  p, , _, _, ω⃗, , ].
We use the second 
order Strang splitting scheme <cit.>, in which within three steps we: 
1) solve the rigid body dynamics equations over a half time step /2, treating the fluid traction as constant in time;
2) solve the IIM/FSI equations over a full time step , treating the configuration of the solid as constant in time; and
3) solve the rigid body dynamics equations over a final half time step /2, treating the fluid traction as constant in time.
The details of the time stepping for the rigid body dynamics equations, the IIM/FSI equations, and the overall algorithm are given below.



  §.§.§ Rigid body time integration scheme

Although we ultimately employ a time step splitting approach that advances the rigid body configuration in two half-steps, 
to simplify the discussion, the approach for a fixed time step size  is detailed here.
Denote by _ the action of the Lagrangian rigid body dynamics solver over the time increment  that
acts on a solution vector Υ⃗. This solution vector includes all of the Eulerian and Lagrangian variables, but only advances 
the volumetric structural variables _, _, ω⃗, , and  while keeping the remaining variables 
fixed.
Briefly, discretizations of Eqs. (<ref>)–(<ref>) are solved over the time increment  to obtain
^n+1_, ^n+1_, 
ω⃗^n+1,k, and ^n+1,k via

    m  ^n+1_ - ^n_/   =∫_Γ^_t^+()  d a,
    ^n+1_- ^n_/   =1/2(^n+1_ + _^n),
    ^n+1,k_0   (^n+1,k)^T  ω⃗^n+1,k-^n_0(^n)^Tω⃗^n/   = ∫_Γ^_t ( - _) ×^+()   d a,
    ^n+1,k - ^n/   =1/4(^n+1,k + ^n)(^n+1,k + ^n).

Eqs. (<ref>) and (<ref>) can be solved through a few subiterations
to obtain the new rotation matrix ^n+1 and 
the angular velocity ω⃗^n+1; starting from ^n+1,k=1≡^n, and ω⃗^n+1,k=1≡ω⃗^n 
at each time step, subiterations are performed by looping over k until ^n+1,k-^n+1,k-1_∞≤ϵ or
 ω⃗^n+1,k-ω⃗^n+1,k-1_∞≤ϵ 
with ϵ=10^-8. These iterations are inexpensive, and between one and three are typically needed to reach the convergence criteria. 
Seventh-order Gaussian quadrature rules are used to approximate integrals on the left side of Eqs. (<ref>) and (<ref>).
The updated configuration of (,t) is then evaluated using Eq. (<ref>). 
Similarly, the moment of inertia is determined using a discretization of Eq. (<ref>),

    ^n+1 =^n+1_0(^n+1)^T.






  §.§.§ IIM time integration scheme



Starting from ^n and u⃗^n at time t^n and p^n-1/2 at 
time t^n-1/2, we compute ^n+1, ^̆n+1, and p^n+1/2. 
Denote by _ the action of the IIM solver over a full time step that
acts on the solution vector Υ⃗ that includes all of the Eulerian and Lagrangian variables but only advances $̆,p, and.
Briefly, using the discrete velocity restriction operator^n=[^n,^n], we first obtain 
initial approximations to the interface position at timet^n+1/2via

    ^n+1 = ^n + /2^nu⃗^n, 
    ^n+1/2=^n+1+^n/2.

Next, we solve for^n+1,^̆n+1, andp^n+1/2via

    ρ(^̆n+1-^̆n/Δ t+A⃗^n+1/2)    =-G⃗  p^n+1/2+ L⃗(^̆n+1 + ^̆n/2)+  ^n+1/2, 
    D⃗·u⃗^n+1   =0, 
    ^n+1-^n/Δ t   = ^n+ = ^n+(^̆n+1+^̆n/2),

in whichA⃗^n+1/2=3/2A⃗^n-1/2A⃗^n-1is obtained from a high-order upwind spatial discretization
of the nonlinear convective term·̆$̆ <cit.>, and
^n+=[^n+,^n+] is the velocity restriction operator at the half time step configuration.
This time stepping scheme requires only linear solvers for the time-dependent incompressible Stokes equations.
We solve this system of equations by the flexible GMRES (FGMRES) algorithm with a preconditioner based on the projection method that uses inexact subdomain solvers <cit.>.
In the initial time step, a two-step predictor-corrector method is used to determine the velocity, deformation, and pressure; see Griffith and Luo <cit.> for further details.






  §.§.§ Fluid-structure interaction time stepping scheme



For the solution vector Υ⃗ over one time step using the Strang splitting scheme we have,

    Υ⃗^n+1 = _Δ t/2_Δ t_Δ t/2Υ⃗^n.

With this time-staggered approach, the overall algorithm to solve the fluid-structure problem 
is
	


	
  *  
		Solve the rigid body dynamics in Eqs. (<ref>)–(<ref>) 
		over a half time step Δ t/2 to advance from step n to n+1/2 and obtain 
		_^n+1/2, _^n+1/2, ω⃗^n+1/2 and ^n+1/2,
		and the new position of the volumetric Lagrangian mesh.
		
	
  *  
		Calculate the penalty force using the most recent position of the volumetric solid mesh and the surface mesh
		that moves with the fluid. 
	
	
  *  
		Solve for the IIM subsystem in Eqs. (<ref>)–(<ref>) over a full time step 
		and obtain the updated Eulerian velocity field u⃗^n+1 and pressure p^n+1/2 as well as the 
		the Lagrangian velocity ^n+1/2 and positions ^n+1 of the surface mesh and the exterior 
		fluid traction forces ^+.
	
	
  *  
		Using the exterior fluid traction force ^+ from Step III, solve the rigid body dynamics
		 in Eqs. (<ref>)–(<ref>) 
		over a half time step Δ t/2 to advance from step n+1/2 to n+1 and obtain 
		_^n+1, _^n+1, ω⃗^n+1 and ^n+1.
	
	
  *  
			Move the Lagrangian mesh of the bulk solid and obtain the new positions ^n+1.	






§ NUMERICAL EXAMPLES


This section presents computational examples to characterize the performance of the present methodology in two and three spatial dimensions. 
As a demonstration of the method's ability to tackle more challenging problems,
applications to two biomedical models are also presented, 
including simulations of the dynamics of a bileaflet mechanical heart valve 
in a pulse duplicator system and the transport of blood clots in a patient-averaged anatomical model of the inferior vena cava.
Where possible, comparisons are made to available experimental or computational results. 
We begin by considering model problems involving a limited number of translational, but not rotational, degrees of freedom (DOF). 
We systematically increase the complexity of the tests by incorporating additional translation and rotational degrees of freedom. We also consider both smooth 
immersed structures as well as structures with sharp corners.
The fluid-solid interface representation that moves with the fluid is discretized 
either by piecewise-linear (P^1) elements for two-dimensional cases, or by piecewise linear (P^1) or piecewise bilinear (Q^1)
elements for three-dimensional cases.
Unless otherwise noted, the structural meshes are constructed so that the ratio of the Lagrangian element size to the Eulerian grid spacing,
denoted by , is ≈ 2 at least along the fluid-structure interface.
In all cases, the Eulerian domain is discretized using an adaptively refined grid.
The Cartesian grid spacing on the finest level of the locally refined grid is h_finest=r^-(N-1)h_coarsest,
in which h_coarsest is the grid spacing on the coarsest level, r is the refinement ratio, and N is the number of refinement levels.
The spring penalty parameters are computationally determined as approximately the largest values allowed by our explicit time stepping algorithm at the time step sizes 
used in the simulations. In each case, the penalty parameter values are tuned using  bisection.
Unless otherwise noted, centimeter-gram-second (CGS) units are used.
For models involving gravitational forcing, gravitational acceleration is set to g = 981  cm·s^-2.
The large scale three-dimensional simulations in Secs. <ref>, <ref>, <ref>, 
and <ref> were performed using the Dogwood cluster provided by the Research Computing Division of University of North Carolina at
Chapel Hill Information Technology Services.
Each node is comprised of 2.4 GHz Intel Xeon E5-2699Av4 processors with the Broadwell-EP micro-architecture, 512 GB memory, and 44 cores per node, and nodes are connected by a high bandwidth Infini-band EDR switching fabric.
The smaller benchmark models are run on a workstation with two 32-core Intel Xeon E5-2680 v3 2.5 GHz processors and 32 GB of memory.





 §.§ Vortex-induced vibration of a cylinder



 
The problem of viscous flow past an elastically mounted two-dimensional cylinder undergoing vortex-induced vibration (VIV) 
has been widely studied both numerically and experimentally because of 
its broad range of engineering applications and its interesting vortex dynamics. 
This problem has also seen substantial use in benchmarking FSI algorithms <cit.>.
Here we consider two separate cases in two spatial dimensions, one using a single vertical degree of freedom, and the second with two degrees of freedom (2-DOF). The governing equations for the 2-DOF cylinder motion are

    d̈^̈ẍ_̈  + ḋ^̇ẋ_̇ + d^x_   = f^x, 
    d̈^̈ÿ_̈ + ḋ^̇ẏ_̇ + d^y_   = f^y,

in which d^x_ and d^y_ are the horizontal and vertical displacements of the center of mass, M_s is the mass
per unit length of the cylinder,  and  are the damping and stiffness constants characterizing the spring, and f^x and f^y are
the instantaneous drag and lift forces, respectively.
To facilitate comparisons to previous work, we define
the non-dimensional horizontal and vertical displacements of the cylinder center in the streamwise and transverse directions, respectively, 
as d^x_=d^x_/D 
and d^y_=d^y_/D, in which D is the diameter of the cylinder.
Taking U_∞ as the free stream flow velocity, 
the mass ratio and reduced velocity are respectively defined as =/ and =U_∞/(f_n D), in which
f_n=√(/)/(2π) is the natural frequency of the structure. The damping ratio is ζ = /(2√()).


 


  §.§.§ 1-DOF transverse oscillation


We first consider
the benchmark problem of a circular cylinder undergoing VIV with a single vertical degree of freedom.
We are interested in capturing the well characterized vortex “lock-in” phenomenon observed in previous studies <cit.>.
Within the lock-in regime, the vortex shedding frequency is close to the natural frequency of the structure, which results
in large amplitude vibrations. Physical parameters are selected to match benchmark results in the literature.
A schematic of the simulation setup is shown in Fig. <ref>(a).
The computational domain is Ω =[-30  cm,45  cm]×[-30  cm,30  cm],
a rectangle of size L_x × L_y = 75  cm× 60  cm.
The cylinder has diameter D=1 cm, is initially at rest and centered at the origin.
A uniform inflow velocity U⃗=(U_∞ = 1  cm·s^-1 ,0  cm·s^-1) is imposed on the left boundary (x=-30 cm), 
and zero normal traction and tangential velocity outflow conditions are imposed at the right boundary (x=45 cm).
Along the bottom (y=-30 cm) and top (y=30 cm) boundaries, zero normal velocity and tangential traction are imposed.
The domain is discretized using N=6 nested grid levels, with coarse grid spacing 
h_coarsest=L_y/64=0.9375  cm and refinement ratio r = 2 between levels,
leading to h_finest≈ 0.029  cm.
With =2, this results in approximately 54 linear elements around the perimeter of the disk.
The time step size is Δ t= (0.1  s·cm^-1)   h_finest with a penalty spring constant κ=(0.00125  g·cm^-2)/^2.
The Reynolds number = U_∞ D/ is fixed at 150, the damping is set to zero (ζ=0), and the mass ratio is m^* = 8/π. 


First, the effect of the reduced velocity on the maximum displacement of the cylinder is studied by systematically varying  within the range 3 ≤≤ 8.
Vortex shedding occurs in all cases. Fig. <ref>(b) shows the maximum displacement with respect to .
These results demonstrate that for ∈ [4,7] there is a large increase in the vibration amplitude.
Fig. <ref>(b) also compares results obtained by our method to previous numerical studies,
including
a geometrically conservative finite volume ALE method <cit.>, a curvilinear immersed boundary method 
<cit.>,
and a finite element based ALE approach <cit.>.
Excellent agreement is obtained over the full range of  values considered here.



In addition, we perform a grid refinement study using the reduced velocity with the largest maximum displacement (i.e. =4)
while fixing the previous values of all other parameters. To achieve this in our AMR framework, we vary the number of refinement levels N between 5 and 7.
Fig. <ref> shows the time-history of the maximum displacement for N=5, 6, and 7 levels of refinement.
The displacement values obtained for N = 6 and 7 closely match each other, whereas the coarser case, using N=5 levels of refinement, under-predicts the maximum displacement 
in the vortex shedding region. 
Notice that these results indicate that using N=6 levels 
of refinement provides essentially grid-converged results for this benchmark problem.

          






  §.§.§ 2-DOF oscillation




The 2-DOF oscillation is studied using the same physical parameters as Blackburn and Karniadakis <cit.>, who used 
a spectral element approach.
The size of the computational domain, the cylinder diameter and initial position, and the physical boundary conditions are all the same as the 1-DOF example in Sec. <ref>.
The Reynolds number is =200, the reduced velocity is U^* = 0.5, the damping ratio is ζ = 0.01, and the mass ratio is m^*=4/π.




We use this example to investigate the order of spatial convergence of the difference in the displacement
between the two Lagrangian 
representations (i.e. the fluid-structure interface representation that moves with the fluid, and the surface of the volume mesh).
As discussed in Sec. <ref>, the two representations are tethered together by
spring-like forces that penalize relative motion between them. 
The Cartesian grid spacing on the coarsest level for all four cases is h_coarsest=60/64 cm. 
We consider four sets of discretizations, ranging from a very coarse composite grid with N=4 levels of refinement to a reasonably fine grid with N=7
levels of refinement, using refinement ratio r=2 between levels.
The penalty spring constant is κ=(0.125  g·s^-2)/h^2, and a small constant time step size = 0.001 s is chosen for all cases.
In our recent work <cit.>, it was shown that
the present IIM algorithm yields second-order convergence in the displacement of the interface
in both Ł2 and  with suitable scalings for the penalty parameter κ.
As in our previous work, for the choice of the κ used here, pointwise second order convergence is also expected for the displacement between the positions of 
the two interface representations. 
Fig. <ref> shows (,t)-(,t)_∞,
the  norm of the discrepancy between the Lagrangian points at the surface mesh and the corresponding points on the boundary of the volume mesh. 
Fig. <ref>(a) indicates that the method converges at second order in the maximum norm at t=55 s.
To investigate the change of the  norm of the error over time, this value is plotted on a semi-log scale in Fig. <ref>(b) for
the four discretizations.
The overall behavior of the error norm remains consistent for all discretizations throughout the entire simulation.
Note that the difference shown in ξ⃗(X⃗,t)-χ⃗(X⃗,t) is proportional to the tethering penalty force that connects the two representations.
Because of the periodic “figure-of-eight”  nature of the cylinder's dynamics, a undulatory pattern in the loading force is to be expected over time.
This pattern is expected to be present in ξ⃗(,t)-χ⃗(,t) difference as well due to its proportionality with the force.
Specifically if we wish to achieve ξ⃗(,t)-χ⃗(,t) = O(h^2), it is necessary 
that the penalty parameter κ also satisfies κ = O(1/h^2), 
so that an applied penalty force of the form = κ (ξ⃗(,t)-χ⃗(,t)) satisfies = O(1) under grid refinement <cit.>.  
 


As another verification, the centerline displacement response of the cylinder is compared to 
the numerical results of Blackburn and Karniadakis <cit.>.
It is expected that the vortex shedding of the oscillating cylinder will lead to a periodic “figure-of-eight”  behavior. 
This is shown in Fig. <ref> 
for three locally refined Cartesian grids with N= 5, 6, and 7 levels along with the solution from 
Blackburn and Karniadakis <cit.>.
The trajectory clearly converges under grid refinement. Moreover, the trajectory of the cylinder
for the finest discretization agrees extremely well with the spectral element solution <cit.>.
The origin of oscillations x_c, which is defined as the intersection point in the periodic figure-of-eight trajectory of the cylinder's
center of mass, and the Strouhal number, which is calculated as =f D/U_∞ with f representing the oscillation frequency, are
reported in Table <ref> for the finest ILE computation (N=7) along with the results 
of Blackburn and Karniadakis and several other studies.
The center of oscillations agrees very well with the original work of Blackburn and Karniadakis <cit.>. Moreover, 
the Strouhal number matches the values reported
from studies by Yang et al. <cit.>, Yang and Stern <cit.>, Kim et al. <cit.>, and Qin et al. <cit.>.
The centerline trajectory and the dimensionless displacement-velocity phases of x/D-u_x/U_∞ and y/D-u_y/U_∞ for this case are shown in
Fig. <ref>.
The phase response obtained from our solution agrees well with previous observations reported 
by Yang and Stern <cit.> and Liu and Hu <cit.>.




  




  §.§.§ 2-DOF oscillation with low mass ratio


We now briefly consider significantly lower mass ratios than the one studied above.
We first consider a case with density ratio of m^*=0.4/π, which is 10 times smaller than the first example.
The remainder of the simulation
parameters are identical to the finest case in the previous example (Sec. <ref>). Fig. <ref>(a) shows that the cylinder
travels more than 7 diameters downstream before undergoing
the same periodic figure-of-eight motion.
Fig. <ref> also reveals that the oscillations for this light cylinder occur
at a higher amplitude and with larger horizontal and, to a lesser extent, vertical velocity magnitudes.    
Fig. <ref> shows the 
instantaneous vorticity contours around this light structure at six time points.
Vortex shedding is observed, with two vortices
shed during one cycle of the quasi-periodic oscillation. The wake footprint and the
associated vorticity patterns
also indicate that the cylinder is undergoing periodic motion.
Notice that previous work has reported severe instabilities
in computing such cases using both weak and strong coupling approaches <cit.>. 
The present method appears to remain stable even for extremely light structures. 
Table <ref> reports
the center of oscillations along the x-axis 
for a wide range of mass ratios along with the results of Kim et al. <cit.>, the only other study that we are aware of to
also consider such low mass ratios.
To the authors' knowledge, results obtained using mass ratios smaller than 0.3 have not been previously reported for this benchmark case.
The present method yields a slightly larger distance than the work of Kim et al. <cit.> for the origin of the oscillations. 
In the work of Kim et al. <cit.>, the mass ratio of m^*=0.3 was reported as the lowest ratio for a stable solution using
a strong coupling approach.
The present method is able to predict the dynamics at substantially smaller mass ratios (nearly two orders of magnitudes smaller) as shown for m^*=0.005 in Table <ref>.
Although stable results were obtained for lower density ratios, the accuracy of the dynamics requires
further validation that can be performed once additional experimental data are available.
Note that to accommodate the extended horizontal movement of the cylinder for m^* < 0.1, the computational domain for these cases was changed to Ω = [-120  cm, 180  cm]×[-120  cm,120  cm],
a rectangle of size L_x × L_y  = 300  cm×  240  cm.
To keep the grid spacing the same as before, N=8 nested grid levels were used, and the Cartesian grid spacing on the coarsest level was set to
h_coarsest=L_y/64.






  §.§.§ Stability and sensitivity studies for κ




Because we use an explicit time stepping approach to couple the fluid and solid degrees of freedom, 
if we fix the spatial and temporal discretization parameters, there will be a largest stable value of 
the penalty spring parameter κ. As mentioned at the beginning of Sec. <ref>, 
in practice, we compute using values of κ close to the stability limit, which we approximately determine using the method of bisection.



To study the sensitivity of the numerical algorithm to the penalty parameter κ and the influence of this parameter on the accuracy of the 
computed solution, we
 reconsider the 2-DOF vortex-induced vibration of a light cylinder with density ratio of m^*=0.4/π (Sec. <ref>). 
 
  
  
Except for the values of κ and the time step size , all other simulation parameters are the same as Sec. <ref>. 
We use a locally refined grid with N=7 levels of refinement, and the final simulation time is t=200 s, which is long after vortex induced vibration has been established. 
  We computationally determine minimum and maximum values of κ that 
  1) satisfy a minimum threshold in displacement between the two Lagrangian representations, 
  and 2) remain stable. 
  Accuracy thresholds are determined to be the minimum values of κ required to achieve tolerances of 0.1   h_finest, 
  which is commonly used in IB models of rigid structures, along with a tighter tolerance of 0.06   h_finest. 
  The maximum stable value of κ determines a stability threshold.
  
Fig. <ref>(a) shows 
the computationally determined accuracy and stability thresholds for this case.



Because it is not practical to determine the precise values of these thresholds across all time step sizes,
error bars indicate the small uncertainty in the lower and upper bounds.
For the smallest time step size considered here, 
the range of values of κ that satisfy both criteria spans nearly three orders of magnitude, 
from the minimum value that achieves the minimum accuracy 
criterion to the maximum value that results in a stable computation with 0.06   h_finest displacement threshold.


For the largest time step, there is at least about three-fold difference between the minimum and maximum values of κ.
Fig. <ref>(b) shows the effect of κ on the ratio between the maximum norm of the two Lagrangian representations
and the finest grid spacing for three stable spring constants of κ = 500, 1000,  and  1400 using the time step =0.001.
Notice that the maximum displacement between the representations scales like κ^-1.
Fig. <ref>(c) shows the trajectory of the center of mass of the cylinder for the same three κ values.
Deviations in the trajectories are small and converge as κ increases.
It is worth noting that the maximum discrepancies for the highest and lowest choices of
κ are less than the spacing on the finest grid level, h_finest.
 




 §.§ Galloping rectangular structure




This example uses a rectangular plate undergoing galloping motion to test the accuracy of the method for models involving only rotational degrees of freedom.
Flow-induced rotational galloping oscillations occur
in many areas of structural <cit.>, aeronautical <cit.>, and mechanical <cit.> engineering applications. 
This problem has also been widely used as a benchmark problem to validate numerical 
algorithms <cit.>.
In this section, only a single rotational degree of freedom is used, and
the translational heave (horizontal) and surge (longitudinal) motions are eliminated.
The governing equation for the mass-spring-damper model with one rotational degree of freedom is

    ^θθ̈  + ^θθ̇ + ^θθ = M^θ,

in which θ is the rotational angle of the body, ^θ is the rotational moment of inertia, ^θ is the torsional damping constant,
K_s^θ is the torsional spring constant, and M^θ is the moment acting on the rigid structure from 
exterior fluid forces.
In our simulations, we consider a rectangular structure with a width-to-thickness ratio of Λ=L/H=4.


To enable comparisons with prior studies, we define the non-dimensional moment of inertia by ^*=^θ/( H^4), 
the non-dimensional damping ratio by ζ_s^θ=^θ /(2√(^θ^θ)), 
and the reduced velocity by U^*=U_∞/(f^θH).
In the latter formula, U_∞ is the free stream velocity and f_s^θ=√(^θ/^θ)/2π is the natural frequency of the body. 
Following the work of Robertson et al. <cit.>, the 
non-dimensional parameters are taken to be I_s^*=400, U^*=40, and ζ_s^θ=0.25.
A schematic of the computational setup and the boundary conditions is given in Fig. <ref>(a).
The rectangular structure is centered at the origin with an initial zero angular velocity and a non-zero angle of θ_0=1^∘.
The computational domain is Ω = [-32  cm,96  cm]×[-32  cm,32  cm],
a rectangle of size L_x × L_y = 128  cm× 64  cm.
The domain is discretized using N=4 nested grid levels, with coarse grid spacing h_coarsest=L_y/32 = 2.0  cm and 
refinement ratio r = 4 between levels, leading to h_finest = 0.0625  cm.
A uniform inflow velocity of U⃗=(U_∞=1 cm·s^-1,0 cm·s^-1) is imposed on the left boundary (x=-32 cm).
Using the free stream velocity, the Reynolds number = U_∞ H/
is set to 250.
A penalty spring constant of κ=400 g·cm^-2·s^-2 is used, and the time step size
 is Δ t=(0.02  s·cm^-1)   h_finest.
Zero normal traction and tangential velocity are imposed at the right boundary (x=96 cm).
Along the bottom (y=-32 cm) and top (y=32 cm) boundaries, the normal velocity and tangential traction are set to zero.
Fig. <ref>(b) shows the time history of the pitch angle of the galloping rectangle
for the damped oscillation of the structure.
Once the vortex shedding state is established, the structure starts to undergo a
periodic rotation with well-characterized frequency and upper bound of the maximum angle.
This behavior is clearly captured by the present method.
Fig. <ref> panels (c)–(e) show
snapshots of the structural rotation along with the vortex structure of the flow at three different times.


To demonstrate the ability of the algorithm in modeling larger rotational angles, an additional case is
considered with ζ_s^θ=0. The initial angle is set to θ_0=-5^∘. All the other parameters, 
including the time step size and penalty spring constant, are the same as before.
As seen in Fig. <ref>(a), the method generates periodic behavior 
with maximum galloping amplitude of approximately θ_max=124^∘.
Fig. <ref> panels (b)–(d) show the flow patterns around the structure at three different times.


To compare the rotational response of the structure with previous work,
the maximum pitch angle and frequency for the above two cases are reported in Table <ref> along with values from other studies. 
The work of
Robertson et al. <cit.> uses a body-fitted spectral element method 
in a non-inertial reference frame.
The methods used by Yang and colleagues <cit.>
are different variations of a strongly-coupled direct forcing approach with 
a field extension strategy for the pressure/velocity derivatives.
Table <ref> demonstrates 
excellent quantitative agreement of both vortex shedding characteristics of θ_max  and the galloping frequency f^θ
for the two cases in comparison with other numerical studies.





 §.§ Freely falling rectangular plate
 


We now consider a model of a falling rigid plate in a water tank that is based on the experiments by Andersen et al. <cit.>.
This example tests the action of the instantaneous fluid forces on an object with sharp corners that leads to extremely complex trajectories.
This benchmark case has also been investigated in the context of fluid-structure interaction
algorithms in prior studies  <cit.>.
We consider two cases from Anderson et al. <cit.>, one undergoing 
fluttering motion at =1147 and the other undergoing tumbling motion at =837.
Different modes of fluttering and tumbling motions 
are captured in both the experiments <cit.> and the simulations.






  §.§.§ Fluttering motion
 



For the case of a freely fluttering plate, the plate thickness is taken to be H = 8.1  × 10^-2 cm with width-to-thickness ratio of Λ = L/H = 14.
  
 
The density of the plate and fluid are = 2.7 g·cm^-3 and = 0.997 g·cm^-3, respectively. To achieve a Reynolds number
of =1147, the fluid viscosity is set to = 8.87 × 10^-3 g·cm^-1·s^-1.
A schematic of the problem setup is shown in Fig. <ref>(a).
The computational domain is Ω = [-40  cm , 40  cm] × [-56  cm, 24  cm], 
a square of size L_x × L_y = 80  cm× 80  cm.
The center of the plate is initially located at (x_0, y_0)=(5  cm,0  cm) 
with an initial angle of θ_0 = 60^∘ with respect to the x-axis.
Zero normal traction and tangential velocity conditions are imposed at the top boundary. 
Along the left and right boundaries zero normal velocity and tangential traction are imposed.
The no-slip condition is imposed at the bottom boundary.
 

The domain is discretized using N=6 nested grid levels, with coarse grid spacing h_coarsest=L_y/16 = 5  cm and refinement ratio r = 4 between levels,
leading to h_finest = 0.00488  cm.
Using =1.8, this leads to the thickness of the plate being discretized by approximately 10 linear elements. 
A constant time step size of Δ t=0.01 ms is used, and the penalty spring constant is κ=7.45 × 10^5 g·cm^-2·s^-2.
Figs. <ref>(b) and  <ref>(c) shows the overall dynamics of the plate during its fluttering free fall. 
Fig. <ref>(b) shows
the vortex structure over the course of the simulation.
Periodic fluttering motion is clearly observed. 
Fig. <ref>(c) shows the trajectory of the center of mass
along with the experimental data of Andersen et al. <cit.>.
Overall, our numerical results are in good agreement with the experimental data for the gliding motion of the plate from side to side as it flutters in its free fall.
Discrepancies around turning points can be, in part, attributed to the interaction of the complex vortex dynamics with sharp corners of the rectangular object.
There are also uncertain differences in the simulation and experimental operating conditions that could explain these differences.
Future studies should further investigate the accuracy of the present IIM coupling strategy for objects with sharp corners, including assessing grid sensitivity.


Fig. <ref> shows close-up views of the fluttering plate at six different time points.
With a low angle of attack, the plate glides from a turning point on one side to reach a new turning point on the opposite side.
Towards the end of its glide and before it changes direction, there is a sharp increase in the magnitude of the angular velocity.
Because of the sharp edge of the structure,
the flow separates on the lower surface (Fig. <ref>(c)) and shortly thereafter flow separation also 
occurs on the upper surface as well (Fig. <ref>(d)). As the plate pitches upward,
it begins to glide in the opposite direction, and this process repeats itself in a periodic manner.
The interaction of vortices in the locations where the plate reverses direction creates a complex pattern of vortices, as shown in the figure.
Fig. <ref> shows
the time history of the horizontal and vertical velocity components as well as the angular velocity
along with a comparison to the experimental results <cit.>. Overall, there is very good agreement between the two results despite the complex dynamics.
Although the algorithm is able to correctly predict the periodic dynamics, 
	small deviations from measured data are observed particularly around troughs and crests of both translational and angular velocity plots. 
	We speculate that these deviations could potentially be due to sensitivity of the results to the complex interaction of the shed vortices 
	with the sharp corners, or possibly other uncertainties that are not accounted for in the comparison.


  


  §.§.§ Tumbling motion
 

For the tumbling case, the plate thickness is kept fixed at H = 8.1 × 10^-2 cm, but the width-to-thickness ratio is changed to Λ = L/H = 8.
The initial location and angle of attack are  (x_0, y_0)=(-20  cm,0  cm) and θ_0 = -45^∘, respectively. 
Note that because of the different width of the plate in this example, the Reynolds number is =837.

The penalty spring constant is κ=5.5 × 10^5 g·cm^-2·s^-2. 
The remaining simulation parameters, including
the computational domain extent and size, spatial resolution, time step size, boundary conditions, and fluid properties, are identical to the fluttering case;
see also the schematic in Fig. <ref>.  



Fig. <ref> 
shows the overall dynamics of the tumbling plate.
As shown in Fig. <ref>(a), the complex vortex structure of the tumbling plate is well resolved by the simulation.
After the plate is released, it begins a gliding motion. 
Shortly after, it pitches upward, similar to the fluttering case. 
Because of the large angular momentum, however,
at the turning point the plate rotates more than 90^∘ clockwise. 
This large rotation, creates a large restoring moment that
causes the plate to rotate slightly counter-clockwise and then continue falling, with an inclination to the 
right side. 
The plate travels a path towards the bottom-right corner of the computational domain by a sequence of descending and accelerating motions.
At the turning points, the plate undergoes a full 360^∘ tumbling rotation, such that the lower surface during the gliding re-configures as upward facing.
Fig. <ref>(b) compares the vortex structure and trajectory of the tumbling plate 
to experimental measurements <cit.>.
As in the fluttering case, the results are in very good agreement with the experimental data, and the trajectory of the plate agrees very well for 
the portion of the trajectory where experimental data are available. 
Fig. <ref> shows enlarged views of the vorticity field at different times.
Flow separation on the lower side of the plate is clearly observed in Fig. <ref>(c). This is followed
by separation of the flow on the opposite side as the plate rotates (Fig. <ref>(d)). The gliding and diving 
towards the bottom left resumes at the end of each tumbling motion, as seen in Fig. <ref>(f).

Fig. <ref> shows the time history of translational and rotational velocities.
Unlike the fluttering case, in which the vertical velocity appeared to have twice as larger frequency than the other two velocities,
here the number of periods demonstrates that approximately the same frequency is observed among all velocities. In addition, the significantly larger angular velocity in this case compared to the fluttering case in Fig. <ref>
indicates faster rotation near the turning points.
The average horizontal and vertical velocity components and the average angular velocity obtained from three full cycles
of the present simulation are found to be u_x=15.87 cm·s^-1, u_y=-11.32 cm·s^-1, 
and ω_z=15.95 rad·s.
For comparison, the experimental measurements of Andersen et al. <cit.> 
for the same average velocities are 
u_x=15.94 cm·s^-1, u_y=-11.5 cm·s^-1, and 
ω_z=14.5 rad·s^-1.
The relative discrepancies in these quantities are 0.19%, 1.57%, and 9.93%, respectively.



 

 


 §.§ Free falling of a dense sphere


This section investigates the dynamics of a freely falling steel bead in 
water using an unconstrained rigid-body structure model.
This problem follows the experimental setup of Mordant and Pinton <cit.>.
The density of the steel bead is =7.85 g·cm^-3, and the diameter is D=0.1 cm.
The Reynolds number is =(V̅ D)/ =430, and a terminal velocity of V̅=38.3 cm·s^-1 is reported in the original work.
With a fluid density of =997 g·cm^-3, this results in a dynamic viscosity of =8.88× 10^-3 g·cm^-1·s^-1. 
The computational domain is Ω = [-16 D,16 D]×[-448   D,64   D]×[-16   D,16   D],
a rectangular cuboid of size L_x × L_y × L_z = 32 D  × 512 D  × 32 D; see Fig. <ref>(a) for a schematic of the problem setup.
Zero normal traction and tangential velocity conditions are imposed at the top boundary.
Along the peripheral walls zero normal velocity and tangential traction are imposed.
The no-slip condition is imposed at the bottom boundary.
The domain is discretized using N=3 nested grid levels, with coarse grid spacing 
h_coarsest=L_x/64 = 0.05  cm and refinement ratio r = 4 between levels,
leading to h_finest=0.003125  cm.
The volumetric mesh of the sphere consists of hexahedral elements leading to a surface mesh composed of  bilinear quadrilateral elements with =2. 
A fixed time step size of =0.01 ms is used, and the penalty spring constant is κ=5 × 10^5 g·cm^-2·s^-2. 
The center of the sphere is initially positioned at the origin and released with zero initial translational and angular velocities. 
Fig. <ref>(b) shows the time history of the vertical velocity. 
After being released, the steel bead starts to accelerate until reaching its terminal velocity. 
The acceleration of the bead is in good agreement
 with the experimental data of Mordant and Pinton <cit.>, and the terminal velocity
 is also in excellent agreement with the experimental measurement.
 Fig. <ref>(c)–(e) show isosurfaces of the Q-criterion <cit.> to visualize the vortex dynamics.  
Note that there is no constraint in the motion of the sphere, and that all of the degrees of freedom including the rotational ones are included in the solution.
Shortly after the sphere reaches its terminal velocity V̅=38.3 cm·s^-1, the flow behind the sphere becomes
irregular due to the moderately large Reynolds number. Although small lateral movements
are observed in the simulation with the deviation of center of mass reaching to about 0.25D = 8 h_finest, no consistent pattern of “zig-zagging” motion is apparent.
Fig. <ref> also compares our results to two other numerical studies <cit.>. 






 §.§ Free falling of a sphere with near-unity density ratio


This section explores the performance of the ILE methodology with respect to two challenging aspects of FSI: the influence of the wall and the ability of the method to handle near-unity 
density ratios in a fully unconstrained rigid-body motion.
This case is based on the experimental data reported by ten Cate et al. <cit.> of a sphere
settling in a confined flow chamber.
The sphere has diameter D=1.5 cm and density =1.120 g·cm^-3, which is close to that of the  surrounding fluid; see Table <ref>.
The Reynolds number is =V̅ D /, in which the terminal velocity V is taken to be the settling velocity in an infinite domain.
Two different cases with Reynolds numbers of =1.5 and =31.9 are considered, which are achieved in the experiments by varying the density and viscosity of the fluid.
The computational domain is Ω = [-5  cm, 5  cm]×[-12  cm, 4  cm]×[-5  cm, 5  cm],
a rectangular cuboid of size L_x × L_y × L_z = 10  cm× 16  cm× 10  cm.
The sphere is initially positioned in the x-z mid-plane 
and towards the top of the box at a height of 12 cm, measured from the bottom of the sphere to the bottom of the box; see Fig. <ref>(a).
Zero normal traction and tangential velocity conditions are applied at the top boundary while no-slip condition is imposed at all other boundaries.
 


The domain is discretized using N=3 nested grid levels, with coarse grid spacing h_coarsest=L_x/32 =0.3125  cm and refinement ratio r = 2 between levels,
leading to h_finest=0.078125  cm .
The volumetric mesh of the sphere consists of hexahedral elements leading to a surface mesh composed of bilinear quadrilateral elements with =2. 
A fixed time step size of =0.5 ms is used, and the penalty spring constant is κ=5 × 10^4 g·cm^-2·s^-2.
Because of the confinement effect of the walls, we expect the settling of the sphere to be different from a free falling sphere in an infinite medium.
Upon release, the particle accelerates until reaching its terminal velocity, which occurs quickly in the lower Reynolds number case.
The sphere then starts to decelerate as it approaches to the bottom the flow chamber. Fig. <ref> panels (b)–(c)
show the time histories of the non-dimensional vertical gap between the sphere and bottom wall (h/D) and the 
vertical velocity of the sphere (V) for both cases. 
The method yields excellent agreement with the 
experimental measurements of ten Cate et al. <cit.> for both the gap size over time and the terminal velocity. Specifically, the present method 
captures the dynamics of the sphere while also remaining stable even as the sphere reaches the bottom of the box. For the lower Reynolds number case, the sphere has 
more time to travel through the box because of its lower speed, and this results in a larger flat region, when
the sphere experiences a nearly constant downward velocity, before it starts to decelerate as it approaches the bottom wall. Indeed, at =1.5, we observe stable 
sedimentary motion of the sphere in the simulation even past the latest time reported in the experiment.
At =31.9, the sphere does not experience an extended period of time at terminal-like velocity as it quickly reaches the bottom of the box. 
The simulation is also stable for this case even after the sphere impacts the bottom; see Fig. <ref>(b).




 §.§ Freely rising sphere
 


 
  
This section considers a rising sphere using a fully unconstrained rigid-body structural model.
Such problems can pose substantial challenges to FSI algorithms using weak coupling schemes, which can become unstable <cit.>.
Experimental results of Horowitz and Williamson <cit.> are used as a benchmark 
to validate the dynamics generated by our numerical method. We consider a challenging case in which the parameters
are chosen such that the sphere oscillates periodically and vigorously in a  “zig-zag” trajectory within a tight vertical plane.
The sphere has diameter D=1.3 cm and mass ratio m^*=0.11.
The experiments are at =  U D/ = 450. To match the experimental Reynolds number,
the dynamic viscosity is = 0.125 g·(cm·s)^-1.
The computational domain is Ω = [-16   D,16   D]×[-112   D,16   D]×[-16   D,16   D]; see the schematic in Fig. <ref>.
The domain is discretized using N=3 nested grid levels, with coarse grid 
spacing h_coarsest=L_x/64=0.65  cm and refinement ratio r = 4 between levels,
leading to h_finest≈ 0.04  cm.
The volumetric mesh of the sphere consists of hexahedral elements leading to a surface mesh composed of bilinear quadrilateral elements with =2. 
The center of the sphere is initially positioned at (0, -96   D, 0) and is released with zero initial translational and angular velocity.
Zero normal traction and tangential velocity conditions are imposed at the top boundary.
The no-slip condition is imposed at the bottom boundary.
Along the peripheral walls zero normal velocity and tangential traction are imposed.
The no-slip condition is imposed at the bottom boundary.
A fixed time step size of =0.01 ms is used, and the penalty spring constant is κ=5.2 × 10^5 g·(cm·s)^-2.
 

Fig. <ref> demonstrates that the sphere motion is approximately planar, and captures the periodic zig-zag trajectory observed experimentally. We emphasize that in our simulation, 
no constraints are imposed on the motion of the sphere. Specifically we do not impose either planar motion or the zig-zag trajectory. 
According to Horowitz and Williamson <cit.>, the periodic motion of the sphere at this particular density ratio
resembles the dynamics at much higher Reynolds number and is always confined to a plane. 
The simulation results 
predict the same kind of in-plane motion as observed in the experiment. As the sphere goes through the zig-zag motion it creates a complex
but organized wake pattern with vortex rings forming at the turning points. The numerical simulation clearly captures the vortices behind the sphere.
Further, there is generally excellent agreement with the trajectory obtained from the experiment.




 §.§ Bileaflet mechanical heart valve at physiological conditions



 
 
Bileaflet mechanical heart valves remain commonly used in heart valve replacement because of their durability.
In this section, we consider the simulation of fluid flow through a geometrically realistic model of a 25 mm St Jude Medical RegentTM bileaflet 
mechanical heart valve in the aortic test section of an experimental pulse duplicator platform <cit.>.
These systems are used in academia, 
industry, and regulatory agencies to assess the performance of prosthetic heart valves.
This simulation aims to demonstrate the ability of the present approach in modeling a realistically complex problem 
that involves multiple moving and stationary parts, 
including the small-scale feature of the hinge geometry.
These parts come in close contact with each other and undergo substantial pressure loading when closed.
In particular, we use experimental pressure and flow data obtained from a real pulse duplicator to establish realistic boundary models for the FSI model. 
These boundary models are calibrated in isolation from the rest of the system in that 
the motion of the leaflets, including the timing of valve opening and closing, is not prescribed and the simulated pressures, flow rates, and 
leaflet kinematics all emerge from integrating these three model components.
The experiments conducted to calibrate the model used saline as a test fluid, 
which we model as a Newtonian fluid with uniform density =1.0 g·cm^-3 and uniform dynamic viscosity μ=1.0 cP.
The computational domain is 
Ω = [0,   L_x]×[0,   L_y]×[0,   L_z] a rectangular cuboid of size
 L_x× L_y× L_z = 7.07  cm× 10.1  cm× 7.07 cm.
Fig. <ref> shows the geometrically detailed three-dimensional model including the valve leaflets and the aortic test section of the pulse duplicator, 
the computational mesh, and relevant boundary conditions.
The Eulerian domain is discretized using N=4 nested grid levels, with coarse grid 
spacing h_coarsest=L_x/34 ≈ 0.208  cm and refinement ratio r = 2 between levels, 
leading to h_finest≈ 0.026  cm. 
The volumetric mesh of each valve leaflet consists of hexahedral elements leading to a surface representation composed of bilinear quadrilateral elements with ≈ 2.
The test section is a stationary surface described by bilinear quadrilateral elements with =3.
The thickness of the valve leaflets is about 0.07 cm and the density is =1.8 g·cm^-3. The gap distance between the two valves is measured to be 0.0275 cm. 
The penalty spring constant associated with the valves is set to κ=5 × 10^6 g·(cm·s)^-2.
The motion of each leaflet is constrained to consist only of rotation about predefined hinge axes.
A restoring spring-damper tortional force is used to keep the right valve leaflet restricted in 
its rotation between θ_min=24^∘ and θ_min=82^∘.
The same rule is applied to keep the left valve leaflet restricted with mirrored angle. 
The test section is kept stationary by means of spring-type penalty forces with κ=1.5 × 10^6 g·cm^-2·s^-2.
A fixed time step size of Δ t=2  μs is used for the simulation.
Three-element Windkessel models establish downstream loading conditions and the upstream driving conditions
for the aortic test section.
A combination of normal traction and zero tangential velocity boundary conditions are used at the inlet and outlet to couple the 
reduced-order models to the detailed description of the flow within the test section. 
The values of the resistances and compliance for the upstream model are C_VIA=0.1  mL·mmHg^-1 and
R_1= R_2 = 0.15  mmHg·mL^-1·s.
The values at the downstream model are R_c=0.0077  mmHg·mL^-1·s ,
 R_p=1.226  mmHg·mL^-1·s, and C =1.039  mL·mmHg^-1; 
see Griffith et al. <cit.> and Lee et al. <cit.> for further discussion on the specification and parameterization of
the Windkessel models.
Solid wall boundary conditions are imposed on the remaining boundaries of the computational domain.


Fig. <ref> shows the velocity magnitudes on the plane bisecting the valves at different time points within the simulated cardiac cycle. At early times when the valves are closed
there exists a large pressure difference across the test section that forces jets of the fluid to escape from small gaps around the valves
with velocity magnitude that reaches about 450  cm·s^-1 at its peak. These hinge gap flows have been 
well characterized for bileaflet mechanical valves result from gaps between the leaflets and the housing and around the hinge areas 
in a fully closed position <cit.>.
These flows create a complex vortical interaction later in the diastole phase of the cardiac cycle, when the valve is closed and supporting a physiological pressure load.
Von Karman like vortex shedding 
clearly occurs during the systolic phase, when the valves are fully open. 
Many prior numerical simulations of mechanical heart valve models have been restricted to imposing either
experimental flow rates, or flow dynamics under prescribed leaflet motion.
Because the flow rate is not imposed in the present model, and because the time-dependent configuration of 
the valve determines the resistance of the aortic test section, this simulation demonstrates a non-trivial test of the numerical method.
 
 



 §.§ Transport of rigid blood clots in the inferior vena cava at exercise flow conditions




 
As a demonstration of the method's capability in modeling the motion of unconstrainted three-dimensional objects in a complex geometry, we simulate transport of rigid blood clots 
through the inferior vena cava (IVC) at a high flow rate of 100 cm^3/s, which corresponds to exercise flow conditions with a maximum Reynolds number of about =1500.
The IVC is a large vein that transports deoxygenated blood from lower extremities of the body back to the right atrium of the heart.
The geometry of the IVC shown in Fig. <ref> is a modified version of the patient-averaged model by Rahbar et al. <cit.> 
that has been recently used in studies of the hemodynamics <cit.>. 
  
The infrarenal IVC downstream of the iliac bifurcation has average hydraulic diameter D_h =2.8 cm. 
The density of the fluid is =1.817 g·cm^-3, and the viscosity is =5.487× 10^-2 g/cm s.
Steady fully-developed parabolic velocity boundary conditions are imposed at the upstream inlets of the iliac veins.
The two inlets are circular with diameter D=2.44 cm, which then transitions to an elliptical shape a short distance downstream.
The surface of the IVC is described using bilinear quadrilateral elements with ≈ 2. 
The IVC is embedded in a rectangular computational domain of size Ω=L_x× L_y × L_z = 50  cm× 25  cm× 50  cm.
The Eulerian domain is discretized using a three-level locally refined grid with a refinement ratio of four between the grid levels, resulting in a grid spacing of 
h_coarsest=25/16 cm=1.5625 cm on the coarsest level and h_finest=25/16× 4^2 cm≈ 0.098 cm on the finest grid level.
At the outlet, the normal traction and tangential velocities are set to zero.
Solid-wall boundary conditions are imposed along the remainder of Ω.
Once steady state condition of the flow has been reached at t=5.0 s, two neutrally buoyant rigid spherical clots (=) of diameter D_i=1 cm are 
released from positions close to the two inlets.
The volumetric mesh of each sphere consists of hexahedral elements leading to a surface mesh composed of  bilinear quadrilateral elements with =2. 
A fixed time step size of =0.1 ms is used. The penalty spring constant associated with both the IVC and the clots is 
κ=1.05 × 10^5 g·(cm·s)^-2. Fig. <ref> shows the transport of 
rigid blood clots under the unconstrained rigid-body model. A few streamlines are also plotted to show the direction of the flow.
Despite the moderately high Reynolds number of the flow and the complex pathway, the clots migrate towards the outlet through the iliac veins while 
being confined within the IVC structure. After passing the confluence of the veins, the clots come in very close contact with each other
in a more complex fluid environment where we expect the emergence of a pair of counter-rotating vortices after the confluence 
as a result of the two streams from iliac veins merging together. The clots continue migrating towards the outlet with a higher speed 
where towards the end, the clot from the right ilict vein tends to get ahead. This can be attributed to the complex nature of the flow in the region
above the confluence.




§ DISCUSSION AND CONCLUSIONS



This work has introduced a numerical approach to simulating 
fluid-structure interaction that we refer to as an immersed Lagrangian-Eulerian (ILE) method.
Results from applying this ILE method to benchmark problems of rigid body fluid-structure interaction with increasing difficulty were presented.
In addition, we also show representative results from applications of this methodology to two biomedical FSI models: 
a bileaflet mechanical heart valve 
under physiological conditions in a model of the aortic test section of a commercial pulse duplicator, and transport of rigid blood clots inside 
a patient-averaged model of the inferior vena cava under exercise conditions. 
Unlike existing partitioned methods for FSI, the ILE formulation uses an immersed approach to couple the fluid and structure subdomains and thereby reduces or even eliminates
the need for grid regeneration during dynamic simulations.
In this formulation, it is crucial to deploy a coupling approach that provides the forces from only the exterior physical fluid region, and in 
this work, we use a coupling scheme
based on the immersed interface method, which enables us to evaluate these external fluid tractions along the fluid-structure interface.
 At least in principle, however, the present method could be used with any coupling strategy that determines the net exterior fluid force acting on the 
fluid-solid interface. Notice that this excludes conventional immersed boundary formulations using regularized delta functions because 
such formulations provide the total force from both the 
the exterior (physical) and interior (nonphysical) fluid regions.
We remark that the jump conditions associated with the singular interfacial force are projected onto continuous Lagrangian basis functions, as described previously <cit.>.
This requires the solution of linear systems of equations. 
We emphasize, however, that these solves involve only surface degrees of freedom, and thus are an order of magnitude 
smaller than the volumetric fluid equations that also are solved in each time step. In addition, the linear systems involve 
discrete L^2 projection equations that can be solved optimally using simple iterative methods (e.g. by a Krylov 
method preconditioned with a diagonal or lumped mass matrix).
Overall, the dominant cost of each time step is in solving the incompressible Navier-Stokes equations, and 
the remaining computations are relatively inexpensive. A strength of the current approach is that it enables the use of 
fast Cartesian grid solvers for the incompressible Navier-Stokes equations.


The FSI coupling strategy allows fluid and solid subproblems to be solved in 
 a partitioned manner as independent, nonconforming discretizations that are coupled through interface conditions. 
In our discretization approach, there exist two Lagrangian representations of the fluid-solid interface, 
including the boundary of the volumetric mesh used in solving the equations of rigid body dynamics and a surface mesh 
that moves with the local fluid velocity.
These two representations are constrained to move together by
a Lagrange multiplier surface force. Exactly imposing the constraint would require the solution of a saddle-point system that couples the 
Eulerian and Lagrangian variables, but we develop a practical numerical scheme
that avoids the complex numerical linear algebra associated with such systems
by relaxing this constraint using a penalty formulation.
In the penalty formulation, the surface mesh moves according to the local fluid velocity but exerts force locally to the
fluid as a weak imposition of the no-slip condition.
Discrepancies in the positions of the boundary representations can be controlled by increasing the penalty parameter. In the present work, 
the maximum relative displacement is always less than 0.1 of the Cartesian grid spacing.



An attractive feature of the present ILE method is that, at least for the specific examples considered herein,
it enables the use of a simple Dirichlet-Neumann coupling scheme 
<cit.> 
without requiring either strong coupling 
or subiterations to maintain stability. In particular, the motion of the solid mesh is driven by the exterior fluid traction, and the motion 
of the solid mesh drives the motion of the fluid-structure interface representation used to impose the no-slip condition.

Although we do not theoretically determine whether the ILE formulation suffers from added mass-related instabilities, 
the computational tests reported herein suggest that it can stabily treat a broad range of density ratios, 
including structures that are less dense than the fluid, more dense than the fluid, and neutrally buoyant. This is demonstrated for multiple benchmark problems. 









For instance, in a 2-DOF model of the oscillation of a cylinder under vortex induced vibration, we obtain stable results for mass ratios up to 40 times smaller than 
smallest value reported in recent prior work  <cit.>. In the literature, added mass instabilities are typically ascribed to the treatment of 
fluid regions that become “uncovered" by the structure 
as it moves through the computational domain.
These include ALE approaches <cit.> and other sharp interface methods <cit.>.
Specialized approaches can be needed with these methods to avoid pressure fluctuations for cases involving low and near-equal density ratios.
Although the total mass and momentum of the fluid are conserved, there are local changes in the fluid mass and momentum in the regions that are both “covered” and “uncovered” by the moving structure.
Evidently, these localized changes in the fluid mass and momentum can induce temporal discontinuities in the velocity.
In the IIM approach introduced by Li and Lai <cit.>, the fluid equations are solved on the entire computational domain, including 
region occupied by the structure, similar to diffuse-interface formulations like those used in Peskin's IB method.










Using the IIM terminology, the continuity of the velocity field across the fluid-structure interface implies (̆,t)  = 0. 
A direct consequence of this
condition for the immersed interface approach is that the jump in the material derivative 
of the velocity is also zero, i.e. D/Dt(,t) =0.
This means that fluid trajectories do not cross the moving interface and 
fluid locations that are “uncovered” by the motion of the structure automatically possess 
velocities that are consistent with the equations of motion.
In contrast, for methods in which the 
fluid domain only exists on one side of the interface, the temporal jump associated with the material derivative 
of the velocity may be non-zero going from one time step to the next as a fixed Eulerian grid point may switch sides between consecutive time steps.
This temporal discontinuity in the fluid acceleration could be important as added mass instabilities 
arise when the inertial effect due to fluid forces are dominant. This could result
in a fluid pressure that is out of phase with the fluid acceleration.
 Although this has yet to be proved rigorously, we hypothesize that the apparent robustness of 
 the method to artificial added mass instabilities results from its
consistent treatment of the momentum of the fluid near the fluid-structure interface. 
Added mass-type instabilities might also be suppressed by our use of rigid-body structure models, which
 are driven by the net fluid force acting on the immersed structure rather than by pointwise forces. 
 Rigorous stability analyses and analytical investigations of the 
methodology similar to those recently developed for overset 
grid methods <cit.> may reveal stability criteria or clarify the absence of added mass-type instabilities in this formulation.
Finally, we note that although the formulation presented here assumes the use of 
a rigid body structural model, it is natural to extend this approach to immersed elastic bodies.


 


§ ACKNOWLEDGEMENTS


We acknowledge research support through NIH Awards HL117063 and HL143336, NSF Awards DMS 1664645, CBET 175193, OAC 1450327, OAC 1652541, and OAC 1931516, 
and the U.S. FDA Center for Devices and Radiological Health (CDRH) Critical Path program.
This research was supported in part by an appointment to the Research Participation Program at the U.S. FDA administered by 
the Oak Ridge Institute for Science and Education through an interagency agreement between the U.S. Department of Energy and FDA.
A.P.S.B also acknowledges support from NSF award OAC 1931368.
Computations were performed using facilities provided by University of North Carolina at Chapel Hill through the 
Research Computing division of UNC Information Technology Services and the high-performance computing clusters at the U.S. FDA.
The findings and conclusions in this article have not been formally disseminated by the FDA and should not be construed to represent
 any agency determination or policy.
The mention of commercial products, their sources, or their use in connection with material reported herein is not to be construed as either an actual or implied endorsement of such products by the Department of Health and Human Services.
We also thank Kenneth Aycock, Saad Qadeer, Jianhua Qin, and Simone Rossi for their constructive comments to improve the manuscript.




	

