On the stability of a loosely-coupled scheme based on a Robin interface condition for fluid-structure interaction

TL;DR

This paper analyzes the stability of a loosely coupled fluid-structure interaction scheme using Robin interface conditions, identifying conditions for stability and optimal parameter choices through theoretical and numerical analysis.

Contribution

It provides a stability analysis of a Robin interface-based fluid-structure interaction scheme, revealing how interface parameters affect stability and accuracy.

Findings

Stability depends on interface parameter and time discretization.

Optimal interface parameter choices improve solution accuracy.

Numerical experiments confirm theoretical stability conditions.

Abstract

We consider a loosely coupled algorithm for fluid-structure interaction based on a Robin interface condition for the fluid problem (explicit Robin-Neumann scheme). We study the dependence of the stability of this method on the interface parameter in the Robin condition. In particular, for a model problem we find sufficient conditions for instability and stability of the method. In the latter case, we found a stability condition relating the time discretization parameter, the interface parameter, and the added mass effect. Numerical experiments confirm the theoretical findings and highlight optimal choices of the interface parameter that guarantee an accurate solution with respect to an implicit one.

Figures

Tables

Table 1. Table 1: Stability of the explicit RN scheme for different values of the parameters.

Mesh	$\Delta t$	$\alpha$	Stability
I	$0.625\cdot {10}^{-4}$	4689	OK
I	$1.25\cdot {10}^{-4}$	4689	NO
I	$1.25\cdot {10}^{-4}$	2000	OK
I	$5\cdot {10}^{-4}$	2000	OK
II	$5\cdot {10}^{-4}$	2000	NO

Equations

$${\rho }_f\frac{\partial \mathbf{u}}{\partial t}+\nabla p=0$$

$$\nabla \cdot \mathbf{u}=0$$

$$\mathbf{u}\cdot \mathbf{n}=\frac{\partial \eta }{\partial t}$$

$${\rho }_s{H}_s\frac{{\partial }^2\eta }{\partial t^2}+\beta \eta -\psi \frac{{\partial }^2\eta }{\partial x^2}=p$$

$${\rho }_f\frac{{\mathbf{u}}^n-{\mathbf{u}}^{n-1}}{\Delta t}+\nabla p^n=0$$

$$\nabla \cdot {\mathbf{u}}^n=0$$

$$u^n=\frac{{\eta }^n-{\eta }^{n-1}}{\Delta t}$$

$${\rho }_s{H}_s{\delta }_{tt}{\eta }^n+\beta {\eta }^n-\psi \frac{{\partial }^2{\eta }^n}{\partial x^2}=p^n$$

$$-\alpha u^n+p^n=-\alpha \frac{{\eta }^n-{\eta }^{n-1}}{\Delta t}+{\rho }_s{H}_s{\delta }_{tt}{\eta }^n+\beta {\eta }^n-\psi \frac{{\partial }^2{\eta }^n}{\partial x^2}\text{on}\Sigma .$$

$${\delta }_{tt}{\eta }^n={\rho }_s{H}_s\frac{{\eta }^n-2{\eta }^{n-1}+{\eta }^{n-2}}{\Delta t^2}.$$

$${\delta }_{tt}{\eta }^n={\rho }_s{H}_s\frac{{\eta }^{n+1}-2{\eta }^n+{\eta }^{n-1}}{\Delta t^2}.$$

$${\rho }_f\frac{{\mathbf{u}}^n-{\mathbf{u}}^{n-1}}{\Delta t}+\nabla p^n=0$$

$$\nabla \cdot {\mathbf{u}}^n=0$$

$$-\alpha u^n+p^n=-\alpha \frac{{\eta }^n-{\eta }^{n-1}}{\Delta t}+{\rho }_s{H}_s\frac{{\eta }^n-2{\eta }^{n-1}+{\eta }^{n-2}}{\Delta t^2}+\beta {\eta }^n-\psi \frac{{\partial }^2{\eta }^n}{\partial x^2}$$

$${\rho }_s{H}_s\frac{{\eta }^{n+1}-2{\eta }^n+{\eta }^{n-1}}{\Delta t^2}+\beta {\eta }^n-\psi \frac{{\partial }^2{\eta }^n}{\partial x^2}=p^n\text{in}\Sigma .$$

$$\beta {\eta }^n-\psi \frac{{\partial }^2{\eta }^n}{\partial x^2}-p^n=-{\rho }_s{H}_s\frac{{\eta }^{n+1}-2{\eta }^n+{\eta }^{n-1}}{\Delta t^2},$$

$$-\alpha u^n=-\alpha \frac{{\eta }^n-{\eta }^{n-1}}{\Delta t}+{\rho }_s{H}_s\frac{{\eta }^n-2{\eta }^{n-1}+{\eta }^{n-2}}{\Delta t^2}-{\rho }_s{H}_s\frac{{\eta }^{n+1}-2{\eta }^n+{\eta }^{n-1}}{\Delta t^2},$$

$$u^n=\frac{{\eta }^n-{\eta }^{n-1}}{\Delta t}+{\rho }_s{H}_s\frac{{\eta }^{n+1}-3{\eta }^n+3{\eta }^{n-1}-{\eta }^{n-2}}{\alpha \Delta t^2}.$$

