A stable added-mass partitioned (AMP) algorithm for elastic solids and incompressible flow

A stable added-mass partitioned (AMP) algorithm
for elastic solids and incompressible flow

D. A. Serino111Research supported by the National Science Foundation Graduate Research Fellowship under Grant No. DGE-1744655.222Research supported by the National Science Foundation under grants DMS-1519934 and DMS-1818926. serind@rpi.edu J. W. Banks333This work was supported by contracts from the U.S. Department of Energy ASCR Applied Math Program.444Research supported by a U.S. Presidential Early Career Award for Scientists and Engineers. banksj3@rpi.edu W. D. Henshaw555This work was supported by contracts from the U.S. Department of Energy ASCR Applied Math Program.666Research supported by the National Science Foundation under grants DMS-1519934 and DMS-1818926. henshw@rpi.edu D. W. Schwendeman777This work was supported by contracts from the U.S. Department of Energy ASCR Applied Math Program.888Research supported by the National Science Foundation under grants DMS-1519934 and DMS-1818926. schwed@rpi.edu Department of Mathematical Sciences, Rensselaer Polytechnic Institute, Troy, NY 12180, USA
Abstract

A stable added-mass partitioned (AMP) algorithm is developed for fluid-structure interaction (FSI) problems involving viscous incompressible flow and compressible elastic solids. The AMP scheme is stable and second-order accurate even when added-mass, and added-damping, effects are large. Deforming composite grids are used to effectively handle the evolving geometry and large deformations. The fluid is updated with an implicit-explicit (IMEX) fractional-step scheme whereby the velocity is advanced in one step, treating the viscous terms implicitly, and the pressure is computed in a second step. The AMP interface conditions for the fluid arise from the outgoing characteristic variables in the solid and are partitioned into a Robin (mixed) interface condition for the pressure, and interface conditions for the velocity. The latter conditions include an impedance-weighted average between fluid and solid velocities using a fluid impedance of a special form. A similar impedance-weighted average is used to define interface values for the solid. The new algorithm is verified for accuracy and stability on a number of useful benchmark problems including a radial-piston problem where exact solutions for radial and azimuthal motions are found and tested. Traveling wave exact solutions are also derived and numerically verified for a solid disk surrounded by an annulus of fluid. Fluid flow in a channel past a deformable solid annulus is computed and errors are estimated from a self-convergence grid refinement study. The AMP scheme is found to be stable and second-order accurate even for very difficult cases of very light solids.

keywords:
fluid-structure interaction, moving overlapping grids, incompressible Navier-Stokes, partitioned schemes, added-mass, added-damping, elastic solids
\biboptions

sort&compress

1 Introduction

Fluid-structure interaction (FSI) problems involving incompressible fluids and elastic solids occur in many areas of engineering and applied science. Examples include modeling flow-induced vibrations of structures (i.e., aircraft wings, undersea cables, wind turbines, and bridges) and simulating blood flow in arteries and veins. These problems are typically modeled using partial differential equations for both the fluid and solid in their respective domains together with coupling conditions involving velocity and stress at the fluid-solid interface. FSI problems are challenging mathematically due, in part, to the evolving fluid and solid domains, and the coupling at the moving interface. Numerical algorithms have been developed to compute solutions of the governing equations, and these schemes can be classified broadly as monolithic or partitioned. Monolithic schemes integrate the equations implicitly as a single large system of coupled equations, while partitioned schemes use separate solvers for the fluid and solid together with some approach to couple the solutions at each time step. Partitioned schemes are often preferred for modularity and computational performance, but they can suffer from instabilities associated with the method used to couple the solutions at the interface. For example, a standard traditional partitioned (TP) scheme employs a coupling approach in which the solid provides a Dirichlet (no-slip) boundary condition for the fluid, and then the fluid supplies a Neumann (traction) boundary condition for the solid. This approach is appealing from a physical point of view in the case when the density of the solid is much greater than that of the fluid, but it is well known that TP schemes suffer from added-mass instabilities in the case when the solid is light relative to the fluid. A method of iteration of the coupling conditions can be included in a TP algorithm in effort to suppress added-mass instabilities, but the effectiveness of this sub-time-step iteration is limited and it increases the cost of the algorithm.

