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

TL;DR

This paper introduces a stable and accurate partitioned algorithm for fluid-structure interaction problems involving incompressible flow and elastic solids, effectively handling large deformations with composite grids.

Contribution

A novel AMP algorithm with impedance-weighted interface conditions for stable, second-order accurate FSI simulations involving viscous incompressible flow and elastic solids.

Findings

Algorithm is stable and second-order accurate.

Verified on benchmark problems with exact solutions.

Effective for very light solids.

Abstract

A stable added-mass partitioned (AMP) algorithm is developed for fluid-structure interaction (FSI) problems involving viscous incompressible flow and compressible elastic solids. 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.

Tables

Table 1. Table 1: Selected solutions for the frequency ω 𝜔 \omega of the dispersion relation in ( 66 ) for δ = 10 − 3 𝛿 superscript 10 3 \delta=10^{-3} , 1 1 1 and 10 3 superscript 10 3 10^{3} . The corresponding values for b / u ¯ 0 𝑏 subscript ¯ 𝑢 0 b/\bar{u}_{0} and b ¯ / u ¯ 0 ¯ 𝑏 subscript ¯ 𝑢 0 \bar{b}/\bar{u}_{0} are also given. Additional parameters chosen are r 0 = 1 2 subscript 𝑟 0 1 2 r_{0}=\frac{1}{2} , R = 1 𝑅 1 R=1 , ρ = 1 𝜌 1 \rho=1 , ν = 0.1 𝜈 0.1 \nu=0.1 and ρ ¯ = μ ¯ = δ ¯ 𝜌 ¯ 𝜇 𝛿 \bar{\rho}=\bar{\mu}=\delta .

$\delta$	$\omega$	$b/{\overline{u}}_0$	$\overline{b}/{\overline{u}}_0$
${10}^{-3}$	$7.664+0.001497i$	$4.938-4.327i$	$-1523+2316i$
$1$	$8.778+0.7854i$	$4.355-6.155i$	$-3.512+1.913i$
${10}^3$	$10.27+0.002055i$	$2.216-5.864i$	$-2.944+0.0005914i$

Table 2. Table 2: Max-norm self-convergence results for the elastic annulus in a channel at t = 1.5 𝑡 1.5 t=1.5 with density ratio δ = 10 𝛿 10 \delta=10 .

Elastic annulus in a channel, $\delta =10$
$h_j$	$E_j^{\left(p\right)}$	$r$	$E_j^{\left(𝐯\right)}$	$r$	$E_j^{\left(\overline{𝐮}\right)}$	$r$	$E_j^{\left(\overline{𝐯}\right)}$	$r$	$E_j^{\left(\overline{𝝈}\right)}$	$r$
1/20	6.3e-3		7.4e-3		1.7e-4		4.1e-4		1.3e-2
1/40	1.7e-3	3.7	1.5e-3	5.0	4.6e-5	3.7	9.2e-5	4.5	3.7e-3	3.6
1/80	4.7e-4	3.6	3.0e-4	4.9	1.2e-5	3.7	2.0e-5	4.5	1.0e-3	3.6
rate	1.86		2.30		1.90		2.18		1.86

Table 3. Table 3: Stability of the traditional-partitioned scheme for the elastic annulus in a channel as a function of the density ratio δ 𝛿 \delta and the grid-spacing. Even for very heavy solids the TP scheme will become unstable when the grid becomes sufficiently fine.

$h$	$\delta =100$	$\delta =200$	$\delta =300$	$\delta =400$	$\delta =500$	$\delta =800$	$\delta =1000$
1/20	unstable	unstable	stable	stable	stable	stable	stable
1/40	unstable	unstable	unstable	unstable	stable	stable	stable
1/80	unstable	unstable	unstable	unstable	unstable	unstable	stable

Table 4. Table 4: Parameters for the solid objects in the channel as shown in Figure  13 .

Solid	$W$	$H$	$w$	$h$	$d$	$(x_c,y_c)$	$r_c$	$(x_0,y_0)$	${\theta }_0$	$\delta$
${\overline{\Omega }}_1$	$1.5$	$0.25$	$-$	$-$	$-$	$(0.25,0)$	$0.05$	$(-1.5,0.85)$	$-{10}^{\circ }$	2
${\overline{\Omega }}_2$	$1.5$	$0.35$	$-$	$-$	$-$	$(-0.25,0)$	$0.05$	$(-1.5,-1)$	${10}^{\circ }$	0.1
${\overline{\Omega }}_3$	$1.5$	$1.5$	$0.30$	$0.30$	$-$	$(0,0)$	$0.1$	$(0,0)$	${20}^{\circ }$	1
${\overline{\Omega }}_4$	$1.0$	$1.0$	$0.35$	$0.25$	$-$	$(0,-0.30)$	$0.05$	$(1.5,0.75)$	${45}^{\circ }$	10
${\overline{\Omega }}_5$	$1.0$	$1.2$	$0.25$	$0.25$	$0.2$	$(-0.35,-0.35)$	$0.05$	$(1.5,-0.85)$	$-{160}^{\circ }$	1000

Equations

$$\nabla \cdot \mathbf{v}$$

$$\rho {\mathbf{v}}_t+\rho (\mathbf{v}\cdot \nabla )\mathbf{v}$$

$$\mathbf{\sigma }$$

$$\mathbf{\tau }$$

$$\Delta p$$

$$\nabla \mathbf{v}:{\left(\nabla \mathbf{v}\right)}^T\equiv \sum _{i=1}^3\sum _{j=1}^3\frac{\partial v_i}{\partial x_j}\frac{\partial v_j}{\partial x_i}.$$

$${\bar{\mathbf{u}}}_t$$

$$\bar{\rho }{\bar{\mathbf{v}}}_t$$

$$\bar{\mathbf{\sigma }}$$

$$\mathbf{x}=\bar{\mathbf{x}}+\bar{\mathbf{u}}(\bar{\mathbf{x}},t),$$

$${\bar{\mathbf{\sigma }}}_t$$

$$\mathbf{v}$$

$$\mathbf{\sigma }\mathbf{n}$$

$${\bar{c}}_p=\sqrt{\frac{\bar{\lambda }+2\bar{\mu }}{\bar{\rho }}},{\bar{c}}_s=\sqrt{\frac{\bar{\mu }}{\bar{\rho }}},$$

$$\mathcal{A}(\bar{\mathbf{\sigma }},\bar{\mathbf{v}})$$

$${\mathcal{A}}_m(\bar{\mathbf{\sigma }},\bar{\mathbf{v}})$$

$$\mathcal{B}(\bar{\mathbf{\sigma }},\bar{\mathbf{v}})$$

$${\mathcal{B}}_m(\bar{\mathbf{\sigma }},\bar{\mathbf{v}})$$

$$\mathcal{A}(\bar{\mathbf{\sigma }},\bar{\mathbf{v}})$$

$$\mathcal{B}(\mathbf{\sigma },\mathbf{v})$$