$$p=-{\rho }_f\mathcal{M}\left(\frac{\partial (\mathbf{u}\cdot \mathbf{n})}{\partial t}\right)\text{in}(0,T)\times \Sigma .$$

$$p^n=-{\rho }_f\mathcal{M}\left(\frac{u^n-u^{n-1}}{\Delta t}\right).$$

$$p^n=-{\rho }_f\mathcal{M}\left(\frac{{\eta }^n-2{\eta }^{n-1}+{\eta }^{n-2}}{\Delta t^2}+\frac{{\rho }_s{H}_s}{\alpha \Delta t^3}\left({\eta }^{n+1}-4{\eta }^n+6{\eta }^{n-1}-4{\eta }^{n-2}+{\eta }^{n-3}\right)\right),$$

$${\eta }^m(x)=\sum _{i=1}^{\infty }{\eta }_i^m{g}_i(x),$$

$${\mu }_i=\frac{L}{i\pi \tanh \left(\frac{i\pi R}{L}\right)},{\lambda }_i={\left(\frac{i\pi }{L}\right)}^2,$$

$${\mu }_{min}\sim h,{\lambda }_{max}\sim h^{-2},{\mu }_{max}\sim h^0,$$

$${\rho }_s{H}_s<\underset{i}{\max }{\gamma }_i,{\gamma }_i=\alpha \Delta t\frac{4{\rho }_f{\mu }_i+\Delta t^2\left(\beta +\psi {\lambda }_i\right)}{16{\rho }_f{\mu }_i+4\alpha \Delta t}.$$

$${\rho }_s{H}_s\frac{{\eta }^{n+1}-2{\eta }^n+{\eta }^{n-1}}{\Delta t^2}+\beta {\eta }^n-\psi \frac{{\partial }^2{\eta }^n}{\partial x^2}+$$

$${\rho }_f\mathcal{M}\left(\frac{{\eta }^n-2{\eta }^{n-1}+{\eta }^{n-2}}{\Delta t^2}+\frac{{\rho }_s{H}_s}{\alpha \Delta t^3}\left({\eta }^{n+1}-4{\eta }^n+6{\eta }^{n-1}-4{\eta }^{n-2}+{\eta }^{n-3}\right)\right)=0.$$

$$\frac{{\rho }_s{H}_s}{\Delta t^2}\left(1+\frac{{\rho }_f{\mu }_i}{\alpha \Delta t}\right){\eta }_i^{n+1}+\left(-\frac{2{\rho }_s{H}_s}{\Delta t^2}+\beta +\psi {\lambda }_i+\frac{{\rho }_f{\mu }_i}{\Delta t^2}-\frac{4{\rho }_s{H}_s{\rho }_f{\mu }_i}{\alpha \Delta t^3}\right){\eta }_i^n+$$

$$\left(\frac{{\rho }_s{H}_s}{\Delta t^2}-2\frac{{\rho }_f{\mu }_i}{\Delta t^2}+\frac{6{\rho }_s{H}_s{\rho }_f{\mu }_i}{\alpha \Delta t^3}\right){\eta }_i^{n-1}+\left(\frac{{\rho }_f{\mu }_i}{\Delta t^2}-\frac{4{\rho }_s{H}_s{\rho }_f{\mu }_i}{\alpha \Delta t^3}\right){\eta }_i^{n-2}+\frac{{\rho }_s{H}_s{\rho }_f{\mu }_i}{\alpha \Delta t^3}{\eta }_i^{n-3}=0.$$

$$\left.\begin{array}[]{ll}\chi(y)=&\displaystyle\frac{\rho_{s}H_{s}}{\Delta t^{2}}\left(1+\frac{\alpha\Delta t}{\rho_{f}\mu_{i}}\right)y^{4}+\left(-\frac{2\alpha\rho_{s}H_{s}}{\rho_{f}\mu_{i}\Delta t}+\frac{\alpha\Delta t}{\rho_{f}\mu_{i}}\left(\beta+\psi\lambda_{i}\right)+\frac{\alpha}{\Delta t}-\frac{4\rho_{s}H_{s}}{\Delta t^{2}}\right)y^{3}\\ &\displaystyle+\left(\frac{\alpha\rho_{s}H_{s}}{\rho_{f}\mu_{i}\Delta t}-2\frac{\alpha}{\Delta t}+\frac{6\rho_{s}H_{s}}{\Delta t^{2}}\right)y^{2}+\left(\frac{\alpha}{\Delta t}-\frac{4\rho_{s}H_{s}}{\Delta t^{2}}\right)y+\frac{\rho_{s}H_{s}}{\Delta t^{2}}.\end{array}\right.$$

$$\chi (-1)$$

$$=\frac{\alpha }{{\rho }_f{\mu }_i\Delta t}\left(4{\rho }_s{H}_s-4{\rho }_f{\mu }_i-\Delta t^2\left(\beta +\psi {\lambda }_i\right)\right)+\frac{16{\rho }_s{H}_s}{\Delta t^2}.$$

$$y\in \left({\phi }_1\left(x\right),{\phi }_2\left(x\right)\right),$$