In recent work fib2014 (); fis2014 (), we developed a new class of Added-Mass Partitioned (AMP) algorithms for FSI problems coupling incompressible flow and elastic solids. The work in fib2014 () for bulk solids and in fis2014 () for thin structures both considered Stokes flow and linear elastic solids, with the fluid-solid interface linearized about a flat surface. The algorithms used a fractional-step approach for the fluid in which the velocity is advanced in one stage followed by the solution of a Poisson problem for the pressure. A key ingredient in the algorithms is the AMP interface conditions, which includes a mixed Robin condition for the pressure. For the case of a bulk solid, this condition includes the acceleration of the solid through a consideration of the outgoing characteristic variables of the hyperbolic system for the elastic solid and the interface conditions. The velocity and stress at the interface are set by impedence-weighted averages so that the AMP conditions behave correctly in the limits of heavy and light solids. Thus, the resulting AMP algorithm is stable for any ratio of the density of the solid to that of the fluid, without sub-time-step iterations, effectively suppressing added-mass instabilities. In beamins2016 (), the approach for thin structures was extended to incompressible flows with nonlinear convection and Euler-Bernoulli beams with finite deflection, and it was found that the AMP scheme remains stable for both heavy and light solids.

The aim of the present work is to develop a stable and second-order accurate AMP algorithm for FSI problems coupling incompressible flow and linear elastic bulk solids. The scheme follows the work in fib2014 (), but with significant modifications to handle finite deformation of the fluid-solid interface, nonlinear advection in the fluid, an IMEX-type scheme for the fluid, and viscous added-damping effects. The Navier-Stokes equations for the fluid are solved in either an Eulerian frame or an Arbitrary-Eulerian-Lagrangian (ALE) frame which is coupled to the motion of the interface, while the equations governing the solid are solved in a static Lagrangian reference frame. Both the fluid and solid domains are covered by composite overlapping grids fsi2012 (), with the deformation of the fluid grid performed using the approach described originally in mog2006 (). For example, Figure 1 shows a composite grid for an FSI problem involving a fluid in a channel flowing past three elastic solids in the shape of the letters “R”, “P” and “I”. The viscous fluid flows from left to right, and it interacts with the three solids whose computed displacements determine the deformed shapes. The image on the far-right of the figure shows an enlarged view of the composite grid near the fluid-solid interface of the right-most deformed solid.

solid domain

fluid domain

interface

Figure 1: Top left: Initial grid is composed of a (blue) background grid and (green) deforming interface-fitted grids representing the fluid domain. The ‘RPI’ shaped solid domain is composed of (red) background grids, (purple) interface-fitted grids, and (pink) annular grids around interior cores where the displacement is set to zero. Bottom left: Solution at described by streamlines in the fluid domain and shaded contours of the magnitude of the solid displacement. Right: Enlarged view of the grid after solid deformation.

The equations of linear elasticity in the solid domain are advanced explicitly using the time-stepping scheme described in smog2012 () for the hyperbolic system of equations in first-order form. The Navier-Stokes equations in the fluid domain are solved using a fractional-step scheme. Unlike the earlier work fib2014 (), the present solver uses an implicit-explicit (IMEX) scheme based on splitStep2003 (). In this scheme, the fluid velocity is advanced first with the viscous terms treated implicitly while the pressure gradient and nonlinear convection terms handled explicitly. The fluid pressure is updated in a subsequent step by solving a Poisson problem. The incorporation of the IMEX scheme has necesitated important changes to the original AMP interface conditions introduced in fib2014 () to handle viscous added-damping effects that arise when the viscous CFL number becomes large. Modifications to the AMP interface conditions are guided by the analysis of carefully chosen FSI model problems, as discussed in the companion manuscript fibrmparXiv (). In particular, the analysis reveals a stable definition of the fluid impedance, which is a critical ingredient in the present AMP algorithm. The focus of the present work is a description of this algorithm for the full FSI problem, and its implementation in two space dimensions using deforming composite grids (DCG). Performance of the algorithm is examined for a set of FSI benchmark problems, including ones in which new exact solutions are derived. These results are used to demonstrate the stability and accuracy of the algorithm, and a comparison to results from a traditional partitioned scheme is also provided.