$$\mathcal{A}(\bar{\mathbf{\sigma }},\bar{\mathbf{v}})=\mathcal{A}({\mathbf{\sigma }}^{(p)},{\mathbf{v}}^{(p)}),{\mathcal{A}}_m(\bar{\mathbf{\sigma }},\bar{\mathbf{v}})={\mathcal{A}}_m({\mathbf{\sigma }}^{(p)},{\mathbf{v}}^{(p)}),\mathbf{x}\in \Gamma (t),$$

$$\mathcal{B}(\mathbf{\sigma },\mathbf{v})=\mathcal{B}({\bar{\mathbf{\sigma }}}^{(p)},{\bar{\mathbf{v}}}^{(p)}),{\mathcal{B}}_m(\mathbf{\sigma },\mathbf{v})={\mathcal{B}}_m({\bar{\mathbf{\sigma }}}^{(p)},{\bar{\mathbf{v}}}^{(p)}),\bar{\mathbf{x}}\in \Gamma (t),$$

$$-p+{\mathbf{n}}^T\mathbf{\tau }\mathbf{n}+{\bar{z}}_p{\mathbf{n}}^T\mathbf{v}={\mathbf{n}}^T\bar{\mathbf{\sigma }}\mathbf{n}+{\bar{z}}_p{\mathbf{n}}^T\bar{\mathbf{v}},\mathbf{x}\in \Gamma (t),$$

$${\mathbf{t}}_m^T\mathbf{\tau }\mathbf{n}+{\bar{z}}_s{\mathbf{t}}_m^T\mathbf{v}={\mathbf{t}}_m^T\bar{\mathbf{\sigma }}\mathbf{n}+{\bar{z}}_s{\mathbf{t}}_m^T\bar{\mathbf{v}},m=1,2,\mathbf{x}\in \Gamma (t),$$

$$\displaystyle-p-{\bar{z}_{p}\Delta t\over\rho}\partial_{n}p=\mathbf{n}^{T}(\bar{\boldsymbol{\sigma}}\mathbf{n}-\boldsymbol{\tau}\mathbf{n})+\bar{z}_{p}\Delta t\mathbf{n}^{T}\bigl{(}\bar{\mathbf{v}}_{t}+\nu\nabla\times\nabla\times\mathbf{v}\bigr{)},\qquad\mathbf{x}\in\Gamma(t),$$

$${\mathbf{n}}^T{\mathbf{v}}^I=$$

$${\mathbf{t}}_m^T{\mathbf{v}}^I=$$

$${\mathbf{n}}^T{\mathbf{\sigma }}^I\mathbf{n}=$$

$${\mathbf{t}}_m^T{\mathbf{\sigma }}^I\mathbf{n}=$$

$$z_f$$

$$\rho v_t+p_x=\mu v_{xx}.$$

$$P=\left(\frac{\rho h}{\Delta t}\right)V.$$

$$P=\left(\frac{\mu }{h}\right)V.$$

$${\mathbf{t}}_m^T\mathbf{v}={\mathbf{t}}_m^T\bar{\mathbf{v}},m=1,2,\mathbf{x}\in \Gamma (t),$$

$${\mathbf{n}}^T\mathbf{v}={\mathbf{n}}^T{\mathbf{v}}^I={\mathbf{n}}^T\bar{\mathbf{v}},\mathbf{x}\in \Gamma (t),$$

$$\displaystyle-{1\over\rho}\partial_{n}p=\mathbf{n}^{T}\bar{\mathbf{v}}_{t}+\nu\mathbf{n}^{T}\bigl{(}\nabla\times\nabla\times\mathbf{v}\bigr{)},\qquad\mathbf{x}\in\Gamma(t).$$

Full text

1. 1 Introduction

2. 2 Governing equations and interface conditions

3. 3 FSI algorithms

4. 3.1 Continuous interface conditions

5. 3.2 Discrete partitioned algorithms

6. 4 Spatial approximations

7. 5 Numerical results

8. 5.1 Radial elastic piston

9. 5.1.1 Radial motion of a circularly symmetric elastic piston

10. 5.1.2 Rotating elastic disk in a fluid annulus

11. 5.2 Radial traveling wave solution

12. 5.3 Elastic annulus in a pressure-driven channel flow

13. 5.4 Multiple elastic solids in a pressure-driven channel flow

14. 6 Conclusions

15. A Radial piston exact solution

16. B Exact solution for the rotating elastic disk in a fluid annulus

17. C Radial traveling wave solution

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 [1, 2], we developed a new class of Added-Mass Partitioned (AMP) algorithms for FSI problems coupling incompressible flow and elastic solids. The work in [1] for bulk solids and in [2] 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 [3], 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 [1], 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 [4], with the deformation of the fluid grid performed using the approach described originally in [5]. 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.