$$\frac{\partial f}{\partial y}\left(x,y\right)\ge b\left(x\right)>0$$

$$\frac{\partial f}{\partial y}\left(x,y\right)\le b\left(x\right)<0,$$

Full text

1 Introduction

Loosely-coupled schemes (also known as explicit) are a very attractive strategy for the numerical solution of the fluid-structure interaction (FSI) problem. Indeed, they are based on the solution of just one fluid and one structure problem at each time step, thus allowing a big improvement in the computational costs in comparison to fully-coupled (implicit) partitioned procedures and monolithic schemes. Another interesting feature of such schemes is that pre-existing fluid and structure solvers could be often employed.

For these reasons, loosely-coupled schemes have been widely used in many engineering applications such as aeroelasticity [27, 28, 12]. However, the stability properties of such schemes deteriorate when the so-called added mass effect becomes relevant. This happens, in particular, when the fluid and structure densities are comparable, as happens in hemodynamics [29]. For example, in [10] it has been proven that the classical explicit Dirichlet-Neumann scheme is unconditionally unstable in the hemodynamic regime, see also [16, 26].

In the recent years, there has been a growing interest in partitioned procedures that are based on Robin interface conditions. The latter are obtained by considering linear combinations of the standard interface conditions owing to the introduction of suitable parameters. The choice of such parameters is crucial for accelerating the convergence of implicit schemes [3, 4, 2, 17, 18]. Some works focused then on the design of stable loosely-coupled schemes for large added mass effect, which are based on Robin interface conditions [25, 19, 13, 15, 8, 7, 21, 6, 9]. These studies proposed specific values of the interface parameters which guarantee good stability properties (possibly in combination with suitable stabilizations).

In this paper, the explicit Robin-Neumann scheme, obtained by equipping the fluid subproblem with a Robin condition with parameter $\alpha$ and the structure one with a Neumann condition, is considered. In particular, it is investigated how the choice of the interface parameter $\alpha$ influences the stability of the numerical solution. To this aim, two analyses on a simplified problem are performed, the first one determining sufficient conditions for instability of the scheme, whereas the second one sufficient conditions for its stability. This will allow us to understand the dependence of stability and instability on the physical and numerical parameters and to properly design stable loosely-coupled schemes which could be easily implemented also by means of available (even commercial) solvers.

To validate the theoretical findings found in the analysis reported in Section 2, in Section 3 we eventually present the results of some numerical experiments, where the issue of accuracy is also discussed by proposing some ”optimal” value of the interface parameters.

2 Position of the problem

2.1 The fluid-structure interaction problem

We introduce in what follows the simplified fluid-structure interaction problem considered in the analyses below, see [10]. We consider the 2D fluid domain $\Omega^{f}$ which is a rectangle $R\times L$, where $R$ is denoted the “radius” and $L$ the length of the domain. Its boundary is given by $\partial\Omega^{f}=\Gamma\bigcup\Sigma$, where $\Sigma$ is the part where the interaction with the structure occurs, see Figure 1. For the structure, we consider a 1D model in the domain $\Omega^{s}=\Sigma$.

For the fluid modeling, we consider a linear incompressible inviscid problem, whereas for the structure the independent rings model [30]. The displacement could happen only in the radial direction. Moreover, we assume small displacements so that the structure deformation is negligible and the fluid domain can be considered fixed. Thus, we have the following FSI problem:

Find fluid velocity $\boldsymbol{u}$, fluid pressure $p$, and structure displacement $\eta$, such that

[TABLE]

where $\boldsymbol{n}$ is the outward normal, $\rho_{f}$ and $\rho_{s}$ the fluid and structure densities, $x$ the axial direction along which $\Sigma$ is located, $H_{s}$ the structure thickness, and $\beta$ and $\psi$ two suitable parameters accounting for the elasticity of the structure. Moreover, we have to equip the fluid problem with boundary conditions (for example of homogeneous type) on $\partial\Omega\setminus\Sigma$ and for the tangential component on $\Sigma$ [10]. Condition $\eqref{eq:fsi-3}$ represents a no-slip condition at the interface $\Sigma$ between fluid and structure (perfect adherence or kinematic condition). Due to the lower space dimension of the structure, the independent rings model (1d) represents also the third Newton law (continuity of the normal stresses or dynamic condition).

2.2 Time discretization and explcit Robin-Neumann scheme

Denoting by $\Delta t$ the time discretization parameter, $v^{n}$ the approximation of $v(t^{n}),\,t^{n}=n\Delta t$, and setting $u^{n}=(\boldsymbol{u}^{n}\cdot\boldsymbol{n})|_{\Sigma}$, we have the following discretized-in-time version of problem (1): Find for each $n$, $\boldsymbol{u}^{n},\,p^{n}$, and $\eta^{n}$, such that

[TABLE]

Notice that we have considered a backward Euler approximation for the fluid problem and we indicated with $\delta_{tt}$ the approximation of the second derivative in time for the structure problem.

Introducing $\alpha>0$, we can substitute in (2) the kinematic condition (2c) with the following linear combination obtained with the dynamic condition (2d):

[TABLE]

Of course, the solution of problem (2a)-(2b)-(3)-(2d) coincides with that of (2).