The overarching technique used to derive Added-Mass Partitioned schemes relies on the use of so-called compatibility boundary conditions, which are derived from the governing PDEs in concert with the interface conditions. This approach has also been used in a variety of other FSI regimes to yield stable and accurate partitioned solvers without the need for sub-iteration. For example in the case of FSI problems involving viscous incompressible flow, AMP schemes have been developed for thin elastic structures beamins2016 () (as noted above) and for rigid solids rbinsmp2017 (); rbins2017 (); rbins3d2018 (). The first AMP schemes were developed for compressible flow banks11a (); sjogreen12 (), and were subsequently extended into a general 2D framework using the DCG approach for FSI problems involving inviscid compressible flows and linear elastic solids fsi2012 (). For the case of inviscid compressible flow, added-mass effects are localized due to finite propagation speeds of disturbances in the fluid vanBrummelen2009 (), and this leads to a somewhat simpler treatment of the coupling conditions in the corresponding AMP scheme relative to the incompressible case considered here. Extensions of the scheme in fsi2012 () to FSI problems coupling inviscid compressible flow and rigid solids is described in lrb2013 (), with the case of nonlinear elastic solids detailed in flunsi2016 ().

In addition to the work cited above, there are other approaches to addressing added-mass-type instabilities in both partitioned and monolithic FSI solvers. For partitioned schemes, a typical strategy is to used under-relaxed subiterations, and previous research has shown that the number of necessary sub-iterations can be reduced by the use of Aitken acceleration or quasi-Newton methods, see Küttler2008 (); MEHL2016869 (). Additionally, sub-iteration schemes based on Robin-Neumann or Robin-Robin coupling, as opposed to the traditional Dirichlet-Neumann coupling, have also been shown to yield performance gains Wang2018 (); BASTING2017312 (); BadiaNobileVergara2008 (); MokWallRamm2001 (); FernandezMullaertVidrascu2014 (); FernandezMullaertVidrascu2013 (); FernandezLandajuela2014 (); Gerardo_GiordaNobileVergara2010 (); BadiaNobileVergara2009 (); NobilePozzoliVergara2014 (). In related research, added-mass effects can also be treated by adding fictitious mass terms in the structure equations BaekKarniadakis2012 (); YuBaekKarniadakis2013 () or artificial compressibility in the fluid equations RabackRuokolainenLyly2001 (); DegrooteSwillensBruggemanHaeltermanSegersVierendeels2010 (); Degroote2011 (). For difficult problems with large added-mass effects, monolithic schemes have been used to eliminate the need for sub-iterations Forti2017 (); FigueroaVignonClementelJansenHughesTaylor2006 (), and the performance of monolithic schemes has improved significantly through the application of multigrid AULISA2018213 (); Heil2004 (); GeeKuttlerWall2011 (); HronTurekMonolithic2006 ().

The remaining sections of the paper are organized as follows. The equations governing the FSI problem are described in Section 2, and the AMP algorithm is described in Section 3. In the latter section, we begin with a discussion of the AMP interface conditions at a continuous level, which is followed by a detailed discussion of the AMP scheme. Section 4 provides a description of the spatial approximations using deforming composite grids (DCG). Numerical results confirming the stability and accuracy of the scheme and discussed in Section 5. Some of the results use new exact solutions of benchmark FSI problems, which are described in the appendices. Conclusions are given in Section 6.

2 Governing equations and interface conditions

We consider the coupled evolution of an incompressible fluid and a linear elastic solid. The fluid occupies the domain , where is a vector of physical coordinates and is time. The equations for the solid are written in terms of the Lagrangian coordinate for a reference configuration at . (An overbar is used here and elsewhere to denote quantities associated with the solid.) The fluid and solid are coupled at an interface described by in physical space and in the corresponding reference space. Figure 1 illustrates these geometric features using a model fluid-structure interaction problem.

The fluid is described by the incompressible Navier-Stokes equations

(1a)
(1b)

where is the velocity, is the stress tensor and is the (constant) density. The fluid stress tensor is given by

(2)

where is the pressure, is the identity matrix and is the viscous stress tensor. The latter quantity is given in terms of the gradient of the fluid velocity by

(3)

where is the dynamic viscosity of the fluid (assumed to be constant).

We consider the governing equations for the fluid in velocity-pressure form. In this form, the continuity equation in (1a) is replaced by a Poisson equation for the pressure given by

(4)

where

(5)

Here, components of the fluid velocity are denoted by , and a similar notation is used for the components of tensors, e.g. the components of are given by with , 2, 3. In the velocity-pressure form of the equations, boundary conditions are required for the the pressure-Poisson equation. A suitable choice that ensures consistency with solutions of the velocity-divergence form is for .

The evolution of the solid is described by the equations of linear elasticity. These equations can be written in the form

(6a)
(6b)

where is the displacement of the solid, its velocity, and its density (assumed constant). The Cauchy stress tensor in (6b) is defined by

(7)

where and are Lamé parameters (taken to be constants). The position of the solid in physical space is determined by the mapping

