Decoupling Techniques for Coupled PDE Models in Fluid Dynamics

Abstract

We review decoupling techniques for coupled PDE models in fluid dynamics. In particular, we are interested in the coupled models for fluid flow interacting with porous media flow and the fluid structure interaction (FSI) models. For coupled models for fluid flow interacting with porous media flow, we present decoupled preconditioning techniques, two-level and multilevel methods, Newton-type linearization-based two-level and multilevel algorithms, and partitioned time-stepping methods. The main theory and some numerical experiments are given to illustrate the effectiveness and efficiency of these methods. For the FSI models, partitioned time-stepping algorithms and a multirate time-stepping algorithm are carefully studied and analyzed. Numerical experiments are presented to highlight the advantages of these methods.

1. Introduction

Coupled PDE models have wide applications in the real world. For example, in fluid dynamics, there are two-phase flow models, fluid structure interaction models, heat transfer models in fluids etc. In this work, we focus on two typical models: coupled models for describing fluid flow interacting with porous media flow, and fluid structure interaction (FSI) models. The first type of models have been validated by experiments [1] and then justified by using homogenization theory [2, 3]. Applications include the environmental engineering problem of groundwater contamination through rivers and the geoscience problem of surface flows filtrating in vuggy porous media. The second type of models come from many practical applications. For example, blood flow interacting with vessel wall, compressible fluids interacting with aircraft wings, as well as slamming and whipping response of ship structure to water flow. These models are typical multidomain coupled PDE models with multiphysics. Due to the heterogenenities in subdomain models, it is very difficult to find a unified approach to solve the different subdomain models simultaneously. Moreover, some coupled models have nonlinearity in either subdomain problems or interface coupling terms. To deal with the difficulties caused by the coupling of different submodels and the nonlinearity [4, 5, 6, 7, 8], we discuss some decoupling techniques [4, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20], which have shown to be very effective and efficient. Among them, we will emphasize the work proposed by our group members and highlight the novelty and importance of these algorithms.

The rest of the paper is organized as follows. In Section 2, we introduce the coupled fluid flow/porous media interacting models. Some decoupling techniques, specifically, decoupled preconditioners, decoupled two-level and multi-level methods, and partitioned time schemes will be presented and analyzed. In Section 3, we present fluid structure interaction models. Partitioned decoupling algorithms include the Robin-Neumann scheme [14, 15], the $\beta$-scheme [13, 21], and the multirate partitioned schemes [22, 23] will be briefly introduced. Concluding remarks are drawn in Section 4.

2. Decoupled algorithms for the coupled models of fluid flow interacting with porous flow

For the linear cases of the coupled models for fluid flow interacting with porous media flow, we refer to the Stokes/Darcy model studied in [7, 11, 18, 20, 24, 25, 26, 27, 28, 29, 30]. For the nonlinear case, we refer to the coupled nonlinear Navier–Stokes/Darcy model [4, 5, 6].

2.1 Coupled models for fluid flow interacting with porous media flow

Let $\mathrm{\Omega }\subset R^d$ be a domain consisting of a fluid region ${\mathrm{\Omega }}_f$ and a porous media region ${\mathrm{\Omega }}_p$ separated by an interface $\mathrm{\Gamma }$, as shown in Figure 1, where $d=2$ or 3, $\mathrm{\Omega }={\mathrm{\Omega }}_f\cup {\mathrm{\Omega }}_p$ and $\mathrm{\Gamma }={\overline{\mathrm{\Omega }}}_f\cap {\overline{\mathrm{\Omega }}}_p$. Let ${\mathit{n}}_f$ and ${\mathit{n}}_p$ denote the unit outward normal directions on $\partial {\mathrm{\Omega }}_f$ and $\partial {\mathrm{\Omega }}_p$. The interface $\mathrm{\Gamma }$ is assumed to be smooth enough as in [6].

Figure 1.

A global domain $\mathrm{\Omega }$ consisting of a fluid region ${\mathrm{\Omega }}_f$ and a porous media region ${\mathrm{\Omega }}_p$ separated by an interface $\mathrm{\Gamma }$.

For incompressible Newtonian fluid flow, Navier–Stokes equations of the stress-divergence form are usually used [31, 32]. $\forall t\ge 0$, $\forall \mathit{x}\in {\mathrm{\Omega }}_f$,