The equations of linear elasticity in the solid domain are advanced explicitly using the time-stepping scheme described in [6] 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 [1], the present solver uses an implicit-explicit (IMEX) scheme based on [7]. 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 [1] 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 [8]. 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 [3] (as noted above) and for rigid solids [9, 10, 11]. The first AMP schemes were developed for compressible flow [12, 13], and were subsequently extended into a general 2D framework using the DCG approach for FSI problems involving inviscid compressible flows and linear elastic solids [4]. For the case of inviscid compressible flow, added-mass effects are localized due to finite propagation speeds of disturbances in the fluid [14], 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 [4] to FSI problems coupling inviscid compressible flow and rigid solids is described in [15], with the case of nonlinear elastic solids detailed in [16].

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 [17, 18]. 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 [19, 20, 21, 22, 23, 24, 25, 26, 27, 28]. In related research, added-mass effects can also be treated by adding fictitious mass terms in the structure equations [29, 30] or artificial compressibility in the fluid equations [31, 32, 33]. For difficult problems with large added-mass effects, monolithic schemes have been used to eliminate the need for sub-iterations [34, 35], and the performance of monolithic schemes has improved significantly through the application of multigrid [36, 37, 38, 39].

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 $\mathbf{x}\in\Omega(t)$, where $\mathbf{x}=(x_{1},x_{2},x_{3})$ is a vector of physical coordinates and $t$ is time. The equations for the solid are written in terms of the Lagrangian coordinate $\bar{\mathbf{x}}=(\bar{x}_{1},\bar{x}_{2},\bar{x}_{3})$ for a reference configuration $\bar{\mathbf{x}}\in{\bar{\Omega}}_{0}$ at $t=0$. (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 $\mathbf{x}\in\Gamma(t)$ in physical space and $\bar{\mathbf{x}}\in{\bar{\Gamma}}_{0}$ 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

[TABLE]

where $\mathbf{v}(\mathbf{x},t)$ is the velocity, $\boldsymbol{\sigma}(\mathbf{x},t)$ is the stress tensor and $\rho$ is the (constant) density. The fluid stress tensor is given by

[TABLE]

where $p(\mathbf{x},t)$ is the pressure, $\mathbf{I}$ is the identity matrix and $\boldsymbol{\tau}(\mathbf{x},t)$ is the viscous stress tensor. The latter quantity is given in terms of the gradient of the fluid velocity by

[TABLE]

where $\mu$ 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

[TABLE]

where

[TABLE]

Here, components of the fluid velocity $\mathbf{v}$ are denoted by $(v_{1},v_{2},v_{3})$, and a similar notation is used for the components of tensors, e.g. the components of $\boldsymbol{\sigma}$ are given by $\sigma_{mn}$ with $m,n=1$, 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 $\nabla\cdot\mathbf{v}=0$ for $\mathbf{x}\in\partial\Omega(t)$.

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

[TABLE]

where $\bar{\mathbf{u}}(\bar{\mathbf{x}},t)$ is the displacement of the solid, $\bar{\mathbf{v}}(\bar{\mathbf{x}},t)$ its velocity, and $\bar{\rho}$ its density (assumed constant). The Cauchy stress tensor $\bar{\boldsymbol{\sigma}}(\bar{\mathbf{x}},t)$ in (6b) is defined by

[TABLE]

where $\bar{\lambda}$ and $\bar{\mu}$ are Lamé parameters (taken to be constants). The position of the solid in physical space is determined by the mapping

[TABLE]

for $\bar{\mathbf{x}}\in{\bar{\Omega}}_{0}$. The corresponding deformation gradient given by $\bar{F}=\nabla_{\bar{\mathbf{x}}}\mathbf{x}=\mathbf{I}+\nabla_{\bar{\mathbf{x}}}\bar{\mathbf{u}}$ 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 [6], 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

[TABLE]

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 $t=0$, and then the stress-strain relation is satisfied approximately for $t>0$ through compatibility conditions applied at the boundary and a stress-strain relaxation term in the numerical scheme (see [16]).

The interface between the fluid and solid, in physical space given by $\mathbf{x}\in\Gamma(t)$, is determined by the mapping in (8) for the boundary of the solid given by $\bar{\mathbf{x}}\in{\bar{\Gamma}}_{0}$. 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:

[TABLE]

where $\mathbf{n}(\mathbf{x},t)$ is the outward unit normal to the fluid domain, i.e. $\mathbf{n}$ 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 $\Gamma(t)$, and initial conditions on $\mathbf{v}$, $\bar{\mathbf{u}}$ and $\bar{\mathbf{v}}$ 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

[TABLE]

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

[TABLE]

corresponding to information propagating out of the solid domain. Here, $\mathbf{t}_{1}(\mathbf{x},t)$ and $\mathbf{t}_{2}(\mathbf{x},t)$ are mutually orthogonal unit vectors tangent to $\Gamma(t)$, and thus orthogonal to $\mathbf{n}$, and $\bar{z}_{p}=\bar{\rho}\bar{c}_{p}$ and $\bar{z}_{s}=\bar{\rho}\bar{c}_{s}$ 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

[TABLE]

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 [1], and uses the characteristic-based matching conditions in (13). At a high level, these are implemented as boundary conditions for the solid solver as

[TABLE]

where the $(p)$ 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 [7]). 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 [6], 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

[TABLE]

with the $(p)$ superscipt on the solid variables suppressed. These conditions, along with $\nabla\cdot\mathbf{v}=0$ 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 [1, 8], 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

[TABLE]

where $\partial_{n}=\mathbf{n}\cdot\nabla$ is the normal derivative and $\nu=\mu/\rho$ is the kinematic viscosity of the fluid. Following [7], we have used the identity, $\Delta\mathbf{v}=-\nabla\times\nabla\times\mathbf{v}$, noting that $\nabla\cdot\mathbf{v}=0$, to replace $\Delta\mathbf{v}$ 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 $\mathbf{x}\in\partial\Omega(t)\backslash\Gamma(t)$, 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 $\mathbf{x}\in\Gamma(t)$ are given by

[TABLE]

while the traction is defined from an inverse impedance average,

[TABLE]

Since the fluid is incompressible there is apparently no natural definition for the fluid impedance $z_{f}$. In our previous paper [1], we used $z_{f}=\rho H/\Delta t$, where $H$ measured the height of the fluid domain, and the corresponding scheme was found to be very insensitive to the choice of $H$. 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 $\nu\Delta t/h^{2}$ may become very large. In this regime, the previous choice of $z_{f}$ does not result in a stable scheme. The analysis of a model problem that accounts for viscous effects performed in the companion paper [8], provides the suitable form for $z_{f}$ which is a key result that leads to a stable scheme.

AMP Fluid Impedence**.**

*The fluid impedance $z_{f}$ accounting for both added-mass and viscous added-damping effects is chosen as *

[TABLE]

where $\rho$ is the fluid density, $h$ measures the fluid grid spacing in the normal direction to the interface, $\Delta t$ is the time-step, $\mu$ is the dynamic viscosity of the fluid, and ${\mathcal{C}}_{\rm AM}=1$ and ${\mathcal{C}}_{\rm AD}=2$ are constants following the analysis in [8].

In agreement with [8], 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 ${\mathcal{C}}_{\rm AD}$ and ${\mathcal{C}}_{\rm AM}$ by a factor of $1/2$ or $2$. 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 $h$ and $\Delta t$, respectively, there is a simple scaling argument leading to the form for $z_{f}$ in (19). Consider a dimensional analysis of the fluid momentum equation in one dimension,

[TABLE]

In addition, let $V$ and $P$ 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

[TABLE]

The implied impedance, $(z_{f})_{{\rm AM}}=\rho h/\Delta t$, 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

[TABLE]

This balance suggests a second expression for the fluid impedance given by $(z_{f})_{{\rm AD}}=\mu/h$. A linear combination of the limiting forms agrees with the choice for $z_{f}$ 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. $\bar{\rho}\rightarrow\infty$ ($\bar{z}_{p}$, $\bar{z}_{s}\rightarrow\infty$), the interface conditions for the fluid should approach those of a (moving) no-slip wall. For $\bar{z}_{s}\rightarrow\infty$, (15b) reduces to

[TABLE]

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

[TABLE]

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 $\bar{z}_{p}\rightarrow\infty$ reduces to

[TABLE]

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. $\bar{\rho}\rightarrow 0$ ($\bar{z}_{p}$, $\bar{z}_{s}\rightarrow 0$), the interface conditions for the fluid should approach those of a free surface. With $\bar{z}_{s}\rightarrow 0$, (15b) becomes

[TABLE]

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 $\mathbf{x}\in\partial\Omega(t)\backslash\Gamma(t)$ which we do not discuss. The interface condition in (15a) or the derived Robin condition in (16) both reduce to

[TABLE]

in the limit $\bar{z}_{p}\rightarrow 0$, 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