(8)

for . The corresponding deformation gradient given by is assumed to be a small perturbation of the identity for a linear elastic solid so that a one-to-one mapping from the reference space to the physical space of the solid exists. Following the approach described in smog2012 (), the solid equations are treated as a first order system of PDEs in time and space. This formulation is achieved by taking the time derivative of (7) to obtain

(9)

which is then appended to previous system of equations in (6). In this form, upwind solvers can be used to advance displacement, velocity and stress of the solid. We note that (7) is enforced at , and then the stress-strain relation is satisfied approximately for through compatibility conditions applied at the boundary and a stress-strain relaxation term in the numerical scheme (see flunsi2016 ()).

The interface between the fluid and solid, in physical space given by , is determined by the mapping in (8) for the boundary of the solid given by . The interface is assumed to be smooth so that a well-defined normal to the interface exists. Along the interface, the following matching conditions hold:

(10a)
(10b)

where is the outward unit normal to the fluid domain, i.e.  points from the fluid domain to the solid domain. Suitable boundary conditions are applied on the boundaries of the fluid and solid domians not included in , and initial conditions on , and are set to close the problem.

Since the equations governing the solid form a hyperbolic system, it is informative to consider the propagation of characteristics. The primary and secondary wave speeds are defined as

(11)

which correspond to compression and sheer waves, respectively. At the interface, the Riemann variables corresponding to incoming and outgoing characteristics can be obtained from the first-order equations in (6b) and (9) projected onto the normal to the interface. These variables for the compression and shear waves are given by

(12a)
(12b)
corresponding to information propagating into the solid domain, and
(12c)
(12d)

corresponding to information propagating out of the solid domain. Here, and are mutually orthogonal unit vectors tangent to , and thus orthogonal to , and and are the solid impedances for compression and shear waves, respectively. The functions given in (12), corresponding to the Riemann variables for the solid, can be used to derive an equivalent set of characteristic-based interface conditions which take into account the direction of propagation. These conditions are

(13a)
(13b)

which are linear combinations of the matching conditions in (10).

This paper is concerned with partitioned algorithms for solving FSI problems involving incompressible fluids and linear elastic solids. Such algorithms employ PDE solvers for the equations governing the fluid and solid separately with the matching conditions at the fluid-solid interface providing boundary conditions for each domain solver (along with conditions on the boundaries not included in the shared interface). For example, in traditional partitioned algorithms, the fluid velocity is determined by the solid, while the solid traction is determined by the fluid. This ordering can often be used for the case of heavy solids, but it suffers from added-mass instabilities when the solid density is light or the solid grid is very fine. In the reverse scheme, referred to as the anti-traditional scheme, the solid imposes a traction condition on the fluid, while the fluid velocity provides a condition for the solid velocity. This reverse ordering works well for light solids or fine grids, but may exhibit instability for heavy solids or coarse grids.

The added-mass partitioned (AMP) scheme described in the next section follows the work in fib2014 (), and uses the characteristic-based matching conditions in (13). At a high level, these are implemented as boundary conditions for the solid solver as

(14a)
and as boundary conditions for the fluid solver as
(14b)

where the superscript denotes suitable predicted values (as described in more detail below) so that the right-hand sides of the boundary conditions in (14) are treated as known. This characteristic-based decomposition of the interface matching conditions is a key step in the derivation of the AMP algorithm, and it leads to a stable scheme for a wide range of fluid and solid densities. In particular, in the limits of heavy solids and light solids, the AMP scheme approaches the traditional and anti-traditional partitioned schemes, respectively.

3 FSI algorithms

In this section, we describe the AMP algorithm, and we also discuss the traditional partitioned (TP) algorithm and the anti-traditional partitioned (ATP) algorithm for comparative purposes. All three partitioned schemes use the same fluid and solid domain solvers. For the fluid domain, we integrate the governing equations in velocity-pressure form using a fractional-step scheme in which the fluid velocity is advanaced in one step and the pressure-Poisson equation is solved to update the fluid pressure in the other step (see splitStep2003 ()). An explicit upwind Godunov method is used to advance the velocity and stress of the solid in its reference domain following the approach described in smog2012 (), while a Lax-Wendroff scheme is used to advance the solid displacement. These two domain solvers, both second-order accurate in space and time, are embedded in the FSI algorithms using a predictor-corrector approach. The first step of the FSI time-stepping loop involves advancing the solid variables to the next time level. The solid displacement predicts the new position of the fluid-solid interface and this information is used to deform the fluid grid and compute its grid velocity. The fluid-domain solver is then called to advance the fluid velocity to the next time level and update the fluid pressure. In the fluid solver, the diffusion terms of the momentum equations are treated implicitly, while the convection and pressure-gradient terms are approximated explicitly. Following this fluid-domain step, the solid interface data is updated using fluid data at the new time level. The steps outlined above are essentially repeated in the corrector step. Depending on the densities of the fluid and solid (and on the grid spacing), the TP and ATP schemes may require addition sub-iterations for stability, in which case the corrector step is repeated.

