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

TL;DR

This paper introduces a stable, second-order accurate partitioned algorithm for fluid-structure interaction problems involving incompressible flow and elastic solids, effectively handling large added-mass effects.

Contribution

A novel AMP scheme that remains stable and accurate for large added-mass effects in FSI problems, with detailed stability analysis and verification.

Findings

The AMP scheme is stable for any solid-to-fluid density ratio.

The algorithm achieves second-order accuracy in complex FSI scenarios.

Validated with exact benchmark solutions including large interface deformations.

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 remains stable and second-order accurate even when added-mass and added-damping effects are large. 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 fluid impedance is defined using material and discretization parameters and follows from a careful analysis of the discretization of the governing equations and coupling conditions near the interface. A normal mode analysis is performed to show that the AMP scheme is stable for a few carefully-selected model problems. Two extensions of the analysis in Banks et al. are considered, including a first-order accurate discretization of a viscous model problem and a second-order accurate discretization of an inviscid model problem. The AMP algorithm is shown to be stable for any ratio of solid and fluid densities, including when added-mass effects are large. The algorithm is verified for accuracy and stability for a set of new exact benchmark solutions where finite interface deformations are permitted. The AMP scheme is found to be stable and second-order accurate even for very difficult cases of very light solids.

Figures

Tables

Table 1. Table 1: Longitudinal motion of an elastic piston: Maximum-norm errors and convergence rates of the numerical solution at t final = 0.6 subscript 𝑡 final 0.6 t_{\rm final}=0.6 computed using the AMP algorithm for ρ ¯ / ρ = δ = 10 3 ¯ 𝜌 𝜌 𝛿 superscript 10 3 \bar{\rho}/\rho=\delta=10^{3} , 1 1 1 and 10 − 3 . superscript 10 3 10^{-3}.

Heavy solid ($\delta ={10}^3$):

Table 2. Table 2: Solutions to the dispersion relation in ( 142 ) for H = 1 𝐻 1 H=1 , H ¯ = 0.5 ¯ 𝐻 0.5 \bar{H}=0.5 , ρ = 1 𝜌 1 \rho=1 , μ = 0.1 𝜇 0.1 \mu=0.1 and δ = ρ ¯ = μ ¯ = λ ¯ 𝛿 ¯ 𝜌 ¯ 𝜇 ¯ 𝜆 \delta=\bar{\rho}=\bar{\mu}=\bar{\lambda} .

$\delta$	$\omega$
${10}^3$	$3.141+7.930\cdot {10}^{-4}$
${10}^0$	$2.351+5.433\cdot {10}^{-1}$
${10}^{-3}$	$6.285+1.784\cdot {10}^{-3}$

Table 3. Table 3: Transverse motion of an elastic piston: Maximum-norm errors and convergence rates of the numerical solution computed using the AMP algorithm for ρ ¯ / ρ = δ = 10 3 ¯ 𝜌 𝜌 𝛿 superscript 10 3 \bar{\rho}/\rho=\delta=10^{3} , 1 1 1 and 10 − 3 . superscript 10 3 10^{-3}. The exact solutions use values of ω 𝜔 \omega from Table  2 and t final = 0.3 subscript 𝑡 final 0.3 t_{\rm final}=0.3 .

Heavy solid ($\delta ={10}^3$):

Equations

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

$$\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}.$$

$$\mathbf{\sigma }=-p\mathbf{I}+\mathbf{\tau },\mathbf{\tau }$$

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

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

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

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

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

$$\mathbf{v}$$

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

$$-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),$$

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

$$\displaystyle\Bigl{.}\mathbf{v}(\mathbf{x},t-\Delta t)\Bigr{|}_{\mathbf{x}={\mathcal{P}}(t-\Delta t)}\approx\Bigl{.}\bigl{(}\mathbf{v}(\mathbf{x},t)-\Delta tD_{t}\mathbf{v}(\mathbf{x},t)\bigr{)}\Bigr{|}_{\mathbf{x}={\mathcal{P}}(t)},$$

$$\displaystyle\Bigl{.}\bar{\mathbf{v}}(\bar{\mathbf{x}},t-\Delta t)\Bigr{|}_{\bar{\mathbf{x}}={\bar{\mathcal{P}}}_{0}}\approx\Bigl{.}\bigl{(}\bar{\mathbf{v}}(\bar{\mathbf{x}},t)-\Delta t\bar{\mathbf{v}}_{t}(\bar{\mathbf{x}},t)\bigr{)}\Bigr{|}_{\bar{\mathbf{x}}={\bar{\mathcal{P}}}_{0}},$$

$$-p+{\mathbf{n}}^T\mathbf{\tau }\mathbf{n}+{\bar{z}}_p\Delta t{\mathbf{n}}^T{D}_t\mathbf{v}={\mathbf{n}}^T\bar{\mathbf{\sigma }}\mathbf{n}+{\bar{z}}_p\Delta t{\mathbf{n}}^T{\bar{\mathbf{v}}}_t,\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}\overset{{\rm def}}{=}{\mathcal{C}}_{\rm AM}\Bigl{(}\frac{\rho h}{\Delta t}\Bigr{)}+{\mathcal{C}}_{\rm AD}\Bigl{(}\frac{\mu}{h}\Bigr{)},$$

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