[TABLE]

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 [1] to moving grids and deforming interfaces. A semi-implicit predictor-corrector scheme is used to solve the fluid equations following the work in [7] and extending the explicit scheme in [1], 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 $\mathbf{x}_{\mathbf{i}}^{n}$ for $\mathbf{i}\in\Omega_{h}$ denote the time-dependent grid points belonging to the fluid domain (including grid points on the boundaries and the interface), and let $\dot{\mathbf{x}}_{\mathbf{i}}^{n}$ denote the velocity of the grid points at the $n^{{\rm th}}$ time step, $t^{n}$. The indices of grid points on the interface of the fluid grid are denoted by $\mathbf{i}\in\Gamma_{h}$. Similarly, let $\bar{\mathbf{x}}_{\bar{\mathbf{i}}}$ for ${\bar{\mathbf{i}}}\in{\bar{\Omega}}_{h}$ denote the grid points belonging to the static overlapping grid for the solid (including boundaries and the interface), while ${\bar{\mathbf{i}}}\in{\bar{\Gamma}}_{h}$ are the indicies associated with grid points on the interface of the solid grid. Discrete approximations for the fluid variables are given by $\mathbf{v}_{\mathbf{i}}^{n}\approx\mathbf{v}(\mathbf{x}_{\mathbf{i}}^{n},t^{n})$ and $p_{\mathbf{i}}^{n}\approx p(\mathbf{x}_{\mathbf{i}}^{n},t^{n})$, while the solid variables are given by $\bar{\mathbf{u}}_{\bar{\mathbf{i}}}^{n}\approx\bar{\mathbf{u}}(\bar{\mathbf{x}}_{\bar{\mathbf{i}}},t^{n})$, $\bar{\mathbf{v}}_{\bar{\mathbf{i}}}^{n}\approx\bar{\mathbf{v}}(\bar{\mathbf{x}}_{\bar{\mathbf{i}}},t^{n})$ and $\bar{\boldsymbol{\sigma}}_{\bar{\mathbf{i}}}^{n}\approx\bar{\boldsymbol{\sigma}}(\bar{\mathbf{x}}_{\bar{\mathbf{i}}},t^{n})$.

In the description of the algorithm, we assume that all quantities belonging to the fluid and solid are known at $t=t^{n}$, and the task is to advance the solution to $t^{n+1}=t^{n}+\Delta t$. 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

[TABLE]

Predicted values for the solid velocity and stress, contained in the vector $\bar{\mathbf{q}}_{\bar{\mathbf{i}}}=(\bar{\mathbf{v}}_{\bar{\mathbf{i}}},\bar{\boldsymbol{\sigma}}_{\bar{\mathbf{i}}})$, are obtained using a second-order accurate Godunov upwind scheme of the form

[TABLE]

where $\bar{\mathbf{F}}_{m,\;{\bar{\mathbf{i}}}}^{\pm}$ represents the upwind fluxes corresponding to the $\bar{x}_{m}$ coordinate direction. In practice, we add a stress-relaxation term to the right-hand side of (28) following the approach in [16] 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, $\bar{\mathbf{u}}_{\bar{\mathbf{i}}}^{(p)}$, is used to predict the position, $\mathbf{x}_{\mathbf{i}}^{(p)}$ for $\mathbf{i}\in\Gamma_{h}$, of the fluid-solid interface at $t^{n+1}$. The overlapping-grid generator is called to move the fluid grid to match the predicted interface position, and to determine the grid velocity $\dot{\mathbf{x}}_{\mathbf{i}}^{(p)}$.

Stage 3 (predict fluid velocity): The fluid velocity at $t^{n+1}$ 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

[TABLE]

where $\mathbf{N}_{h}$ and $\mathbf{L}_{h}$ represent grid operators associated with the explicit and implicit terms in (29), respectively, given by

[TABLE]

Since $\mathbf{v}_{\mathbf{i}}^{(p)}$ is determined implicitly in (29), boundary conditions are required for $\mathbf{i}\in\partial\Omega_{h}$. 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

[TABLE]

where the interface normal, $\mathbf{n}$, and tangent vectors, $\mathbf{t}_{m}$, are evaluated at grids points $\mathbf{i}\in\Gamma_{h}$ along the predicted interface at $t^{n+1}$. 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:

[TABLE]

where $\mathbf{V}^{p}_{h}$ is defined by

[TABLE]

Note that (32), with $\mathbf{V}^{p}_{h}$ in (33), forms a set of equations for $\mathbf{v}_{\mathbf{i}}^{(p)}$, $\mathbf{i}\in\Gamma_{h}$, 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 $\mathbf{v}_{\mathbf{i}}^{(p)}$.

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

[TABLE]

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 $\bar{z}_{s}=0$), and then replacing the interface conditions in (32) with $\mathbf{v}_{\mathbf{i}}^{(p)}=\mathbf{V}^{p}_{h}(\mathbf{v}_{\mathbf{i}}^{(p)})$ for $\mathbf{i}\in\Gamma_{h}$. 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

[TABLE]

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 $\alpha_{\mathbf{i}}$ 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

[TABLE]

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

[TABLE]

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

[TABLE]

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 ${\bar{\mathbf{i}}}\in{\bar{\Gamma}}_{h}$ ($\mathbf{i}\in\Gamma_{h}$) as

[TABLE]

and

[TABLE]

respectively. These projected values are used to specify $\bar{\mathbf{v}}_{\bar{\mathbf{i}}}^{(p)}$ and $\bar{\boldsymbol{\sigma}}^{(p)}_{{\bar{\mathbf{i}}}}\mathbf{n}$ on the interface, and then the corresponding values in the ghost points are set using extrapolation.

In the TP scheme, we use $\bar{\boldsymbol{\sigma}}^{(p)}_{{\bar{\mathbf{i}}}}\mathbf{n}={\boldsymbol{\sigma}}^{(p)}_{\mathbf{i}}\mathbf{n}$ to set the predicted solid traction on the interface, but then use the compatibility condition

[TABLE]

which can be interpreted as a set of discrete Neumann conditions for the components of $\bar{\mathbf{v}}_{\bar{\mathbf{i}}}^{(p)}$. The time-derivative of the solid stress on the right-hand side of (39) is obtained using a second-order accurate approximation involving $\bar{\boldsymbol{\sigma}}_{{\bar{\mathbf{i}}}}^{(p)}$ and the solid stress at $t^{n}$ and $t^{n-1}$. For the ATP scheme, we set the solid velocity on the interface using $\bar{\mathbf{v}}_{\bar{\mathbf{i}}}^{(p)}=\mathbf{v}_{\mathbf{i}}^{(p)}$, and use the compability condition

[TABLE]