The AMP algorithm differs from that of the TP and ATP schemes in the implementation of the fluid-solid interface conditions. To this end, we begin in Section 3.1 by describing the various treatments of the interface conditions, including compatability conditions, at a continuous level. This discussion is then followed in Section 3.2 by a detailed description of the various FSI algorithms at a discrete level.

3.1 Continuous interface conditions

Using the formulas for the fluid stress tensor in (2) and (3), the AMP interface conditions in (14b) for the fluid become

(15a)
(15b)

with the superscipt on the solid variables suppressed. These conditions, along with applied on the boundary, are a sufficient set of conditions for the fluid equations in velocity-pressure form. The left-hand side of the interface conditions in (15b) only involve the fluid velocity (and its derivatives), and these conditions with the divergence-free constraint are used as boundary conditions for the momentum equation in (1b) to advance the fluid velocity. The stress and velocity of the solid appearing on the right-hand side of (15b) are obtained from an explicit time step of the hyperbolic system of equations for the solid. Once the fluid velocity is advanced to the next time step, a Poisson problem is solved to update the fluid pressure. Following the analysis in fib2014 (); fibrmparXiv (), the interface condition in (15a) is used with the momentum equation in (1b) to derive a Robin condition for the pressure-Poisson problem that takes the form

(16)

where is the normal derivative and is the kinematic viscosity of the fluid. Following splitStep2003 (), we have used the identity, , noting that , to replace on the right-hand side of (16) in favor of the curl-curl operator. This is done for improved stability of the fractional-step scheme. The condition in (16), along with suitable conditions for , is used for the Poisson equation in (4) for the pressure.

The AMP algorithm requires an additional interface projection to ensure the fluid and solid interface values match at the end of the time-step. This projection defines common interface values for the velocity and traction using an impedance-weighted average. In particular, the normal and tangential components of the velocity on the interface are given by

(17a)
(17b)

while the traction is defined from an inverse impedance average,

(18a)
(18b)

Since the fluid is incompressible there is apparently no natural definition for the fluid impedance . In our previous paper fib2014 (), we used , where measured the height of the fluid domain, and the corresponding scheme was found to be very insensitive to the choice of . This previous scaling was sufficient when using explicit time-stepping in the fluid. In the current paper, however, we use an IMEX scheme in the fluid where viscous effects are more important since the viscous CFL parameter may become very large. In this regime, the previous choice of does not result in a stable scheme. The analysis of a model problem that accounts for viscous effects performed in the companion paper fibrmparXiv (), provides the suitable form for which is a key result that leads to a stable scheme.

AMP Fluid Impedence.

The fluid impedance accounting for both added-mass and viscous added-damping effects is chosen as

(19)

where is the fluid density, measures the fluid grid spacing in the normal direction to the interface, is the time-step, is the dynamic viscosity of the fluid, and and are constants following the analysis in fibrmparXiv ().

In agreement with fibrmparXiv (), the choice (19) is found to yield a stable AMP algorithm for all FSI problems considered, including those covering a wide range of solid densities from very heavy to light. On the other hand, stable AMP algorithms also result when scaling the constants and by a factor of or . This indicates that there is some robustness in the choice of the constants in (19). However, it appears to be important to keep both terms in (19), since failing to do so can lead to subtle instabilities at late times.

We note that with appropriate scales for length and time from the analysis given by  and , respectively, there is a simple scaling argument leading to the form for in (19). Consider a dimensional analysis of the fluid momentum equation in one dimension,

(20)

In addition, let and be reference scales for velocity and pressure (stress), respectively. In the inertia-dominated limit, the two terms on the left-hand side of (20) set the scale for pressure leading to the balance

The implied impedance, , is expected to be dominant for added-mass effects. In the viscous-dominated limit when added-damping effects are important, the pressure gradient term in (20) balances the viscous stress term leading to