Now, the idea is to consider a partitioned method where condition (3) is given to the fluid problem, whereas (2d) is in fact the structure problem. Since the fluid problem has been discretized with an implicit method, in (3) we use the following approximation of the second derivative

[TABLE]

Instead, for the structure problem (2d) we use the explicit leap-frog approximation

[TABLE]

We observe that, due to the explicit time discretization of the structure problem, fluid and structure are in fact decoupled and, accordingly, we can introduce the following algorithm.

Algorithm 1**.**

Explicit Robin-Neumann algorithm.*

Given $\boldsymbol{u}^{0},\,\eta^{1},\,\eta^{0}$, for $n\geq 1$, at time step $t^{n}$:*

Solve the fluid problem with a Robin condition at the interface $\Sigma$:

[TABLE] 2. 2.

Then, solve the structure problem (Neumann condition at the interface):

[TABLE]

In the next section, we study how the stability of the previous algorithm is affected by the choice of the parameter $\alpha$.

3 Stability analysis

3.1 Preliminaries

First, we notice that using the previous algorithm, the discrete kinematic condition (2c) is not satisfied anymore. Indeed, from (5) we have

[TABLE]

where it is understood from now on that the equalities we derive hold true at the interface $\Sigma$. By introducing the latter expression in (4c), we obtain

[TABLE]

which leads to

[TABLE]

The latter equality provides a ”correction” of the discrete kinematic condition (2c) as a consequence of the explicit treatment.

Following then [10], we consider the added mass operator $\mathcal{M}:H^{-1/2}(\Sigma)\to H^{1/2}(\Sigma)$, which allows us to write the following relation between fluid pressure and velocity at the interface, under the assumption of null external pressure:

[TABLE]

At the time discrete level, we have

[TABLE]

Inserting (6) for both $u^{n}$ and $u^{n-1}$, we obtain

[TABLE]

which gives a relation between pressure and displacement at the interface.

We can write $\eta^{m}$ for any $m$ as a linear combination of the $L^{2}$ orthonormal basis functions $\left\{g_{i}(x)=\sqrt{2/L}\sin\left(\frac{i\pi x}{L}\right)\right\}$:

[TABLE]

for suitable coefficients $\eta^{m}_{i}$, see [10, 3]. Notice that $g_{i}$ are eigenfunctions of both the added mass operator $\mathcal{M}$ and of the Laplace operator $\mathcal{L}=-b\,\partial_{xx}|_{\Sigma}$, with eigenvalues given by

[TABLE]

respectively.

It is useful for later purposes to highlight that the eigenvalues $\mu_{i}$ and $\lambda_{i}$ of the discrete versions of the operators $\mathcal{M}$ and $\mathcal{L}$ obtained with finite elements, feature the following properties [10, 3]:

[TABLE]

where $h$ is the space discretization parameter.

3.2 Sufficient conditions for instability

We present in what follows a first result that provides sufficient conditions that guarantee conditional instability of the explicit Robin-Neumann scheme. This results generalizes the one proven in [10] about the unconditional instability of the Dirichlet-Neumann scheme ($\alpha\to+\infty$, see Proposition 3 in [10]).

Proposition 1**.**

The explicit Robin-Neumann scheme is unstable if

[TABLE]

Proof.

We start by inserting in the interface condition (5) the expression of $p^{n}$ given by (8), obtaining

[TABLE]

[TABLE]

Notice that the previous is a relation in the discrete structure displacement solely. Multiplying it by the basis function $g_{i}$ and integrating over the interface $\Sigma$, we obtain

[TABLE]

By multiplying the last identity by $\frac{\alpha\Delta t}{\rho_{f}\mu_{i}}$, we obtain the following characteristic polynomial corresponding to the previous difference equation:

[TABLE]

Now, we compute the value of the previous polynomial for $y=-1$:

[TABLE]

It follows that $\chi(-1)<0$ under condition (10). Since $\lim_{y\to-\infty}\chi(y)=+\infty$, it follows that in this case there exists at least one real root $\bar{y}<-1$ of the polynomial associated to the difference equation, implying that the method is unstable. ∎

3.3 Sufficient conditions for stability

We discuss in the following result some sufficient conditions that guarantee that the explicit Robin-Neumann scheme is conditionally stable. The idea is to start again from the polynomial (11) and discuss when its four roots have all modulus less than 1.

To this aim, we first introduce the following version of the implicit function theorem.

Theorem 1**.**

Let $f\in\mathcal{C}^{1}\left(\mathbb{R}^{2}\right)$ and suppose that for all $x\in\Omega$, an open interval, and for all

[TABLE]

where $\varphi_{1},\varphi_{2}:\Omega\rightarrow\mathbb{R}$ are continuous functions, we have either

[TABLE]

or

[TABLE]

for some continuous function $b:\Omega\rightarrow\mathbb{R}$. Let $g:\Omega\rightarrow\mathbb{R}$ be such that for all $x\in\Omega$

[TABLE]

Then, there is a unique function $\xi:\Omega\rightarrow\mathbb{R}$ such that, for all $x\in\Omega$, $\xi\left(x\right)\in\left(\varphi_{1}\left(x\right),\varphi_{2}\left(x\right)\right)$ and $f\left(x,\xi\left(x\right)\right)=0.$ Furthermore, for all $x\in\Omega$