which can be interpreted as a set of discrete Neumann conditions for the components of $\bar{\boldsymbol{\sigma}}_{\bar{\mathbf{i}}}^{(p)}\mathbf{n}$. 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 ${\mathcal{G}}$, which is defined by a set of component grids $\{G_{g}\}$, $g=1,\,\ldots,\,{\mathcal{N}}$, 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 $\mathbf{x}$ to a unit computational space $\mathbf{r}$ (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

[TABLE]

where ${n_{d}}$ is the number of dimensions. Each component grid for the solid is defined by a mapping from the (fixed) reference space $\bar{\mathbf{x}}$ to a unit computational space, and is given by

[TABLE]

The corresponding physical space $\mathbf{x}$ for the solid is then determined by the mapping in (8) using the computed displacement $\bar{\mathbf{u}}(\bar{\mathbf{x}},t)$. 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 [40] is used to regenerate the fluid interface-fitted grids, and then the Ogen grid generator [41] is called to regenerate the overlapping grid connectivity between the (moving) interface-fitted grids of the fluid and the (static) background grids.

The Navier-Stokes equations for the fluid in velocity-pressure form are transformed to the unit-square coordinates, $\mathbf{r}$, using standard chain-rule differentiation formulas. On moving grids, the equations are transformed to a moving coordinate system which introduces the grid velocity, $\dot{\mathbf{x}}$, 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 $\mathbf{r}$ following the approach in [42, 7, 5]. 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 [6, 4, 16]. 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, $\Delta t$, 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 $\Delta t$ 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 [6, 16]).

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 $\bar{r}=r_{I}(0)$ is surrounded by an incompressible fluid which occupies the annular region $r_{I}(t)<r<R$. 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.

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, $\bar{u}_{r}(\bar{r},t)$, having the form

[TABLE]

where $\omega$ is the frequency of the oscillation, $\beta$ determines the amplitude, and $J_{1}(z)$ is the Bessel function of the first kind of order one. The corresponding time-dependent interface radius, $r_{I}(t)$, is given by

[TABLE]

where $r_{0}=r_{I}(0)$ is the initial interface radius. The radial component of the fluid velocity and its pressure are given by

[TABLE]

where

[TABLE]

Here, $V(t)$ and $P(t)$ are the radial velocity and pressure of the fluid at $r=R$, respectively, and $\dot{V}(t)$ is the radial acceleration. These functions are specified in the exact solution based on the choice for the interface motion given in (44). We take $\beta=0.05$ to be the specified amplitude of the motion and $\omega=\pi$ as the specified frequency. Additionally, we choose $r_{0}=0.5$, $R=1$, $\rho=1$, and $\bar{\rho}=\bar{\mu}=\bar{\lambda}=\delta$, where $\delta$ is a parameter that determines the ratio of the density of the solid to that of the fluid.

Numerical solutions are computed using a composite grid, denoted by ${\mathcal{G}}_{\rm rep}^{(j)}$, and shown in Figure 4. The grid ${\mathcal{G}}_{\rm rep}^{(j)}$ has a target grid spacing of $\Delta s=1/(10j)$ for a chosen resolution parameter $j$. The solid disk is represented by an inner Cartesian (red) grid and an interface-fitted annular (pink) grid. These two grids are fixed in the reference domian, ${\bar{\Omega}}_{0}$, of the solid. The fluid domain is covered by an interface-fitted annular (green) grid that deforms over time with the radial motion of the interface, and an outer annular (blue) grid that is fixed to the outer boundary of the fluid.

Figure 4 also shows the behavior of the solution at times $t=0.3$ and $0.9$ computed using the AMP algorithm and the grid ${\mathcal{G}}_{\rm rep}^{(8)}$. The density ratio for this solution is $\delta=1$. We note that while the viscous terms in the fluid momentum equation vanish identically in the derivation of the exact solution, the numerical scheme requires a choice for the visosity which is taken to be $\nu=0.05$. The interface between the fluid and solid expands and contracts over time. At $t=0.3$ the interface is near its largest radius, while at $t=0.9$ it is near its smallest.

Grid convergence studies of the max-norm errors are given in the left column of plots of Figure 5 for density ratios of $\delta=10^{3}$, $\delta=1$ and $\delta=10^{-3}$. The computations are stable for all values of $\delta$ considered, and we observe that the max-norm errors converge at close to second-order accuracy as indicated by the (black) reference line in each plot.

5.1.2 Rotating elastic disk in a fluid annulus

The second FSI benchmark problem has a geometry similar to that shown in Figure 3, but instead of the radial motion shown we consider a circumferential motion of the disk. The exact solution, derived in B, for the circumferential components of the fluid velocity and solid displacement takes the form

[TABLE]

where

[TABLE]

Here,

[TABLE]

and $Y_{1}(z)$ is the Bessel function of the second kind of order one. The parameter $\bar{u}_{0}$ determines the amplitude of the rotational motion. Solutions for the circumferential components of the fluid velocity and the solid displacement are taken from the real parts of $v_{\theta}(r,t)$ and $\bar{u}_{\theta}(r,t)$ defined in (45). The pressure is then determined by

[TABLE]

The frequency $\omega$ is obtained from the solution of the dispersion relation in (66). For our computations, we choose $r_{0}=\frac{1}{2}$, $R=1$, $\rho=1$, $\nu=0.1$, $\bar{\rho}=\bar{\mu}=\delta$ and $\bar{u}_{0}=10^{-5}$. Selected values of $\omega$ and the corresponding values for $b/\bar{u}_{0}$ and $\bar{b}/\bar{u}_{0}$ are given in Table 1 for different values of $\delta$.

Numerical solutions are computed using the AMP algorithm and the composite grid ${\mathcal{G}}_{\rm rep}^{(j)}$ as described for the previous problem. Figure 6 shows the composite grid, and the computed solution at $t=0.1$ and $0.3$ for the case $\delta=1$. The Lamé parameter $\bar{\lambda}$ does not appear in the exact solution for this problem, but a value is required for the numerical scheme and we set $\bar{\lambda}=\delta$ in the numerical calculations.

Grid convergence studies of the max-norm errors in the right column of plots of Figure 5 for density ratios of $\delta=10^{3}$ , $\delta=1$ and $\delta=10^{-3}$. The computations are stable, as before, and the max-norm errors all converge at close to second-order accuracy.

5.2 Radial traveling wave solution

We next consider traveling wave solutions involving both radial and circumferential motion of an incompressible fluid in a circular region $0<r<r_{0}$ coupled to a surrounding linearly elastic solid in the annular region $r_{0}<r<R$. The geometry of the FSI problem is shown in Figure 7. An exact solution is derived for this problem in C assuming a Stokes fluid and an linearized interface about $r=r_{0}$. In this solution, the fluid velocity and pressure take the form

[TABLE]

where $\omega$ is a frequency, $n$ is a positive integer, and $\hat{v}_{r}(r)$, $\hat{v}_{\theta}(r)$ and $\hat{p}(r)$ are coefficient functions. Similarly, the displacement in the solid is given by

[TABLE]

with coefficient functions $\hat{\bar{u}}_{r}(r)$ and $\hat{\bar{u}}_{\theta}(r)$. Following the analysis in C, the coefficient functions are found to be

[TABLE]

and

[TABLE]

where $\lambda^{2}=i\omega/\nu$, $k_{s}=\omega/\bar{c}_{s}$, $k_{p}=\omega/\bar{c}_{p}$, $(p_{I},d,\bar{d}_{1},\bar{d}_{2},\bar{d}_{3},\bar{d}_{4})$ are constants of integration, and $J_{n}$ and $Y_{n}$ are Bessel functions of order $n$. An application of traction-free boundary conditions at $r=R$ and matching conditions involving velocity and stress at the interface $r=r_{0}$ leads to a homogeneous linear system of equations for the six constants. Nontrivial solutions are found if the system is singular, which leads to the dispersion relation in (87) involving the frequency $\omega$ and circumferential wave number $n$.