This balance suggests a second expression for the fluid impedance given by . A linear combination of the limiting forms agrees with the choice for in (19).

The impedance-weighted averages (17) and (18) are defined to ensure that the interface conditions behave correctly in the limits of a heavy and light solids. In the limit of a heavy solid, i.e.  (, ), the interface conditions for the fluid should approach those of a (moving) no-slip wall. For , (15b) reduces to

(21)

so that the tangential components of the fluid velocity match those of the solid wall. We note that the tangential components of the interface velocity defined in (17b) agree with those given in (21) in this limit. The normal component of the fluid is given by the normal component of the interface velocity, defined in (17a), which becomes

(22)

so that the normal component of the velocity of the fluid also matches that of the solid wall. In this heavy-solid limit, the divergence-free constraint applied on the boundary becomes a compatibility condition used to specify values of the normal component of the fluid velocity at ghost points. The Robin condition for the pressure given in (16) in the limit reduces to

(23)

This Neumann condition for the pressure is often used for the pressure-Poisson problem at a no-slip wall.

In the limit of a light solid, i.e.  (, ), the interface conditions for the fluid should approach those of a free surface. With , (15b) becomes

(24)

and thus the tangential components of the fluid stress match those of the light solid. In this limit, the conditions in (24) and the divergence-free constraint on the fluid velocity provide the boundary conditions needed to advance the velocity using (1b), along with suitable boundary conditions for which we do not discuss. The interface condition in (15a) or the derived Robin condition in (16) both reduce to

(25)

in the limit , and this Dirichlet condition on pressure is used for the pressure-Poisson problem. Essentially, the solid traction sets the fluid traction along the interface in the light-solid limit, and the definitions of the interface traction in (18) only confirms this assignment. Finally, the impedance-weighted averages in (17) reduce to

(26)

so that the velocity of the interface is specified by the fluid velocity, as expected, which then prescribes the velocity of the solid as a boundary condition.

3.2 Discrete partitioned algorithms

In this section, the details of the AMP algorithm are presented, along with the special cases of the TP and ATP schemes for comparison. The description of the AMP algorithm here extends the framework of the original scheme given in fib2014 () to moving grids and deforming interfaces. A semi-implicit predictor-corrector scheme is used to solve the fluid equations following the work in splitStep2003 () and extending the explicit scheme in fib2014 (), while a second-order accurate upwind scheme is used to advance the solid equations. The evolving fluid domain is discretized on a moving overlapping grid, which is discussed in more detail in Section 4. The solid domain is discretized using a static overlapping grid in its reference configuration.

Algorithm 1 outlines the steps for the AMP scheme. Let for denote the time-dependent grid points belonging to the fluid domain (including grid points on the boundaries and the interface), and let denote the velocity of the grid points at the time step, . The indices of grid points on the interface of the fluid grid are denoted by . Similarly, let for denote the grid points belonging to the static overlapping grid for the solid (including boundaries and the interface), while are the indicies associated with grid points on the interface of the solid grid. Discrete approximations for the fluid variables are given by and , while the solid variables are given by , and .

Algorithm 1 Added-mass partitioned (AMP) scheme

In the description of the algorithm, we assume that all quantities belonging to the fluid and solid are known at , and the task is to advance the solution to . Algorithm 1 gives a summary of all of the stages in the AMP scheme. The first five stages of the algorithm form the predictor step.

Stage 1 (predict solid): Predicted values for the solid displacement are obtained by integrating (6a) using the explicit Lax-Wendroff-type scheme

(27)

Predicted values for the solid velocity and stress, contained in the vector , are obtained using a second-order accurate Godunov upwind scheme of the form

(28)

where represents the upwind fluxes corresponding to the coordinate direction. In practice, we add a stress-relaxation term to the right-hand side of (28) following the approach in flunsi2016 () so that the stress-strain relation in (7) is satisfied approximately throughout the solid domain.

Stage 2 (predict fluid grid motion): The computed solid displacement, , is used to predict the position, for , of the fluid-solid interface at . The overlapping-grid generator is called to move the fluid grid to match the predicted interface position, and to determine the grid velocity .

Stage 3 (predict fluid velocity): The fluid velocity at is predicted using an explicit Adams-Bashforth scheme for the convection and pressure-gradient terms and an implicit Crank-Nicholson scheme for the viscous terms. The scheme takes the form

(29)

where and represent grid operators associated with the explicit and implicit terms in (29), respectively, given by

(30)