[TABLE]

Proof.

Let us consider the case $\frac{\partial f}{\partial y}\left(x,y\right)\geq b\left(x\right)>0$. The other case follows from it after replacing $f\left(x,y\right)$ with $-f\left(x,y\right)$ and $b\left(x\right)$ with $-b\left(x\right)$. Fix $x\in\Omega$. By strict monotonicity of the function $y\mapsto f\left(x,y\right)$ in the interval $\left(\varphi_{1}\left(x\right),\varphi_{2}\left(x\right)\right)$, there is at most one value $\xi\left(x\right)\in\left(\varphi_{1}\left(x\right),\varphi_{2}\left(x\right)\right)$ for which $f\left(x,\xi\left(x\right)\right)=0$. This proves uniqueness.

Now, the function $y\mapsto f\left(x,y\right)$ takes the value $f\left(x,g\left(x\right)\right)$ at $y_{0}=g\left(x\right)$. Consider now the point

[TABLE]

and assume without loss of generality that $y_{0}\leq y_{1}$, that is $f\left(x,y_{0}\right)\leq 0$. By the hypotheses, we have $\left[y_{0},y_{1}\right]\subseteq\left(\varphi_{1}\left(x\right),\varphi_{2}\left(x\right)\right)$, so that for all $y\in\left[y_{0},y_{1}\right]$ we have $\frac{\partial f}{\partial y}\left(x,y\right)\geq b\left(x\right)>0$ and

[TABLE]

By the intermediate value theorem, there is a point $\xi\left(x\right)\in\left[y_{0},y_{1}\right]$ such that $f\left(x,\xi\left(x\right)\right)=0$. This prove existence.

Finally,

[TABLE]

means

[TABLE]

so that

[TABLE]

proving the last part of the theorem. ∎

Now, we are ready to introduce the main result of this work.

Theorem 2**.**

For all positive $\rho_{s}$, $\rho_{f}$, $\alpha$, $H_{s}$ and for all finite sequences of positive eigenvalues $\left\{\mu_{i}\right\}_{i=1}^{N}$, $\left\{\lambda_{i}\right\}_{i=1}^{N}$, there exists a positive $\delta$ such that, if $0<\Delta t<\delta$, then, for all $i=1,\ldots,N$, the polynomial (11) has four simple roots in the open unit disc in the complex plane.

Proof.

It is convenient to normalize $\chi\left(y\right)$, setting

[TABLE]

where

[TABLE]

For technical reasons, we will equivalently study the polynomial

[TABLE]

showing that all its four roots have modulus greater than $1$. Observe first that when $z=0$, $P\left(x\right)$ reduces to

[TABLE]

By classical results, this implies that for sufficiently small $z$, $P\left(x\right)$ has four simple roots as close as desired to $x=1$ in the complex plane (see [20], page 122). Unfortunately, this is not sufficient for our purposes, and for this reason we need a deeper analysis.

Now set $x=1+U$. This simplifies our formulas, since we look for roots that are close to $1$. This gives

[TABLE]

Recalling that $U$ is a complex variable, we write it as $U=u+iv$ where $u$ and $v$ are its real and imaginary parts, respectively. Thus, the equation

[TABLE]

reduces to the system

[TABLE]

Notice that the solution of the second equation $v=0$ reduces the first equation to

[TABLE]

which does not have any real solution $u>-1$ (all summands are positive). Since we look for roots $u+iv$ close to [math], we disregard the solution of the second equation $v=0$ and focus on the other solution

[TABLE]

which reduces the first equation to

[TABLE]

where

[TABLE]

Since we look for roots that go to [math] with $z$, we may assume that $u$ is $O\left(z\right)$ as $z\rightarrow 0$. The original equation (18) can therefore be approximated with

[TABLE]

further, disregarding all higher order terms in $z$, it can be approximated with

[TABLE]

If $u=O\left(z\right)$, then the third term in the above equation can be neglected and the same equation can be approximated with

[TABLE]

which gives the approximate solution

[TABLE]

If $u=O\left(z^{2}\right)$, then the first term in equation (19) can be neglected and (19) can be approximated with

[TABLE]

which gives the approximate solution

[TABLE]

We now need to estimate the derivative with respect to $u$ of

[TABLE]

that is

[TABLE]

We are now ready to apply Theorem 1, with the two following choices for $g$ suggested by the previous approximate solutions:

[TABLE]

Assume first $u$ close to $g_{1}\left(z\right).$ Precisely, assume that for some $K_{1}>0$, $\varphi_{1}\left(z\right)=g_{1}\left(z\right)-K_{1}z^{2}\leq u\leq\varphi_{2}\left(z\right)=g_{1}\left(z\right)+K_{1}z^{2}$, so that hypothesis (13a) in Theorem 1 is satisfied. In particular, $u=g_{1}\left(z\right)+O\left(z^{2}\right)$ and, from (20),

[TABLE]

Thus, for some small positive $\delta_{1},$ if $0<z<\delta_{1}$ and $g_{1}\left(z\right)-K_{1}z^{2}\leq u\leq g_{1}\left(z\right)+K_{1}z^{2}$, we have