$$\begin{array}[]{l}\hbox{// {Predictor steps}}\vskip 3.0pt plus 1.0pt minus 1.0pt\\ 1.\text{ Predict solid:}\\ \hskip 10.79999pt\begin{cases}\bar{\mathbf{u}}_{\bar{\mathbf{i}}}^{(p)}=\bar{\mathbf{u}}_{\bar{\mathbf{i}}}^{n}+\Delta t\bar{\mathbf{v}}_{\bar{\mathbf{i}}}^{n}+\frac{\Delta t^{2}}{2\bar{\rho}}\bar{\nabla}_{h}\cdot\bar{\boldsymbol{\sigma}}_{\bar{\mathbf{i}}}^{n},&\qquad{\bar{\mathbf{i}}}\in{\bar{\Omega}}_{h}.\\ \bar{\mathbf{q}}_{\bar{\mathbf{i}}}^{(p)}=\bar{\mathbf{q}}_{\bar{\mathbf{i}}}^{n}-\Delta t\sum_{m=1}^{3}\frac{1}{{\Delta\bar{x}}_{m}}\Bigl{(}\bar{\mathbf{F}}_{m,\;{\bar{\mathbf{i}}}}^{+}-\bar{\mathbf{F}}_{m,\;{\bar{\mathbf{i}}}}^{-}\Bigr{)},&\qquad{\bar{\mathbf{i}}}\in{\bar{\Omega}}_{h},\end{cases}\vskip 6.0pt plus 2.0pt minus 2.0pt\\ 2.\text{ Predict fluid grid: advance fluid grid to $t^{n+1}$ using $\bar{\mathbf{u}}_{\bar{\mathbf{i}}}^{(p)}$ for ${\bar{\mathbf{i}}}\in{\bar{\Gamma}}_{h}$, and compute grid velocity.}\vskip 6.0pt plus 2.0pt minus 2.0pt\\ 3.\text{ Predict fluid velocity:}\\ \hskip 10.79999pt\begin{cases}\mathbf{v}_{\mathbf{i}}^{(p)}=\mathbf{v}_{\mathbf{i}}^{n}+\frac{\Delta t}{2}\bigl{(}3\mathbf{N}_{h}(\mathbf{v}_{\mathbf{i}}^{n},p_{\mathbf{i}}^{n})-\mathbf{N}_{h}(\mathbf{v}_{\mathbf{i}}^{n-1},p_{\mathbf{i}}^{n-1})\bigr{)}+\frac{\Delta t}{2}\bigl{(}\mathbf{L}_{h}(\mathbf{v}_{\mathbf{i}}^{(p)})+\mathbf{L}_{h}(\mathbf{v}_{\mathbf{i}}^{n})\bigr{)},&\qquad\mathbf{i}\in\Omega_{h}\backslash\Gamma_{h},\\ \mathbf{t}_{m}^{T}\boldsymbol{\tau}_{\mathbf{i}}^{(p)}\mathbf{n}+\bar{z}_{s}\mathbf{t}_{m}^{T}\mathbf{v}_{\mathbf{i}}^{(p)}=\mathbf{t}_{m}^{T}\bar{\boldsymbol{\sigma}}_{\bar{\mathbf{i}}}^{(p)}\mathbf{n}+\bar{z}_{s}\mathbf{t}_{m}^{T}\bar{\mathbf{v}}_{{\bar{\mathbf{i}}}}^{(p)},&\qquad\mathbf{i}\in\Gamma_{h},\;{\bar{\mathbf{i}}}\in{\bar{\Gamma}}_{h},\\ \nabla_{h}\cdot\mathbf{v}_{\mathbf{i}}^{(p)}=0,&\qquad\mathbf{i}\in\Gamma_{h},\\ \mathbf{n}^{T}\mathbf{v}_{\mathbf{i}}^{(p)}=\frac{z_{f}}{z_{f}+\bar{z}_{p}}\mathbf{n}^{T}\mathbf{V}^{p}_{h}(\mathbf{v}_{\mathbf{i}}^{(p)})+\frac{\bar{z}_{p}}{z_{f}+\bar{z}_{p}}\mathbf{n}^{T}\bar{\mathbf{v}}_{\bar{\mathbf{i}}}^{(p)},\qquad\mathbf{t}_{m}^{T}\mathbf{v}_{\mathbf{i}}^{(p)}=\mathbf{t}_{m}^{T}\mathbf{V}^{p}_{h}(\mathbf{v}_{\mathbf{i}}^{(p)}),&\qquad\mathbf{i}\in\Gamma_{h},\;{\bar{\mathbf{i}}}\in{\bar{\Gamma}}_{h},\\ \text{Velocity boundary conditions on $\partial\Omega_{h}\backslash\Gamma_{h}$.}\end{cases}\vskip 6.0pt plus 2.0pt minus 2.0pt\\ 4.\text{ Predict fluid pressure:}\\ \hskip 10.79999pt\begin{cases}\Delta_{h}p_{\mathbf{i}}^{(p)}=-\rho\nabla_{h}\mathbf{v}_{\mathbf{i}}^{(p)}:\bigl{(}\nabla_{h}\mathbf{v}_{\mathbf{i}}^{(p)}\bigr{)}^{T}+\alpha_{\mathbf{i}}\nabla_{h}\cdot\mathbf{v}_{\mathbf{i}}^{(p)},&\quad\mathbf{i}\in\Omega_{h},\\ -p_{\mathbf{i}}^{(p)}-\frac{\bar{z}_{p}\Delta t}{\rho}(\mathbf{n}\cdot\nabla_{h})p_{\mathbf{i}}^{(p)}=\mathbf{n}^{T}\bigl{(}\bar{\boldsymbol{\sigma}}_{\bar{\mathbf{i}}}^{(p)}\mathbf{n}-\boldsymbol{\tau}_{\mathbf{i}}^{(p)}\mathbf{n}\bigr{)}+\bar{z}_{p}\Delta t\mathbf{n}^{T}\bigl{(}(\bar{\mathbf{v}}_{t})_{\bar{\mathbf{i}}}^{(p)}+\nu\nabla_{h}\times\nabla_{h}\times\mathbf{v}_{\mathbf{i}}^{(p)}\bigr{)},&\quad\mathbf{i}\in\Gamma_{h},\;{\bar{\mathbf{i}}}\in{\bar{\Gamma}}_{h},\\ \text{Pressure boundary conditions on $\partial\Omega_{h}\backslash\Gamma_{h}$.}\end{cases}\vskip 6.0pt plus 2.0pt minus 2.0pt\\ 5.\text{ Project solid interface:}\\ \hskip 10.79999pt\begin{cases}\mathbf{n}^{T}\bar{\mathbf{v}}_{\bar{\mathbf{i}}}^{I}=\frac{z_{f}}{z_{f}+\bar{z}_{p}}\mathbf{n}^{T}\mathbf{v}_{\mathbf{i}}^{(p)}+\frac{\bar{z}_{p}}{z_{f}+\bar{z}_{p}}\mathbf{n}^{T}\bar{\mathbf{v}}_{\bar{\mathbf{i}}}^{(p)}+{\frac{1}{z_{f}+\bar{z}_{p}}\bigl{(}\mathbf{n}^{T}\bar{\boldsymbol{\sigma}}^{(p)}_{\mathbf{i}}\mathbf{n}-\mathbf{n}^{T}{\boldsymbol{\sigma}}^{(p)}_{\mathbf{i}}\mathbf{n}\bigr{)}},&\qquad{\bar{\mathbf{i}}}\in{\bar{\Gamma}}_{h},\mathbf{i}\in\Gamma_{h},\\ \mathbf{t}_{m}^{T}\bar{\mathbf{v}}_{\bar{\mathbf{i}}}^{I}=\frac{z_{f}}{z_{f}+\bar{z}_{s}}\mathbf{t}_{m}^{T}\mathbf{v}_{\mathbf{i}}^{(p)}+\frac{\bar{z}_{s}}{z_{f}+\bar{z}_{s}}\mathbf{t}_{m}^{T}\bar{\mathbf{v}}_{\bar{\mathbf{i}}}^{(p)}+{\frac{1}{z_{f}+\bar{z}_{s}}\bigl{(}\mathbf{t}_{m}^{T}\bar{\boldsymbol{\sigma}}^{(p)}_{\mathbf{i}}\mathbf{n}-\mathbf{t}_{m}^{T}{\boldsymbol{\sigma}}^{(p)}_{\mathbf{i}}\mathbf{n}\bigr{)}},&\qquad{\bar{\mathbf{i}}}\in{\bar{\Gamma}}_{h},\mathbf{i}\in\Gamma_{h},\\ \mathbf{n}^{T}\bar{\boldsymbol{\sigma}}^{I}_{{\bar{\mathbf{i}}}}\mathbf{n}=\frac{z_{f}^{-1}}{z_{f}^{-1}+\bar{z}_{p}^{-1}}\mathbf{n}^{T}{\boldsymbol{\sigma}}^{(p)}_{\mathbf{i}}\mathbf{n}+\frac{\bar{z}_{p}^{-1}}{z_{f}^{-1}+\bar{z}_{p}^{-1}}\mathbf{n}^{T}\bar{\boldsymbol{\sigma}}_{\bar{\mathbf{i}}}^{(p)}\mathbf{n}+{\frac{1}{z_{f}^{-1}+\bar{z}_{p}^{-1}}\bigl{(}\mathbf{n}^{T}\bar{\mathbf{v}}^{(p)}_{\bar{\mathbf{i}}}-\mathbf{n}^{T}\mathbf{v}^{(p)}_{\mathbf{i}}\bigr{)}},&\qquad{\bar{\mathbf{i}}}\in{\bar{\Gamma}}_{h},\mathbf{i}\in\Gamma_{h}.\\ \mathbf{t}_{m}^{T}\bar{\boldsymbol{\sigma}}^{I}_{{\bar{\mathbf{i}}}}\mathbf{n}=\frac{z_{f}^{-1}}{z_{f}^{-1}+\bar{z}_{s}^{-1}}\mathbf{t}_{m}^{T}{\boldsymbol{\sigma}}^{(p)}_{\mathbf{i}}\mathbf{n}+\frac{\bar{z}_{s}^{-1}}{z_{f}^{-1}+\bar{z}_{s}^{-1}}\mathbf{t}_{m}^{T}\bar{\boldsymbol{\sigma}}_{\bar{\mathbf{i}}}^{(p)}\mathbf{n}+{\frac{1}{z_{f}^{-1}+\bar{z}_{s}^{-1}}\bigl{(}\mathbf{t}_{m}^{T}\bar{\mathbf{v}}^{(p)}_{\bar{\mathbf{i}}}-\mathbf{t}_{m}^{T}\mathbf{v}^{(p)}_{\mathbf{i}}\bigr{)}},&\qquad{\bar{\mathbf{i}}}\in{\bar{\Gamma}}_{h},\mathbf{i}\in\Gamma_{h},\\ {\bar{\mathbf{v}}_{\bar{\mathbf{i}}}^{(p)}\leftarrow\bar{\mathbf{v}}_{\bar{\mathbf{i}}}^{I},\quad\bar{\boldsymbol{\sigma}}_{\bar{\mathbf{i}}}^{(p)}\mathbf{n}\leftarrow\bar{\boldsymbol{\sigma}}_{\bar{\mathbf{i}}}^{I}\mathbf{n},}&\qquad{\bar{\mathbf{i}}}\in{\bar{\Gamma}}_{h},\mathbf{i}\in\Gamma_{h},\\ \text{Apply solid boundary conditions and set all ghost points.}\end{cases}\vskip 6.0pt plus 2.0pt minus 2.0pt\\ \hbox{// {Corrector steps}}\vskip 3.0pt plus 1.0pt minus 1.0pt\\ 6.\text{ Correct fluid grid: recompute grid velocity using $\bar{\mathbf{v}}_{\bar{\mathbf{i}}}^{I}$ for ${\bar{\mathbf{i}}}\in{\bar{\Gamma}}_{h}$.}\vskip 6.0pt plus 2.0pt minus 2.0pt\\ 7.\text{ Correct fluid velocity:}\\ \hskip 10.79999pt\begin{cases}\mathbf{v}_{\mathbf{i}}^{n+1}=\mathbf{v}_{\mathbf{i}}^{n}+\frac{\Delta t}{2}\bigl{(}\mathbf{N}_{h}(\mathbf{v}_{\mathbf{i}}^{(p)},p_{\mathbf{i}}^{(p)})+\mathbf{N}_{h}(\mathbf{v}_{\mathbf{i}}^{n},p_{\mathbf{i}}^{n})\bigr{)}+\frac{\Delta t}{2}\bigl{(}\mathbf{L}_{h}(\mathbf{v}_{\mathbf{i}}^{n+1})+\mathbf{L}_{h}(\mathbf{v}_{\mathbf{i}}^{n})\bigr{)},&\qquad\mathbf{i}\in\Omega_{h}\backslash\Gamma_{h},\\ \mathbf{t}_{m}^{T}\boldsymbol{\tau}_{\mathbf{i}}^{n+1}\mathbf{n}+\bar{z}_{s}\mathbf{t}_{m}^{T}\mathbf{v}_{\mathbf{i}}^{n+1}=\mathbf{t}_{m}^{T}\bar{\boldsymbol{\sigma}}_{\bar{\mathbf{i}}}^{I}\mathbf{n}+\bar{z}_{s}\mathbf{t}_{m}^{T}\bar{\mathbf{v}}_{{\bar{\mathbf{i}}}}^{I},&\qquad\mathbf{i}\in\Gamma_{h},\;{\bar{\mathbf{i}}}\in{\bar{\Gamma}}_{h},\\ \nabla_{h}\cdot\mathbf{v}_{\mathbf{i}}^{n+1}=0,&\qquad\mathbf{i}\in\Gamma_{h},\\ \mathbf{n}^{T}\mathbf{v}_{\mathbf{i}}^{n+1}=\frac{z_{f}}{z_{f}+\bar{z}_{p}}\mathbf{n}^{T}\mathbf{V}_{h}(\mathbf{v}_{\mathbf{i}}^{n+1})+\frac{\bar{z}_{p}}{z_{f}+\bar{z}_{p}}\mathbf{n}^{T}\bar{\mathbf{v}}_{\bar{\mathbf{i}}}^{I},\qquad\mathbf{t}_{m}^{T}\mathbf{v}_{\mathbf{i}}^{n+1}=\mathbf{t}_{m}^{T}\mathbf{v}_{\mathbf{i}}^{(e)},\qquad m=1,2,&\qquad\mathbf{i}\in\Gamma_{h},\;{\bar{\mathbf{i}}}\in{\bar{\Gamma}}_{h},\\ \text{Velocity boundary conditions on $\partial\Omega_{h}\backslash\Gamma_{h}$.}\end{cases}\vskip 6.0pt plus 2.0pt minus 2.0pt\\ 8.\text{ Correct fluid pressure.}\\ \hskip 10.79999pt\begin{cases}\Delta_{h}p_{\mathbf{i}}^{n+1}=-\rho\nabla_{h}\mathbf{v}_{\mathbf{i}}^{n+1}:\bigl{(}\nabla_{h}\mathbf{v}_{\mathbf{i}}^{n+1}\bigr{)}^{T}+\alpha_{\mathbf{i}}\nabla_{h}\cdot\mathbf{v}_{\mathbf{i}}^{n+1},&\quad\mathbf{i}\in\Omega_{h},\\ -p_{\mathbf{i}}^{n+1}-\frac{\bar{z}_{p}\Delta t}{\rho}(\mathbf{n}\cdot\nabla_{h})p_{\mathbf{i}}^{n+1}=\mathbf{n}^{T}\bigl{(}\bar{\boldsymbol{\sigma}}_{\bar{\mathbf{i}}}^{I}\mathbf{n}-\boldsymbol{\tau}_{\mathbf{i}}^{n+1}\mathbf{n}\bigr{)}+\bar{z}_{p}\Delta t\mathbf{n}^{T}\bigl{(}(\bar{\mathbf{v}}_{t})_{\bar{\mathbf{i}}}^{I}+\nu\nabla_{h}\times\nabla_{h}\times\mathbf{v}_{\mathbf{i}}^{n+1}\bigr{)},&\quad\mathbf{i}\in\Gamma_{h},\;{\bar{\mathbf{i}}}\in{\bar{\Gamma}}_{h},\\ \text{Pressure boundary conditions on $\partial\Omega_{h}\backslash\Gamma_{h}$.}\end{cases}\vskip 6.0pt plus 2.0pt minus 2.0pt\\ 9.\text{ Correct solid interface.}\\ \hskip 10.79999pt\begin{cases}{\bar{\mathbf{v}}_{\bar{\mathbf{i}}}^{n+1}=\mathbf{v}_{\mathbf{i}}^{n+1}},&\qquad{\bar{\mathbf{i}}}\in{\bar{\Gamma}}_{h},\mathbf{i}\in\Gamma_{h},\\ {\bar{\boldsymbol{\sigma}}^{n+1}_{{\bar{\mathbf{i}}}}\mathbf{n}={\boldsymbol{\sigma}}^{n+1}_{\mathbf{i}}\mathbf{n},}&\qquad{\bar{\mathbf{i}}}\in{\bar{\Gamma}}_{h},\mathbf{i}\in\Gamma_{h},\\ \text{Reset ghost points corresponding to ${\bar{\mathbf{i}}}\in{\bar{\Gamma}}_{h}$.}\end{cases}\\ \end{array}$$

$$\displaystyle\mathbf{N}_{h}(\mathbf{v}_{\mathbf{i}},p_{\mathbf{i}})=-\bigl{(}(\mathbf{v}_{\mathbf{i}}-\dot{\mathbf{x}}_{\mathbf{i}})\cdot\nabla_{h}\bigr{)}\mathbf{v}_{\mathbf{i}}-\frac{1}{\rho}\nabla_{h}\,p_{\mathbf{i}},\qquad\mathbf{L}_{h}(\mathbf{v}_{\mathbf{i}})=\nu\Delta_{h}\mathbf{v}_{\mathbf{i}},$$

$$\displaystyle\mathbf{V}^{p}_{h}(\mathbf{v}_{\mathbf{i}}^{(p)})\overset{{\rm def}}{=}\mathbf{v}_{\mathbf{i}}^{n}+\frac{\Delta t}{2}\left(3\mathbf{N}_{h}(\mathbf{v}_{\mathbf{i}}^{n},p_{\mathbf{i}}^{n})-\mathbf{N}_{h}(\mathbf{v}_{\mathbf{i}}^{n-1},p_{\mathbf{i}}^{n-1})\right)+\frac{\Delta t}{2}\Big{(}\mathbf{L}_{h}(\mathbf{v}_{\mathbf{i}}^{(p)})+\mathbf{L}_{h}(\mathbf{v}_{\mathbf{i}}^{n})\Big{)}.$$

$${\mathbf{n}}^T{\mathbf{v}}_{\mathbf{i}}^{(p)}=\frac{z_f}{z_f+{\bar{z}}_p}{\mathbf{n}}^T{\mathbf{V}}_h^p({\mathbf{v}}_{\mathbf{i}}^{(p)})+\frac{{\bar{z}}_p}{z_f+{\bar{z}}_p}{\mathbf{n}}^T{\bar{\mathbf{v}}}_{\bar{\mathbf{i}}}^{(p)},\mathbf{i}\in {\Gamma }_h,\bar{\mathbf{i}}\in {\bar{\Gamma }}_h,$$

$${\mathbf{n}}^T{\mathbf{v}}_{\mathbf{i}}^{(p)}={\mathbf{n}}^T{\bar{\mathbf{v}}}_{\bar{\mathbf{i}}}^{(p)},\mathbf{i}\in {\Gamma }_h,\bar{\mathbf{i}}\in {\bar{\Gamma }}_h,$$

$$\{\begin{array}{ll}\rho {\partial }_t\mathbf{v}+\nabla p=\mu \Delta \mathbf{v},& \mathbf{x}\in \Omega ,\\ \Delta p=0,& \mathbf{x}\in \Omega ,\\ \nabla \cdot \mathbf{v}=0,& \mathbf{x}\in \Gamma ,\end{array}$$

$$\{\begin{array}{ll}\bar{\rho }{\partial }_t\bar{\mathbf{v}}=\nabla \cdot \bar{\mathbf{\sigma }},& \mathbf{x}\in {\bar{\Omega }}_0,\\ {\partial }_t\bar{\mathbf{\sigma }}=\bar{\lambda }(\nabla \cdot \mathbf{v})\mathbf{I}+\bar{\mu }(\nabla \mathbf{v}+(\nabla \mathbf{v}{)}^T),& \mathbf{x}\in {\bar{\Omega }}_0,\end{array}$$

$$\{\begin{array}{ll}\mathbf{v}=\bar{\mathbf{v}},& \mathbf{x}\in \Gamma ,\\ \mathbf{\sigma }\mathbf{n}=\bar{\mathbf{\sigma }}\mathbf{n},& \mathbf{x}\in \Gamma .\end{array}$$

$$q(x,y,t)=\sum _{k=-N_x/2}^{N_x/2}{e}^{2\pi ikx/L}{\hat{q}}_k(y,t),x\in [0,L],$$

$$\rho {\partial }_t{v}_1+i{k}_{x}p$$

$$\rho {\partial }_t{v}_2+{\partial }_{y}p$$

Full text

{AMS}

65M12, 74F10, 74S10, 76M20, 76D99

1 Introduction

In this paper, we consider fluid-structure interaction (FSI) of incompressible fluids and bulk elastic solids. FSI arises in many scientific and engineering applications including flow-induced vibrations of structures (i.e., aircraft wings, undersea cables, wind turbines, and bridges) and simulating blood flow in arteries and veins. FSI is typically modeled by solving suitable partial differential equations in the fluid and solid domains together with coupling conditions involving velocity and stress at the fluid-solid interface. FSI algorithms can either be categorized as monolithic schemes, where the solutions for the fluid and solid are advanced implicitly as one large system, or partitioned schemes, where the evolution of the fluid and solid solutions are decoupled from the interface conditions. While partitioned schemes are generally more computationally efficient than monolithic schemes, they may suffer from instabilities when added-mass effects are present. The goal of current work is to develop a robust partitioned scheme.

In the traditional partitioned (TP) scheme, 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. Without sub-iterations, it is well known that the TP scheme is unstable for light solids. On the other hand, we can consider the anti-traditional partitioned (ATP) scheme, which reverses the role of the solid and fluid. In this scheme, the solid provides a Neumann (traction) boundary condition for the fluid and the fluid supplies a Dirichlet (no-slip) boundary condition for the solid. This scheme is only stable when the solid is sufficiently light, see [1]. Stability of the aforementioned schemes is typically addressed by using an iteration of the interface conditions. Research on acceleration techniques such as Aitken acceleration or quasi-Newton methods has proven that the number of sub-iterations can be greatly reduced, see [2, 3]. Reduction of sub-iterations can also be achieved by considering Robin-Neumann or Robin-Robin coupling instead of the traditional Dirichlet-Neumann coupling [4, 5, 6, 7, 8, 9, 10, 11, 12, 13].

In recent work [1], we developed a new class of Added-Mass Partitioned (AMP) algorithms for FSI problems coupling incompressible flow and elastic solids. The fluid AMP interface conditions are Robin conditions motivated from the propagation of characteristics out of the solid domain. One important property of the AMP scheme is that it approaches the TP scheme in the heavy solid limit and the ATP scheme in the light solid limit. The algorithms use 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. During the pressure step, a Robin condition which involves both the interface traction and acceleration is applied. This condition is derived by manipulating the characteristic conditions using the governing equations of the fluid. The central aim of this paper is to extend the AMP algorithm in [1] for IMEX fluid solvers, where the viscous terms are treated implicitly and the convection and pressure gradient terms are treated explicitly. The adaptation to IMEX solvers necessitated important modifications for stability. An impedance-based projection of the interface velocity and stress (motivated from [14] for the case of compressible fluids and linearly-elastic solids) is applied after the fluid and solid updates to ensure the heavy and light solid limits are achieved. An implicit version of this projection is also used during the stage for the fluid velocity. These projections rely on defining the fluid impedance, which is a key ingredient which leads to a stable scheme.

A fluid impedance was defined for incompressible flows and bulk solids in [1]. This definition leads to a stable scheme for explicit time step restrictions which are typically limited by viscous terms, but led to instabilities for IMEX schemes where the time step is less restrictive and chosen by the convective terms. In this paper, a suitable definition for the fluid impedance is derived through an analysis of a discretization of the fluid equations near the interface. In this analysis, it was revealed that the impedance has an inertial component for treating added-mass instabilities and a viscous component for treating added-damping instabilities. For this choice of fluid impedance, an analysis of the AMP scheme is performed to show stability for any ratio of the mass of the solid to that of the fluid. The stability analysis extends [1] in two different directions. Previously in [1], stability was analyzed for a first-order accurate scheme applied to an inviscid incompressible fluid and acoustic solid in a Cartesian geometry. First we consider an extension of the model problem to include viscous fluids. Due to the introduction of viscosity, there are both pressure and shear forces at the interface and instabilities can arise from both added-mass and added-damping effects. This new analysis is complicated by the introduction of more dimensionless parameters, namely the viscous CFL number in the fluid. Despite the added complication, the AMP scheme is shown to be stable for all possible mass ratios and viscous CFL numbers when applied to this new model problem. Secondly, we consider the same inviscid model problem in [1] and analyze a second-order accurate scheme. In both extensions, we consider a detailed study of the traditional and anti-traditional schemes. The stability and accuracy of the AMP algorithm is tested on new exact solutions for a wide range of fluid and solid densities. The stability and accuracy of the AMP algorithm is verified numerically for new exact solutions. In our companion paper [15], the AMP algorithm is implemented using deforming composite grids for curvilinear geometries.

The remaining sections of the paper are organized as follows. The equations governing the FSI problem are described in Section 2. The AMP algorithm is summarized in Section 3. In Section 4, the fluid impedance is derived following an analysis of an FSI discretization near the interface. The stability of the AMP algorithm is analyzed in Section 5 for a viscous model problem and in Section 6 for an inviscid model problem. Section 7 provides numerical results confirming the stability and accuracy of the scheme. Some of the results use new exact solutions of benchmark FSI problems. Conclusions are given in Section 8.

2 Governing equations

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 velocity-pressure form of the incompressible Navier-Stokes equations is given by

[TABLE]

where

[TABLE]

Here, $\mathbf{v}(\mathbf{x},t)$ is the velocity, $p(\mathbf{x},t)$ is the pressure, $\rho$ is the (constant) density, and $\mu$ is the (constant) dynamic viscosity. The fluid stress tensor is given by

[TABLE]

where $\mathbf{I}$ is the identity matrix and $\boldsymbol{\tau}(\mathbf{x},t)$ is the viscous stress tensor. In the velocity-pressure form of the equations, an extra boundary condition is required and a suitable choice is to impose $\nabla\cdot\mathbf{v}=0$ for $\mathbf{x}\in\partial\Omega(t)$, see [16].

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 belonging to the solid.) The position of the solid in physical space is determined by the mapping

[TABLE]

where $\bar{\mathbf{u}}(\bar{\mathbf{x}},t)$ is the displacement of the solid. The Cauchy stress tensor $\bar{\boldsymbol{\sigma}}(\bar{\mathbf{x}},t)$ for a linearly-elastic solid is defined by

[TABLE]

where $\bar{\lambda}$ and $\bar{\mu}$ are Lamé parameters (taken to be constants). The solid equations are considered as a first-order system of PDEs in time and space, following [17], and are given by

[TABLE]

where $\bar{\mathbf{v}}(\bar{\mathbf{x}},t)$ is the velocity of the solid, and $\bar{\rho}$ is its density (assumed constant). In this form, upwind solvers can be used to advance displacement, velocity and stress of the solid. We note that (5) is enforced at $t=0$.

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. At the interface, the conversion between reference and physical coordinates is determined by the mapping in (4). 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.

3 AMP interface conditions and algorithm

In this section, we derive the AMP interface conditions at a continuous level and discuss their implementation in the AMP algorithm. The derivation follows the work in [1], but there are important modifications discussed to accommodate the IMEX fractional-step scheme used in the AMP algorithm to solve the equations in the fluid domain. These modifications are guided by a consideration of the behavior of the AMP interface conditions in the limits of very light and very heavy solids.

3.1 AMP interface conditions

The starting point for the derivation is the matching conditions involving velocity and stress in (7). Following [1], a linear combination of these conditions are expressed in terms of the outgoing characteristic variables of the solid, i.e.

[TABLE]

where $\mathbf{n}$ is the unit normal, $\mathbf{t}_{m}$, $m=1,2$, are mutually orthogonal unit vectors tangent to the interface, and $\bar{z}_{p}$ and $\bar{z}_{s}$ are impedances involving the characteristic velocities of the solid given by

[TABLE]

In the AMP algorithm, the conditions in (8) are interpreted as providing boundary conditions for the fluid at the interface with the outgoing characteristic quantities of the solid on the left-hand side considered to be known from a previous stage of the algorithm. While these conditions, along with $\nabla\cdot\mathbf{v}=0$ for $\mathbf{x}\in\Gamma(t)$, are sufficient conditions for the fluid equations in velocity-pressure form, a further manipulation is required to obtain suitable conditions to be used for the fractional-step solver. The objective is to separate the conditions in (8) to obtain a condition to be used in the IMEX time-stepping scheme for the fluid velocity and a condition for the subsequent update for the pressure.

For the Poisson problem for the fluid pressure, the interface condition in (8a) is used with the momentum equation in (1a) to derive a Robin condition for the pressure. The momentum equation involves the acceleration of the fluid, and this quantity may be obtained on the moving fluid-solid interface using the Taylor approximation

[TABLE]

where $D_{t}=\partial_{t}+\mathbf{v}\cdot\nabla$ is the material derivative, ${\mathcal{P}}(t)$ is a point on the moving interface and $\Delta t$ is a time-step. The corresponding approximation for the solid is

[TABLE]

where ${\bar{\mathcal{P}}}_{0}$ is the Lagrangian position associated with ${\mathcal{P}}(t)$. Using (9) and (10) in (8a), and assuming the fluid and solid velocities match on the interface at times $t-\Delta t$ and $t$, we obtain

[TABLE]

We may now eliminate the fluid acceleration using (1a) to obtain the following Robin condition for the fluid pressure:

[TABLE]

where $\partial_{n}=\mathbf{n}\cdot\nabla$ is the normal derivative and $\nu=\mu/\rho$ is the kinematic viscosity of the fluid. Following [16], 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 (12) in favor of the curl-curl operator. This is done for improved stability of the fractional-step scheme. The condition in (12), along with suitable conditions for $\mathbf{x}\in\partial\Omega(t)\backslash\Gamma(t)$, is used for the Poisson equation in (1b) for the pressure.

As was noted in [1], the remaining interface conditions in (8b), together with the continuity equation, can be used as boundary conditions to advance the fluid velocity. This was found to be an effective approach for an explicit integration of the momentum equations. To ensure that the fluid velocity and tractions match at the end of the time step, an interface projection is performed to give a common interface velocity $\mathbf{v}^{I}$ and interface traction $\boldsymbol{\sigma}^{I}\mathbf{n}$. In analogy to the interface projection used for compressible fluids in [14, 18], which is based on a characteristic analysis, the projection for incompressible fluids is also proposed to be of the form of an impedance-weighted average. For the velocity, the projection has the form

[TABLE]

while for the traction, the projection is based on an inverse impedance-weighted average of the form

[TABLE]

These projections introduce a fluid impedance, $z_{f}$, which is well defined for compressible fluids, but has no obvious definition for incompressible fluids. However, an analysis of a discrete approximation to the governing equations given in Section 4 suggests a form for $z_{f}$ given by

[TABLE]

where $h$ is an appropriate mesh spacing and $({\mathcal{C}}_{\rm AM},{\mathcal{C}}_{\rm AD})$ are constants whose approximate values are provided by the analysis. The projections in (13) and (14) can be used to set values of the fluid and solid velocity and traction at the interface.

For the IMEX scheme considered here, a modification of the previous approach in [1] is required in the implementation of the interface conditions for the fluid velocity. The issue is informed by considering the limits of very light and heavy solids. In the limit of a very light solid ($\bar{z}_{p},\bar{z}_{s}\rightarrow 0$), for example, the Robin condition in (12) becomes a Dirichlet condition for the pressure, while the interface conditions in (8b) reduce to matching conditions involving the shear stress of the fluid. The latter conditions, along with the continuity constraint, provide Neumann conditions on the fluid velocity. These conditions for the fluid pressure and velocity correspond to those for a free surface, and the latter are suitable for the implicit solution of the fluid velocity in the IMEX fractional-step scheme.

The difficulty is revealed in the limit of a very heavy solid ($\bar{z}_{p},\bar{z}_{s}\rightarrow\infty$). In this limit, the Robin condition in (12) becomes a Neumann condition for the fluid pressure balancing the acceleration of the interface as determined by the solid. This condition is analogous to the usual Neumann boundary condition for the pressure at a rigid boundary obtained from the fluid momentum equations as a compatibility condition (see [16] for example). The interface conditions in (8b) reduce to matching conditions involving the tangential components of velocity. However, the matching condition on the normal component of velocity,

[TABLE]

implied by (8a) in the limit of a heavy solid has been lost in the derivation of (12). A remedy can be obtained by using the interface projection for the normal component of the velocity in (13a) as a boundary condition for the implicit solution of the fluid velocity in the IMEX fractional-step scheme. The implementation of this approach is described next in the discussion of the AMP algorithm.

3.2 AMP algorithm

Algorithm 1 provides a consise description of the AMP time-stepping scheme (see [15] for additional details of the implementation of the algorithm). The algorithm advances the solution from a time $t^{n}$ to $t^{n+1}=t^{n}+\Delta t$. It is assumed that the fluid domain is represented by a grid consisting of interior points $\mathbf{i}\in\Omega_{h}$, boundary points $\mathbf{i}\in\partial\Omega_{h}$ and interface points $\mathbf{i}\in\Gamma_{h}$, where $\mathbf{i}=(i_{1},i_{2},i_{3})$ is a multi-index. Similarly, the solid reference domain is covered by a grid with interior points ${\bar{\mathbf{i}}}\in{\bar{\Omega}}_{h}$, boundary points ${\bar{\mathbf{i}}}\in\partial{\bar{\Omega}}_{h}$ and interface points ${\bar{\mathbf{i}}}\in{\bar{\Gamma}}_{h}$. Discrete operators, such as $\nabla_{\!h}$ and $\Delta_{h}$, denote approximations of the corresponding differential operators on the grid.

The time-stepping scheme uses a predictor-corrector approach. Steps 1–5 of Algorithm 1 describe the preditor steps. Predicted values for the solid displacement $\bar{\mathbf{u}}_{\bar{\mathbf{i}}}$ are obtained in Step 1 using a Lax-Wendroff-type scheme, while the solid velocity and stress $\bar{\mathbf{q}}_{\bar{\mathbf{i}}}=(\bar{\mathbf{v}}_{\bar{\mathbf{i}}},\bar{\boldsymbol{\sigma}}_{{\bar{\mathbf{i}}}})$ are advanced using a Godunov type scheme with numerical fluxes $\bar{\mathbf{F}}_{m,\;{\bar{\mathbf{i}}}}^{\pm}$ corresponding to the $\bar{x}_{m}$ coordinate direction. In Step 2, the solid displacement is used to compute the deformed fluid grid at time $t^{n+1}$.

The fluid velocity is predicted in Step 3. Here, $\mathbf{N}_{h}$ and $\mathbf{L}_{h}$ represent grid operators associated with the explicit and implicit terms in the velocity update, respectively, given by

[TABLE]

where $\dot{\mathbf{x}}_{\mathbf{i}}$ is the velocity of the grid. The explicit terms are advanced using an Adams-Bashforth scheme, while the implicit terms use Crank-Nicholson. The boundary conditions on the interface makes use of a predicted velocity, coming from the interior equation applied on the boundary, and defined by

[TABLE]

In particular, this velocity is used in the impedance-weighted average condition

[TABLE]

which is obtained from the projection in (13a). Here the the term involving the jump in the stress is dropped (because a suitable approximation for the fluid stress is unavailable at the predictor stage). Notice that (19) is an implicit condition on $\mathbf{v}_{\mathbf{i}}^{(p)}$ which appears on the left and right-hand sides. In the light-solid limit ($\bar{z}_{p}\rightarrow 0$), the boundary condition in (19) reduces to $\mathbf{n}^{T}\mathbf{v}_{\mathbf{i}}^{(p)}=\mathbf{n}^{T}\mathbf{V}^{p}_{h}(\mathbf{v}_{\mathbf{i}}^{(p)})$, which simply sets the normal component of the fluid velocity to be equal to that given by the interior time-stepping scheme applied on the boundary. In the heavy-solid limit ($\bar{z}_{p}\rightarrow\infty$), (19) becomes

[TABLE]

which recovers the matching condition in (16). Our later analysis of a viscous model problem (Section 5) and an inviscid model problem (Section 6), and subsequent numerical results (Section 7), verify that the boundary conditions used to advance the fluid velocity in the fractional-step scheme lead to stable and accurate results for a wide range of solid densities.

Steps 4 and 5 complete the set of steps belonging to the predictor stage of the algorithm. The predicted fluid pressure is computed in Step 4 by solving a discrete Poisson problem. This elliptic problem uses a discrete approximation of the Robin condition in (12). Finally, interface values for the solid velocity and traction are obtained in Step 5 using the impedance-weighted projections in (13) and (14). These interface values overwrite the corresponding predicted values of the solid on the boundary.

The set of corrector steps consisting of Steps 6–9 essentially mirror those of the predictor. The fluid grid is recomputed in Step 6 using an updated grid velocity obtained from solid velocity. In Step 7, the discrete fluid velocity at $t^{n+1}$ is determined, now using an Adams-Moulton corrector. The fluid pressure is updated in Step 8 according to the solution of a discrete Poisson problem. Lastly, the solid velocity and traction are set equal to the corrected fluid values in Step 9, so that the fluid and solid interface values agree at the end of the corrector stage of the AMP time-stepping scheme.

4 Derivation of the fluid impedance

The focus of this section is an analysis of a FSI problem leading to an expression for the fluid impedance that guides the choice for $z_{f}$ introduced in (15). A value for the fluid impedance is required in the formulas for the interface projections in (13) and (14), and also for the implementation of the IMEX fractional-step scheme for the fluid velocity. The interface projections were also used in our earlier paper [1], as these provide formulas for the interface velocity and traction that ensure the fluid and solid velocities and tractions match at the interface and that smoothly accommodate the limiting cases of light and heavy solids. The earlier paper introduced a fluid impedance given by $z_{f}=\rho H/\Delta t$, where $H$ was a measure of the depth of the fluid layer, and it was found that this choice led to a stable AMP algorithm when using an explicit fractional-step scheme for the fluid. It was also noted that the scheme was quite insensitive to the choice of $H$. For the present AMP scheme, we use an IMEX fractional-step scheme. Since the viscous terms in the fluid equations are now treated implicitly, the viscous CFL number, $\Lambda=\nu\Delta t/h^{2}$, can be large in which case the choice for $z_{f}$ used in [1] is no longer sufficient for stability. We have found that the difficulty can be resolved by considering an analysis of a more general FSI model problem for which the viscous terms in the fluid equations contribute.

Consider a FSI model problem in which the fluid occupies the two-dimensional domain, $\Omega$, given by $0<x<L$, $y>0$, while the solid exists on the domain, ${\bar{\Omega}}_{0}$, for $0<x<L$, $y<0$, see Figure 1. The fluid-solid interface, $\Gamma$, of length $L$ is linearized about a flat surface, $y=0$. The equations governing the model problem are

[TABLE]

Solutions of the model problem are assumed to be periodic in $x$ with period equal to $L$, and bounded as $y\rightarrow\pm\infty$.

The equations governing the fluid and solid are discretized in the $x$-direction on a uniform grid, $x_{\ell}=\ell\Delta x$ for $\ell=0,1,\ldots N_{x}$, with grid spacing $\Delta x=L/N_{x}$. Since the problem is periodic, each variable can be represented as a discrete Fourier series

[TABLE]

where $\hat{q}_{k}(y,t)$ are the Fourier coefficient functions and $N_{x}$ is an integer, assumed to be even for convenience. Taking a finite Fourier transform of the fluid equations in (20a) gives

[TABLE]

where $k_{x}=2\pi k/L$. The hats and $k$ subscripts on the coefficient functions in (22) have been dropped for notational convenience. The equations for the Fourier coefficient functions are now discretized in time. Define the grid functions $\mathbf{v}^{n}(y)\approx\mathbf{v}(y,t^{n})$ and $p^{n}\approx p(y,t^{n})$, where $t^{n}=n\Delta t$ for a (fixed) time step $\Delta t$. An implicit scheme to advance the solution from $t^{n}$ to $t^{n+1}$, based on backward-Euler time-stepping, is given by

[TABLE]

Assume that the coefficient functions for the solid variables have been advanced to $t=t^{n+1}$ using an upwind scheme, for example, and that $b_{p}^{n+1}$ and $b_{s}^{n+1}$ are, respectively, the normal and tangential components of the outgoing characteristic variables of the solid at $t^{n+1}$. Using (8), the boundary conditions for the fluid at $y=0$ take the form

[TABLE]

where the components of the fluid shear stress in (24) are given by

[TABLE]

The implicit scheme in (23) with boundary conditions in (24) at $y=0$ and boundedness as $y\rightarrow\infty$ determine the grid functions for the fluid at $t^{n+1}$ in terms the fluid velocity at $t^{n}$ and the outgoing solid data $(b_{p}^{n+1},b_{s}^{n+1})$.

Consider perturbations in the grid functions of the fluid at $t^{n+1}$ for $y>0$ subject to perturbations in the interface data $b_{p}^{n+1}$ and $b_{s}^{n+1}$ at $y=0$. Define

[TABLE]

and

[TABLE]

where $(\delta V_{1},\delta V_{2},\delta P)$ and $(\delta B_{p},\delta B_{s})$ are small perturbations. Assuming the fluid velocity at $t^{n}$ is fixed, the variational equations corresponding to (23) are

[TABLE]

Solution to these equations that remain bounded as $y\rightarrow\infty$ are

[TABLE]

where

[TABLE]

Substituting the solution for the perturbations of the fluid variables into the variational equations (24) for the interface conditions leads to the linear system

[TABLE]

where

[TABLE]

The coefficients $a_{ij}$ in the linear system are defined in terms of the dimensionless parameters

[TABLE]

The solution of the linear system

[TABLE]

determines the variation in the interface values of the fluid velocity and pressure in terms of the variations in the outgoing characteristic variables of the solid.

The AMP algorithm uses impedance-weighted averages to set values for the velocity and pressure at the interface. For example, the normal component of velocity at the interface is given by (13a). In terms of the variational problem, (13a) reduces to

[TABLE]

assuming that the fluid velocity and stress on the right-hand side are held fixed. In view of the solution in (35), we have

[TABLE]

which, after some manipulation, gives

[TABLE]

Of particular interest are the limiting cases when the viscous CFL number, $\Lambda$, is small and large. A straightforward analysis of the dimensionless parameter $R$ in (38) gives

[TABLE]

to leading order. In view of (39), we define the relatively simple approximation

[TABLE]

The plots in Figure 2 show that $R/\tilde{R}\approx 1$ over a wide range of values for $\theta_{s}$, so that the fluid impedance given by

[TABLE]

is a good approximation of the more complicated form given in (38). The model problem analysis of Section 5 confirms that this choice leads to a stable scheme.

Formula (40) provides the generic form of the fluid impedance we use, but it remains to make a choice for $k_{x}$ so that the approximation can be used for a discrete approximation in physical coordinates (as opposed to the Fourier transformed space). Note that in a discrete approximation, the possible wave numbers $k_{x}$ appearing in (40) are bounded by approximately $1/h$, where $h$ is a measure of the grid spacing in the tangential direction. For the present model problem with the the pseudo-spectral approximation (21), for example, we have $|k_{x}|\leq\pi/\Delta x$, while a second-order difference approximation would roughly imply $|k_{x}|\leq 2/\Delta x$. Experience [19] shows that added-damping instabilities are generally caused by relatively high-frequency modes on the grid, and this suggests taking $k_{x}=1/h$ which leads to a definition of the fluid impedance of the form

[TABLE]

as was done in (15). The extensive numerical results in Section 7 and [15] confirm that this is an appropriate choice, and furthermore that the scheme is rather insensitive to the choice of $h$, ${\mathcal{C}}_{\rm AM}$ and ${\mathcal{C}}_{\rm AD}$.

5 Stability analysis of a viscous model problem

The stability of the AMP algorithm is explored in the context of two FSI model problems. The first model problem, discussed in this section, involves a viscous incompressible (Stokes) fluid in contact with a simplified elastic solid. This analysis extends the work in [1] to the case of a viscous fluid where both added-mass and added-damping effects are important, and for an IMEX-type scheme in the fluid. The stability analysis of a second FSI model problem involving an inviscid incompressible fluid is carried out in Section 6.

We will compare the stability of the AMP scheme to that of the traditional partitioned (TP) scheme and the anti-traditional partitioned (ATP) scheme. In the TP scheme, 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. The ATP scheme reverses the role of the solid and fluid. In this scheme, the solid provides a Neumann (traction) boundary condition for the fluid and the fluid supplies a Dirichlet (no-slip) boundary condition for the solid.

5.1 Viscous model problem

The viscous model problem analyzed here is similar to the one discussed in Section 4 and illustrated in Figure 1. An incompressible Stokes fluid satisfies the system of equations in (20a) for $\mathbf{x}\in\Omega$. The solid satisfies (20b) for $\mathbf{x}\in{\bar{\Omega}}_{0}$, but with $\bar{\lambda}$ set equal to $-\bar{\mu}$. This choice simplifies the system of equations for the solid somewhat since the compressive and shear wave speeds both equal $\bar{c}=\sqrt{\bar{\mu}/\bar{\rho}}$, i.e. the shear wave speed which is particularly relevant for the viscous model problem. It is convenient to consider the hyperbolic equations for the solid in characteristic form. These equations are

[TABLE]

where

[TABLE]

are the variables associated with the incoming $(a_{1},a_{2})$ and outgoing $(b_{1},b_{2})$ characteristics at the interface, and

[TABLE]

The interface is linearized about a flat surface $\Gamma$ given by $y=0$, and the matching conditions between the fluid and the solid are given in (20c).

5.2 Discretization

The discretization of the equations in the $x$-direction follows the approach used previously in Section 4. The equations for the fluid are transformed using the finite Fourier series in (21), which results in a system of equations for the corresponding Fourier coefficient functions given in (20a). These equations are then discretized in time using an IMEX-type scheme given by

[TABLE]

Here, $v_{1}^{n}(y)$, $v_{2}^{n}(y)$ and $p^{n}(y)$ approximate $v_{1}(y,t^{n})$, $v_{2}(y,t^{n})$ and $p(y,t^{n})$, respectively, at $t^{n}=n\Delta t$ for a fixed time step $\Delta t$. Recall that $k_{x}=2\pi k/L$, and that the hats and $k$ subscripts on the Fourier coefficients have been suppressed. Note that the components of the fluid velocity are advanced in time using (implicit) backward Euler for the viscous terms and (explicit) forward Euler for the pressure gradient terms. An elliptic equation is solved at each time step to update the pressure. It is convenient to keep the discrete equations for the fluid variables continuous in $y$, and we assume that solutions are bounded as $y\rightarrow\infty$.

The characteristic equations for the solid in (42) are similarly transformed using the finite Fourier series in (21), and then the resulting equations are discretized in time and space using an upwind-type scheme given by

[TABLE]

where, for example, $a_{1,j}^{n+1}\approx a_{1}(y_{j},t_{n})$ with $y_{j}=j\Delta y$ and $t^{n}=n\Delta t$, and where $\lambda_{x}=\bar{c}k_{x}\Delta t$ and $\lambda_{y}=\bar{c}\Delta t/\Delta y$. The grid in the $y$-direction is collocated about the interface at $y=0$. The terms involving transverse derivatives are treated implicitly to stabilize the pseudo-spectral approximation. For reference, the solid velocity and stress are related to the characteristic variables by

[TABLE]

We assume bounded solutions of (46) as $y_{j}\rightarrow-\infty$.

5.3 Interface coupling

We explore the stability of partitioned schemes for the model problem that use different interface coupling approaches. For any of the approahes, corresponding to the AMP, TP, and ATP schemes, the discrete equations require a certain number of boundary conditions at the interface. For example, the evolution of the fluid equations in (45) require three boundary conditions on the interface, $y=0$, to determine the interface velocity and pressure. Similarly, the evolution of the solid equations in (46) require two boundary conditions at $y=0$ corresponding to the two incoming characteristic variables.

We first describe the coupling based on the AMP interface conditions given in Section 3. We assume the fluid and solid solutions are known at time $t^{n}$. The solid variables are advanced first to $t^{n+1}$ on grid points $j=0,-1,-2,\ldots$ using the evolution equations in (46). The solid interface velocity and stress are computed using

[TABLE]

The fluid velocity is advanced using (45a)–(45b). Two boundary conditions are required at $y=0$ to obtain the fluid velocity at $t^{n+1}$. The condition on the outgoing solid tangential characteristic in (8b) becomes

[TABLE]

The normal component of the velocity is projected to obtain the proper limiting behaviors for heavy and light solids. This condition, taken from (19), reduces to

[TABLE]

where the fluid impedance is given by

[TABLE]

according to the derivation in Section 4. The predicted velocity, $V^{p}(v_{2}^{n+1})$, in (50) is given by

[TABLE]

This definition is analogous to that in (18), but with the substitution $\partial_{y}v_{2}^{n+1}=-ik_{x}v_{1}^{n+1}$ noting that $\nabla\cdot\mathbf{v}=0$ on the boundary. The pressure is updated using (45c) along with the AMP pressure condition described in (12). For the present scheme, this condition reduces to

[TABLE]

again using $\partial_{y}v_{2}^{n+1}=-ik_{x}v_{1}^{n+1}$. The acceleration of the solid on the interface, denoted by $\dot{\bar{v}}_{2,0}^{n+1}$ in (53), is taken to be

[TABLE]

After solving for the fluid velocity and pressure, interface quantities from the fluid are obtained using

[TABLE]

where all fluid quantities on the right-hand side are evaluated at $y=0$. The interface velocity and traction are projected from fluid and solid values using (13) and (14). These equations reduce to

[TABLE]

where $m=1,2$. Finally, the ghost points at $j=1$ for the incoming solid characteristics are set using

[TABLE]

which is a first-order accurate approximation (consistent with the order of accuracy of the upwind scheme).

We next consider the coupling conditions for the TP and ATP schemes. These conditions can be obtained from the coupling conditions for the AMP scheme in the limits of heavy ($\bar{z}\rightarrow\infty$) and light ($\bar{z}\rightarrow 0$) solids. For the AMP algorithm, the fluid velocity and pressure conditions are given in (49), (50) and (53), while the final interface values are defined by the projections in (56). For the TP algorithm ($\bar{z}\rightarrow\infty$), the AMP conditions in (49) and (50) reduce to Dirichlet conditions on the fluid velocity given by

[TABLE]

The pressure condition in equation (53) becomes a Neumann condition given by

[TABLE]

For the TP scheme, the interface velocity is taken to be the solid velocity, $v_{m}^{I}=\bar{v}_{m,0}^{n+1}$, and the interface traction is taken to be the fluid traction, $\sigma_{m2}^{I}=\sigma_{m2,f}^{n+1}$, $m=1.2$.

For the ATP scheme, we consider the light-solid limit ($\bar{z}\rightarrow 0$) of the AMP conditions. In this limit, the condition on the outgoing solid tangential characteristic in (49) reduces to a Neumann condition for the velocity given by

[TABLE]

while the condition in (50) becomes

[TABLE]

where $V^{p}(v_{2}^{n+1})$ is given by (52). Using (45b), it can be shown that the condition in (61) can be replaced by

[TABLE]

which is equivalent to setting the fluid velocity to be divergence-free on the interface. For the ATP scheme, the pressure condition in (53) reduces to

[TABLE]

For the ATP scheme, the interface velocity is taken to be the fluid velocity, $v_{m}^{I}=v_{m,f}^{n+1}$, and the interface traction is taken to be the solid traction, $\sigma_{m2}^{I}=\bar{\sigma}_{m2,0}^{n+1}$, $m=1,2$.

5.4 Stability analysis

In order to assess the stability of the AMP, TP and ATP schemes, we search for normal mode solutions to the discrete evolution equations. In the fluid, solutions are of the form

[TABLE]

where $A$ is an amplification factor. Substituting (64) into (45) and integrating gives

[TABLE]

where

[TABLE]

Here, $\Lambda$ represents the viscous CFL number and we have imposed boundedness of the solution in (65) as $y\rightarrow\infty$. The constants, $v_{m,f}^{0}$ and $p_{f}^{0}$, are obtained by imposing the appropriate boundary conditions at $y=0$, namely, (49),(50) and (53) for the case of the AMP scheme. For the TP scheme, the three constraints are the two boundary conditions for the components of the velocity in (60) and the condition on the pressure in (59), while the ATP scheme uses the boundary conditions in (60), (62) and (63).

Having found solutions for the velocity and pressure of the fluid, these solutions can be used (along with the appropriate boundary conditions at $y=0$ for the AMP, TP or ATP coupling) to eliminate the fluid variables on the boundary in favor of the solid variables. The issue of stability, then, reduces to examining the behavior of the evolution equations for the solid with the appropriate boundary conditions. Solutions of these evolution equations are sought in the form

[TABLE]

where $\phi$ is a spatial eigenvalue and $\tilde{\mathbf{r}}$ is a constant eigenvector. The scheme is said to be weakly stable if there are no non-trivial solutions with $|A|>1$. Our strategy for determining regions of stability will be to search for unstable modes with $|A|>1$, and then identify regions of the parameter space where no non-trivial solutions exist. To do this, we begin by finding the general solution for the spatial grid functions satisfying the discrete equations and regularity condition as $j\rightarrow-\infty$, assuming $|A|>1$. We then apply the conditions at the interface to determine whether non-trivial solutions exist.

When the normal mode ansatz in (67) is substituted into the evolution equations for the solid in (46), a homogeneous system arises of the form

[TABLE]

The matrix $\mathcal{F}$ is given by

[TABLE]

where

[TABLE]

The determinant of $\mathcal{F}$ is given by

[TABLE]

The system is singular when $f(\phi)=0$, which occurs when $A=1$ or $\eta(\phi)\eta(1/\phi)+(A\lambda_{x})^{2}=0$. Since we are searching for unstable modes with $|A|>1$, we are interested in the latter case which leads to roots given by

[TABLE]

The product of the roots is equal to one ($\phi_{+}\phi_{-}=1$). Since we are searching for solutions that are bounded as $j\rightarrow-\infty,$ we are only interested in the root with modulus greater than one.

Lemma: If $|A|>1$ and if $\lambda_{x}$ and $\lambda_{y}$ are chosen to satisfy a CFL condition, then there is precisely one root, either $\phi_{+}$ or $\phi_{-}$, denoted by $\phi_{*}$ that has magnitude strictly greater than one, i.e. $|\phi_{*}|>1.$

This result follows from a similar argument to that given in [20]. We first consider the scheme applied to the pure initial-value problem (Cauchy problem). Setting $\phi=e^{i\vartheta}$ in (70), we determine a region of the $(\lambda_{x},\lambda_{y})$ plane for which $|A|\leq 1$ for all $\vartheta\in[0,2\pi]$. This region is found numerically as discussed in Section A.1 and shown in Figure 14, and it corresponds to a region where the Cauchy problem is (Fourier) stable. It is also shown that this stability region includes a region of $(\lambda_{x},\lambda_{y})$ satisfying a reasonable CFL restriction, namely $\lambda_{x}^{2}+\lambda_{y}^{2}\leq 1$. Next, since $|A|\leq 1$ when $|\phi|=1$, we have that $|\phi|\neq 1$ when $|A|>1$. Thus, if $|A|>1$ and if $(\lambda_{x},\lambda_{y})$ remains within the CFL restriction, then $\phi$ cannot cross the unit circle, $|\phi|=1$, as $(\lambda_{x},\lambda_{y})$ vary. It is therefore only necessary to prove that the lemma holds for one set of parameters. For $\lambda_{x}=0$, the discretization reduces to four uncoupled upwind schemes for linear advection. In this case, equation (70) is equivalent to $\eta(\phi)\eta(1/\phi)=0,$ which has solutions $\phi_{+}=(A-1+\lambda_{y})/\lambda_{y}$ and $\phi_{-}=1/\phi_{+}$. When $|A|>1$ and $\lambda_{y}\in(0,1]$, $|\phi_{+}|>1$ and therefore $\phi_{*}=\phi_{+}.$ Thus, the condition holds for all $(\lambda_{x},\lambda_{y})$ provided the CFL condition is satisfied.

When $\phi=\phi_{*},$ there are two eigenvectors, namely

[TABLE]

The solution which remains bounded as $j\rightarrow-\infty$ is given by

[TABLE]

where $k_{1}$ and $k_{2}$ are constants to be determined by the two interface conditions in (57). The application of these interface conditions leads to another homogeneous system of equations given by

[TABLE]

The coefficients of the matrix $\mathcal{G}$ are provided in Section A.2 for the AMP, TP, and ATP schemes. Solutions for the amplification factor $A$ are roots of the transcendental equation given by

[TABLE]

These roots dependent on the choice of interface coupling and four dimensionless parameters ($\Lambda,Z,\lambda_{x},\lambda_{y}$), where

[TABLE]

Proving stability of the partitioned scheme for a choice of the interface coupling and dimensional parameters is equivalent to showing that no roots of (78) exist such that $|A|>1$. The number of roots with $|A|>1$ can be assessed using the argument principle. Define

[TABLE]

There are branch points of ${G}(\zeta)$ in the region $|\zeta|>1$, and a single-valued branch of ${G}(\zeta)$ can be defined so that its branch cuts lie outside the unit disk. The only singularity of ${G}(\zeta)$ in the region $|\zeta|\leq 1$ is a pole of order $2$ at the origin, and thus $\mathcal{P}=N-2$, where $N$ corresponds to the number of roots of $g(A)$ with $|A|>1$.

An analytic evaluation of the integral for $\mathcal{P}$ is unavailable, and so we consider a numerical evaluation. The four-dimensional parameter space ($Z$,$\Lambda$,$\lambda_{x}$,$\lambda_{y}$) is discretized on a 31$\times$31$\times$20$\times$20 array. The parameters $Z$ and $\Lambda$ are equally spaced on a logarithmic scale on the interval $[10^{-3},10^{3}]$, while $\lambda_{x}$ and $\lambda_{y}$ are equally spaced on the interval $[0.05,0.95]$. At each grid point, $\mathcal{P}$ is computed numerically with $|\mathcal{P}+2|\leq\delta$ corresponding to stability, where $\delta$ is a small parameter taken to be $10^{-5}$. The results of this computation are shown in Figure 3 for the AMP, TP and ATP schemes. A grid point in the $(\Lambda,Z)$ plane is marked as stable if the computations of $\mathcal{P}$ for all values of $\lambda_{x}$ and $\lambda_{y}$ in the search region are stable. The point is marked as unstable otherwise. The results shown in the left plot indicate that the AMP scheme applied to the viscous model problem is stable for all points in the $(\Lambda,Z)$ plane, whereas the results shown in the right plot indicate that the TP and ATP schemes have large regions of instability. For example, the region in red shows the stable region for the TP scheme, which occurs for heavy solids ($Z$ small) and coarser meshes ($\Lambda$ large). The stability region for the ATP scheme shown in blue corresponds to light solids ($Z$ large). The following theorem summarizes the results for the AMP scheme.

Theorem: The AMP scheme applied to viscous model problem is weakly stable $|A|\leq 1$ provided $\lambda_{x}^{2}+\lambda_{y}^{2}\leq 1$, which gives the usual CFL-type time-step restriction

[TABLE]

This is a sufficient but not a necessary condition. The proof follows from the argument principle and a numerical evaluation of $\mathcal{P}$.

6 Stability analysis of an inviscid model problem

In this section, the stability of the AMP, TP and ATP algorithms are explored for a two-dimensional FSI model problem involving an inviscid incompressible fluid. A similar stability analysis was performed in [1] for a first-order accurate scheme on a staggered grid. The algorithms considered here are second-order accurate, as in the full AMP algorithm discussed in Section 3, and thus the present analysis extends the previous results for less dissipative partitioned schemes.

6.1 Inviscid model problem

The fluid is assumed to be inviscid and incompressible, and it occupies the region $\Omega$ given by $0<x<L$ and $0<y<H$. The solid lies below the fluid in the semi-infinite region ${\bar{\Omega}}_{0}$ given by $0<x<L$ and $y<0$, see Figure 4. The solid is taken to be an “acoustic” solid that supports a motion in the vertical direction only. The fluid-solid interface $\Gamma$ is linearized about the flat surface, $y=0$, and it is assumed that the vertical components of the velocity and stress between the fluid and solid match along $\Gamma$. The fluid pressure is taken to be zero along the top boundary, $y=H$, and the solutions for the fluid and solid are assumed to be periodic in the $x$-direction with period equal to $L$. The governing equations for this FSI model problem are

[TABLE]

We note that the horizontal component of the fluid velocity is omitted from the equations for the fluid in (82a) since it decouples in the simplified FSI problem and can be determined from the remaining fluid variables once the solution is known. The solid variables are assumed to be bounded as $y\rightarrow-\infty$. The model problem described here corresponds to the MP-IA model problem considered in [1].

6.2 Discretization

The equations for the fluid in (82a) are transformed to Fourier space using (21). The resulting equations for the Fourier coefficient functions are

[TABLE]

The solution of the elliptic equation for the pressure in (83b), with $p=0$ at $y=H$, is

[TABLE]

where $p_{I}(t)$ is an interface pressure which is specified later by the chosen coupling at the interface. For the purpose of the analysis, we only require the fluid velocity on the interface, which we define to be $v_{I}(t)$. The evolution of the interface velocity involves the gradient of pressure at $y=0$, which is given by

[TABLE]

The effective fluid height, $H_{\text{eff}}$, in the formula for the pressure gradient depends on the normalized wave number, and it takes the limiting values of $H_{\text{eff}}=H$ as $k_{x}\rightarrow 0$ and $H_{\text{eff}}=0$ as $k_{x}\rightarrow\infty$.

The equations for the solid in (82b) are also transformed to Fourier space, and then expressed in the characteristic form

[TABLE]

where $a(y,t)$, $b(y,t)$ and $d(y,t)$ are Fourier coefficients of the characteristic variables

[TABLE]

Equations (86a) and (86b) are advection equations for the incoming and outgoing characteristic variables $a$ and $b$, respectively, while (86c) describes the evolution of the tangential component of the solid stress.

The solid characteristic variables are discretized in the $y$-direction using a uniform collocated grid, $y_{j}=j\Delta y$ for $j=0,-1,-2\ldots$, with grid spacing $\Delta y$. Although a staggered grid is used for the analysis in [1], and for the stability analyses of other FSI algorithms in [21, 22], it is cleaner to use a collocated grid for the second-order accurate scheme examined here. Define the grid functions $b^{n}_{j}\approx b(y_{j},t^{n})$, $a^{n}_{j}\approx a(y_{j},t^{n})$ and $d^{n}_{j}\approx d(y_{j},t^{n})$, where $t^{n}=n\Delta t$ for a (fixed) time step $\Delta t$. The advection equations in (86a) and (86b) are approximated using a second-order accurate Lax-Wendroff-type scheme having the form

[TABLE]

where $\lambda_{x}=k_{x}\bar{c}_{p}\Delta t$ and $\lambda_{y}=\bar{c}_{p}\Delta t/\Delta y$. The centered approximations of the spatial derivatives in (88) are defined in terms of the undivided difference operators $\delta_{+}u^{n}_{j}=u^{n}_{j+1}-u^{n}_{j}$, $\delta_{-}u^{n}_{j}=u^{n}_{j}-u^{n}_{j-1}$ and $\delta_{0}u^{n}_{j}=u^{n}_{j+1}-u^{n}_{j-1}$ for an arbitrary grid function $u^{n}_{j}$. The evolution equation in (86c) is approximated using a second-order accurate BDF scheme of the form

[TABLE]

For reference, the Fourier coefficients of the solid velocity and normal stress are related to the incoming and outgoing characteristic variables by

[TABLE]

and the boundedness condition at infinity implies that

[TABLE]

6.3 Interface coupling

The evolution of the solid requires a boundary condition on the interface, $y=0$, corresponding to the incoming characteristic variable. Similarly, the evolution of the fluid requires a boundary condition at $y=0$, which can be interpreted as an interface pressure needed to complete the solution of the elliptic problem for the fluid pressure in (84). These two boundary conditions come from the matching conditions at the interface involoving the vertical velocity and normal stress. In terms of the discretization of the Fourier coefficents for the model problem, the coupling at the interface takes on different forms depending on the choice of the partitioned scheme. Since the behavior of the fluid for the model problem is determined by the solution of an elliptic problem for the pressure and an evolution equation for the vertical velocity on the interface, the coupling of the fluid to the behavior of the solid can be reduced to a modified boundary condition for the solid depending on the partitioned scheme employed as we discuss below.

Let us first consider the AMP algorithm and its interface coupling scheme. The general algorithm is described in Section 3 and we follow this basic description but with suitable modifications for the choices of the discretizations made for the present model problem. Let us assume that the discrete characteristic variables for the solid given by $a^{n}_{j}$, $b^{n}_{j}$ and $d^{n}_{j}$ are known at time levels $t^{n}$ and $t^{n-1}$ for all grid points, $y_{j}=j\Delta y$, $j=0,-1,\ldots$, including a ghost point at $j=1$. We also assume that the interface velocity, $v_{I}(t)$, is known at $t=t^{n}$ and $t^{n-1}$. The evolution equations for the solid variables in (88) and (89) are used to advance the variables to $t^{n+1}$, and from these we compute the solid velocity and normal stress on the interface using

[TABLE]

We next compute an extrapolated interface velocity, $v_{I}^{(e)}$, using a second-order accurate BDF-type integration of the evolution equation for the fluid velocity. This gives the formula

[TABLE]

where the pressure gradient on the interface is given in (85). This extrapolated value for the fluid is averaged with the interface velocity for the solid from (92) to give

[TABLE]

where we take $z_{f}=\rho H/\Delta t$. The choice for $z_{f}$ here is the same one used in [1], and it agrees with (15) for ${\mathcal{C}}_{\rm AM}h=H$. We have found that the results are insensitive to the choice of $z_{f}$ for this inviscid model problem. For later convenience, define the weights

[TABLE]

and we note that $\theta_{f}+{\bar{\theta}}_{p}=1$. The AMP condition for the fluid pressure given in (12) requires an acceleration of the interface, and for this we use the backward difference formula

[TABLE]

Using (93) in (94) to eliminate $v_{I}^{(e)}$, and then using the computed result in (96) to eliminate $v_{I}^{n+1}$ gives an alternate expression for the acceleration having the form

[TABLE]

which shows that the acceleration of the interface is also a weighted average of the acceleration from the fluid and that computed from the solid. Finally, the AMP condition for the fluid pressure in (12) reduces to

[TABLE]

Using (97) in (98) to eliminate $\dot{v}_{I}^{n+1}$, and then solving for the interface pressure gives

[TABLE]

where $M$ is a mass ratio given by

[TABLE]

We may now complete the time-step for the solid variables by specifying values for the incoming and outgoing characteristic variables in the ghost points using

[TABLE]

where $a_{I}^{n+1}$ is an interface value for the incoming characteristic variable determined by

[TABLE]

Thus, the boundary condition for the incoming characteristic variable is specified by $a_{I}^{n+1}$ (using a second-order accurate average) and the outgoing characteristic variable at the ghost point is given by a simple extrapolation.

For the AMP algorithm, the interface velocity and pressure used in (102) to define the interface value for the incoming characteristic variable are given by the formulas in (94) and (99), respectively, while the corresponding interface values for the TP and ATP algorithms use difference formulas. For the TP algorithm, the interface velocity is given by the solid, while the interface stress is determined by the fluid. This specification leads to

[TABLE]

where the solid velocity is given in (92) and the interface acceleration is given in (96). For the ATP algorithm, the interface velocity is given by the fluid, while the interface stress is determined by the solid, and this gives

[TABLE]

where the fluid velocity is given in (93) and the solid stress is given in (92). It is worth noting that values specified for the TP and ATP algorithms in (103) and (104), respectively, agree with values given by the AMP algorithm in the limiting cases of a heavy solid ($\bar{z}_{p}\rightarrow\infty$) and a light solid ($\bar{z}_{p}\rightarrow 0$).

6.4 Stability analysis

In order to assess the stability of the AMP, TP and ATP schemes, we consider normal mode solutions of the discrete equations for the solid with the appropriate interface conditions. For this analysis, we let

[TABLE]

where $A$ is an amplification factor and $(\tilde{a}_{j},\tilde{b}_{j},\tilde{d}_{j})$ are spatial grid functions. As in stability analysis of Section 5, showing the scheme is weakly stable is equivalent to showing there are no non-trivial solutions where $|A|>1$ for a given set of parameters. We first make the assumption $|A|>1$ and solve the spatial problem by finding general solutions of the discretization and boundedness condition as $j\rightarrow-\infty.$ The interface conditions are then used to determine if the amplification factor satisfies the assumption $|A|>1.$

The discrete equations for the solid in (88) and (89), with $\tilde{\mathbf{b}}_{j}=[\tilde{a}_{j},\tilde{b}_{j},\tilde{d}_{j}]^{T}$, can be written as a system of linear, second-order difference equations of the form

[TABLE]

where

[TABLE]

Solutions of the difference equations in (106) are sought in the form $\tilde{\mathbf{b}}_{j}=\tilde{\mathbf{r}}\phi^{j}$, where $\tilde{\mathbf{r}}$ is a constant vector and $\phi$ is a scalar. This implies the homogeneous linear system

[TABLE]

Non-trivial vectors $\tilde{\mathbf{r}}$ exist for values of $\phi$ satisfying

[TABLE]

It can be shown that the algebraic condition in (108) is equivalent to a quartic polynomial in $\phi$ for which there are four roots, $\phi_{\alpha}$, $\alpha=1,2,3,4$. Since we seek spatial grid functions that are bounded as $j\rightarrow-\infty$, we are only interested in the roots satisfying $|\phi_{\alpha}|\geq 1$.

Lemma: If $|A|>1$ and if $\lambda_{x}$ and $\lambda_{y}$ are chosen to satisfy a CFL condition, then the roots of the quartic polynomial implied by (108) satisfy

[TABLE]

The analogous lemma for the viscous stability analysis was proven in Section 5. The following proof is structured similarly. First, the stability for the scheme applied to the pure initial-value problem (Cauchy problem) is considered. This corresponds to setting $\phi=e^{i\vartheta}$ in (108) and finding the region of stability $|A|\leq 1$ in the $\lambda_{x}$-$\lambda_{y}$ plane for all $\vartheta\in[0,2\pi].$ This region is found in Section A.1 and shown in Figure 15. This region contains the quarter circle defined by $\lambda_{x}^{2}+\lambda_{y}^{2}\leq 1$ and $\lambda_{x},\lambda_{y}\geq 0.$ When $|\phi|=1,$ we must have $|A|\leq 1.$ Therefore $|A|>1$ implies $|\phi|\neq 1.$ Due to the continuity of $\phi,$ if $|A|>1$ and if $(\lambda_{x},\lambda_{y})$ remains within the CFL restriction, $\phi$ cannot cross the unit circle, $|\phi|=1$, as $(\lambda_{x},\lambda_{y})$ vary. Therefore, it is only necessary to prove that the conditions on the roots in (109) hold for one set of parameters. For $\lambda_{x}=0$, the difference scheme reduces to a pair of decoupled Lax-Wendroff schemes (and a BDF scheme applied to $\partial_{t}d=0$), where it is well known (see [23] for example) that the conditions in (109) hold provided $|\lambda_{y}|\leq 1$. Thus, the condition holds for all $(\lambda_{x},\lambda_{y})$ provided the CFL condition is satisfied.

In view of the above Lemma, we can now write the general solution for $\mathbf{b}^{n}_{j}=[a^{n}_{j},b^{n}_{j},d^{n}_{j}]^{T}$, satisfying the regularity condition and assuming $|A|>1$, in the form

[TABLE]

where $k_{1}$ and $k_{2}$ are scalar constants. Let $\tilde{\mathbf{r}}_{\alpha}=[\tilde{q}_{\alpha},\tilde{r}_{\alpha},\tilde{s}_{\alpha}]^{T}$ for $\alpha=1,2$. We note that the third equation in the matrix system (107) implies

[TABLE]

so that the third component of $\mathbf{b}^{n}_{j}$ can be eliminated in terms of the first two components. The final two constraints are given by the interface conditions in (101), and this leads to a homogeneous system for $\mathbf{k}=[k_{1},k_{2}]^{T}$ of the form

[TABLE]

The coefficients of the matrix $\mathcal{H}$ involve the amplitude factor $A$, and these are given in Section A.3 for the AMP, TP and ATP schemes. For non-trivial solutions, we seek solutions for $A$ such that

[TABLE]

We first proceed with an analysis of the 1D case assuming $\lambda_{x}=0$. After finding solutions for $A$ when $\lambda_{x}=0$, we use a continuation procedure to determine solutions for $A$ numerically as $\lambda_{x}$ increases from zero.

6.4.1 1D schemes ($\lambda_{x}=0$)

As noted earlier, the discrete evolution equations in (88) for the solid reduce to a pair of standard Lax-Wendroff schemes (without source terms) for the variables associated with incoming and outgoing characteristics when $\lambda_{x}=0$. Accordingly, the first two components of (106) become

[TABLE]

Setting $\tilde{a}_{j}=\tilde{a}_{0}\phi_{a}^{j}$ and $\tilde{b}_{j}=\tilde{b}_{0}\phi_{b}^{j}$ yields the eigenvalues

[TABLE]

If $|A|>1$ and if $\lambda_{y}$ satisfies the CFL constraint $|\lambda_{y}|\leq 1$, then one eigenvalue from each plus-minus pair in (119) has magnitude greater than one, while the other has magnitude less than one in agreement with (109). Let $\phi_{1}$ to be the eigenvalue from $\phi_{a,\pm}$ with magnitude greater than one and let $\phi_{2}$ to be the corresponding eigenvalue $\phi_{b,\pm}$. Assuming the principal branch of the square-root function in (119), we set

[TABLE]

Since the evolution equations for $a_{j}^{n}$ and $b_{j}^{n}$ are decoupled, the corresponding eigenvectors in (110) are unit vectors, $\tilde{\mathbf{r}}_{1}=[1,0,0]^{T}$ and $\tilde{\mathbf{r}}_{2}=[0,1,0]^{T}$, i.e. $\tilde{q}_{2}=\tilde{r}_{1}=\tilde{s}_{1}=\tilde{s}_{2}=0$, and the matrix $\mathcal{H}$ in (116) becomes diagonal. Thus, the determinant condition in (117) reduces to

[TABLE]

The component, $\mathcal{H}_{22}$, is associated with the extrapolation of the outgoing characteristic variable from (101), and this factor is zero only if $\phi_{2}=1$. However, $|\phi_{2}|>1$, so that $\mathcal{H}_{22}\neq 0$ and thus $k_{2}=0$ in (110). The remaining component, $\mathcal{H}_{11}$, is associated with the interface condition on the incoming characteristic variable, and this term takes different forms depending on the choice of the interface coupling (see Section A.3). Manipulation of the constraint, $\mathcal{H}_{11}=0$, for the AMP, TP and ATP schemes leads to polynomials in $A$ of degrees 6, 6 and 4, respectively, whose coefficients depend on $\lambda_{y}$ and $M_{0}=\rho H/(\bar{z}_{p}\Delta t)$. Due to the algebraic complexity, the software package Maple is used to generate the polynomials and compute the roots for given values of $(\lambda_{y},M_{0})$. We note that some roots of the polynomials may not be solutions of the determinant constraint and so these spurious roots are discarded. If no valid roots are found satisfying $|A|>1$ for a given value of $(\lambda_{y},M_{0})$, then we conclude that $k_{1}$ in (110) must also be zero so that only the trivial solution exists and thus the scheme is (weakly) stable.

We find it convenient to explore the roots of the three polynomials in terms of the parameters $\lambda_{y}$ and $\mathcal{M}$, where

[TABLE]

The latter parameter is the ratio of the mass contained in a fluid column of height $H$ to that of the solid with the same cross section and of height $\Delta y$. Starting with the AMP scheme, we compute valid roots of its polynomial for an $800\times 800$ array of parameter values on a logarithmic scale for the range $10^{-6}\leq\lambda_{y}\leq 1$ and $10^{-6}\leq\mathcal{M}\leq 10^{7}$. For these values, no roots with $|A|>1$ are found, which provides a strong indication of stability of the 1D AMP scheme for $0<\lambda_{y}\leq 1$ and $\mathcal{M}>0$. Shaded contours of the maximum values of $|A|$ are displayed in the left plot of Figure 5 for a range of values of $\mathcal{M}$ and $\lambda_{y}$. The roots of the polynomials associated with the 1D TP and ATP schemes are also computed, and parameter values for which $|A|>1$ are found indicating that regions of instability exist for these two schemes. The right plot of Figure 5 shows regions of stability of the TP and ATP schemes where no roots with $|A|>1$ are found. We note that for a fixed value of $\lambda_{y}=\bar{c}_{p}\Delta t/\Delta y\in(0,1]$, decreasing the mesh spacings in both time and space corresponds to increasing the mass ratio $\mathcal{M}$ defined in (122). The right plot in the figure indicates that the TP scheme becomes unstable as the mesh is refined, while the ATP scheme ultimately becomes stable. In contrast, the AMP scheme is always stable, and these observations agree with the analysis of the first-order accurate difference schemes performed in [1].

The AMP, TP and ATP schemes can be implemented for the one-dimensional FSI model problem, and then run numerically for different values of the parameters to check the results of the stability analysis. Initial conditions are specified so that the exact solution decays exponentially in time, see Section A.4 for details. Since the interface velocity is expected to decay in time, we can assess the stability of the scheme by checking the magnitude of the interface velocity. If this velocity grows in time, becoming larger than the initial interface velocity, then we consider the scheme to be unstable for the choice of parameter values. For the AMP scheme, no parameter values are found for which the numerical results are unstable. For the TP and ATP schemes, the colored marks in Figure 6 indicate stability or instability of the numerical scheme for selected parameter values. The black curves in the figure correspond to the stability boundaries obtained from the stability analysis. We note that there is good agreement between the numerical results of the schemes and that given by the analysis.

6.4.2 2D schemes ($|\lambda_{x}|\neq 0$)

The strategy for finding solutions of the determinant condition in (117) for $|\lambda_{x}|\neq 0$ is to use the solutions for $A$ of the 1D determinant condition as a starting point for a numerical root solver that performs a continuation as $\lambda_{x}$ varies. In addition to the principal continuation parameter $\lambda_{x}$, the roots also depend on the parameters $\lambda_{y}$ and $M=\rho H_{\text{eff}}/(\bar{z}_{p}\Delta t)$, where $H_{\text{eff}}$ was defined previously in (85). As before, we find it convenient to examine the roots in terms of the mass ratio $\mathcal{M}$ introduced in (122) and the height ratio, $\eta=H_{\text{eff}}/H$, which varies on $(0,1]$, and we note that $M=\mathcal{M}\eta$. Using our continuation procedure, we are able to compute the roots of the determinant condition as a function of the parameters $\lambda_{x}$, $\lambda_{y}$ and $M=\mathcal{M}\eta$. Though there are multiple roots, $A_{z}$, ${z}=1,2,\ldots,{\mathcal{Z}}$, for a given set of parameters, we are only interested in the amplification factor with the maximum modulus. Define

[TABLE]

to be the maximum amplification factor in magnitude for a given set of parameters. The pseudo-spectral approximation of the model problem given in (82b) depends on the normalized wave number $k_{x}$. The corresponding discrete evolution equations for each Fourier mode of the approximation must all be stable for the chosen partitioned scheme to be stable (either the AMP, TP or ATP schemes). Since $\lambda_{x}$ and $\eta$ both depend on $k_{x}$, we are interested in the values of $A_{\max}$ for all possible values for $\lambda_{x}$ and $\eta$. Define

[TABLE]

where $\mathcal{D}(\lambda_{y})$ gives the interval of $\lambda_{x}$ satisfying the Cauchy stability bound for a given value of $\lambda_{y}$ (see Section A.1). We take $\mathcal{D}(\lambda_{y})=[0,(1-\lambda_{y}^{2})^{1/2}]$. Since the maximum amplification factor for the 1D case (with $\lambda_{x}=0$) must be less than or equal to $\mathcal{A}_{\max}$, we are effectively examining whether the maximum amplification can increase when $|\lambda_{x}|>0$. Figure 7 shows shaded contours of $\mathcal{A}_{\max}(\mathcal{M},\lambda_{y})$ on the left, and the stability regions of the TP and ATP schemes on the right. These plots are similar to the ones for the 1D analysis in Figure 5, indicating that the stability results provided by the 1D analysis are essentially unchanged when variations in the transfer direction are considered. In fact, Figure 8 shows the stability regions for the TP and ATP schemes given by the 1D and 2D analyses overlayed. Here, we observe that the stability regions given by the 2D analysis are only slightly smaller than those given by the 1D analysis.

We can also search for a maximum amplification factor over all mass ratios, $\mathcal{A}_{\text{CFL}},$ to see the CFL region for the AMP scheme in $(\lambda_{x},\lambda_{y})$ plane. Define

[TABLE]

This definition allows us to find a time step restriction regardless of the mass ratio of the problem. Figure 9 shows a surface plot of $\mathcal{A}_{\text{CFL}}$ and the stability region, $|A|\leq 1$ in the $\lambda_{x}$-$\lambda_{y}$ plane.

Theorem: By numerical evaluation of $\mathcal{A}_{\text{CFL}}(\lambda_{x},\lambda_{y})$, we have found that the AMP algorithm is weakly stable, i.e. $|\mathcal{A}_{\text{CFL}}|\leq 1$, provided that $\lambda_{x}^{2}+\lambda_{y}^{2}\leq 1$. This implies the time-step restriction

[TABLE]

This is a sufficient but not a necessary condition.

7 Numerical results for an elastic piston

We now present numerical results for two FSI problems to verify the accuracy and stability of the AMP scheme. The two FSI problems considered involve the interaction of a fluid column with an elastic piston. In the first problem, we examine longitudinal motion of the piston, while transverse motion of the piston is considered in the second problem. Exact solutions are found for both FSI problems, and these are used to verify the accuracy and stability of the AMP algorithm for a range of the problem parameters.

7.1 Longitudinal motion of an elastic piston

The geometry of the elastic piston problem is shown in Figure 10. The plot on the left shows the configuration at $t=0$. The fluid occupies the physical domain between $y=0$ and $y=H$ initially, while the solid lies in its reference domain between ${\bar{y}}=-\bar{H}$ and ${\bar{y}}=0$. It is assumed that there is no dependence in the $x$-direction so that the fluid-solid interface remains flat at a position $y=y_{I}(t)$ as shown in the plot on the right. In the fluid domain, $\Omega(t)$, it is assumed that the horizontal component of velocity $v_{1}$ is zero, and thus the vertical component $v_{2}$ is a function of $t$ alone according to the continuity equation. The fluid pressure is a linear function $y$, and is given by

[TABLE]

where $p_{I}(t)$ is the pressure on the interface and $p_{H}(t)$ is a specified fluid pressure at $y=H$. The momentum equation for the fluid in the vertical direction reduces to

[TABLE]

where $\rho$ is the (constant) fluid density.

With the displacement of the solid in the horizontal direction assumed to be zero, the equation for the vertical component of the displacement in the reference domain, $\bar{\Omega}(0)$, becomes

[TABLE]

where $\bar{c}_{p}$ is the longitudinal wave speed. The general solution of (129), assuming a zero-displacement condition at ${\bar{y}}=-\bar{H}$, is

[TABLE]

where $\bar{F}(\bar{\tau})$ is any smooth function. The vertical position of the fluid-solid interface is determined by the solid displacement evaluated at ${\bar{y}}=0$. Thus,

[TABLE]

and the matching conditions on velocity and stress imply

[TABLE]

The interface conditions in (131) and (132) can be used in (128) to eliminate the fluid variables, which leads to a boundary condition involving the solid displacement and its derivatives. This boundary condition can be written in the form

[TABLE]

We are interested in finding a function $\bar{F}(\bar{\tau})$ in (130) so that $\bar{u}_{2}({\bar{y}},t)$ satisfies initial conditions on $\bar{u}_{2}$ and $\bar{u}_{2,t}$, and also satisfies the nonlinear boundary condition in (133) for a specified fluid pressure $p_{H}(t)$. While this could be done in principle, we choose instead to construct an exact solution by making a choice for $\bar{F}(\bar{\tau})$, and then backing out the corresponding initial conditions on $\bar{u}_{2}({\bar{y}},t)$ and the fluid pressure $p_{H}(t)$. This is done in our verification tests using the choice

[TABLE]

where $\alpha$ and $\omega$ are parameters. With $\bar{u}_{2}({\bar{y}},t)$ and $p_{H}(t)$ known from (130) and (133), respectively, the remaining fluid variables can be obtained from (127), (131) and (132). In particular, we note that for this choice, the position of the fluid-solid interface oscillates with frequency $\omega$ and an amplitude $a$ given by

[TABLE]

Numerical results are obtained for the case $H=1,\rho=1$ and $\mu=0.01$ for the fluid, and using $\bar{H}=0.5$ and $\bar{\mu}=\bar{\lambda}=\bar{\rho}=\delta$ for the solid. The interface position is specified by $a=0.1$ and $\omega=2\pi$. The density ratio, $\bar{\rho}/\rho=\delta$, is taken to be $10^{-3}$, 1 and $10^{3}$, representing FSI problems with light, moderate and heavy solids, respectively. Numerical solutions are computed using the AMP algorithm on a two-dimensional rectangular configuration (as shown in Figure 10) with periodic boundary conditions taken at $x=0$ and $x=L$ consistent with a one-dimensional solution. Table 1 gives the maximum-norm errors for solutions computed using the AMP algorithm at $t_{{\rm final}}=0.6$ with grid resolutions $h=1/(20j)$ for $j=1,2,4,8.$ The errors in the table indicate that the solution is converging at second-order accuracy. The solution is plotted in Figure 11 for the case $\delta=1$.

7.2 Transverse motion of an elastic piston

Exact solutions can also be constructed for an FSI problem involving transverse motion of an elastic piston, see Figure 12. For this case, the vertical components of the fluid velocity and solid displacement are taken to be zero, while the corresponding horizontal components are assumed to be functions of $y$ and $t$ alone. As a result, the interface remains in its initial position, $y_{I}(t)=y_{I}(0)=0$, and so the solid reference coordinate ${\bar{y}}$ is equivalent to the physical coordinate $y$.

For this problem, the equations governing the horizontal components of the fluid velocity and solid displacement reduce to

[TABLE]

and

[TABLE]

respectively, where $\nu=\mu/\rho$ is the kinematic viscosity of the fluid and $\bar{c}_{s}$ is shear wave speed of the solid. Solutions of (134) and (135) are sought in the form

[TABLE]

where $\omega$ is a parameter (ultimately an eigenvalue). The coefficient functions, $\hat{v}_{1}(y)$ and $\hat{\bar{u}}_{1}(y)$, in (136) satisfy

[TABLE]

where

[TABLE]

Solutions of the second-order ODEs in (137) satisfying a no-slip condition in the fluid at $y=H$ and a zero-displacement condition in the solid at $y=-\bar{H}$ are

[TABLE]

where $(b,\bar{b})$ are constants. The interface matching conditions on velocity and shear stress imply

[TABLE]

Using (139) in (140) leads to the system

[TABLE]

Nontrivial solutions for $(b,\bar{b})$ exist if the determinant of the coefficient matrix in (141) is zero, which implies

[TABLE]

Assuming a value for $\omega$ is found satisfying (142), the constants $b$ and $\bar{b}$ are given by

[TABLE]

where $\bar{u}_{0}$ is the interface displacement at $t=0$. The fluid pressure is constant and equal to zero due to the matching of the normal component of stress at the interface ($\sigma_{22}=\bar{\sigma}_{22}=0$). The real part of (136) is used as the solution.

Values of $\omega$ with $\Re(\omega)>0$ satisfying $\mathcal{D}(\omega)=0$ correspond to solutions of the transverse elastic piston problem that decay in time. Selected values of $\omega$ are listed in Table 2 for $H=1$, $\bar{H}=0.5$, $\rho=1$ and $\mu=0.1$, and for different values of $\delta=\bar{\rho}=\bar{\mu}=\bar{\lambda}$. Table 3 gives the maximum-norm errors for solutions computed using the AMP algorithm for $\bar{u}_{0}=0.1$. The results are presents for solutions at $t_{{\rm final}}=0.3$ using grid resolutions $h=1/(20j)$, for $j=1,2,4,8.$ The errors in the table indicate that the solution is converging at second-order accuracy. The solution is plotted in Figure 13 for the case $\delta=1$.

8 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. The fluid is advanced using a fractional-step IMEX scheme with the viscous terms treated implicitly. 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 fluid impedance is derived from an analysis of a solution to an FSI problem near the interface. Stability of the AMP scheme is analyzed for two model problems, one involving a first-order accurate discretization of a viscous model problem and a second-order accurate discretization of an inviscid model problem. The elastic-piston, a set of new benchmark problems, was developed to verify stability and accuracy of the AMP scheme. These solutions are exact and include finite interface deformations either normal or tangent to the surface.

Appendix A Results supporting the stability analysis of the model problems

The subsections below give additional results used for the stability analysis of the viscous and inviscid model problems of Sections 5 and 6, respectively.

A.1 CFL region for Cauchy problem

The CFL stability regions are describedfor the discretizations of the model problems discussed in Sections 5 and 6 applied to the pure initial-value problem (Cauchy problem). For both discretization, we assume that the spatial eigenfunction is given by $\phi=e^{i\vartheta}$ for all $\vartheta\in[0,2\pi]$, i.e. $|\phi|=1$. For the viscous model problem, the amplification factor $A$ can be found from (70) for a given value of $\phi$, while (108) gives $A$ in terms of $\phi$ for the inviscid model problem. The CFL stability region is defined as the region in the $(\lambda_{x},\lambda_{y})$ plane for which the amplification factor satisfies $|A|\leq 1$. This region can be found numerically by computing the $A_{{\rm max}}=\max_{|\phi|=1}|A|$ for a range of values in the $(\lambda_{x},\lambda_{y})$ plane. The left plot of Figure 14 shows shaded countours of $A_{{\rm max}}$, while the plot on the right indicates the region for which $A_{{\rm max}}\leq 1$. Note that the constraint $\lambda_{x}^{2}+\lambda_{y}^{2}\leq 1$ is sufficient for CFL stability. Similar results are shown in Figure 15 for the analysis of the inviscid model problem.

A.2 Matrix system in the stability analysis of the viscous model problem

The stability analysis of the viscous model problem in Section 5 led to the matrix system

[TABLE]

where the components $\mathcal{G}_{ij}$ depend on the choice of coupling at the interface for the partitioned scheme. For the AMP scheme, the components can be written in the form

[TABLE]

where

[TABLE]

For the TP scheme, the components of the linear system have the form

[TABLE]

The components of the linear system for the ATP scheme are

[TABLE]

A.3 Matrix system in the stability analysis of the inviscid model problem

The stability analysis of the inviscid model problem in Section 6 led to the matrix system

[TABLE]

where the components $\mathcal{H}_{ij}$ depend on the choice of coupling at the interface for the partitioned scheme as before. The first row of $\mathcal{G}$ is associated with boundary condition on the incoming characteristic. For the AMP scheme, the components of this row are given by

[TABLE]

where

[TABLE]

For the TP and ATP schemes, the components of the first row are given by

[TABLE]

where $\alpha=1,2$. The second row of $\mathcal{G}$ is common to all three schemes as it represents an extrapolation of the outgoing characteristic into the interface ghost point. This extrapolation condition gives

[TABLE]

A.4 Exact solution for the 1D inviscid model problem

The verification of the stability analysis of the inviscid model problem in Section 6.4.1 requires the exact solution of the model problem. Assuming that the solution has no variation in the $x$-direction, the ODE for the pressure with the boundary condition at $y=H$ gives

[TABLE]

where $v_{I}(t)$ is the interface velocity (to be determined). In terms of the characteristic variables of the solid, the solution is

[TABLE]

where $a_{0}(y)$ and $b_{0}(y)$, $y<0$, are given by the initial conditions for $\bar{v}_{2}$ and $\bar{\sigma}_{22}$, and $a_{I}(t)$ specifies the incoming characteristic variable of the solid at the interface (to be determined). Evaluating the outgoing characteristic of the solid at the interface gives $b_{I}(t)=b_{0}(-\bar{c}_{p}t)$, and this can be matched with the fluid to give

[TABLE]

The solution of the ODE in (153) for the interface velocity has the form

[TABLE]

where $\lambda=\bar{z}_{p}/(\rho H)$ and $v_{I}(0)$ is assumed to be given as an initial condition. The fluid pressure in (151) is now specified, and the incoming characteristic of the solid at the interface can now be obtained from

[TABLE]

This completes the characteristic description of the exact solution for the solid in (152), which gives

[TABLE]

The initial conditions used to specify the exact solution in Section 6.4.1 are

[TABLE]

and the parameters of the model problem are taken as $\rho=H=\bar{c}_{p}=1$, while $\bar{\rho}$ is varied as needed to generate the plots in Figure 6.