Since is determined implicitly in (29), boundary conditions are required for . Standard conditions are used for the portion of the boundary away from the interface, while care must be used in the application of the AMP interface conditions for the velocity as described in Section 3.1 so that the correct limiting behaviors are obtained. On the interface, the characteristic condition in (15b) and the divergence-free constraint give

(31a)
(31b)

where the interface normal, , and tangent vectors, , are evaluated at grids points along the predicted interface at . Since (29) excludes the grids points on the interface, another set of equations are required to complete the system for the implicit time step. We use the following impedance-weighted averages to set the velocity on the interface:

(32a)
(32b)

where is defined by

(33)

Note that (32), with in (33), forms a set of equations for , , which couples to the evolution equations in (29) and the boundary conditions in (31). Together, this linear system is solved to obtain the predicted fluid velocity .

The TP scheme sets the fluid velocity on the interface to equal that of the solid, so that the interface conditions in (32) are replaced by

(34)

The divergence-free constraint in (31b) is used to specify the normal component of the fluid velocity in the ghost points, while the tangential components of the velocity in the ghost points are extrapolated. The ATP scheme, on the other hand, applies a tangential stress boundary condition on the fluid along with the divergence-free constraint. These conditions are equivalent to the conditions in (31), but with the terms in (31a) involving the fluid and solid velocity dropped (effectively setting ), and then replacing the interface conditions in (32) with for . We note that these limiting cases for TP and ATP schemes correspond to the limits of the continuous conditions discussed above in Section 3.1.

Stage 4 (predict fluid pressure): Predicted values for the fluid pressure are obtained by solving the discrete Poisson problem

(35)

with the predicted fluid velocity from the previous stage used to determine the right-hand side of (35). We have included a divergence damping term to the right-hand side of (35) with coefficient  to suppress the growth of the divergence of the velocity in the numerical solution that may occur due to discretization errors. The discrete approximation of the Robin condition in (16) is

(36)

where the acceleration of the solid along the interface appearing on the right-hand side of (36) is obtained from the solid velocity using a second-order accurate finite difference approximation in time.

In the TP scheme, the Robin condition in (36) is replaced by a discrete approximation of (23) given by

whereas the ATP scheme uses an approximation of (25) given by

As before, the both cases are consistent with the limiting cases of a heavy (TP) and light (ATP) solid.

Stage 5 (predict solid interface): The projections given by (17) and (18) are used in this stage to define values for the solid velocity and traction on the interface () as

(37a)
(37b)

and

(38a)
(38b)

respectively. These projected values are used to specify and on the interface, and then the corresponding values in the ghost points are set using extrapolation.

In the TP scheme, we use to set the predicted solid traction on the interface, but then use the compatibility condition

(39)

which can be interpreted as a set of discrete Neumann conditions for the components of . The time-derivative of the solid stress on the right-hand side of (39) is obtained using a second-order accurate approximation involving and the solid stress at and . For the ATP scheme, we set the solid velocity on the interface using , and use the compability condition

(40)

which can be interpreted as a set of discrete Neumann conditions for the components of . The acceleration of the solid on the interface is computed using a finite-difference approximation similar to that described to compute the stress for the TP scheme.

This completes the set of stages used for the predictor step of Algorithm 1. The remaining four stages form the corrector step. Stages 7 and 8 in the corrector step are similar in form to the corresponding two stages in the predictor step. In Stage 9 the solid velocity and traction are set equal to the current fluid values. This ensures that the interface jump conditions are exactly matched at the end of the correction step. The boundary conditions for the TP and ATP scheme in the corrector stages are handled in an analogous fashion to the predictor stages.

4 Spatial approximations

The domain for the FSI problem is discretized using an overlapping grid , which is defined by a set of component grids , , that cover both the fluid and solid domains. Static Cartesian background grids are used for the bulk of each domain, while boundary and interface-fitted grids are used to resolve curved boundaries and moving fluid-solid interfaces (see Figure 2). In the fluid, each component grid is defined by a mapping from the physical space to a unit computational space (a square in two dimensions or a box in three dimensions). The mapping is time dependent in general to account for the moving fluid-solid interface, and is given by

(41)

where is the number of dimensions. Each component grid for the solid is defined by a mapping from the (fixed) reference space to a unit computational space, and is given by

(42)