[TABLE]

which shows that, with the above choices, hypothesis (12) in Theorem 1 is satisfied.

Now we need to check hypothesis (13b) in Theorem 1, that is

[TABLE]

or, in other words,

[TABLE]

We have:

[TABLE]

since the first two terms at numerator vanish. Thus, it follows from Theorem 1 that there exists a unique function

[TABLE]

such that $f\left(z,u_{1}\left(z\right)\right)=0$ for all $z\in\left(0,\delta_{1}\right).$ Furthermore, for all $z\in\left(0,\delta_{1}\right)$

[TABLE]

which implies that

[TABLE]

as $z\rightarrow 0^{+}.$

Similarly, let us assume now $u$ close to $g_{2}\left(z\right).$ Precisely, assume that for some $K_{2}>0$, $\psi_{1}\left(z\right)=g_{2}\left(z\right)-K_{2}z^{3}\leq u\leq\psi_{2}\left(z\right)=g_{2}\left(z\right)+K_{1}z^{3}$. Then, we have $u=g_{2}(z)+O\left(z^{3}\right)$ and

[TABLE]

Thus, for some small positive $\delta_{2},$ if $0<z<\delta_{2}$ and $g_{2}\left(z\right)-K_{2}z^{3}\leq u\leq g_{2}\left(z\right)+K_{2}z^{3}$, we have

[TABLE]

Now we need to check if

[TABLE]

or, in other words, if

[TABLE]

We have:

[TABLE]

Thus, it follows from Theorem 1 that there is a unique function

[TABLE]

such that $f\left(z,u_{2}\left(z\right)\right)=0$ for all $z\in\left(0,\delta_{2}\right).$ Furthermore, for all $z\in\left(0,\delta_{2}\right)$

[TABLE]

which implies that

[TABLE]

as $z\rightarrow 0^{+}.$

Let us now compute $v^{2}$ by means of (17) for these two choices. Setting $u=u_{1}\left(z\right)=\frac{A}{2}z+O\left(z^{2}\right)$, we obtain

[TABLE]

Setting now $u=u_{2}\left(z\right)=-\frac{AC_{i}B_{i}}{2\left(A+B_{i}\right)^{2}}z^{2}+O\left(z^{3}\right)$, we obtain

[TABLE]

Thus, $(u_{1},v_{1})$ and $(u_{2},v_{2})$ are solutions of system (16). This means that $u_{1}\pm iv_{1}$ and $u_{2}\pm iv_{2}$ are the 4 roots of (15) and $1+u_{1}\pm iv_{1}$ and $1+u_{2}\pm iv_{2}$ are the four roots of (14).

The final step of the proof is to show that these four roots for $P(x)$ have modulus greater than 1. Accordingly, we have

[TABLE]

Thus, there exists a positive $\delta_{3}$ such that $\left(1+u_{1}\right)^{2}+v_{1}^{2}>1$ for $z<\delta_{3}$. Moreover, we have

[TABLE]

Thus, there exists a positive $\delta_{4}$ such that $\left(1+u_{2}\right)^{2}+v_{2}^{2}>1$ for $z<\delta_{4}$.

This proves that the four roots of $\chi(y)$ in (11) have all modulus strictly less than 1, provided that $\Delta t<\delta=\min\{\delta_{1},\delta_{2},\delta_{3},\delta_{4}\}$.

∎

3.4 Discussion

The results proven above give us sufficient conditions for instability or stability of the explicit Robin-Neumann scheme. In particular, in the last case, we found that, given $\alpha>0$, for $\Delta t$ small enough this scheme is absolute stable.

In view of the numerical experiments, we discuss hereafter more in detail the hypotheses of the previous results. To this aim, in what follows we write three conditions that imply (10):

• i)

[TABLE]

This is obtained by taking the greatest possible value of $i$ in (10);

• ii)

[TABLE]

This is obtained by solving (21) in the variable $\alpha$;

• iii)

[TABLE]

This is obtained by taking $i=1$ in (10), deleting the term $\Delta t^{2}\left(\beta+\psi\lambda_{1}\right)$ and solving in $\alpha$.

Dependence on $\rho_{s}H_{s}$. By looking at conditions (22) and (23), we see that when $\rho_{s}H_{s}<\max\left(\eta_{1},\eta_{2}\right)$, the explicit Robin-Neumann scheme is unstable if $\alpha$ is large enough. In particular, for decreasing values of $\rho_{s}H_{s}<\eta_{1}$ (i.e. for increasing added mass effect), the value of $\alpha_{1}$ decreases, enlarging the range of $\alpha$ such that the scheme is unstable. The same argument holds true for $\alpha_{2}$ when $\rho_{s}H_{s}<\eta_{2}$. 2. 2.

Dependence on $\Delta t$. From Theorem 2, we have that for any fixed $\alpha>0$, the explicit Robin-Neumann scheme is stable provided that $\Delta t$ is small enough. This result is consistent with Proposition 1. Indeed, for all indices $i$, we have

[TABLE]

Accordingly, we have also from (22)

[TABLE]