Numerical solutions are computed for the case $r_{0}=1$ and $R=1.2$. The fluid has $\rho=1$ and $\nu=0.1$, while the solid has $\bar{\mu}=\bar{\lambda}=\bar{\rho}=\delta$. The domain for the FSI problem is represented by the composite grid, ${\mathcal{G}}_{\rm rtw}^{(j)}$, with target grid spacing $\Delta s=1/(10j)$ as shown in Figure 7. The circular fluid domain is covered by a fixed background Cartesian grid (blue) and interface-fitted grid (green). The solid domain is represented by an annular grid (red). While the exact solution assumes a linearization about a fixed interface position, $r=r_{0}$, the numerical solution does not make this assumption so that the interface-fitted grid of the fluid moves in time with the (very small) deformation of the interface. Also, the numerical solution does not assume a Stokes fluid. However, since the velocity of the fluid is very small for this problem, the nonlinear convective terms in the momentum equation are negligible. Initial conditions for the AMP time-stepping scheme are taken from the exact (linearized) solution at $t=0$ for an $(n,\omega)$-pair satisfying the dispersion relation in (87) and the corresponding set of constants $(p_{I},d,\bar{d}_{1},\bar{d}_{2},\bar{d}_{3},\bar{d}_{4})$ with a normalization taken to be ${\bar{u}}_{0}=10^{-7}$. The exact solution is described in terms of complex-valued functions, but we use the real part to set the initial conditions and to compare with the numerical solution.

Figure 7 shows representative solutions at $t=0.1$ and $0.3$ computed using the grid ${\mathcal{G}}_{\rm rtw}^{(8)}$. The density ratio is chosen to be $\delta=1$, and a $(n,\omega)$-pair satisfying the dispersion relation is given by $n=3$ and $\omega=3.491-1.154\,i$. The shaded contours of the fluid pressure and the norm of the solid displacement both show three-fold symmetry in the circumferential direction in agreement with the choice of $n=3$. Grid convergence studies of the max-norm errors are given in Figure 8 for density ratios of $\delta=10$, $\delta=1$ and $\delta=10^{-1}$. The computations are stable and the max-norm errors all converge at close to second-order accuracy.

5.3 Elastic annulus in a pressure-driven channel flow

In this section, we consider the flow of an incompressible fluid in a channel past an embedded annular elastic solid which is anchored by a fixed inner boundary. The fluid with density $\rho=1$ and viscosity $\nu=0.05$ exists in the rectangular region $[-3,3]\times[-1.75,1.75]$. The solid annulus has center $\mathbf{x}=(0,0)$ with outer radius $R_{d}=0.8$ and inner radius $R_{i}=0.4$. The solid parameters are $\bar{\rho}=\bar{\lambda}=\bar{\mu}=\delta$, where the density ratio $\delta$ is a parameter defining light and heavy solids as before. Initially the fluid is at rest and the displacement of the solid is zero. A flow in the channel is driven by a pressure gradient determined by a variable pressure assigned at the inflow boundary on the left. To define a smooth solution for testing convergence, the pressure at the left inflow boundary is smoothly ramped from zero to a final value of $p_{0}=1$ at $t=1$. Specifically, we set $p=p_{0}{\eta}(t)$ at $x=-3$, where ${\eta}(t)$ is a ramp function given by

[TABLE]

The ramp function satisfies ${\eta}={\eta}^{\prime}={\eta}^{\prime\prime}={\eta}^{\prime\prime\prime}=0$ at $t=0$, and it has three continuous derivatives at $t=1$. We also set the tangential component of the velocity to be zero at the left inflow boundary. No-slip boundary conditions are applied at the top and bottom walls of the channel, while the pressure and normal derivatives of the velocity are set to zero at the right outflow boundary. The usual matching conditions involving velocity and stress are taken at the fluid-solid interface, and the inner solid boundary is taken as a zero displacement condition. While an exact solution is not available for this problem, the set-up is very clean so that it makes a good benchmark problem for testing various FSI algorithms as is done below.

The composite grid used for this problem is denoted by ${\mathcal{G}}_{d}^{(j)}$ with target grid spacing equal to $\Delta s=1/(10j)$. An enlarged view of the central portion of the grid is shown in Figure 9. The solid is represented by a single annular (red) grid. The fluid domain is covered by a background Cartesian (blue) grid and an annular (green) grid of radial width equal to $0.2$ initially. This latter grid is attached to the interface and deforms over time.

Figures 10 and 11 describe the evolution of the fluid-solid interface over several snapshots in time. Figure 10 shows the shape of the interface in physical space at different times by indicating the displacement of a set of marker points on the interface from their reference positions at $t=0$ to their current positions at the times given. Figure 11, on the other hand, shows the radial and circumferential components of the interface displacement as a function of the polar angle, $\theta$, for the same set of times. The angle $\theta=0$ corresponds to the trailing edge of the solid and $\theta$ increases in the counter-clockwise direction. At time $t=0$ the interface is circular. The pressure at the inflow boundary ramps up from $p=0$ at $t=0$ to $p=1$ at $t=1$, and by $t=2$ the interface has deformed with the left side being pushed to the right. Following the initial deformation of the solid to the right by the drag of the fluid, the solid undergoes elastic oscillations in time about this initially deformed shape. The oscillations in the solid are damped over time due to energy transfer to the fluid and the viscous dissipation that occurs there. We observe that the outer boundary of the solid achieves its maximum circumferential displacement (in magnitude) at approximately $\theta=90^{\circ}$ and $\theta=270^{\circ}$, while its maximum radial displacement occurs at the leading edge $\theta=180^{\circ}$.

To estimate the accuracy of the computed solution, a self-convergence test is performed. Solutions are computed on the composite grids ${\mathcal{G}}_{d}^{(j)}$, $j=2$, 4 and 8, and the errors and convergence rates are estimated using a Richardson extrapolation procedure following the approach described in [43]. Table 2 presents the max-norm error results for this self-convergence test. It is seen that the max-norm errors are converging at close to second-order.

Numerical solutions are also computed for this FSI problem using the traditional partitioned (TP) scheme. Table 3 shows the stability of the computed solution as a function of the density ratio $\delta=\bar{\rho}/\rho$ and the target grid spacing $h$ for the composite grid ${\mathcal{G}}_{d}^{(j)}$. We observe that the TP scheme is unstable on even the coarsest grid used for $\delta=200$ (and smaller) indicating that this value is already a very light solid case for the TP scheme. Larger values of $\delta$ are stable on the coarsest grid, but the results indicate that the TP scheme becomes unstable for a fixed value of $\delta$ as the grid is refined. This result is in agreement with the stability results discussed in [1] and [8].

5.4 Multiple elastic solids in a pressure-driven channel flow