The corresponding physical space for the solid is then determined by the mapping in (8) using the computed displacement . Within the fluid and solid domains, the component grids may overlap and interpolation formulas are used to match the discrete solutions across the component grids. The position of the fluid-solid interface is determined by the interface-fitted component grids of the solid, and then the interface-fitted fluid grids deform to match the moving interface position in physical space. A hyperbolic grid generator HyperbolicGuide () is used to regenerate the fluid interface-fitted grids, and then the Ogen grid generator OGEN () is called to regenerate the overlapping grid connectivity between the (moving) interface-fitted grids of the fluid and the (static) background grids.

interpolation

ghost

unused

Figure 2: Left: composite of a background grid (, blue) and a boundary-fitted grid (, green) in physical space for the domain defined by the interior of the red boundary. The grid points on with green dots interpolate from and the grid points on with blue dots interpolate from Middle: Plot of showing interpolation points, ghost points (grid points which exist outside the physical boundary), and unused points (grid points which do not affect the computation). Right: The green boundary fitted grid, is mapped to a unit square. The plot shows interpolation points and ghost points.

The Navier-Stokes equations for the fluid in velocity-pressure form are transformed to the unit-square coordinates, , using standard chain-rule differentiation formulas. On moving grids, the equations are transformed to a moving coordinate system which introduces the grid velocity, , into the advection terms as indicated in (29) and (30). The resulting equations are then discretized using standard finite-difference approximations for the derivatives with respect to following the approach in ICNS (); splitStep2003 (); mog2006 (). The equations of linear elasticity for the solid are also transformed to the unit-square coordinates and then integrated in time using a (Godunov-based) upwind scheme similar to smog2012 (); fsi2012 (); flunsi2016 (). The approximations of the equations governing the fluid and solid are second-order accurate in the present implementation of the AMP time-stepping scheme.

The time step, , is taken as the minimum value between the stable time steps for the fluid and solid domain solvers independently. There is no restriction on the time step imposed by the AMP interface conditions. For the fluid solver, the time step is generally determined by a CFL-type stability constraint based on the advection terms in the fluid momentum equations. The viscous terms are integrated implicitly so that there is no stability constraint on arising from these terms. The time step for the solid solver is determined based on a CFL-type stability constraint involving the compression wave speed, but there also may be an adjustment imposed by the stress-relaxation term (see smog2012 (); flunsi2016 ()).

5 Numerical results

We now present numerical results to verify the stability and accuracy of the AMP scheme for a variety of FSI problems. We start with fundamental problems for which exact solutions are known. The first problem involves an annular fluid region surrounding an elastic solid disk with the motion of the solid constrained to the radial direction normal to the fluid-solid interface. The second problem uses the same geometry, but assumes circumferential motion of the solid. Exact solutions exist for both problems, as discussed in A and B, while the results are given in Sections 5.1.1 and 5.1.2 below. We next test the AMP scheme for a more general FSI problem involving a radial traveling wave. Results for this problem are given in Section 5.2 with the corresponding exact solution for a linearized traveling wave derived in C. For each of these test cases, we examine the convergence of the numerical solution using the exact solution for different fluid-solid density ratios to confirm that the AMP scheme is both stable and second-order accurate. The remaining two tests involve more complex FSI problems for which exact solutions are unavailable. The first problem, discussed in Section 5.3, considers a pressure-driven flow in a channel about an elastic solid. A self-convergence study is performed for this clean benchmark problem which further verifies the accuracy and stability of the AMP scheme. We also examine the behavior of the TP scheme for this problem. The final test problem involves the flow past multiple solids with a wide range of densities. This test is discussed in Section 5.4 and is designed to illustrate the robustness of the AMP scheme.

5.1 Radial elastic piston

Let us begin by considering two benchmark FSI problems for which exact solutions are available. The geometry of the two problems, shown in Figure 3, is similar. A solid elastic disk with initial radius is surrounded by an incompressible fluid which occupies the annular region . In the first problem, it is assumed that the motion of the solid and fluid is confined to the radial direction, while circumferential motion is assumed in the second problem.

solid:

interface:

fluid:

solid:

interface:

fluid:

Figure 3: The geometry for the radial elastic piston solution. The left figure shows the initial domain without deformation and the right figure shows the deformed domain for

5.1.1 Radial motion of a circularly symmetric elastic piston

For the case of radial motion, the solid disk behaves as an elastic piston that drives a radial flow in the fluid surrounding it (as illustrated in Figure 3). The exact solution for this problem is derived in A. It is convenient to define the solution in terms of a specified radial component of the solid displacement, , having the form

(43)

where is the frequency of the oscillation, determines the amplitude, and is the Bessel function of the first kind of order one. The corresponding time-dependent interface radius, , is given by