Observe also that, by (23), for $\rho_{s}H_{s}<\eta_{2}$ instability of the scheme follows if

[TABLE]

This means that in order to have stability according to Theorem 2 we need at least $\delta\leq c\alpha^{-1}$. The dependence of $\delta$ on $\alpha$ is still under study and will be hopefully discussed in future studies. 3. 3.

Dependence on $h$. By exploiting the behaviour of $\mu_{min}$ and $\lambda_{max}$ with respect to $h$, see (9), we find from (21) that $\lim_{h\rightarrow 0}\overline{\eta}=+\infty$ for $\Delta t$ fixed. This means that the stability properties of the method deteriorates when $h$ decreases. This is justified by the fact that when $h$ becomes small, the fluid and structure solutions should match a large number of d.o.f. at the interface and, due to the implicit treatment for the fluid time discretization and the explicit one for the structure, this matching becomes more and more difficult for smaller values of $h$. 4. 4.

$\boldsymbol{\Delta}\boldsymbol{t},\boldsymbol{h}\to\boldsymbol{0}$. From (9), we have that if $\Delta t\sim h$, then $\lim_{\Delta t,h\to 0}\alpha_{1}\simeq\frac{16\rho_{s}H_{s}\rho_{f}}{\psi+4\rho_{f}-4\rho_{s}H_{s}}$. When the added mass effect is large enough ($\rho_{s}H_{s}$ small enough), this limit is positive and bounded, unlike the case $\Delta t\to 0,h$ fixed (cf. point 2). This means that in this case we have instability for a wide range of values of $\alpha$ even as $\Delta t\to 0$. This suggests that the value of $\delta$ in Theorem 2 should be smaller (up to a constant) than $h$. Indeed, we have the following result.

Lemma 1**.**

The value of $\delta$ in Theorem 2 satisfies the relation

[TABLE]

for a suitable constant $c$.

Proof.

Fix $\alpha>0$. Call $\delta=f(h)$ the relationship between $\delta$ and $h$. Let assume that the thesis is not true, that is that there exists a sequence of values $h_{j}$ and correspondingly $\delta_{j}=f(h_{j})$, such that

[TABLE]

From Theorem 2, we have that stability is guaranteed if $\Delta t=\Delta t_{j}=\delta_{j}/2$. On the other side, from (24) and the choice of $\Delta t$ above, we have that $h_{j}=o(\delta_{j})=o(\Delta t_{j})$. From (21), we have that in the range $h=o\left(\Delta t\right)$, $\lim_{\Delta t,h\rightarrow 0}\overline{\eta}=+\infty$, obtaining unconditional instability. This contradicts the previous finding about stability. This means that the thesis is true. ∎ 5. 5.

The cases $\alpha=+\infty$ and $\alpha=0$. Theorem 2 holds true for any $\alpha\in(0,+\infty)$. The case $\alpha=+\infty$ corresponds to the explicit Dirichlet-Neumann scheme. In this case, the polynomial $\chi(y)$ in (11) corresponds to the one found in [10] (see Proposition 3 therein), where it is shown that at least one root has modulus greater than one independently of $\Delta t$ (cf. also (21)).

Regarding $\alpha=0$, we obtain $\chi(y)=\frac{\rho_{s}H_{s}}{\Delta t^{2}}(y-1)^{4}$. This means that the solution does not blow up, even if it is not strictly absolute stable. This is in accordance with the fact that in this case the numerical solution does not evolve during the time evolution, being always equal to the initial condition (the same Neumann datum is transferred at the interface). Thus, accuracy is completely lost. From this observation, we can argue that too small values of $\alpha$, even though give stability, do not lead to an accurate solution.

4 Numerical results

4.1 Preliminaries

In this section we present some numerical results with the aim of validating the theoretical findings reported in the previous section. In particular, we studied the effectiveness of the analyses obtained for the simplified models, when applied to complete three-dimensional fluid and structure models. All the simulations are inspired from hemodynamics. This problem is of great interest for our purposes since it is characterized by a high added mass effect, so that the stability of explicit methods is a challenging issue. Moreover, the thickness of the structure is small enough to make acceptable the use of a membrane model in the stability analysis of the previous section.

We considered the coupling between the 3D incompressible Navier-Stokes equations written in the Arbitrary Lagrangian-Eulerian formulation [11] and the 3D linear infinitesimal elasticity, see [23, 18] for all the details. For the time discretization, we used the BDF scheme of order 1 for both the fluid and structure, with a semi-implicit treatment of the fluid convective term, whereas for the numerical discretization we employed $P1bubble-P1$ finite elements for the fluid and $P1$ finite elements for the structure. The fluid domain at each time step is obtained by extrapolation of previous time steps (semi-implicit approach [14, 5, 24]). We used the following data: fluid viscosity $\mu=0.035\,g/(cm\,s)$, fluid density $\rho_{f}=1\,g/cm^{3}$, structure density $\rho_{s}=1.1\,g/cm^{3}$, Poisson ratio $\nu=0.49$, Young modulus $E=3\cdot 10^{5}\,Pa$, surrounding tissue parameter for the structure problem [22] $\gamma_{ST}=1.5\cdot 10^{5}\,Pa/cm$.