To test the robustness of the AMP algorithm, we consider a final example of an incompressible fluid in a channel flowing past five embedded solids with different shapes and densities. The fluid has density $\rho=1$ and viscosity $\nu=0.02$, and it occupies the rectangular region $[-3.5,5]\times[-2,2]$. The form of the boundary conditions at the inlet and outlet of the channel, and along its top and bottom walls, are the same as that described for the pressure-driven flow in Section 5.3. For the present problem, the pressure is ramped smoothly from $p=0$ at $t=0$ to $p=p_{0}=0.1$ at $t=0.5$ by prescribing $p(-3,y,t)=p_{0}\eta(2t)$ for $-2<y<2$ with the ramp function $\eta$ given by (49).

Figure 12 describes the baseline reference geometries for the elastic solids ${\bar{\Omega}}_{k}$, $k=1,\ldots,5$, in the problem. These solids are positioned in the channel as shown in Figure 13. The geometric parameters that describe each solid are given in Table 4, along with the position $(x_{0},y_{0})$ of the geometric center of each solid in the channel and the angle $\theta_{0}$ of the solid (in degrees) relative to the baseline configurations shown in Figure 12. Each solid is anchored to the channel at a small circular hole of radius $r_{c}$ and center $(x_{c},y_{c})$ where the displacement is set to zero. We set $\bar{\rho}=\bar{\lambda}=\bar{\mu}=\delta$ for each solid with the different values for $\delta$ also listed in the table. The solid densities, which range from very light, $\delta=0.1$, to very heavy, $\delta=1000$, provide a significant test for a partitioned scheme. Also note that even large values for $\delta$ become difficult for the traditional scheme as the grids are refined as noted from Table 3.

The composite grid, denoted by $\mathcal{G}_{\rm mb}^{(j)}$ with target grid spacing $\Delta s=1/(10j)$, covers the fluid domain of the channel and reference domains of the five solids as shown in Figure 13. The fluid domain is composed of a (blue) background grid and a (green) deforming interface-fitted grid about each of the five solids. The (static) reference domain of each solid is composed of a (red) background grid, a (purple) interface-fitted grid, and a (pink) annular grid about the circular cut-out. An enlarged view of the composite grid showing more detail of the overlapping grid structure is also shown in the figure.

Figure 14 shows computed results at two times, $t=8$ and $t=10$, with a comparison between two different grid resolutions, a ‘coarse’ grid $\mathcal{G}_{\rm mb}^{(4)}$ (left column) and a ‘fine’ grid $\mathcal{G}_{\rm mb}^{(8)}$ (right column). The plots in the figure show (instantaneous) streamlines in the fluid domain and shaded contours of the magnitude of displacement in each solid domain. The shaded contours are colored on different scales for each solid as noted in the figure caption. The streamlines are colored by the magnitude of the flow speed. The solids are seen to divert the flow with the fastest flow found above the top of the I-shaped solid $\bar{\Omega}_{4}$. The flow is time-dependent and the solids undergo oscillations over time as elastic waves propagate through the bodies. The lighter solids are seen to deform more than the heavier ones. The rectangular beam $\bar{\Omega}_{2}$ at the bottom near the inlet of the channel, for example, is very light and is easily deformed by the fluid. As a result, it also undergoes a fairly large compression and expansion. Three arms of the plus-shaped solid $\bar{\Omega}_{3}$, are seen to be significantly pushed by the fluid while the top of the I-shaped solid is deformed to the right. Note that since the deformation of the bodies is time-dependent, the matching fluid velocity can be non-zero on the fluid-solid interfaces and the fluid streamlines appear to enter or exit the bodies at some locations. The results from the coarse and fine grid are very similar indicating a measure of grid convergence. Importantly, the AMP scheme is shown to be robust and stable for solid bodies with a wide range of densities.

6 Conclusions

A stable added-mass partitioned (AMP) algorithm was developed for fluid-structure interaction problems involving viscous incompressible flow and compressible elastic solids. The new algorithm is stable, without sub-time-step iterations, for both heavy and very light solids and effectively suppresses both added-mass and added-damping effects. Key elements of the new AMP scheme are a Robin interface condition for the pressure and an impedance-weighted interface projection based on a new form for the fluid impedance. The scheme was implemented in two dimensions using deforming composite grids to treat complex geometry with multiple deforming bodies. The fluid is advanced using a fractional-step IMEX scheme with the viscous terms treated implicitly, while the hyperbolic system for the solid is integrated using an explicit upwind scheme. Both domain solvers are individually second-order accurate in both time and space. A number of new benchmark problems were developed including a radial-piston problem where exact solutions for radial and azimuthal motions were found. Traveling-wave exact solutions for a solid disk surrounded by an annulus of fluid were also derived. These exact solutions were used to verify the stability and second-order accuracy of the scheme. A clean benchmark problem involving a pressure-driven flow in a channel past a deformable solid annulus was developed and simulated; a self-convergence grid-refinement study was used to demonstrate second-order accuracy. A final example of flow past five solid bodies with densities ranging from $10^{-1}$ to $10^{3}$ illustrated the flexibility and robustness of the new scheme.

Appendix A Radial piston exact solution

The geometry of a FSI problem involving radial motion of a circularly symmetric elastic piston is shown in Figure 3. An incompressible fluid occupies the annular region between the fluid-solid interface at $r=r_{I}(t)$ and an outer boundary at $r=R$. The fluid surrounds an elastic solid disk for $0<r<r_{I}(t)$. It is assumed that the solid deforms in the radial direction only (as shown in the figure), and that the corresponding fluid motion is also confined to the radial direction. Under these assumptions, the fluid continuity equation (1a) in cylindrical coordinates reduces to

[TABLE]

so that the radial component of the fluid velocity is given by

[TABLE]

where $V(t)$ is the radial velocity of the fluid at $r=R$. The radial component of the momentum equation (1b) reduces to

[TABLE]

Integration with respect to $r$ gives

[TABLE]

where $P(t)$ is the fluid pressure at $r=R$. The evolution of the fluid velocity and pressure is completely determined by the interface position, $r_{I}(t)$, and the velocity and pressure at the outer radius, $V(t)$ and $P(t)$. These functions are found by matching the evolution of the fluid with that of the solid, which is considered next.

In the solid reference domain, the radial component of the momentum equation reduces to

[TABLE]

where $\bar{u}_{r}(\bar{r},t)$ is the radial component of the solid displacement and $r_{0}=r_{I}(0)$ is the initial radius of the interface. Bounded solutions of (54) as $\bar{r}\rightarrow 0$ are sought in the form

[TABLE]

where $J_{1}$ is the Bessel function of the first kind of order one, $\omega$ is the frequency of the oscillating displacement, and $\beta$ determines the magnitude of the oscillation. In terms of the radial displacement in (55), the interface position in physical space is given by

[TABLE]

Matching the radial velocity of the fluid given in (51) to that of the solid determined from (55) on the interface gives

[TABLE]

which specifies $V(t)$. The matching condition on the radial component of stress is

[TABLE]

Using (53) and (55) in the stress condition (58) gives

[TABLE]