$$\left\{\begin{array}{l}{\rho }_f\left(\frac{\partial \mathit{u}}{\partial t}+\left(\mathit{u}\cdot \nabla \right)\mathit{u}\right)-\mathit{div}\mathit{T}\left(\mathit{u},p\right)={\mathit{f}}_f,\\ \mathit{div}\mathit{u}=0,\end{array}\right.$$E1

where ${\rho }_f$ is the fluid density, $\mathit{u}$ is the velocity vector, $p$ is the pressure, ${\mathit{f}}_f$ is the external force,

$$\mathit{T}\left(\mathit{u},p\right)=2\nu \mathit{D}\left(\mathit{u}\right)-p\mathit{I},\text{with}\mathit{D}\left(\mathit{u}\right)=\frac{1}{2}\left[\nabla \mathit{u}+{\left(\nabla \mathit{u}\right)}^T\right],$$E2

is the stress tensor with $\nu >0$ being the kinematic viscosity. By dropping the term $\frac{\partial \mathit{u}}{\partial t}$ in (1), the steady state Navier–Stokes equations read as: $\forall \mathit{x}\in {\mathrm{\Omega }}_f$,

$$\left\{\begin{array}{l}{\rho }_f\left(\mathit{u}\cdot \nabla \right)\mathit{u}-div\mathit{T}\left(\mathit{u},p\right)={\mathit{f}}_f,\\ div\mathit{u}=0.\end{array}\right.$$E3

In strong form,

$$-div\mathit{T}\left(\mathit{u},p\right)=-\nu \mathrm{\Delta }\mathit{u}+\nabla p,$$E4

because the fluid flow is assumed to be divergence free and $div\left({\nabla }^T\mathit{u}\right)=0$ holds.

Among various porous media flow models, Darcy’s law is the most favored. The governing variable in ${\mathrm{\Omega }}_p$ is the so-called piezometric head or pressure head,

$$\varphi =z+\frac{p_p}{{\rho }_{f}g}.$$E5

Here, $z$ is the elevation from a reference level (for simplicity, $z$ is assumed to be $0$). Darcy’s law states that the velocity ${\mathit{u}}_p$ (also called seepage velocity) in the porous media region is proportional to the gradient of $\varphi$ [27, 33].

$${\mathit{u}}_p=-\mathit{K}\nabla \varphi .$$E6

We assume that

$${\alpha }_1\left(\mathit{x},\mathit{x}\right)\le \left(\mathit{Kx},\mathit{x}\right)\le {\alpha }_2\left(\mathit{x},\mathit{x}\right),\forall \mathit{x}\in {\mathrm{\Omega }}_p.$$E7

Moreover, the divergence of the seepage velocity equals to the source term. This leads to the following steady state equation:

$$-div\left(\mathit{K}\nabla \varphi \right)=f_p.$$E8

In the time-dependent case, the governing equations in ${\mathrm{\Omega }}_p$ reads as:

$$S_0{\varphi }_t-div\left(\mathit{K}\nabla \varphi \right)=f_p,$$E9

where $S_0$ is a specific storage and $f_p$ is a source term.

No matter time-dependent or steady state, the key part of the coupled model is a set of interface conditions, which describe the interaction mechanism of the two different types of flows. The following interface conditions have been extensively used and studied in the literature [1, 2, 3, 7, 27, 28]:

$$\left\{\begin{array}{l}\mathit{u}\cdot {\mathit{n}}_f={\mathit{u}}_p\cdot {\mathit{n}}_f=-\mathit{K}\nabla \varphi \cdot {\mathit{n}}_f,\\ -\nu \left(\nabla {\mathit{un}}_f\right)\cdot {\mathit{n}}_f+p=\mathit{\rho g\varphi },\\ -\nu \left(\nabla {\mathit{un}}_f\right)\cdot {\tau }_i=\frac{{\mathit{\nu \alpha }}_{\mathit{BJS}}}{\sqrt{{\mathit{\nu \tau }}_i\cdot \mathit{K}{\tau }_i}}\mathit{u}\cdot {\tau }_i,i=1,\dots ,d-1.\end{array}\right.$$E10

Here, ${\left\{{\tau }_i\right\}}_{i=1}^{d-1}$ is the unit tangent vector on $\mathrm{\Gamma }$, ${\alpha }_{\mathit{BJS}}$ is a positive parameter depending on the properties of the porous medium. The first interface condition ensures mass conservation across $\mathrm{\Gamma }$. The second one is the balance of normal forces across the interface. The third condition is well known as Beavers-Joseph-Saffman’s law [1, 2], which states that the slip velocity is proportional to the shear stress along $\mathrm{\Gamma }$. Without loss of generality, we impose homogeneous Dirichlet boundary conditions on both of the external boundaries:

$$\left\{\begin{array}{ll}\mathit{u}=\mathbf{0}& \mathrm{on}{\mathrm{\Gamma }}_f,\\ \varphi =0& \mathrm{on}{\mathrm{\Gamma }}_p.\end{array}\right.$$E11

The proper functional spaces for $\mathit{u}$, $p$ and $\varphi$ are

$$\begin{array}{c}{\mathit{X}}_f=\left\{\mathit{v}\in {\mathit{H}}^1\left({\mathrm{\Omega }}_f\right)={\left(H^1\left({\mathrm{\Omega }}_f\right)\right)}^d|\mathit{v}=\mathbf{0}\mathrm{on}{\mathrm{\Gamma }}_f\right\},\\ Q=L^2\left({\mathrm{\Omega }}_f\right),X_p=\left\{\psi \in H^1\left({\mathrm{\Omega }}_p\right)|\psi =0\mathrm{on}{\mathrm{\Gamma }}_p\right\}.\end{array}$$E12

Moreover, we denote $\underset{¯}{X}={\mathit{X}}_f\times X_p$ for ease of presentation. Multiplying test functions to (3) and (8), integrating by parts and plugging in the interface boundary conditions (10)–(11), we have the weak form of the coupled Navier–Stokes/Darcy model: find $\underset{¯}{u}=\left(\mathit{u},\varphi \right)\in \underset{¯}{X},p\in Q$ such that

$$\left\{\begin{array}{ll}a\left(\underset{¯}{u},\underset{¯}{v}\right)+c\left(\mathit{u},\mathit{u},\mathit{v}\right)+b\left(\mathit{v},p\right)=f\left(\underset{¯}{v}\right)& \forall \underset{¯}{v}=\left(\mathit{v},\psi \right)\in \underset{¯}{X},\\ b\left(\mathit{u},q\right)=0& \forall q\in Q,\end{array}\right.$$E13

where

$$\begin{array}{l}a\left(\underset{¯}{u},\underset{¯}{v}\right)=a_f\left(\mathit{u},\mathit{v}\right)+a_p\left(\varphi ,\psi \right)+a_{\mathrm{\Gamma }}\left(\underset{¯}{u},\underset{¯}{v}\right),b\left(\mathit{v},p\right)=-{\int }_{{\mathrm{\Omega }}_f}p\nabla \cdot \mathit{v},\\ c\left(\mathit{u},\mathit{v},\mathit{w}\right)=\rho {\int }_{{\mathrm{\Omega }}_f}\left(\mathit{u}\cdot \nabla \right)\mathit{v}\cdot \mathit{w},f\left(\underset{¯}{v}\right)={\int }_{{\mathrm{\Omega }}_f}{\mathit{f}}_f\cdot \mathit{v}+\mathit{\rho g}{\int }_{{\mathrm{\Omega }}_p}{f}_p\psi \end{array}$$E14

with

$$\begin{array}{l}{a}_f\left(\mathit{u},\mathit{v}\right)=\nu {\int }_{{\mathrm{\Omega }}_f}\nabla \mathit{u}:\nabla \mathit{v}+\sum _{i=1}^{d-1}\frac{{\mathit{\nu \alpha }}_{\mathit{BJS}}}{\sqrt{{\mathit{\nu \tau }}_i\cdot \mathit{K}{\tau }_i}}{\int }_{\mathrm{\Gamma }}\left(\mathit{u}\cdot {\tau }_i\right)\left(\mathit{v}\cdot {\tau }_i\right),\\ a_p\left(\varphi ,\psi \right)=\mathit{\rho g}{\int }_{{\mathrm{\Omega }}_p}\nabla \psi \cdot \mathit{K}\nabla \varphi ,a_{\mathrm{\Gamma }}\left(\underset{¯}{u},\underset{¯}{v}\right)=\mathit{\rho g}{\int }_{\mathrm{\Gamma }}\left(\varphi \mathit{v}-\psi \mathit{u}\right)\cdot {\mathit{n}}_f.\end{array}$$E15

For the wellposedness of the coupled Navier–Stokes/Darcy model, we refer to the recent results in [5, 6, 34]. It is shown in [34] that if the viscosity is sufficiently large and the normal velocity across the interface is sufficiently small, then the problem (14) is wellposed. We follow the assumptions in [34]. In addition to these assumptions on model parameters and variables, we shall frequently use these properties: $a\left(\cdot ,\cdot \right)$ is bounded and coercive; $b\left(\cdot ,\cdot \right)$ is bounded and satisfies the inf-sup condition [27, 35]; and the nonlinear term can be bounded by using the $H^1$ norm of the three components [31, 36].

We partition ${\mathrm{\Omega }}_f$ and ${\mathrm{\Omega }}_p$ by quasi-uniform triangulations ${\mathcal{T}}_{f,h}$ and ${\mathcal{T}}_{p,h}$ with a characteristic meshsize $h$. Here, we require that the two subdomain triangulations coincide at $\mathrm{\Gamma }$. The corresponding finite element spaces are denoted by ${\mathit{X}}_{f,h}\times Q_h\subset {\mathit{X}}_f\times Q$ and $X_{p,h}\subset X_p$, respectively. Moreover, ${\mathit{X}}_{f,h}\times Q_h$ needs to be stable, i.e., there exists a positive constant $\beta$ such that

$$\underset{{\mathit{v}}_h\in {\mathit{X}}_{f,h}}{sup}\frac{b\left({\mathit{v}}_h,q_h\right)}{{\left|{\mathit{v}}_h\right|}_{1,{\mathrm{\Omega }}_f}}\ge \beta \parallel q_h{\parallel }_{0,{\mathrm{\Omega }}_f}\forall q_h\in Q_h.$$E16

We assume that the solution of (13) is smooth enough and the finite element spaces have the following typical approximation properties: for all $\left(\mathit{u},p\right)\in {\mathit{H}}^{k+1}\left({\mathrm{\Omega }}_f\right)\cap {\mathit{X}}_f\times H^k\left({\mathrm{\Omega }}_f\right)$ and $\varphi \in H^{k+1}\left({\mathrm{\Omega }}_p\right)\cap X_p$,

$$\underset{{\mathit{v}}_h\in {\mathit{X}}_{f,h},q_h\in Q_h}{inf}\left\{h{\left|\mathit{u}-{\mathit{v}}_h\right|}_{1,{\mathrm{\Omega }}_f}+\parallel \mathit{u}-{\mathit{v}}_h{\parallel }_{0,{\mathrm{\Omega }}_f}+h\parallel p-q_h{\parallel }_{0,{\mathrm{\Omega }}_f}\right\}\lesssim h^{k+1}\left({\left|\mathit{u}\right|}_{k+1,{\mathrm{\Omega }}_f}+{\left|p\right|}_{k,{\mathrm{\Omega }}_f}\right),$$E17

$$\underset{{\psi }_h\in X_{p,h}}{inf}\left\{h{\left|\varphi -{\psi }_h\right|}_{1,{\mathrm{\Omega }}_p}+\parallel \varphi -{\psi }_h{\parallel }_{0,{\mathrm{\Omega }}_p}\right\}\lesssim h^{k+1}{\left|\varphi \right|}_{k+1,{\mathrm{\Omega }}_p}.$$E18

To satisfy the discrete inf-sup condition and the approximation properties (17)–(18), if $k=1$, one may apply the Mini elements [31, 37] in ${\mathrm{\Omega }}_f$ and the piecewise linear elements in ${\mathrm{\Omega }}_p$; if $k\ge 2$, the $k$-th order Taylor-Hood elements [31, 37, 38] can be applied in ${\mathrm{\Omega }}_f$ and the piecewise $k$-th order elements can be adopted in ${\mathrm{\Omega }}_p$.

Coupled Algorithm: A conventional finite element discretization applied to the model problem (13) leads to the discrete problem: Find ${\underset{¯}{u}}_h=\left({\mathit{u}}_h,{\varphi }_h\right)\in {\underset{¯}{X}}_h={\mathit{X}}_{f,h}\times X_{p,h},p_h\in Q_h$ such that

$$\left\{\begin{array}{ll}a\left({\underset{¯}{u}}_h,{\underset{¯}{v}}_h\right)+c\left({\mathit{u}}_h,{\mathit{u}}_h,{\mathit{v}}_h\right)+b\left({\mathit{v}}_h,p_h\right)=f\left({\underset{¯}{v}}_h\right)& \forall {\underset{¯}{v}}_h=\left({\mathit{v}}_h,{\psi }_h\right)\in {\underset{¯}{X}}_h,\\ b\left({\mathit{u}}_h,q_h\right)=0,& \forall q_h\in Q_h.\end{array}\right.$$E19

2.2 Decoupled algorithms in the preconditioning steps

As an illustration of decoupled preconditioning techniques, we will consider the linear Stokes/Darcy model, whose weak form is (19) while the nonlinear term is dropped. We note that the discrete model in the operator form is

$$\left[\begin{array}{ccc}{A}_p& A_{\mathrm{\Gamma }}^T& \mathbf{0}\\ -A_{\mathrm{\Gamma }}& A_f& B_f^T\\ \mathbf{0}& B_f& \mathbf{0}\end{array}\right]\left[\begin{array}{l}{\varphi }_h\\ {\mathit{u}}_h\\ p_h\end{array}\right]=\left[\begin{array}{c}{\mathit{f}}_{f,h}\\ {\mathit{f}}_{p,h}\\ g_h\end{array}\right].$$E20

Here $A_f,A_p,A_{\mathrm{\Gamma }}$, and $B=\left(\mathbf{0},B_f\right)$ are the corresponding linear operators induced by the corresponding bilinear forms in (15) and (16). We denote

$$A=\left[\begin{array}{cc}{A}_p& A_{\mathrm{\Gamma }}^T\\ -A_{\mathrm{\Gamma }}& A_f\end{array}\right],\text{and}M=\left[\begin{array}{cc}\mathit{A}& B^T\\ B& \mathbf{0}\end{array}\right].$$E21

By discarding the coupling interface terms and plugging in $\frac{1}{\nu }\mathit{I}$ (which leads to pressure mass matrix) in the (2, 2) block, we have the following block-diagonal decoupled preconditioner:

$$P_M=\left[\begin{array}{cc}{\mathit{A}}_0& \mathbf{0}\\ \mathbf{0}& \frac{1}{\nu }{\mathit{I}}_h\end{array}\right],\text{with}{A}_0=\left[\begin{array}{cc}{A}_p& \mathbf{0}\\ \mathbf{0}& A_f\end{array}\right].$$E22

By keeping one of $B$ and $B^T$, one can easily construct block-triangular decoupled preconditioner. GMRES method is used as the outer iterative method. Block diagonal or block-triangular preconditioners are used in the inner iteration. The effectiveness and the efficiency of the preconditioners have been verified in [9, 39, 40]. In the implementation, $A_0^{-1}$ and the inverse for $\frac{1}{\nu }\mathit{I}$ should be realized by applying a Multigrid algorithm or domain decomposition methods. Particularly, when Krylov subspace methods are used, these inverses should be applied inexactly (for example, using one V-cycle Multigrid algorithm to provide approximate inverses).

Let us denote

$$P_{\pm }=\left(\begin{array}{cc}{A}_0& 0\\ 0& \pm \frac{1}{\nu }{\mathit{I}}_h\end{array}\right),$$E23

We also propose another preconditioner of the block triangular type:

$$P_{T_1}=\left(\begin{array}{cc}{A}_0& 0\\ B& -\frac{1}{\nu }{\mathit{I}}_h\end{array}\right),$$E24

by retaining the divergence operator. This preconditioner is still decoupled as the computation can be carried out in a block forward substitution manner. As an illustration, we illustrate the importance of using preconditioners in Table 1.

$h$	$\mathit{DOF}$	$N\left(I\right)$	$N\left(P_{-}\right)$	$N\left(P_{T_1}\right)$
$2^{-2}$	268	186	41	20
$2^{-3}$	948	432	45	22
$2^{-4}$	3556	$\ast \ast$	48	22
$2^{-5}$	13,764	$\ast \ast$	46	22
$2^{-6}$	54,148	$\ast \ast$	46	22

Table 1.

Number of iterations for the GMRES method without and with the two preconditioners $P_{-}$ and $P_{T_1}$.

We use $\ast \ast$ to indicate that the iteration does not converge within the prescribed maximum number of iterations. $N\left(P\right)$ refers to the number of iterations with a preconditioner $P$, and no preconditioner is applied when $P$ is simply $I$. From the table, it is clear that both $P_{-}$ and $P_{T_1}$ accelerate the convergence of the GMRES method. The number of iterations based on the two preconditioners is independent of the mesh refinement. More numerical experiments for testing the robustness with respect to the physical parameters can be found in [9, 39].

2.3 Decoupling and linearization by two-level and multi-level algorithms

For the mixed Stokes/Darcy model, Mu and Xu in [11] propose a two-grid method in which the coarse grid solution is used to supplement the boundary conditions at the interface for both of the two subproblems. The Two-grid Algorithm proposed in [11] is composed by the following two steps.

1. Solve the linear part of problem ${\left(2.7\right)}_H$ on a coarse grid: find ${\underset{¯}{u}}_H=\left({\mathit{u}}_H,{\varphi }_H\right)\in {\underset{¯}{X}}_H\subset {\underset{¯}{X}}_h,p_H\in Q_H\subset Q_h$ such that

   $$\left\{\begin{array}{l}a\left({\underset{¯}{u}}_H,{\underset{¯}{v}}_H\right)+b\left({\underset{¯}{v}}_H,p_H\right)=\left(f,{\underset{¯}{v}}_H\right),\forall {\underset{¯}{v}}_H=\left({\mathit{v}}_H,{\psi }_H\right)\in {\underset{¯}{X}}_H,\\ b\left({\underset{¯}{u}}_H,q_H\right)=0,\forall q_H\in Q_H;\end{array}\right.$$E25

2. Solve a modified fine grid problem: find ${\underset{¯}{u}}^H=\left({\mathit{u}}^h,{\varphi }^h\right)\in {\underset{¯}{X}}_h,p^h\in Q_h$ such that

$$\left\{\begin{array}{l}a\left({\underset{¯}{u}}^H,v_h\right)+b\left(v_h,p^h\right)=\left(f,v_h\right)-a_{\mathrm{\Gamma }}\left({\underset{¯}{u}}_H,v_h\right),\forall v_h\in {\underset{¯}{X}}_h,\\ b\left({\underset{¯}{u}}^H,q_h\right)=0,\forall q_h\in Q_h.\end{array}\right.$$E26

The main theoretical results for the two grid algorithm are as follows.

Theorem 1.Let$\left({\underset{¯}{u}}_h,p_h\right)$be the solution of coupled Stokes/Darcy model, and$\left({\underset{¯}{u}}^H,p^h\right)$be defined by and(27)on the fine grid. The following error estimates hold:

$${‖{\varphi }_h-{\varphi }^h‖}_{H_p}\lesssim H^2,$$E27

$${‖{\mathit{u}}_h-{\mathit{u}}^h‖}_{H_f}\lesssim H^{3/2},$$E28

$${‖p_h-p^h‖}_Q\lesssim H^{3/2}.$$E29

Based on this algorithm, some other improvements have been made. For example, in [18, 19, 20], by sequentially solving the Stokes submodel and the Darcy submodel, the authors can make ${‖{\varphi }_h-{\varphi }^h‖}_{H_p}$, ${‖{\mathit{u}}_h-{\mathit{u}}^h‖}_{H_f}$, and ${‖p_h-p^h‖}_Q$ are all of order $H^2$. Furthermore, Hou constructed a new auxiliary problem [16] for the Darcy submodel, and proved that Mu and Xu’s two-grid algorithm can retain ${‖{\mathit{u}}_h-{\mathit{u}}^h‖}_{H_f}$ and ${‖p_h-p^h‖}_Q$ order of $H^2$. It is remarkable that Mu and Xu’s two-grid algorithm is naturally parallel and of optimal order, if $h$ is of order $H^2$.

The extension to a multilevel decoupled algorithm can be found in [26]. The Multilevel Algorithm is as follows:

1. Solve the linear part of problem ${\left(2.7\right)}_H$ on a coarse grid: find ${\underset{¯}{u}}_H=\left({\mathit{u}}_H,{\varphi }_H\right)\in {\underset{¯}{X}}_H\subset {\underset{¯}{X}}_h,p_H\in Q_H\subset Q_h$ such that

   $$\left\{\begin{array}{l}a\left({\underset{¯}{u}}_H,{\underset{¯}{v}}_H\right)+b\left({\underset{¯}{v}}_H,p_H\right)=f\left({\underset{¯}{v}}_H\right),\forall {\underset{¯}{v}}_H=\left({\mathit{v}}_H,{\psi }_H\right)\in {\underset{¯}{X}}_H,\\ b\left({\underset{¯}{u}}_H,q_H\right)=0,\forall q_H\in Q_H;\end{array}\right.$$E30

2. Set $h_0=H$, for j = 1 to $L$,

   find $u^{h_j}=\left({\mathit{u}}^{h_j},{\varphi }^{h_j}\right)\in {\underset{¯}{X}}_{h_j},p^{h_j}\in Q_{h_j}$ such that

   $$\left\{\begin{array}{ll}a\left(u^{h_j},v_{h_j}\right)+b\left(v_{h_j},p^{h_j}\right)=f\left(v_{h_j}\right)-a_{\mathrm{\Gamma }}\left(u^{h_{j-1}},v_{h_j}\right),& \forall v_{h_j}\in {\underset{¯}{X}}_{h_j},\\ b\left(u^{h_j},q_{h_j}\right)=0,& \forall q_{h_j}\in Q_{h_j}.\end{array}\right.$$E31

   end.

In our multilevel algorithm, we refine the grid step by step, the coupled problem is only solved on the coarsest mesh, and linear decoupled subproblems are solved in parallel on successively refined meshes. We see that the algorithm is very effective and efficient. Moreover, the theory of the two grid algorithm guarantees that the approximation properties are good. As an illustration of the effectiveness of the multilevel algorithm [41], we present numerical results in Table 2.

$h$	${\left|{\varphi }^h-\varphi \right|}_1$	${\left|u^h-u\right|}_1$	${\left|v^h-v\right|}_1$	${‖p^h-p‖}_0$
$2^{-1}$	$4.592\times {10}^{-2}$	$1.550\times {10}^{-1}$	$1.066\times {10}^{-1}$	$8.410\times {10}^{-2}$
$2^{-2}$	$1.152\times {10}^{-2}$	$3.958\times {10}^{-2}$	$2.664\times {10}^{-2}$	$1.752\times {10}^{-2}$
$2^{-4}$	$7.280\times {10}^{-4}$	$2.466\times {10}^{-3}$	$1.652\times {10}^{-3}$	$1.040\times {10}^{-3}$
$2^{-8}$	$5.296\times {10}^{-6}$	$9.981\times {10}^{-6}$	$7.922\times {10}^{-6}$	$1.694\times {10}^{-5}$

Table 2.

Errors between the solutions of multilevel algorithm and the exact solutions (second order discretization).

In the following, we use steady state NS/Darcy model to illustrate how to apply two-level and multilevel methods to decouple the coupled nonlinear PDE models. The algorithm combines the two-level algorithms and the Newton-type linearization [4, 36, 42]. Our Newton Type Linearization Based Two-level Algorithm consists of the following three steps [17, 43].

1. Solve the coupled problem (19) on a coarse grid triangulation with the meshsize $H$: Find ${\underset{¯}{u}}_H=\left({\mathit{u}}_H,{\varphi }_H\right)\in {\underset{¯}{X}}_H$ and $p_H\in Q_H$ such that

   $$\left\{\begin{array}{ll}a\left({\underset{¯}{u}}_H,{\underset{¯}{v}}_H\right)+c\left({\mathit{u}}_H,{\mathit{u}}_H,{\mathit{v}}_H\right)+b\left({\mathit{v}}_H,p_H\right)=f\left({\underset{¯}{v}}_H\right),& \forall {\underset{¯}{v}}_H=\left({\mathit{v}}_H,{\psi }_H\right)\in {\underset{¯}{X}}_H,\\ b\left({\mathit{u}}_H,q_H\right)=0,& \forall q_H\in Q_H.\end{array}\right.$$E32

2. On a fine grid triangulation with the meshsize $h\le H$, sequentially solve two decoupled and linearized local subproblems:

   Step a. Solve a discrete Darcy problem in ${\mathrm{\Omega }}_p$: Find ${\varphi }_h^{\ast }\in X_{p,h}$ such that

   $$a_p\left({\varphi }_h^{\ast },{\psi }_h\right)=\mathit{\rho g}{\left(f_p,{\psi }_h\right)}_{{\mathrm{\Omega }}_p}+\mathit{\rho g}{\left({\mathit{u}}_H\cdot {\mathit{n}}_f,{\psi }_h\right)}_{\mathrm{\Gamma }}\forall {\psi }_h\in X_{p,h}.$$E33

   Step b. Solve a modified Navier–Stokes model using the Newton type linearization: Find ${\mathit{u}}_h^{\ast }\in {\mathit{X}}_{f,h}$ and $p_h^{\ast }\in Q_h$ such that

   $$\left\{\begin{array}{ll}{a}_f\left({\mathit{u}}_h^{\ast },{\mathit{v}}_h\right)+c\left({\mathit{u}}_H,{\mathit{u}}_h^{\ast },{\mathit{v}}_h\right)+c\left({\mathit{u}}_h^{\ast },{\mathit{u}}_H,{\mathit{v}}_h\right)+b\left({\mathit{v}}_h,p_h^{\ast }\right)={\left({\mathit{f}}_f,{\mathit{v}}_h\right)}_{{\mathrm{\Omega }}_f}& \\ -\mathit{\rho g}{\left({\varphi }_h^{\ast },{\mathit{v}}_h\cdot {\mathit{n}}_f\right)}_{\mathrm{\Gamma }}+c\left({\mathit{u}}_H,{\mathit{u}}_H,{\mathit{v}}_h\right)& \forall {\mathit{v}}_h\in {\mathit{X}}_{f,h},\\ b\left({\mathit{u}}_h^{\ast },q_h\right)=0& \forall q_h\in Q_h.\end{array}\right.$$E34

3. On the same fine grid triangulation, solve two subproblems by using the newly obtained solution.

   Step a. Solve a discrete Darcy problem in ${\mathrm{\Omega }}_p$: Find ${\varphi }^h\in X_{p,h}$ such that

   $$a_p\left({\varphi }^h,{\psi }_h\right)=\mathit{\rho g}{\left(f_p,{\psi }_h\right)}_{{\mathrm{\Omega }}_p}+\mathit{\rho g}{\left({\mathit{u}}_h^{\ast }\cdot {\mathit{n}}_f,{\psi }_h\right)}_{\mathrm{\Gamma }}\forall {\psi }_h\in X_{p,h}.$$E35

   Step b. Correct the solution of the fluid flow model: Find ${\mathit{u}}^h\in {\mathit{X}}_{f,h}$ and $p^h\in Q_h$ such that

   $$\left\{\begin{array}{ll}{a}_f\left({\mathit{u}}^h,{\mathit{v}}_h\right)+c\left({\mathit{u}}_H,{\mathit{u}}^h,{\mathit{v}}_h\right)+c\left({\mathit{u}}^h,{\mathit{u}}_H,{\mathit{v}}_h\right)+b\left({\mathit{v}}_h,p^h\right)={\left({\mathit{f}}_f,{\mathit{v}}_h\right)}_{{\mathrm{\Omega }}_f}& \\ -\mathit{\rho g}{\left({\varphi }^h,{\mathit{v}}_h\cdot {\mathit{n}}_f\right)}_{\mathrm{\Gamma }}+c\left({\mathit{u}}_H,{\mathit{u}}_h^{\ast },{\mathit{v}}_h\right)+c\left({\mathit{u}}_h^{\ast },{\mathit{u}}_H-{\mathit{u}}_h^{\ast },{\mathit{v}}_h\right)& \forall {\mathit{v}}_h\in {\mathit{X}}_{f,h},\\ b\left({\mathit{u}}^h,q_h\right)=0& \forall q_h\in Q_h.\end{array}\right.$$E36

We remark here that the problem (36) and the problem (34) differ only in the right hand side. Similarly, the problem (35) and the problem (33) have the same stiffness matrix. In sum, the advantages of our work exist in that the scaling between the two meshsizes is better, the algorithm is decoupled and linear on the fine grid level and the two submodels in the last two steps share the same stiffness matrices.

For the coupled problem (19), by using the properties (16)–(18), the error estimates in the energy norm can be derived by using a fixed-point framework [4, 31]. Moreover, the Aubin-Nitsche duality argument can result in the $L^2$ error analysis of the problem (19). In summary, we have.

Lemma 1.Let$\left(\mathit{u},\varphi ,p\right)\in {\mathit{H}}^{k+1}\left({\mathrm{\Omega }}_f\right)\times H^{k+1}\left({\mathrm{\Omega }}_p\right)\times H^k\left({\mathrm{\Omega }}_f\right)$be the solution of the Navier–Stokes/Darcy model(13)and$\left({\mathit{u}}_h,{\varphi }_h,p_h\right)$be the FE solution of(19). We assume that$\nu$is sufficiently large and$h$is sufficiently small. There holds the following energy norm estimate for the problem(19).

$${\left|\mathit{u}-{\mathit{u}}_h\right|}_{1,{\mathrm{\Omega }}_f}+{\left|\varphi -{\varphi }_h\right|}_{1,{\mathrm{\Omega }}_p}+\parallel p-p_h{\parallel }_{0,{\mathrm{\Omega }}_f}\lesssim h^k.$$E37

Moreover, we have the following $L^2$ error estimate:

$$\parallel \mathit{u}-{\mathit{u}}_h{\parallel }_{0,{\mathrm{\Omega }}_f}+\parallel \varphi -{\varphi }_h{\parallel }_{0,{\mathrm{\Omega }}_p}\lesssim h^{k+1}.$$E38

The energy norm estimate (37) is the Lemma 2 in [4]. The $L^2$ error estimate (38) corresponds to the Lemma 3 in [4]. Detailed proofs of these results can be found in [4].

The following lemma concludes the error estimate of $\left({\varphi }_h^{\ast },{\mathit{u}}_h^{\ast },p_h^{\ast }\right)$ in the energy norm.

Lemma 2.(Error analysis of the intermediate step two-level solution) Let$\left(\varphi ,\mathit{u},p\right)$and$\left({\varphi }_h^{\ast },{\mathit{u}}_h^{\ast },p_h^{\ast }\right)$be defined by the problems(13)and(33)–(34), respectively. Under the assumptions of Lemma 1, there holds

$${\left|\varphi -{\varphi }_h^{\ast }\right|}_{1,{\mathrm{\Omega }}_p}\lesssim H^{k+1}+h^k,$$E39

$${\left|\mathit{u}-{\mathit{u}}_h^{\ast }\right|}_{1,{\mathrm{\Omega }}_f}+\parallel p-p_h^{\ast }{\parallel }_{0,{\mathrm{\Omega }}_f}\lesssim H^{k+1}+h^k.$$E40

$$\parallel \mathit{u}-{\mathit{u}}_h^{\ast }{\parallel }_{0,{\mathrm{\Omega }}_f}\lesssim H^{2k+1}+H^{k+1}h+h^{k+1}.$$E41

From Lemma 1, we note that the optimal finite element solution error in the energy norm is of order $O\left(h^k\right)$. Combining the conclusions of this Lemma, we see that the intermediate step two-level solution error is still optimal in the energy norm if the scaling between the two meshsizes is taken to be $h=H^{\frac{k+1}{k}}$. The $L^2$ error analysis is extended from [4, 20, 31, 37]. Let $\mathbf{u}$ and $\varphi$ be the nonsingular solution of (13). From Lemma 1, the optimal $L^2$ error for the finite element solution is of order $O\left(h^{k+1}\right)$. To make sure $\parallel \mathit{u}-{\mathit{u}}_h^{\ast }{\parallel }_{0,{\mathrm{\Omega }}_f}$ is also of order $O\left(h^{k+1}\right)$, the scaling between the two grids has to be taken as $h=max\left\{H^{\frac{k+1}{k}},H^{\frac{2k+1}{k+1}}\right\}$. For instance, if $k=1$, we have to set $h=H^{\frac{3}{2}}$ to make sure the $L^2$ error of ${\mathit{u}}_h^{\ast }$ is optimal. We now show that the final step two-level solution is indeed a good approximation to the solution of problem (13) .

Theorem 2.(Error analysis of the final step two-level solution) Let$\left(\varphi ,\mathit{u},p\right)$and$\left({\varphi }^h,{\mathit{u}}^h,p^h\right)$be the solutions of(13)and(35)–(36)respectively. Under the assumptions of Lemma 1, the following error estimates hold:

$${\left|\varphi -{\varphi }^h\right|}_{1,{\mathrm{\Omega }}_p}+{\left|\mathit{u}-{\mathit{u}}^h\right|}_{1,{\mathrm{\Omega }}_f}+\parallel p-p^h{\parallel }_{0,{\mathrm{\Omega }}_f}\lesssim H^{2k+1}+H^{k+1}h+h^k,$$E42

Proposition 1.Let$\left(\varphi ,\mathit{u},p\right)$and$\left({\varphi }^h,{\mathit{u}}^h,p^h\right)$be the solutions of(13)and(35)–(36)respectively. If we take$h=H^{\frac{2k+1}{k}}$for$k=1,2$and$h=H^{\frac{k+1}{k-1}}$for$k\ge 3$, then there holds the following error estimate.

$${\left|\varphi -{\varphi }^h\right|}_{1,{\mathrm{\Omega }}_p}+{\left|\mathit{u}-{\mathit{u}}^h\right|}_{1,{\mathrm{\Omega }}_f}+\parallel p-p^h{\parallel }_{0,{\mathrm{\Omega }}_f}\lesssim h^k.$$E43

Finally, we would like to make some comments on the mixed Stokes/Darcy model. We note that by dropping those trilinear terms (32), (34) and (36), our two-level algorithm can be naturally applied to the coupled Stokes/Darcy model. We note that the above algorithms can be naturally extended to multi-level algorithms by recursively calling the above two-level algorithms [17, 43]. The extension and the corresponding analysis can be found in [26, 43].

2.4 Decoupled algorithms by partitioned time schemes

The fully evolutionary Stokes/Darcy equations will be used as the model problem in this subsection to illustrate the partitioned time schemes. We neglect the fluid density and the porosity effects in this subsection. We review some decoupled methods that converge within a reasonable amount of time, and are stable when the physical parameters are small. More precisely, partitioned time methods can efficiently solve the surface subproblem and the subsurface subproblem separately.

For the estimate of the stability, we assume that the solution of the Stokes/Darcy problem is long-time regular [44]:

$$\begin{array}{l}\mathit{u}\in {\mathit{W}}^{2,\infty }\left(0,\infty ;{\mathit{L}}^2\left({\mathrm{\Omega }}_f\right)\right)\cap {\mathit{W}}^{1,\infty }\left(0,\infty ;{\mathit{H}}^2\left({\mathrm{\Omega }}_f\right)\right),\\ \varphi \in W^{2,\infty }\left(0,\infty ;L^2\left({\mathrm{\Omega }}_p\right)\right)\cap W^{1,\infty }\left(0,\infty ;H^2\left({\mathrm{\Omega }}_p\right)\right),\\ p\in L^{2,\infty }\left(0,\infty ;H^1\left({\mathrm{\Omega }}_f\right)\right).\end{array}$$E44

The simplest time scheme for the evolutionary coupled Stokes/Darcy model is the Backward Euler Algorithm, which reads as: Given $\left({\mathit{u}}_h^n,p_h^n,{\varphi }_h^n\right)\in {\mathit{X}}_{f,h}\times Q_h\times X_{p,h}$, find $\left({\mathit{u}}_h^{n+1},p_h^{n+1},{\varphi }_h^{n+1}\right)\in {\mathit{X}}_{f,h}\times Q_h\times X_{p,h}$, such that for all ${\mathit{v}}_h\in {\mathit{X}}_{f,h},q_h\in Q_h,{\psi }_h\in X_{p,h}$,

$$\left\{\begin{array}{l}\left(\frac{{\mathit{u}}_h^{n+1}-{\mathit{u}}_h^n}{\mathrm{\Delta }t},{\mathit{v}}_h\right)+a_f\left({\mathit{u}}_h^{n+1},{\mathit{v}}_h\right)-\left(p_h^{n+1},\nabla \cdot {\mathit{v}}_h\right)+g{\left({\mathit{v}}_h\cdot {\mathit{n}}_f,{\varphi }_h^{n+1}\right)}_{\mathrm{\Gamma }}=\left({\mathit{f}}_f^{n+1},{\psi }_h\right),\\ \left(q_h,\nabla \cdot {\mathit{u}}_h^{n+1}\right)=0,\end{array}\right.$$E45

$${\mathit{gS}}_0\left(\frac{{\varphi }_h^{n+1}-{\varphi }_h^{n-1}}{\mathrm{\Delta }t},{\mathit{v}}_h\right)+a_p\left({\varphi }_h^{n+1},{\psi }_h\right)-g{\left({\mathit{u}}_h^{n+1}\cdot {\mathit{n}}_f,{\psi }_h\right)}_{\mathrm{\Gamma }}=g\left(f_p^{n+1},{\psi }_h\right).$$E46

However, this scheme is fully coupled and each time step one has to solve a coupled system including both (45) and (46), although, on the other hand, this scheme enjoys the desirable strong stability and convergence properties. In [12], Mu and Zhu propose the following backward Euler forward Euler scheme and combine it with the two-level spatial discretization. We neglect the two-level spatial discretization in this presentation. Here, the Forward Euler means it discretizes the coupling term explicitly. Backward Euler Forward Euler Scheme (BEFE): given $\left({\mathit{u}}_h^n,p_h^n,{\varphi }_h^n\right)\in {\mathit{X}}_{f,h}\times Q_h\times X_{p,h}$, find $\left({\mathit{u}}_h^{n+1},p_h^{n+1},{\varphi }_h^{n+1}\right)\in {\mathit{X}}_{f,h}\times Q_h\times X_{p,h}$, such that for all ${\mathit{v}}_h\in {\mathit{X}}_{f,h},q_h\in Q_h,{\psi }_h\in X_{p,h}$,

$$\left\{\begin{array}{l}\left(\frac{{\mathit{u}}_h^{n+1}-{\mathit{u}}_h^n}{\mathrm{\Delta }t},{\mathit{v}}_h\right)+a_f\left({\mathit{u}}_h^{n+1},{\mathit{v}}_h\right)-\left(p_h^{n+1},\nabla \cdot {\mathit{v}}_h\right)+g{\left({\mathit{v}}_h\cdot {\mathit{n}}_f,{\varphi }_h^n\right)}_{\mathrm{\Gamma }}=\left({\mathit{f}}_f^{n+1},{\psi }_h\right),\\ \left(q_h,\nabla \cdot {\mathit{u}}_h^{n+1}\right)=0,\end{array}\right.$$E47

$$\begin{array}{l}{\mathit{gS}}_0\left(\frac{{\varphi }_h^{n+1}-{\varphi }_h^n}{\mathrm{\Delta }t},{\mathit{v}}_h\right)+a_p\left({\varphi }_h^{n+1},{\psi }_h\right)-g{\left({\mathit{u}}_h^n\cdot {\mathit{n}}_f,{\psi }_h\right)}_{\mathrm{\Gamma }}=g\left(f_p^{n+1},{\psi }_h\right).\end{array}$$E48

The analysis of this can be found in [12, 45]. In particular, the longtime stability (cf. (51)) of BEFE method was proved in [45] in the sense that no form of Gronwall’s inequality was used.

Theorem 3.Assume the following time step condition is satisfied

$$\mathrm{\Delta }t\lesssim \mathit{min}\left\{\nu k_{\mathit{min}}^2,S_0{\nu }^2{k}_{\mathit{min}}\right\},$$E49

Then, BELF algorithm achieves the optimal convergence rate uniformly in time. The solution of the BEFE method satisfies the uniform in time error estimates:

1. If ${\mathit{f}}_f\in L^{\infty }\left(0,+\infty ;{\mathit{L}}^2\left({\mathrm{\Omega }}_f\right)\right),f_p\in L^{\infty }\left(0,+\infty ;L^2\left({\mathrm{\Omega }}_p\right)\right)$, then

   $${‖{\mathit{u}}_h^n‖}^2+{‖{\varphi }_h^n‖}^2\le C,\forall n\ge 0.$$E50

2. If ${‖{\mathit{f}}_f‖}_{L^{\infty }\left(0,+\infty ;{\mathbf{L}}^2\left({\mathrm{\Omega }}_f\right)\right)},{‖f_p‖}_{L^{\infty }\left(0,+\infty ;L^2\left({\mathrm{\Omega }}_p\right)\right)}$ are uniformly bounded in $\mathrm{\Delta }t$, then

   $${‖{\mathit{u}}_h^n‖}^2+{‖{\varphi }_h^n‖}^2+\mathrm{\Delta }t\sum _{l=0}^n\left({‖\nabla {\mathit{u}}_h^l‖}^2+{‖\nabla {\varphi }_h^l‖}^2\right)\le C,\forall n\ge 0.$$E51

The advantage of this scheme is that it is parallel. As revealed in the time step restriction (50), the disadvantage of this method is that it may become highly unstable when the parameters $S_0$ and $k_{\mathit{min}}$ are small. Another way for uncoupling surface/subsurface flow models is using splitting schemes which require sequential sub-problem solves at each time step [46]. As an example, we note that in solving (49), one can replace ${\mathit{u}}_h^n$ by using the most updated solution obtained in the Stokes step. This will lead to Backward Euler time-split scheme [46]. We skip the details of this time-split method, interested readers can refer to [46]. By this way, one can design different sequential splitting schemes. Noting that the BEFE method is only of first order, in some other decoupled Implicit-explicit (IMEX) methods, one can combine of the three level implicit method with the coupling terms treated by the explicit method to achieve high order. For example, Crank–Nicolson Leap-Frog method [47, 48], second-order backward-differentiation with Gear’s extrapolation, Adam-Moulton-Bashforth [49]. We present one of them: the Crank–Nicolson Leap-Frog Method for the evolutionary Stokes/Darcy model: given $\left({\mathit{u}}_h^{n-1},p_h^{n-1},{\varphi }_h^{n-1}\right),\left({\mathit{u}}_h^n,P_h^n,{\varphi }_h^n\right)\in {\mathit{X}}_{f,h}\times Q_h\times X_{p,h}$, find $\left({\mathit{u}}_h^{n+1},p_h^{n+1},{\varphi }_h^{n+1}\right)\in {\mathit{X}}_{f,h}\times Q_h\times X_{p,h}$, such that for all ${\mathit{v}}_h\in {\mathit{X}}_{f,h},q_h\in Q_h,{\psi }_h\in X_{p,h}$,

$$\left\{\begin{array}{l}\left(\frac{{\mathit{u}}_h^{n+1}-{\mathit{u}}_h^{n-1}}{2\mathrm{\Delta }t},{\mathit{v}}_h\right)+\left(\nabla \cdot \left(\frac{{\mathit{u}}_h^{n+1}-{\mathit{u}}_h^{n-1}}{2\mathrm{\Delta }t},{\mathit{v}}_h\right),\nabla \cdot {\mathit{v}}_h\right)+a_f\left(\frac{{\mathit{u}}_h^{n+1}+{\mathit{u}}_h^{n-1}}{2},{\mathit{v}}_h\right)\\ -\left(\frac{p_h^{n+1}-p_h^{n-1}}{2},\nabla \cdot {\mathit{v}}_h\right)+g{\left({\mathit{v}}_h,{\varphi }_h^n\right)}_{\mathrm{\Gamma }}=\left({\mathit{f}}_f^n,{\mathit{v}}_h\right),\\ \left(q_h,\nabla \cdot \frac{{\mathit{u}}_h^{n+1}-{\mathit{u}}_h^{n-1}}{2}\right)=0,\end{array}\right.$$E52

$$\begin{array}{l}{\mathit{gS}}_0\left(\frac{{\varphi }_h^{n+1}-{\varphi }_h^{n-1}}{2\mathrm{\Delta }t},{\varphi }_h\right)+a_p{\left(\frac{{\varphi }_h^{n+1}+{\varphi }_h^{n-1}}{2},{\psi }_h\right)}_h-g{\left({\mathit{u}}_h^n,{\psi }_h\right)}_{\mathrm{\Gamma }}=g\left(f_p^n,{\psi }_h\right).\end{array}$$E53

The Crank–Nicolson-Leap-Frog possesses strong stability and convergence properties [47, 48]. Most importantly, the time-step condition for the scheme does not depend on ${\kappa }_{\mathit{min}}$.

For the numerical experiments of these partitioned time schemes, we refer the readers to [12].

3. FSI models and decoupled algorithms

FSI models include a fluid model whose general form is (1), a structure model, plus certain interface conditions that describe the interaction mechanism (see Figure 2 for an illustration of FSI models in the reference configuration and the deformed configuration). To differentiate the notations in different subdomains, we will use a subscript “$f$” to denote the variables in the fluid domain, and a subscript “$s$” to denote the variables in the structure domain.

Figure 2.

An illustration of fluid structure interaction: the reference configuration (left) and the deformed configuration (right).

In general, the structure model reads as:

$$\left\{\begin{array}{l}{\rho }_s\left(\frac{\partial {\mathit{u}}_s}{\partial t}+\left({\mathit{u}}_s\cdot \nabla \right){\mathit{u}}_s\right)-\mathit{d}\mathit{iv}\mathit{T}\left({\mathit{u}}_s,p_s\right)={\mathit{f}}_s,\forall \mathit{x}\in {\mathrm{\Omega }}_s,\\ div{\mathit{u}}_s=0,\forall \mathit{x}\in {\mathrm{\Omega }}_s.\end{array}\right.$$E54

Here, ${\rho }_s$ is the density of the structure, ${\mathit{u}}_s$ is the structure velocity, $p_s$ is the structure pressure. In structure mechanics, the displacement $\mathit{d}$ is usually used as a primary variable ($\dot{\mathit{d}}={\mathit{u}}_s$), and the stress term in the linear case can be described by using Hooke’s law. As the structure model is usually based on Lagrangian coordinates, researchers usually introduce the so-called Arbitrary Lagrangian Eulerian description for FSI models. In some special cases of FSI models, one can apply the simplified structure model such as 1D structure model or linear elasticity model for structure part, and the simplified fluid model such as linear Stokes or inviscid flow model for fluid part.

The fluid motion and the structure motion are coupled through certain interface conditions that describe the compatibility of the kinematics and transactions at the fluid–structure interface. For applications with non-slip interface conditions, both velocity and normal stress are continuous across the interface $\mathrm{\Gamma }$, which may be described as

$$\left\{\begin{array}{l}{\mathit{u}}_f={\mathit{u}}_s,\mathrm{on}\mathrm{\Gamma },\\ \mathit{T}\left({\mathit{u}}_f,p_f\right)\mathit{n}=\mathit{T}\left({\mathit{u}}_s,p_s\right)\mathit{n},\mathrm{on}\mathrm{\Gamma }.\end{array}\right.$$E55

Here, $\mathit{n}$ denotes one of ${\mathit{n}}_f$ and ${\mathit{n}}_s$. With suitable initial conditions and boundary conditions such as fluid flux boundary condition and structure Dirichlet boundary condition, the FSI models are complete.

There have been many advanced numerical methods for various FSI models. Our focus here is a linear fluid model coupled with a thin wall structure. The reason we choose this model is that this kind of models, if calculated using the standard (Dirichlet-Neumann) explicit decoupling schemes will lead to the so-called added mass effect [50]. Moreover, the algorithms dealing with the added-mass difficulties are the most exciting development in this direction. Therefore, we consider dropping the nonlinear term in (1), and consider the following 1D structure model:

$$\left\{\begin{array}{l}{\rho }_s\epsilon {\partial }_t\dot{\mathit{d}}+{\mathit{L}}_s\left(\mathit{d},\dot{\mathit{d}}\right)=-\mathit{T}\left({\mathit{u}}_f,p_f\right)\mathit{n},\mathrm{on}\mathrm{\Gamma },\\ {\mathit{u}}_s=\dot{\mathit{d}}\mathrm{on}\mathrm{\Gamma }.\end{array}\right.$$E56

Here, $\epsilon$ denotes the structure thickness, ${\mathit{L}}_s\left(\mathit{d},\dot{\mathit{d}}\right)$ denotes the operator in the structure model, which may include both the elastic term and the damping term [14, 15].

3.1 Partitioned algorithms for FSI models

First of all, we comment here that the decoupled preconditioning techniques can also be naturally applied to FSI models. In the preconditioning step, one can apply either the Multigrid approach [51] or domain decomposition methods [52].

In this presentation, we focus on the two most recent approaches for the linear Stokes model coupled with thin wall structure. The first approach is called partitioned Robin-Neumann scheme, in which the fluid subproblem is imposed with Robin boundary condition at the interface while the structure subproblem is imposed with Neumann boundary conditioner at the interface [14, 15]. The second approach is called kinetically coupled $\beta$-scheme, which is actually decoupled in the sense that computations are realized subdomain by subdomain [13, 21]. The derivation of the $\beta$-scheme is based on operator splitting.

Partitioned Robin-Neumann scheme:

1. Given the initial solution ${\mathit{u}}_f^0$,$p_f^0$ and ${\mathit{d}}^0$.

2. For $m=\mathrm{0,1,2,3},\dots ,N-1$,

   • Fluid step: find ${\mathit{u}}_f^{m+1}$ and $p_f^{m+1}$ such that

     $$\left\{\begin{array}{lll}\frac{{\rho }_f}{\mathrm{\Delta }t}\left({\mathit{u}}_f^{m+1}-{\mathit{u}}_f^m\right)-div\mathit{T}\left({\mathit{u}}_f^{m+1},p_f^{m+1}\right)=\mathbf{0}& \text{in}& {\mathrm{\Omega }}_f,\\ div{\mathit{u}}_f^{m+1}=0& \text{in}& {\mathrm{\Omega }}_f,\\ \mathit{T}\left({\mathit{u}}_f^{m+1},p_f^{m+1}\right)\mathit{n}+\frac{{\rho }_s\epsilon }{\mathrm{\Delta }t}{\mathit{u}}_f^{m+1}=\frac{{\rho }_s\epsilon }{\mathrm{\Delta }t}{\dot{\mathit{d}}}^m+\mathit{T}\left({\mathit{u}}_f^m,p^m\right)\mathit{n}& \mathrm{on}& \mathrm{\Gamma }.\end{array}\right.$$E57

   • Structure step: find ${\mathbf{d}}_s^{m+1}$ such that

     $$\left\{\begin{array}{lll}\frac{{\rho }_s\epsilon }{\mathrm{\Delta }t}\left({\dot{\mathit{d}}}^{m+1}-{\dot{\mathit{d}}}^m\right)+{\mathit{L}}_s\left({\mathit{d}}^{m+1},{\dot{\mathit{d}}}^{m+1}\right)=-\mathit{T}\left({\mathit{u}}_f^{m+1},p_f^{m+1}\right)\mathit{n}& \mathrm{on}& \mathrm{\Gamma },\\ {\dot{\mathit{d}}}^{m+1}=\frac{1}{\mathrm{\Delta }t}\left({\mathit{d}}^{m+1}-{\mathit{d}}^m\right)& \mathrm{on}& \mathrm{\Gamma },.\end{array}\right.$$E58

3. End

Here, $\frac{{\rho }_s\epsilon }{\mathrm{\Delta }t}$ is treated as a Robin coefficient. These partitioned iterative methods were firstly introduced in [53], as added-mass free alternatives to the standard Dirichlet-Neumann scheme. Some extensions and generalizations can be found in [14, 15].

Different from the partitioned Robin-Neumann scheme. The kinematically coupled $\beta$-scheme for the time-discrete problem is given as follows. The stability and the convergence rate of this scheme are analyzed in [13, 21].

1. Given the initial solution ${\mathit{u}}_f^0$,$p_f^0$ and ${\mathit{d}}^0$.

2. For $m=\mathrm{0,1,2,3},\dots ,N-1$,

   • Structure step: find ${\stackrel{\sim }{\mathit{u}}}_s^{m+1}$ such that

     $$\left\{\begin{array}{ll}{\rho }_s\epsilon \frac{{\stackrel{\sim }{\mathit{u}}}_s^{m+1}-{\mathit{u}}_s^m}{\mathrm{\Delta }t}+{\mathbf{L}}_s\left({\mathit{d}}^{m+1},{\dot{\mathit{d}}}^{m+1}\right)=-{\mathit{\beta \sigma }}_f\left({\mathit{u}}_f^m,p_f^m\right)\mathit{n}& \mathrm{on}\mathrm{\Gamma },\\ {\dot{\mathit{d}}}^{m+1}={\stackrel{\sim }{\mathit{u}}}_s^{m+1},{\mathit{d}}^{m+1}={\mathit{d}}^m+\mathrm{\Delta }t{\stackrel{\sim }{\mathit{u}}}_s^{m+1}& \mathrm{on}\mathrm{\Gamma }.\end{array}\right.$$E59

   • Fluid step: find ${\mathit{u}}_f^{m+1}$, $p_f^{m+1}$ and ${\mathit{u}}_s^{m+1}$ such that

     $$\left\{\begin{array}{ll}\frac{{\rho }_f}{\mathrm{\Delta }t}\left({\mathit{u}}_f^{m+1}-{\mathit{u}}_f^m\right)-div{\sigma }_f\left({\mathit{u}}_f^{m+1},p_f^{m+1}\right)=\mathbf{0}& \text{in}{\mathrm{\Omega }}_f,\\ div{\mathit{u}}_f^{m+1}=0& \text{in}{\mathrm{\Omega }}_f,\\ {\rho }_s\epsilon \frac{{\mathit{u}}_s^{m+1}-{\stackrel{\sim }{\mathit{u}}}_s^{m+1}}{\mathrm{\Delta }t}=-{\sigma }_f\left({\mathit{u}}_f^{m+1},p_f^{m+1}\right)\mathit{n}+{\mathit{\beta \sigma }}_f\left({\mathit{u}}_f^m,p_f^m\right)\mathit{n}& \mathrm{on}\mathrm{\Gamma },\\ {\mathit{u}}_f^{m+1}={\mathit{u}}_s^{m+1}& \mathrm{on}\mathrm{\Gamma }.\end{array}\right.$$E60

3. End.

3.2 Multirate time step approach for FSI models

Due to different time scales in many FSI problems, it is natural and essential to develop multirate time-stepping schemes [22, 23] that mimic the physical phenomena. For illustration, we will examine the application of the multirate technique to the $\beta$-scheme, since similar performance is observed for both the Robin-Neumann scheme and the $\beta$-scheme in numerical experiments. Furthermore, the decoupled multirate $\beta$-scheme can be extended to more general FSI problems involving nonlinearity, irregular domains, and large structural deformations.

In order to apply multirate time-stepping scheme to FSI problems, a question is in which region the larger time step size should be used. Numerical experiments suggest that the version with a larger time step size in the fluid solver (cf. Figure 3) results in a better accuracy. The corresponding method is described in the following. We comment here that the multirate $\beta$-scheme is nothing else but the $\beta$-scheme itself when the time-step ratio $r=1$.

Figure 3.

An illustration of a multirate time stepping technique.

Multirate time-stepping $\beta$-scheme:

1. Given the initial solution ${\mathit{u}}_f^0$,$p_f^0$ and ${\mathit{d}}^0$.

2. For $k=\mathrm{0,1,2,3},\dots ,N-1$, $\mathrm{set}{m}_k=r\cdot k$.

   • Structure steps:$\text{for}m=m_k,m_k+1,m_k+2,\dots ,m_{k+1}-1$,

     $$\left\{\begin{array}{ll}{\rho }_s\epsilon \frac{{\stackrel{\sim }{\mathit{u}}}_s^{m+1}-{\mathit{u}}_s^m}{\mathrm{\Delta }{t}_s}+{\mathit{L}}_s\left({\mathit{d}}^{m+1},{\dot{\mathit{d}}}^{m+1}\right)=-{\mathit{\beta \sigma }}_f\left({\mathit{u}}_f^{m_k},p_f^{m_k}\right)\mathit{n}& \mathrm{on}\mathrm{\Gamma },\\ {\dot{\mathit{d}}}^{m+1}={\stackrel{\sim }{\mathit{u}}}_s^{m+1},{\mathit{d}}^{m+1}={\mathit{d}}^m+\mathrm{\Delta }{t}_s{\stackrel{\sim }{\mathit{u}}}_s^{m+1}& \mathrm{on}\mathrm{\Gamma }.\end{array}\right.$$E61

   • Fluid step:

     $$\left\{\begin{array}{ll}\frac{{\rho }_f}{\mathrm{\Delta }{t}_f}\left({\mathit{u}}_f^{m_{k+1}}-{\mathit{u}}_f^{m_k}\right)-div{\sigma }_f\left({\mathit{u}}_f^{m_{k+1}},p_f^{m_{k+1}}\right)=\mathbf{0}& \text{in}{\mathrm{\Omega }}_f,\\ div{\mathit{u}}_f^{m_{k+1}}=0& \text{in}{\mathrm{\Omega }}_f,\\ {\rho }_s\epsilon \frac{{\mathit{u}}_s^{m_{k+1}}-{\stackrel{\sim }{\mathit{u}}}_s^{m_{k+1}}}{\mathrm{\Delta }{t}_f}=-{\sigma }_f\left({\mathit{u}}_f^{m_{k+1}},p_f^{m_{k+1}}\right)\mathit{n}+{\mathit{\beta \sigma }}_f\left({\mathit{u}}_f^{m_k},p_f^{m_k}\right)\mathit{n}& \mathrm{on}\mathrm{\Gamma },\\ {\mathit{u}}_f^{m_{k+1}}={\mathit{u}}_s^{m_{k+1}}& \mathrm{on}\mathrm{\Gamma }.\end{array}\right.$$E62

3. End.

3.3 Numerical experiment

In this subsection, we present numerical experiments to demonstrate the convergence and stability performance of the multirate $\beta$-scheme (60)–(61) for coupling a Stokes flow with a thin-walled structure by using a benchmark model. The model consists of a 2-D rectangular fluid domain ${\mathrm{\Omega }}_f=\left[0,L\right]\times \left[0,R\right]$ with $L=6$ and $R=0.5$ and a 1-D structure domain $\mathrm{\Gamma }=\left[0,L\right]\times R$ that meanwhile also plays the role of the fluid–solid interface, as shown in Figure 2. The displacement of the interface is assumed to be infinitesimal and the Reynolds number in the fluid is assumed to be small (Figure 4). All the quantities will be given in terms of the centimeter-gram-second (CGS) system of units.

Figure 4.

Geometrical configuration.

The physical parameters are set as follows: ${\rho }_f=1.0$, ${\rho }_s=1.1$, $\mu =0.035$; ${\mathit{L}}_s\left(\mathit{d},\dot{\mathit{d}}\right)=c_1{\partial }_x^2\mathit{d}+c_0\mathit{d}$, where $c_1=\frac{\mathit{E\epsilon }}{2\left(1+\nu \right)}$, $c_0=\frac{\mathit{E\epsilon }}{R^2\left(1-{\nu }^2\right)}$ with $\epsilon =0.1$, the Poisson ratio $\nu =0.5$ and the Young’s modulus $E=0.75\cdot {10}^6$. A pressure-wave

$$P\left(t\right)=P_{\mathit{max}}\left(1-cos\left(2\mathit{t\pi }/T^{\ast }\right)\right)/2\text{with}{P}_{\mathit{max}}=2\cdot {10}^4,$$E63

is prescribed on the fluid inlet boundary for $T^{\ast }=5\cdot {10}^{-3}$ (seconds). Zero traction is enforced on the fluid outlet boundary and no-slip condition is imposed on the lower boundary $y=0$. For the solid, the two endpoints are fixed with $\mathit{d}=\mathbf{0}$ at $x=0$ and $x=6$.

The first experiment is set up to compare the Robin-Neumann scheme with the $\beta$-scheme, the two stable decoupled methods recently developed for the benchmark model. Figure 5 displays the displacements computed by the Robin-Neumann scheme and the $\beta$-scheme, together with the coupled implicit scheme for reference, where the mesh size and the time step size are $h=0.05$ and $\mathrm{\Delta }t={10}^{-4}$. It is clearly seen that both decoupled schemes converge as well as the coupled implicit scheme. Moreover, little difference is observed between the two decoupled schemes numerically. This suggests we focus on the $\beta$-scheme for investigating the multirate time-step technique.

Figure 5.

Comparisons of the numerical results obtained by the coupled implicit scheme, the RN scheme, and the $\beta$-scheme under the setting: $h=0.05$ and $\mathrm{\Delta }{t}_s={10}^{-4}$.

We then conduct numerical experiments to investigate the effects of the time-step ratio $r$. Figure 6 illustrates that a larger time step size in the fluid part results in a more accurate numerical solution than that obtained by using a larger step size in the structure part. In the test, we fix $h=0.1$ and $\mathrm{\Delta }t={10}^{-5}$. In addition, we observe that, when $\frac{\mathrm{\Delta }{t}_s}{\mathrm{\Delta }{t}_f}$ is further increased to be $\frac{\mathrm{\Delta }{t}_s}{\mathrm{\Delta }{t}_f}=5$ or $\frac{\mathrm{\Delta }{t}_s}{\mathrm{\Delta }{t}_f}=10$, there are substantial numerical instability. This screens out the possibility of using a larger time step size in the structure part.

Figure 6.

Comparison of the $\beta$-scheme and the multirate $\beta$-schemes with two different time-step ratios ($h=0.1$ and $\mathrm{\Delta }t={10}^{-5}$ are fixed).

To examine how the stability and approximation are affected when the time step in the fluid region is too large, we fix the time step $\mathrm{\Delta }{t}_s$ and $h$ while varying the time-step ratio $r=\mathrm{1,5,10,20,50}$. Figure 7 displays the computed displacements at $t=0.015$ with the structure time step size $\mathrm{\Delta }{t}_s={10}^{-5}$, the mesh size $h=0.1$ (left) and $h=0.01$ (right). In the left figure, we observe that the structure displacements computed by using $r=1,2,5,10$ approximate very well to that by using the coupled implicit scheme. To further investigate the stability and the convergence of the multirate $\beta$-scheme, in the right part of Figure 7, we reduce the mesh size to be $h=0.01$ while fixing the time step size. The numerical results confirm that the multirate $\beta$-scheme is still stable even the time-step ratio is reasonably large.

Figure 7.

Numerical displacements under the settings: $h=0.1$ (left) $h=0.01$ (right) and $\mathrm{\Delta }{t}_s={10}^{-5}$.

In Figure 8, we present the numerical results of the fluid pressure distribution at $t=0.005,0.01,0.015$. From the top to the bottom, numerical results are: a reference solution by the coupled implicit scheme, the numerical solution by the $\beta$-scheme, and the solution by the multirate $\beta$-scheme with $r=10$. By comparing the results, we see that the multirate $\beta$-scheme provides a very good approximation.

Figure 8.

Fluid pressure distribution at $t=0.005,0.010,0.015$ obtained by the coupled implicit scheme (top), the multirate $\beta$-scheme with $r=1$ (middle) and $r=10$ (bottom) with $h=0.01$ and $\mathrm{\Delta }{t}_s=0.00001$. (a) t = 0.005, (b) t = 0.010, (c) t = 0.015, (d) t = 0.005, (e) t = 0.010, (f) t = 0.015, (g) t = 0.005, (h) t = 0.010, (i) t = 0.015.

In order to examine the order of convergence, we start with $h=0.1$ and $\mathrm{\Delta }{t}_s=0.0001$, and then refine the mesh size by a factor of $2$ and the time step size by a factor of $4$. The space–time size settings are:

$${\left\{h,\mathrm{\Delta }{t}_s\right\}}^i=\left\{0.1\cdot {\left(0.5\right)}^i\mathrm{,0.0001}\cdot {\left(0.25\right)}^i\right\},i=\mathrm{0,1,2,3,4}.$$E64

We compare the numerical solutions of the multirate $\beta$-scheme with the reference solution. The reference solution is computed by using the coupled implicit scheme with a high space–time grid resolution $\left(h=3.125\times {10}^{-3},\mathrm{\Delta }t={10}^{-6}\right)$ as that in [14]. In the multirate scheme, $r=1$ and $r=10$. The relative errors of the primary variables (${\mathit{u}}_f$, $p_f$ and $\mathit{d}$) at $t=0.015$ are displayed in Figure 9. From the comparisons, we see that the numerical errors are approximately reduced by a factor of $4$ as the mesh size and the time step size are refined once. Therefore, the multirate $\beta$-scheme is of a second order in $h$ and a first order in $t$.

Figure 9.

Relative errors of primary variables with the spacing $h$ and time step size $\mathrm{\Delta }{t}_s$ in (64). (a) Relative error of uf, (b) Relative error of pf, (c) Relative error of d.

Finally, in order to demonstrate the advantage of the multirate $\beta$-scheme, we compare in Table 3 the CPU times of the concerned numerical algorithms under various settings. We fixed $\mathrm{\Delta }{t}_s={10}^{-5}$ and varying the mesh sizes as $h=\frac{1}{10},\frac{1}{20},\frac{1}{40},\frac{1}{80},\frac{1}{160}$. From the table, it is observed that the multirate $\beta$-scheme takes much less computational cost than that of the coupled implicit scheme, particularly when $r$ is large.

	Implicit scheme	Multirate $\beta$-scheme r = 1	Multirate $\beta$-scheme r = 10
$h=\frac{1}{10}$	14.90	4.02	0.74
$h=\frac{1}{20}$	48.64	16.00	2.82
$h=\frac{1}{40}$	179.83	66.67	11.6
$h=\frac{1}{80}$	797.76	297.96	49.23
$h=\frac{1}{160}$	3165.26	1270.30	206.32

Table 3.

CPU times (in seconds) for the coupled implicit scheme and the multirate $\beta$-scheme (with $r=1$ or $10$) under different settings of mesh sizes ($\mathrm{\Delta }{t}_s={10}^{-5}$ is fixed).

4. Concluding remarks

In numerical methods for coupled multidomain PDE models, there are usually two types of methods: the monolithic methods and the partitioned (or decoupled) methods. The monolithic methods usually require a code developed for the coupled problem being solved globally. In contrast, the partitioned approach preserves software modularity because one can use existing subdomain solvers. The advantages of using monolithic methods exist in that they provide better approximation accuracy and usually have better stability than the decoupled methods. The drawback is that the computation costs based on the monolithic approaches are usually high. In comparison, the partitioned approaches allow reusing existing software which is an attractive advantage. However, the accuracy and stability of the partitioned method need to be taken into consideration. Nevertheless, from our research works, we also note that one can apply decoupling techniques in monolithic methods, for example, decoupled preconditioners and schemes which are parallel in time. On the other hand, in partitioned algorithms, one can also try to apply coupling numerical techniques such as extrapolations using previous time-step solutions, interpolation using the coarse-grid solution, or extrapolations between subdomain solutions. In this work, we review the decoupling algorithms for the coupled models of fluid flow interacting with porous media flow and FSI models. We show how to decouple the coupled PDE models using preconditioning, two-level and multi-level algorithms, and partitioned time schemes. We attach importance to the decoupling numerical techniques while also emphasizing how to preserve the coupled multidomain physics features. This review provides a general framework for designing decoupled algorithms for coupled PDE models and exhibits the philosophy of the interplays between PDE models and the corresponding numerical methods.

Funding

Mingchao Cai’s work was supported in part by the NSF Awards (170023, 1831950), NIH-BUILD grant through UL1GM118973, NIH-RCMI grant through U54MD013376, and XSEDE HPC resource TG-MCH200022. Mo Mu’s work was supported in part by the HongKong RGC CERG HKUST16301218.