The fluid domain is a cylinder with length $L=5\,cm$ and radius $R=0.5\,cm$, whereas the structure domain is the external cylindrical crown with thickness $H_{s}=0.1\,cm$. We considered a couple of meshes with 4680 tetrahedra and 1050 vertices for the fluid and 1260 vertices for the structure (mesh I). Another couple of meshes (mesh II) has been obtained by halving the values of the space discretization parameter.

At the inlet we prescribed a Neumann condition given by the following pressure function

[TABLE]

with absorbing resistance conditions at the outlets [25, 23].

All the numerical results have been obtained with the parallel Finite Element library LIFEV [1].

4.2 On the stability of the numerical solution

In the first set of numerical experiments, we study the stability of the solution. The time discretization parameter is $\Delta t=0.0005\,s$. In Figure 2 we report the mean pressure over the middle cross section of the cylinder obtained for different values of $\alpha$ and with an implicit method.

From these results we observe stability of the numerical solution obtained with the explicit Robin-Neumann method for some values of the interface parameter $\alpha$. The accuracy seems to deteriorate for decreasing values of $\alpha$ (cf. point 5 in the Discussion of Sect. 3.4). Notice that with $\alpha=2500$ the numerical solution (not reported here) blows up. The same happens for bigger values of $\alpha$. This is consistent with the result proven in Proposition 1, for which an unstable solution is obtained for $\alpha$ greater than a threshold when the added mass effect is large enough, see (22)-(23) (cf. also point 1 in the Discussion of Sect. 3.4).

In Table 1 we indicate if stability is achieved for different space and time discretizations parameters.

From the first two rows, we observe that, given a value of $\alpha$, stability is achieved only if $\Delta t$ is small enough. This was expected from the theoretical findings, see Theorem 2. From the first three lines, we observe that the value of $\Delta t$ needed to have stability could be increased when $\alpha$ is decreased. This is in accordance with point 2 in the Discussion of Sect. 3.4. Finally, from the last two rows, we find that stability is achieved for given values of $\alpha$ and $\Delta t$ if the mesh is not too fine, thus confirming the observation made in point 3 of the Discussion of Sect. 3.4.

4.3 On the accuracy of the numerical solution

Two very interesting topics that are not yet discussed are: i) how to select not empirically a reasonable value of $\alpha$ that should fall in the range of stability and ii) which is the dependence of the accuracy of explicit Robin-Neumann methods on $\alpha$. In what follows we provide some preliminary hints about these two points.

The idea we propose is to use an optimal value of $\alpha$ which should guarantee efficient convergence of the implicit Robin-Neumann scheme. Since such a value makes the convergence factor of the method (thus the error at each iteraton) small, we expect that the use of the same value of $\alpha$ for the explicit counterpart of the method should reduce the error accumulated at each time step.

In particular, we propose here to use the optimality result proven in [18] for the implicit Robin-Robin method in the case of cylindrical interfaces. This leads to two optimal values $\alpha^{opt}_{f}$ and $\alpha^{opt}_{s}$ in the Robin interface conditions. Since in the hemodynamic regime the convergence properties of the implicit Robin-Neumann scheme with the optimal value $\alpha^{opt}_{f}$ are very similar to that of the ”optimal” implicit Robin-Robin scheme [17], we propose here to use $\alpha=\alpha^{opt}_{f}$ as an effective value that should provide a stable and accurate solution for the explicit Robin-Neumann scheme.

We consider the same numerical experiment as above, with time discretization parameter $\Delta t=0.001\,s$. In Figure 3, we report the mean pressure over the middle cross section for $\alpha=\alpha_{f}^{opt}=2250$ and $\alpha=1500,3000$.

For $\alpha>3500$ the solution blows up, whereas for $\alpha<1500$ the accuracy deteriorates. We observe that the solution obtained with the value $\alpha_{f}^{opt}$ proposed a priori is very close to the optimal one ($\alpha=3000$) found empirically. This result highlights that an effective choice of $\alpha$ that guarantees stability and accuracy is possible, at least in the case of a cylindrical domain.

Remark 1**.**

Notice that, from empirical observations running the implicit Robin-Robin scheme, we found that $\alpha_{f}^{opt}$ is greater than zero and does not blow up when $\Delta t\to 0$. Thus, from point 2 of the Discussion of Sect. 3.4, we have that $\delta\leq c(\alpha^{opt})^{-1}$ still makes sense also when $\Delta t\to 0$, since $\delta>0$.

4.4 Conclusions

In this work we have proposed an analysis of stability of a loosely coupled scheme of Robin-Neumann type for fluid-structure interaction, possibly featuring a large added mass effect. In order to make the results found in this paper reliable for concrete applications, we need now to apply them to realistic geometries and boundary data. This would make the explicit Robin-Neumann scheme an effective strategy to be considered for example in hemodynamics, where the added mass effect is elevated. This is actually outside the aims of the paper which is focused on the analysis of sufficient conditions for conditional stability and instability. For this reason, we are currently studying this topic for a future development of this work. What in our opinion is promising is that our numerical experiments in the hemodynamic regime highlighted an excellent agreement with the theoretical findings and lead to accurate solutions in the range of stability.