An exact solution of the radial elastic piston problem is available based on the choice for the solid displacement in (55), which leads to the radial position of the fluid-solid interface given in (56). Matching the radial velocity and stress at the interface specifies $V(t)$ and $P(t)$ from (57) and (59), respectively, and these functions determine the fluid velocity and pressure given in (51) and (53).

Appendix B Exact solution for the rotating elastic disk in a fluid annulus

An exact solution can be constructed for the case of circumferential motion of a circularly symmetric elastic disk. It is assumed that the radial components of the fluid velocity and the solid displacement are zero, and thus the interface is located at $r=r_{0}$ for all time. It is also assumed that the circumferential components of the fluid velocity, $v_{\theta}$, and the solid displacement, $\bar{u}_{\theta}$, depend on $r$ and $t$ alone, so that the equations governing these quantities reduce to

[TABLE]

and

[TABLE]

respectively. Since there is no displacement in the radial direction, the reference coordinate, $\bar{r}$, is equivalent to the physical coordinate, $r$. Setting

[TABLE]

and substituting these into (60) and (61) it is found that the coefficient functions, $\hat{v}_{\theta}(r)$ and $\hat{\bar{u}}_{\theta}(r)$, both satisfy Bessel equations of order one. Assuming a no-slip boundary condition on the fluid velocity at $r=R$ and a bounded solution for the solid displacement as $r\rightarrow 0$, gives solutions of the form

[TABLE]

where $J_{1}$ and $Y_{1}$ are Bessel functions of the first and second kind, respectively. Here, $\lambda^{2}=i\omega/\nu,$ $k_{s}=\omega/\bar{c}_{s},$ and $(b,\bar{b})$ are constants. Matching conditions on velocity and stress imply

[TABLE]

Substituting (63) into the matching conditions in (64), and using the identity

[TABLE]

leads to the system

[TABLE]

Nontrivial solutions for $(b,\bar{b})$ exist if the following determinant condition is met:

[TABLE]

Values of $\omega$ with $\Re(\omega)>0$ satisfying $\mathcal{D}_{2}(\omega)=0$ imply solutions of the rotating elastic disk problem that decay in time. For such a value of $\omega$, a unique choice for $(b,\bar{b})$, satisfying (65), can be expressed as

[TABLE]

where $\bar{u}_{0}=\bar{u}_{\theta}(r_{0},0)$ is the initial displacement of the interface. Solutions for the circumferential components of the fluid velocity and the solid displacement are taken from the real parts of $v_{\theta}(r,t)$ and $\bar{u}_{\theta}(r,t)$ given in (62), while the fluid pressure is given by

[TABLE]

assuming $p=0$ on the interface $r=r_{I}(0)=r_{0}$.

Appendix C Radial traveling wave solution

We consider traveling wave solutions of an FSI problem involving a Stokes fluid in the circular region $0<r<r_{0}$ initially surrounded by a linearly elastic solid in the annular region $r_{0}<r<R$. The solution is linearized about small interface deformations, therefore it is assumed that the interface is located at $r=r_{0}$ for all time.

In the fluid, solutions are sought of the form

[TABLE]

where $\omega$ is a frequency, $n$ is a positive integer, and $\hat{v}_{r}(r)$, $\hat{v}_{\theta}(r)$ and $\hat{p}(r)$ are coefficient functions. Upon substitution, the equations governing the Stokes fluid become

[TABLE]

An elliptic equation for the pressure can be obtained from (68)–(70), and it takes the form

[TABLE]

Bounded solutions of the Cauchy-Euler equation in (71) at $r=0$ are

[TABLE]

where $p_{I}$ is a constant corresponding to the pressure at the (linearized) fluid-solid interface. Using (68) to eliminate $\hat{v}_{\theta}$ from (69) gives

[TABLE]

where $\lambda^{2}=i\omega/\nu$. Bounded solutions at $r=0$ have the form

[TABLE]

where $d$ is a constant and $J_{n}(z)$ is the Bessel function of the first kind of order $n$. The corresponding solution for $\hat{v}_{\theta}$, obtained from (68), is given by

[TABLE]

In the solid, we consider the solid displacement in the form

[TABLE]

where $\phi$ is a potential and $\mathbf{H}$ is a vector field. Let $\mathbf{e}_{r}=[\cos\theta,\sin\theta,0]^{T},$ $\mathbf{e}_{\theta}=[-\sin\theta,\cos\theta,0]^{T},$ and $\mathbf{e}_{z}=[0,0,1]^{T}$ represent the polar basis. Since we are interested in displacements in the two-dimensional $(r,\theta)$ plane, we take $\mathbf{H}=H\mathbf{e}_{z}$. With these assumptions, the equations governing linear elasticity reduce to a longitudinal and shear wave equations for $\phi(r,\theta,t)$ and $H(r,\theta,t)$, respectively, having the form

[TABLE]

As before, we seek solutions of the form

[TABLE]

Upon substitution into (77), we find

[TABLE]

where $k_{p}=\omega/\bar{c}_{p}$ and $k_{s}=\omega/\bar{c}_{s}$. General solutions of the two equations in (79) involve Bessel functions of the first and second kinds, and are given by

[TABLE]

where $(\bar{d}_{1},\bar{d}_{2},\bar{d}_{3},\bar{d}_{4})$ are constants. Using (76), the coefficient functions for the solid displacements in the radial and circumferential directions are given by

[TABLE]

respectively. For later purposes, we also require coefficient functions for the stress components $(\bar{\sigma}_{rr},\bar{\sigma}_{r\theta})$ in terms of $\hat{\phi}$ and $\hat{H}$, which are given by

[TABLE]

The solution for the fluid involves the two constants $(p_{I},d)$, while the solution for the solid involves the four constants $(\bar{d}_{1},\bar{d}_{2},\bar{d}_{3},\bar{d}_{4})$. Constraints for these six constants can be obtained by imposing boundary conditions on the solid at $r=R$ and the matching conditions at the fluid-solid interface, $r=r_{0}$. At the outer boundary of the solid, we assume traction-free boundary conditions of the form

[TABLE]

while matching conditions on the velocity and stress at the interface give

[TABLE]

The two constraints in (84) with the four constraints in (85) lead to a homogeneous system of equations of the form

[TABLE]

for the vector of constants $\mathbf{d}=[p_{I},d,\bar{d}_{1},\bar{d}_{2},\bar{d}_{3},\bar{d}_{4}]^{T}$. Nontrivial solutions for $\mathbf{d}$ exist if $\det(M)=0$, which leads to a dispersion relation between the frequency $\omega$ and the circumferential wave number $n$ of the form

[TABLE]

For a given $(n,\omega)$-pair satisfying (87), a nontrivial $\mathbf{d}$ can be found subject to the normalization

[TABLE]

where $\bar{u}_{0}$ is an initial displacement of the interface. Solutions of (87) for $(n,\omega)$ and the corresponding constants in $\mathbf{d}$ are given in Section 5.2 for the case $\bar{u}_{0}=10^{-7}$.