Decoupling Techniques for Coupled PDE Models in Fluid Dynamics

*Mingchao Cai, Mo Mu and Lian Zhang*

## **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 timestepping 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.

**Keywords:** decoupling, linearization, Stokes/Darcy model, FSI model, finite element, two-level method, Robin-Neumann scheme, β-scheme

## **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–8], we discuss some decoupling techniques [4, 9–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 *β*-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–30]. For the nonlinear case, we refer to the coupled nonlinear Navier–Stokes/Darcy model [4–6].

#### **2.1 Coupled models for fluid flow interacting with porous media flow**

Let Ω ⊂*R<sup>d</sup>* be a domain consisting of a fluid region Ω*<sup>f</sup>* and a porous media region Ω*<sup>p</sup>* separated by an interface Γ, as shown in **Figure 1**, where *d* ¼ 2 or 3, Ω ¼ Ω*<sup>f</sup>* ∪ Ω*<sup>p</sup>* and <sup>Γ</sup> <sup>¼</sup> <sup>Ω</sup>*f*∩Ω*p*. Let *<sup>n</sup><sup>f</sup>* and *<sup>n</sup><sup>p</sup>* denote the unit outward normal directions on *<sup>∂</sup>*Ω*<sup>f</sup>* and *∂*Ω*p*. The interface Γ is assumed to be smooth enough as in [6].

For incompressible Newtonian fluid flow, Navier–Stokes equations of the stressdivergence form are usually used [31, 32]. ∀*t*≥0, ∀*x* ∈ Ω*<sup>f</sup>* ,

**Figure 1.**

*A global domain* Ω *consisting of a fluid region* Ω*<sup>f</sup> and a porous media region* Ω*<sup>p</sup> separated by an interface* Γ*.*

*Decoupling Techniques for Coupled PDE Models in Fluid Dynamics DOI: http://dx.doi.org/10.5772/intechopen.105997*

$$\begin{cases} \rho\_f \left( \frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla) \mathbf{u} \right) - \operatorname{div} \mathbf{T}(\mathbf{u}, p) = \mathbf{f}\_f, \\\ \operatorname{div} \mathbf{u} = \mathbf{0}, \end{cases} \tag{1}$$

where *ρ<sup>f</sup>* is the fluid density, *u* is the velocity vector, *p* is the pressure, *f<sup>f</sup>* is the external force,

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

is the stress tensor with *ν*>0 being the kinematic viscosity. By dropping the term *∂u <sup>∂</sup><sup>t</sup>* in (1), the steady state Navier–Stokes equations read as: ∀*x*∈ Ω*<sup>f</sup>* ,

$$\begin{cases} \rho\_f(\mathfrak{u}\cdot\nabla)\mathfrak{u} - \mathrm{div}\,\boldsymbol{T}(\mathfrak{u},\boldsymbol{p}) = \mathfrak{f}\_f, \\ \mathrm{div}\,\mathfrak{u} = \mathbf{0}. \end{cases} \tag{3}$$

In strong form,

$$-\text{div}\,T(\mathfrak{u},p) = -\nu\Delta\mathfrak{u} + \nabla p,\tag{4}$$

because the fluid flow is assumed to be divergence free and div <sup>∇</sup>*<sup>T</sup><sup>u</sup>* � � <sup>¼</sup> 0 holds. Among various porous media flow models, Darcy's law is the most favored. The governing variable in Ω*<sup>p</sup>* is the so-called *piezometric head* or *pressure head*,

$$
\phi = \mathbf{z} + \frac{p\_p}{\rho\_f \mathbf{g}}.\tag{5}
$$

Here, *z* is the elevation from a reference level (for simplicity, *z* is assumed to be 0). Darcy's law states that the velocity *u<sup>p</sup>* (also called seepage velocity) in the porous media region is proportional to the gradient of *ϕ* [27, 33].

$$\mathfrak{u}\_p = -\mathbf{K} \nabla \phi. \tag{6}$$

We assume that

$$a\_1(\mathfrak{x}, \mathfrak{x}) \le (\mathbf{K}\mathfrak{x}, \mathfrak{x}) \le a\_2(\mathfrak{x}, \mathfrak{x}), \quad \forall \mathfrak{x} \in \Omega\_{\mathfrak{p}}.\tag{7}$$

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

$$-\text{div}\,(\mathbf{K}\nabla\phi) = f\_p.\tag{8}$$

In the time-dependent case, the governing equations in Ω*<sup>p</sup>* reads as:

$$\mathbf{S}\_0 \phi\_t - \text{div } (\mathbf{K} \nabla \phi) = f\_p,\tag{9}$$

where *S*<sup>0</sup> is a specific storage and *f <sup>p</sup>* 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–3, 7, 27, 28]:

$$\begin{cases} \mathbf{u} \cdot \mathbf{n}\_f = \mathbf{u}\_p \cdot \mathbf{n}\_f = -\mathbf{K} \nabla \phi \cdot \mathbf{n}\_f, \\\\ -\nu \left( \nabla \mathbf{u} \mathbf{n}\_f \right) \cdot \mathbf{n}\_f + p = \rho \mathbf{g} \phi, \\\\ -\nu \left( \nabla \mathbf{u} \mathbf{n}\_f \right) \cdot \boldsymbol{\tau}\_i = \frac{\nu \alpha\_{B \mid \mathbf{S}}}{\sqrt{\nu \boldsymbol{\tau}\_i \cdot \mathbf{K} \boldsymbol{\tau}\_i}} \mathbf{u} \cdot \boldsymbol{\tau}\_i, \quad i = \mathbf{1}, \ldots, d - \mathbf{1}. \end{cases} \tag{10}$$

Here, f g *τ<sup>i</sup> d*�1 *<sup>i</sup>*¼<sup>1</sup> is the unit tangent vector on <sup>Γ</sup>, *<sup>α</sup>BJS* is a positive parameter depending on the properties of the porous medium. The first interface condition ensures mass conservation across Γ. 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 Γ. Without loss of generality, we impose homogeneous Dirichlet boundary conditions on both of the external boundaries:

$$\begin{cases} \mathfrak{u} = \mathbf{0} & \text{on} \quad \Gamma\_f, \\\\ \phi = \mathbf{0} & \text{on} \quad \Gamma\_p. \end{cases} \tag{11}$$

The proper functional spaces for *u*, *p* and *ϕ* are

$$\mathbf{X}\_f = \left\{ \mathbf{v} \in H^1(\mathfrak{Q}\_f) = \left( H^1(\mathfrak{Q}\_f) \right)^d \, \middle| \, \mathbf{v} = \mathbf{0} \text{ on } \Gamma\_f \right\},$$

$$\mathbf{Q} = L^2(\mathfrak{Q}\_f), \quad X\_p = \left\{ \boldsymbol{\psi} \in H^1(\mathfrak{Q}\_p) \, \middle| \, \boldsymbol{\psi} = \mathbf{0} \text{ on } \Gamma\_p \right\}.\tag{12}$$

Moreover, we denote *X* ¼ *X<sup>f</sup>* � *Xp* 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 *u* ¼ ð Þ *u*, *ϕ* ∈*X*, *p* ∈ *Q* such that

$$\begin{cases} a(\underline{u}, \underline{v}) + c(\mathfrak{u}, \mathfrak{u}, \mathfrak{v}) + b(\mathfrak{v}, p) = f(\underline{v}) & \forall \underline{v} = (\mathfrak{v}, \psi) \in \underline{X}, \\ b(\mathfrak{u}, q) = \mathbf{0} & \forall q \in Q, \end{cases} \tag{13}$$

where

$$\begin{split} a(\underline{u}, \underline{v}) &= a\_f(\boldsymbol{u}, \boldsymbol{\nu}) + a\_p(\boldsymbol{\phi}, \boldsymbol{\nu}) + a\_\Gamma(\underline{u}, \underline{v}), \quad b(\boldsymbol{\nu}, p) = -\int\_{\Omega\_f} p \nabla \cdot \boldsymbol{\nu}, \\ c(\boldsymbol{u}, \boldsymbol{\nu}, \boldsymbol{w}) &= \rho \int\_{\Omega\_f} (\boldsymbol{u} \cdot \nabla) \boldsymbol{\nu} \cdot \boldsymbol{w}, \quad f(\underline{v}) = \int\_{\Omega\_f} \boldsymbol{f}\_f \cdot \boldsymbol{\nu} + \rho \mathbf{g} \Big|\_{\Omega\_p} \boldsymbol{f}\_p \boldsymbol{\nu} \end{split} \tag{14}$$

with

$$\begin{split} a\_{\boldsymbol{f}}(\boldsymbol{\upmu},\boldsymbol{\upnu}) &= \nu \int\_{\Omega\_{\boldsymbol{f}}} \nabla \boldsymbol{\upmu} : \nabla \boldsymbol{\upnu} + \sum\_{i=1}^{d-1} \frac{\nu \alpha\_{\boldsymbol{\upbeta} \mid \boldsymbol{\upbeta}}}{\sqrt{\nu \boldsymbol{\uppi}\_{i} \cdot \mathbf{K} \boldsymbol{\uppi}\_{i}}} \int\_{\Gamma} (\boldsymbol{\upmu} \cdot \boldsymbol{\uppi}\_{i}) (\boldsymbol{\upnu} \cdot \boldsymbol{\uppi}\_{i}), \\ a\_{\boldsymbol{\upmu}}(\boldsymbol{\upmu},\boldsymbol{\upmu}) &= \rho \mathbf{g} \int\_{\Omega\_{\boldsymbol{\upmu}}} \nabla \boldsymbol{\upmu} \cdot \mathbf{K} \nabla \boldsymbol{\upphi}, \quad a\_{\boldsymbol{\upalpha}}(\underline{\boldsymbol{\upmu}},\underline{\boldsymbol{\upnu}}) = \rho \mathbf{g} \int\_{\Gamma} (\boldsymbol{\upphi} \boldsymbol{\upnu} - \boldsymbol{\upmu} \boldsymbol{\upmu}) \cdot \boldsymbol{\uppi}\_{\boldsymbol{\upmu}}. \end{split} \tag{15}$$

**148**

*Decoupling Techniques for Coupled PDE Models in Fluid Dynamics DOI: http://dx.doi.org/10.5772/intechopen.105997*

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*ð Þ �, � is bounded and coercive; *b*ð Þ �, � is bounded and satisfies the inf-sup condition [27, 35]; and the nonlinear term can be bounded by using the *H*<sup>1</sup> norm of the three components [31, 36].

We partition Ω*<sup>f</sup>* and Ω*<sup>p</sup>* by quasi-uniform triangulations T *<sup>f</sup>*,*<sup>h</sup>* and T *<sup>p</sup>*,*<sup>h</sup>* with a characteristic meshsize *h*. Here, we require that the two subdomain triangulations coincide at Γ. The corresponding finite element spaces are denoted by *X<sup>f</sup>*,*<sup>h</sup>* � *Qh* ⊂ *X<sup>f</sup>* � *Q* and *Xp*,*<sup>h</sup>* ⊂*Xp*, respectively. Moreover, *X<sup>f</sup>*,*<sup>h</sup>* � *Qh* needs to be stable, i.e., there exists a positive constant *β* such that

$$\sup\_{\boldsymbol{\sigma}\_h \in \underline{\mathbf{X}}\_{f,h}} \frac{b\left(\boldsymbol{\sigma}\_h, q\_h\right)}{\left|\boldsymbol{\sigma}\_h\right|\_{1, \Omega\_f}} \ge \beta \|q\_h\|\_{0, \Omega\_f} \qquad \forall q\_h \in \mathcal{Q}\_h. \tag{16}$$

We assume that the solution of (13) is smooth enough and the finite element spaces have the following typical approximation properties: for all ð Þ *<sup>u</sup>*, *<sup>p</sup>* <sup>∈</sup> *<sup>H</sup><sup>k</sup>*þ<sup>1</sup> <sup>Ω</sup>*<sup>f</sup>* � �∩*X<sup>f</sup>* � *<sup>H</sup><sup>k</sup>* <sup>Ω</sup>*<sup>f</sup>* � � and *ϕ*∈ *H<sup>k</sup>*þ<sup>1</sup> Ω*<sup>p</sup>* � �∩*Xp*,

$$\inf\_{\boldsymbol{\sigma}\_{k}\in\mathcal{K}\_{f,k},q\_{k}\in\mathcal{Q}\_{k}}\left\{h|\boldsymbol{\mu}-\boldsymbol{\sigma}\_{h}|\_{1,\Omega\_{f}}+\|\boldsymbol{\mu}-\boldsymbol{\sigma}\_{h}\|\_{0,\Omega\_{\mathcal{I}}}+h\|\boldsymbol{p}-q\_{h}\|\_{0,\Omega\_{\mathcal{I}}}\right\}\lesssim h^{k+1}\left(|\boldsymbol{\mu}|\_{k+1,\Omega\_{\mathcal{I}}}+|\boldsymbol{p}|\_{k,\Omega\_{\mathcal{I}}}\right),\tag{17}$$

$$\inf\_{\boldsymbol{\mu}\_{\mathcal{Y}h}\in X\_{p,k}} \left\{ h|\phi - \boldsymbol{\mu}\_{h}|\_{\mathbf{1},\boldsymbol{\Omega}\_{\mathcal{Y}}} + \|\phi - \boldsymbol{\mu}\_{h}\|\_{0,\boldsymbol{\Omega}\_{\mathcal{Y}}} \right\} \lesssim h^{k+1}|\phi|\_{k+1,\boldsymbol{\Omega}\_{\mathcal{Y}}}.\tag{18}$$

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 Ω*<sup>f</sup>* and the piecewise linear elements in Ω*p*; if *k*≥2, the *k*-th order Taylor-Hood elements [31, 37, 38] can be applied in Ω*<sup>f</sup>* and the piecewise *k*-th order elements can be adopted in Ω*p*.

**Coupled Algorithm**: A conventional finite element discretization applied to the model problem (13) leads to the discrete problem: Find *uh* ¼ *uh*, *ϕ<sup>h</sup>* ð Þ∈*Xh* ¼ *X<sup>f</sup>*,*<sup>h</sup>* � *Xp*,*<sup>h</sup>*, *ph* ∈ *Qh* such that

$$\begin{cases} a(\underline{u}\_h, \underline{v}\_h) + c(\underline{u}\_h, \underline{u}\_h, \underline{v}\_h) + b(\underline{v}\_h, \underline{p}\_h) = f(\underline{v}\_h) \quad \forall \underline{p}\_h = (\underline{v}\_h, \underline{\mu}\_h) \in \underline{X}\_h, \\ b(\underline{u}\_h, q\_h) = 0, \quad \forall q\_h \in Q\_h. \end{cases} \tag{19}$$

#### **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

$$
\begin{bmatrix} A\_p & A\_\Gamma^T & \mathbf{0} \\ -A\_\Gamma & A\_f & B\_f^T \\ \mathbf{0} & B\_f & \mathbf{0} \end{bmatrix} \begin{bmatrix} \phi\_h \\ \mathbf{u}\_h \\ p\_h \end{bmatrix} = \begin{bmatrix} f\_{f,h} \\ f\_{p,h} \\ g\_h \end{bmatrix} . \tag{20}
$$

Here *Af* , *Ap*, *A*Γ, and *B* ¼ **0**, *Bf* � � are the corresponding linear operators induced by the corresponding bilinear forms in (15) and (16). We denote

$$A = \begin{bmatrix} A\_p & A\_\Gamma^T \\ -A\_\Gamma & A\_f \end{bmatrix}, \quad \text{and} \quad M = \begin{bmatrix} \mathbf{A} & \mathbf{B}^T \\ \mathbf{B} & \mathbf{0} \end{bmatrix}. \tag{21}$$

By discarding the coupling interface terms and plugging in <sup>1</sup> *<sup>ν</sup> I* (which leads to pressure mass matrix) in the (2, 2) block, we have the following block-diagonal decoupled preconditioner:

$$P\_M = \begin{bmatrix} A\_0 & \mathbf{0} \\ \mathbf{0} & \mathbf{1}\_f \\ \mathbf{0} & \nu \end{bmatrix}, \quad \text{with} \quad A\_0 = \begin{bmatrix} A\_p & \mathbf{0} \\ \mathbf{0} & A\_f \end{bmatrix}. \tag{22}$$

By keeping one of *B* and *BT*, 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*�<sup>1</sup> <sup>0</sup> and the inverse for <sup>1</sup> *<sup>ν</sup> 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} = \begin{pmatrix} A\_0 & 0\\ 0 & \pm \frac{1}{\nu} I\_h \end{pmatrix},\tag{23}$$

We also propose another preconditioner of the block triangular type:

$$P\_{T\_1} = \begin{pmatrix} A\_0 & 0 \\ B & -\frac{1}{\nu} I\_h \end{pmatrix},\tag{24}$$

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**.

We use ∗ ∗ to indicate that the iteration does not converge within the prescribed maximum number of iterations. *N P*ð Þ 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 *PT*<sup>1</sup> 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 *Decoupling Techniques for Coupled PDE Models in Fluid Dynamics DOI: http://dx.doi.org/10.5772/intechopen.105997*


**Table 1.**

*Number of iterations for the GMRES method without and with the two preconditioners P*� *and PT*<sup>1</sup> *.*

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 2ð Þ *:*7 *<sup>H</sup>* on a coarse grid: find *uH* ¼ ð Þ *uH*, *ϕ<sup>H</sup>* ∈*XH* ⊂*Xh*, *pH* ∈ *QH* ⊂ *Qh* such that

$$\begin{cases} a(\underline{u}\_H, \underline{v}\_H) + b(\underline{v}\_H, p\_H) = (f, \underline{v}\_H), \quad \forall \underline{v}\_H = (\underline{v}\_H, \underline{\nu}\_H) \in \underline{X}\_H, \\ b(\underline{u}\_H, q\_H) = 0, \quad \forall q\_H \in Q\_H; \end{cases} \tag{25}$$

2. Solve a modified fine grid problem: find *uH* <sup>¼</sup> *<sup>u</sup><sup>h</sup>*, *<sup>ϕ</sup><sup>h</sup>* � �∈*Xh*, *ph* <sup>∈</sup> *Qh* such that

$$\begin{cases} a(\underline{\boldsymbol{u}}^{H}, \boldsymbol{v}\_{h}) + b(\boldsymbol{v}\_{h}, \boldsymbol{p}^{h}) = (\boldsymbol{f}, \boldsymbol{v}\_{h}) - a\_{\Gamma}(\underline{\boldsymbol{u}}\_{H}, \boldsymbol{v}\_{h}), \quad \forall \boldsymbol{v}\_{h} \in \underline{\mathbf{X}}\_{h}, \\ b(\underline{\boldsymbol{u}}^{H}, \boldsymbol{q}\_{h}) = \mathbf{0}, \quad \forall \boldsymbol{q}\_{h} \in \mathbf{Q}\_{h}. \end{cases} \tag{26}$$

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

**Theorem 1.** *Let uh*, *ph* � � *be the solution of coupled Stokes/Darcy model, and uH*, *ph* � � *be defined by and (27) on the fine grid. The following error estimates hold:*

$$\left. \right| \left| \phi\_h - \phi^h \right| \big|\_{H\_p} \lesssim H^2,\tag{27}$$

$$\left||\mathfrak{u}\_{h} - \mathfrak{u}^{h}\right||\_{H\_{\tilde{f}}} \lesssim H^{3/2},\tag{28}$$

$$\left. \right| \left| p\_h - p^h \right| \right|\_{Q} \lesssim H^{3/2}. \tag{29}$$

Based on this algorithm, some other improvements have been made. For example, in [18–20], by sequentially solving the Stokes submodel and the Darcy submodel, the authors can make *<sup>ϕ</sup><sup>h</sup>* � *<sup>ϕ</sup><sup>h</sup>* � � � � *Hp* , *<sup>u</sup><sup>h</sup>* � *<sup>u</sup><sup>h</sup>* � � � � *Hf* , and *ph* � *ph* � � � � *Q* are all of order *H*<sup>2</sup> . Furthermore, Hou constructed a new auxiliary problem [16] for the Darcy submodel, and proved that Mu and Xu's two-grid algorithm can retain *<sup>u</sup><sup>h</sup>* � *<sup>u</sup><sup>h</sup>* � � � � *Hf* and *ph* � *ph* � � � � *<sup>Q</sup>* order of *<sup>H</sup>*<sup>2</sup> . It is remarkable that Mu and Xu's two-grid algorithm is naturally parallel and of optimal order, if *h* is of order *H*<sup>2</sup> .

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 2ð Þ *:*7 *<sup>H</sup>* on a coarse grid: find *uH* ¼ ð Þ *uH*, *ϕ<sup>H</sup>* ∈*XH* ⊂*Xh*, *pH* ∈ *QH* ⊂ *Qh* such that

$$\begin{cases} a(\underline{\boldsymbol{u}}\_{H}, \underline{\boldsymbol{v}}\_{H}) + b(\underline{\boldsymbol{v}}\_{H}, p\_{H}) = f(\underline{\boldsymbol{v}}\_{H}), \forall \underline{\boldsymbol{v}}\_{H} = (\boldsymbol{v}\_{H}, \boldsymbol{\mu}\_{H}) \in \underline{\boldsymbol{X}}\_{H}, \\\ b(\underline{\boldsymbol{u}}\_{H}, q\_{H}) = \boldsymbol{0}, \forall q\_{H} \in \boldsymbol{Q}\_{H}; \end{cases} \tag{30}$$

2. Set *h*<sup>0</sup> ¼ *H*, for j = 1 to *L*,

find *<sup>u</sup>hj* <sup>¼</sup> *<sup>u</sup>hj* , *ϕhj* � �∈*Xhj* , *phj* ∈ *Qhj* such that

$$\begin{cases} a\left(u^{h\_i}, v\_{h\_j}\right) + b\left(v\_{h\_j}, p^{h\_j}\right) = f\left(v\_{h\_j}\right) - a\_\Gamma \left(u^{h\_{j-1}}, v\_{h\_j}\right), \quad \forall v\_{h\_j} \in \underline{\mathbf{X}}\_{h\_j}, \\\ b\left(u^{h\_j}, q\_{h\_j}\right) = \mathbf{0}, \quad \forall q\_{h\_j} \in \mathbf{Q}\_{h\_j}. \end{cases} \tag{31}$$

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**.

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 *uH* ¼ ð Þ *uH*, *ϕ<sup>H</sup>* ∈*XH* and *pH* ∈ *QH* such that

$$\begin{cases} a(\underline{\boldsymbol{u}}\_{H}, \underline{\boldsymbol{v}}\_{H}) + c(\boldsymbol{u}\_{H}, \boldsymbol{u}\_{H}, \boldsymbol{\nu}\_{H}) + b(\boldsymbol{\nu}\_{H}, p\_{H}) = f(\underline{\boldsymbol{v}}\_{H}), \quad \forall \underline{\boldsymbol{v}}\_{H} = (\boldsymbol{\nu}\_{H}, \boldsymbol{\mu}\_{H}) \in \underline{\mathbf{X}}\_{H},\\ b\left(\boldsymbol{u}\_{H}, q\_{H}\right) = \mathbf{0}, \quad \forall q\_{H} \in \mathbf{Q}\_{H}.\end{cases} \tag{32}$$

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


**Table 2.**

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

**Step a**. Solve a discrete Darcy problem in Ω*p*: Find *ϕ*<sup>∗</sup> *<sup>h</sup>* ∈ *Xp*,*<sup>h</sup>* such that

$$\mathfrak{a}\_p(\phi\_h^\*, \psi\_h) = \rho \mathfrak{g}\left(f\_p, \psi\_h\right)\_{\mathfrak{Q}\_p} + \rho \mathfrak{g}\left(\mathfrak{u}\_H \cdot \mathfrak{u}\_f, \psi\_h\right)\_\Gamma \quad \forall \psi\_h \in \mathcal{X}\_{p,h}.\tag{33}$$

**Step b**. Solve a modified Navier–Stokes model using the Newton type linearization: Find *u*<sup>∗</sup> *<sup>h</sup>* ∈ *X<sup>f</sup>*,*<sup>h</sup>* and *p*<sup>∗</sup> *<sup>h</sup>* ∈ *Qh* such that

$$\begin{cases} \begin{aligned} a\_f \left( \mathfrak{u}\_h^\*, \mathfrak{v}\_h \right) + c \left( \mathfrak{u}\_H, \mathfrak{u}\_h^\*, \mathfrak{v}\_h \right) + c \left( \mathfrak{u}\_h^\*, \mathfrak{u}\_H, \mathfrak{v}\_h \right) + b \left( \mathfrak{v}\_h, p\_h^\* \right) = \left( \mathfrak{f}\_f, \mathfrak{v}\_h \right)\_{\Omega\_f} \\ -\rho \lg \left( \phi\_h^\*, \mathfrak{v}\_h \cdot \mathfrak{u}\_f \right)\_\Gamma + c \left( \mathfrak{u}\_H, \mathfrak{u}\_H, \mathfrak{v}\_h \right) \quad \forall \mathfrak{v}\_h \in \mathbf{X}\_{f,h}, \\ b \left( \mathfrak{u}\_h^\*, q\_h \right) = 0 \quad \forall q\_h \in \mathcal{Q}\_h. \end{aligned} \tag{34}$$


$$\mathfrak{a}\_p(\phi^h, \psi\_h) = \rho \mathfrak{g}\left(f\_p, \psi\_h\right)\_{\Omega\_p} + \rho \mathfrak{g}\left(\mathfrak{u}\_h^\* \cdot \mathfrak{n}\_f, \psi\_h\right)\_\Gamma \qquad \forall \psi\_h \in X\_{p,h}.\tag{35}$$

**Step b**. Correct the solution of the fluid flow model: Find *u<sup>h</sup>* ∈ *X<sup>f</sup>*,*<sup>h</sup>* and *p<sup>h</sup>* ∈ *Qh* such that

$$\begin{cases} a\_f(\boldsymbol{\mu}^h, \boldsymbol{\nu}\_h) + c(\boldsymbol{\mu}\_H, \boldsymbol{\mu}^h, \boldsymbol{\nu}\_h) + c(\boldsymbol{\mu}^h, \boldsymbol{\mu}\_H, \boldsymbol{\nu}\_h) + b(\boldsymbol{\nu}\_h, p^h) = \left(\boldsymbol{\mathcal{f}}\_f, \boldsymbol{\nu}\_h\right)\_{\Omega\_f} \\ \quad - \rho \boldsymbol{g}\left(\boldsymbol{\phi}^h, \boldsymbol{\nu}\_h \cdot \boldsymbol{\nu}\_f\right)\_{\Gamma} + c\left(\boldsymbol{\mu}\_H, \boldsymbol{\mu}^\*\_h, \boldsymbol{\nu}\_h\right) + c\left(\boldsymbol{\mu}^\*\_h, \boldsymbol{\mu}\_H - \boldsymbol{\mu}^\*\_h, \boldsymbol{\nu}\_h\right) \quad \forall \boldsymbol{\nu}\_h \in \mathbf{X}\_{f,h} \\ b\left(\boldsymbol{\mu}^h, q\_h\right) = 0 \quad \forall q\_h \in \mathbf{Q}\_h. \end{cases} \tag{36}$$

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*<sup>2</sup> error analysis of the problem (19). In summary, we have.

**Lemma 1.** *Let* ð Þ *<sup>u</sup>*, *<sup>ϕ</sup>*, *<sup>p</sup>* <sup>∈</sup> *<sup>H</sup><sup>k</sup>*þ<sup>1</sup> <sup>Ω</sup>*<sup>f</sup>* � � � *<sup>H</sup><sup>k</sup>*þ<sup>1</sup> <sup>Ω</sup>*<sup>p</sup>* � � � *<sup>H</sup><sup>k</sup>* <sup>Ω</sup>*<sup>f</sup>* � � *be the solution of the Navier–Stokes/Darcy model (13) and uh*, *ϕh*, *ph* � � *be the FE solution of (19). We assume that ν is sufficiently large and h is sufficiently small. There holds the following energy norm estimate for the problem (19)*.

$$\|\mathfrak{u} - \mathfrak{u}\_h\|\_{1,\Omega\_\sharp} + |\phi - \phi\_h|\_{1,\Omega\_\mathbb{P}} + \|\mathfrak{p} - \mathfrak{p}\_h\|\_{0,\Omega\_\sharp} \lesssim h^k. \tag{37}$$

Moreover, we have the following *L*<sup>2</sup> error estimate:

$$\|\mathfrak{u} - \mathfrak{u}\_h\|\_{0,\Omega\_\sharp} + \|\phi - \phi\_h\|\_{0,\Omega\_p} \lesssim h^{k+1}.\tag{38}$$

The energy norm estimate (37) is the Lemma 2 in [4]. The *L*<sup>2</sup> 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 *ϕ*<sup>∗</sup> *<sup>h</sup>* , *u*<sup>∗</sup> *<sup>h</sup>* , *p*<sup>∗</sup> *h* � � in the energy norm.

**Lemma 2.** *(***Error analysis of the intermediate step two-level solution***) Let* ð Þ *<sup>ϕ</sup>*, *<sup>u</sup>*, *<sup>p</sup> and <sup>ϕ</sup>*<sup>∗</sup> *<sup>h</sup>* , *u*<sup>∗</sup> *<sup>h</sup>* , *p*<sup>∗</sup> *h* � � *be defined by the problems (13) and (33)–(34), respectively. Under the assumptions of Lemma 1, there holds*

$$\left|\phi - \phi\_h^\*\right|\_{1,\Omega\_\mathbb{P}} \lesssim H^{k+1} + h^k,\tag{39}$$

$$\left|\mathfrak{u} - \mathfrak{u}\_h^\*\right|\_{1,\Omega\_\ell} + \|p - p\_h^\*\|\_{0,\Omega\_\ell} \lesssim H^{k+1} + h^k. \tag{40}$$

$$\|\|\boldsymbol{u} - \boldsymbol{u}\_h^\*\|\|\_{0,\Omega\_\sharp} \lesssim H^{2k+1} + H^{k+1}h + h^{k+1}.\tag{41}$$

From Lemma 1, we note that the optimal finite element solution error in the energy norm is of order *O h<sup>k</sup>* � �. 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 *<sup>h</sup>* <sup>¼</sup> *<sup>H</sup><sup>k</sup>*þ<sup>1</sup> *<sup>k</sup>* . The *L*<sup>2</sup> error analysis is extended from [4, 20, 31, 37]. Let **u** and *ϕ* be the nonsingular solution of (13). From Lemma 1, the optimal *<sup>L</sup>*<sup>2</sup> error for the finite element solution is of order *O h<sup>k</sup>*þ<sup>1</sup> � �. To make sure <sup>∥</sup>*<sup>u</sup>* � *<sup>u</sup>*<sup>∗</sup> *<sup>h</sup>* <sup>∥</sup>0,Ω*<sup>f</sup>* is also of order *O h<sup>k</sup>*þ<sup>1</sup> � �, the scaling between the two grids has to be taken as *<sup>h</sup>* <sup>¼</sup> max *<sup>H</sup><sup>k</sup>*þ<sup>1</sup> *<sup>k</sup>* , *H*<sup>2</sup>*k*þ<sup>1</sup> *k*þ1 n o. For instance, if *<sup>k</sup>* <sup>¼</sup> 1, we have to set *<sup>h</sup>* <sup>¼</sup> *<sup>H</sup>*<sup>3</sup> <sup>2</sup> to make sure the *L*<sup>2</sup> error of *u*<sup>∗</sup> *<sup>h</sup>* 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* ð Þ *ϕ*, *u*, *p and ϕ<sup>h</sup>*, *u<sup>h</sup>*, *ph* � � *be the solutions of (13) and (35)–(36) respectively. Under the assumptions of Lemma 1, the following error estimates hold:*

$$\left|\phi - \phi^{h}\right|\_{1,\Omega\_{p}} + \left|\mathfrak{u} - \mathfrak{u}^{h}\right|\_{1,\Omega\_{\tilde{\mathcal{I}}}} + \left\|p - p^{h}\right\|\_{0,\Omega\_{\tilde{\mathcal{I}}}} \lesssim H^{2k+1} + H^{k+1}h + h^{k},\tag{42}$$

**Proposition 1.** *Let* ð Þ *<sup>ϕ</sup>*, *<sup>u</sup>*, *<sup>p</sup> and <sup>ϕ</sup><sup>h</sup>*, *<sup>u</sup><sup>h</sup>*, *ph* � � *be the solutions of (13) and (35)–(36) respectively. If we take h* <sup>¼</sup> *<sup>H</sup>*<sup>2</sup>*k*þ<sup>1</sup> *<sup>k</sup> for k* <sup>¼</sup> 1, 2 *and h* <sup>¼</sup> *<sup>H</sup><sup>k</sup>*þ<sup>1</sup> *<sup>k</sup>*�<sup>1</sup> *for k*≥3*, then there holds the following error estimate.*

$$\left|\phi - \phi^{h}\right|\_{1,\Omega\_{p}} + \left|\mathfrak{u} - \mathfrak{u}^{h}\right|\_{1,\Omega\_{\xi}} + \left\|p - p^{h}\right\|\_{0,\Omega\_{\xi}} \lesssim h^{k}.\tag{43}$$

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 *Decoupling Techniques for Coupled PDE Models in Fluid Dynamics DOI: http://dx.doi.org/10.5772/intechopen.105997*

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{split} & \mathfrak{u} \in \mathcal{W}^{2,\infty}(0,\infty; \boldsymbol{L}^{2}(\boldsymbol{\Omega}\_{\boldsymbol{f}})) \cap \mathcal{W}^{1,\infty}(0,\infty; \boldsymbol{H}^{2}(\boldsymbol{\Omega}\_{\boldsymbol{f}})), \\ & \phi \in \mathcal{W}^{2,\infty}(0,\infty; \boldsymbol{L}^{2}(\boldsymbol{\Omega}\_{\boldsymbol{p}})) \cap \mathcal{W}^{1,\infty}(0,\infty; \boldsymbol{H}^{2}(\boldsymbol{\Omega}\_{\boldsymbol{p}})), \\ & \boldsymbol{p} \in \boldsymbol{L}^{2,\infty}(0,\infty; \boldsymbol{H}^{1}(\boldsymbol{\Omega}\_{\boldsymbol{f}})). \end{split} \tag{44}$$

The simplest time scheme for the evolutionary coupled Stokes/Darcy model is the **Backward Euler Algorithm**, which reads as: Given *u<sup>n</sup> <sup>h</sup>*, *pn <sup>h</sup>*, *ϕ<sup>n</sup> h* � �<sup>∈</sup> *<sup>X</sup><sup>f</sup>*,*<sup>h</sup>* � *Qh* � *Xp*,*<sup>h</sup>*, find *u<sup>n</sup>*þ<sup>1</sup> *<sup>h</sup>* , *p<sup>n</sup>*þ<sup>1</sup> *<sup>h</sup>* , *ϕ<sup>n</sup>*þ<sup>1</sup> *h* � �<sup>∈</sup> *<sup>X</sup><sup>f</sup>*,*<sup>h</sup>* � *Qh* � *Xp*,*<sup>h</sup>*, such that for all *<sup>v</sup><sup>h</sup>* <sup>∈</sup> *<sup>X</sup><sup>f</sup>*,*<sup>h</sup>*, *qh* <sup>∈</sup> *Qh*, *<sup>ψ</sup><sup>h</sup>* <sup>∈</sup>*Xp*,*<sup>h</sup>*,

$$\begin{cases} \left(\frac{\mathbf{u}\_h^{n+1} - \mathbf{u}\_h^n}{\Delta t}, \boldsymbol{\nu}\_h\right) + a\_f(\boldsymbol{u}\_h^{n+1}, \boldsymbol{\nu}\_h) - \left(p\_h^{n+1}, \nabla \cdot \boldsymbol{\nu}\_h\right) + \mathbf{g}\left(\boldsymbol{\nu}\_h \cdot \boldsymbol{\nu}\_f, \phi\_h^{n+1}\right)\_\Gamma = \left(\boldsymbol{f}\_f^{n+1}, \boldsymbol{\nu}\_h\right),\\ \left(q\_h, \nabla \cdot \boldsymbol{\mu}\_h^{n+1}\right) = 0,\end{cases} \tag{45}$$

$$\mathrm{gS}\_{0}\left(\frac{\phi\_{h}^{n+1} - \phi\_{h}^{n-1}}{\Delta t}, \boldsymbol{\nu}\_{h}\right) + a\_{p}\left(\phi\_{h}^{n+1}, \boldsymbol{\nu}\_{h}\right) - \mathrm{g}\left(\mathfrak{u}\_{h}^{n+1} \cdot \boldsymbol{\nu}\_{f}, \boldsymbol{\nu}\_{h}\right)\_{\Gamma} = \mathrm{g}\left(\boldsymbol{f}\_{p}^{n+1}, \boldsymbol{\nu}\_{h}\right). \tag{46}$$

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 *un <sup>h</sup>*, *p<sup>n</sup> <sup>h</sup>*, *ϕ<sup>n</sup> h* � �<sup>∈</sup> *<sup>X</sup><sup>f</sup>*,*<sup>h</sup>* � *Qh* � *Xp*,*<sup>h</sup>*, find *<sup>u</sup><sup>n</sup>*þ<sup>1</sup> *<sup>h</sup>* , *pn*þ<sup>1</sup> *<sup>h</sup>* , *ϕ<sup>n</sup>*þ<sup>1</sup> *h* � �<sup>∈</sup> *<sup>X</sup><sup>f</sup>*,*<sup>h</sup>* � *Qh* � *Xp*,*<sup>h</sup>*, such that for all *v<sup>h</sup>* ∈ *X<sup>f</sup>*,*<sup>h</sup>*, *qh* ∈ *Qh*, *ψ<sup>h</sup>* ∈ *Xp*,*<sup>h</sup>*,

$$\begin{cases} \left(\frac{\mathbf{u}\_h^{n+1} - \mathbf{u}\_h^n}{\Delta t}, \boldsymbol{\nu}\_h\right) + a\_f \left(\mathbf{u}\_h^{n+1}, \boldsymbol{\nu}\_h\right) - \left(p\_h^{n+1}, \nabla \cdot \boldsymbol{\nu}\_h\right) + \mathbf{g}\left(\boldsymbol{\nu}\_h \cdot \boldsymbol{\nu}\_f, \boldsymbol{\phi}\_h^n\right)\_\Gamma = \left(\boldsymbol{f}\_f^{n+1}, \boldsymbol{\nu}\_h\right), \\\left(q\_h, \nabla \cdot \mathbf{u}\_h^{n+1}\right) = 0, \\\ \mathbf{g}\boldsymbol{\mathcal{S}}\_0 \left(\frac{\boldsymbol{\phi}\_h^{n+1} - \boldsymbol{\phi}\_h^n}{\Delta t}, \boldsymbol{\nu}\_h\right) + a\_p \left(\boldsymbol{\phi}\_h^{n+1}, \boldsymbol{\nu}\_h\right) - \mathbf{g}\left(\boldsymbol{\mu}\_h^n \cdot \boldsymbol{\nu}\_f, \boldsymbol{\nu}\_h\right)\_\Gamma = \mathbf{g}\left(\boldsymbol{f}\_p^{n+1}, \boldsymbol{\nu}\_h\right). \end{cases} \tag{48}$$

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*

$$
\Delta t \lesssim \min \left\{ \nu k\_{\min}^2, \mathcal{S}\_0 \nu^2 k\_{\min} \right\}, \tag{49}
$$

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

$$\text{i.i.} \,\text{If} \, f\_f \in L^\infty(\mathbf{0}, +\infty; L^2(\mathfrak{Q}\_f)), f\_p \in L^\infty(\mathbf{0}, +\infty; L^2(\mathfrak{Q}\_p)), \text{ then}$$

$$\left\| \left| \mathfrak{u}\_h^n \right| \right\|^2 + \left\| \left| \mathfrak{d}\_h^n \right| \right\|^2 \le \mathbf{C}, \quad \forall n \ge \mathbf{0}. \tag{50}$$

$$\text{iii. If } \left\| f\_f \right\|\_{L^{\infty}\left(0, +\infty; L^2\left(\mathfrak{d}\_f\right)\right)}, \left\| f\_p \right\|\_{L^{\infty}\left(0, +\infty; L^2\left(\mathfrak{d}\_p\right)\right)} \text{ are uniformly bounded in } \Delta t \text{, then}$$

$$\left(\left\|\boldsymbol{\boldsymbol{u}}\_{h}^{n}\right\|\right)^{2} + \left\|\boldsymbol{\phi}\_{h}^{n}\right\|^{2} + \Delta t \sum\_{l=0}^{n} \left(\left\|\nabla \boldsymbol{u}\_{h}^{l}\right\|\right^{2} + \left\|\nabla \boldsymbol{\phi}\_{h}^{l}\right\|^{2}\right) \leq \mathsf{C}, \quad \forall n \geq 0. \tag{51}$$

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*<sup>0</sup> and *k min* are small. Another way for uncoupling surface/ subsurface flow models is using splitting schemes which require sequential subproblem solves at each time step [46]. As an example, we note that in solving (49), one can replace *u<sup>n</sup> <sup>h</sup>* 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 *u<sup>n</sup>*�<sup>1</sup> *<sup>h</sup>* , *pn*�<sup>1</sup> *<sup>h</sup>* , *ϕ<sup>n</sup>*�<sup>1</sup> *h* � �, *u<sup>n</sup> <sup>h</sup>*, *P<sup>n</sup> <sup>h</sup>*, *ϕ<sup>n</sup> h* � �<sup>∈</sup> *<sup>X</sup><sup>f</sup>*,*<sup>h</sup>* � *Qh* � *Xp*,*<sup>h</sup>*, find *<sup>u</sup><sup>n</sup>*þ<sup>1</sup> *<sup>h</sup>* , *pn*þ<sup>1</sup> *<sup>h</sup>* , *<sup>ϕ</sup><sup>n</sup>*þ<sup>1</sup> *h* � �<sup>∈</sup> *<sup>X</sup><sup>f</sup>*,*<sup>h</sup>* � *Qh* � *Xp*,*<sup>h</sup>*, such that for all *v<sup>h</sup>* ∈ *X<sup>f</sup>*,*<sup>h</sup>*, *qh* ∈ *Qh*, *ψ<sup>h</sup>* ∈*Xp*,*<sup>h</sup>*,

$$\begin{cases} \left(\frac{\mathbf{u}\_h^{n+1} - \mathbf{u}\_h^{n-1}}{2\Delta t}, \boldsymbol{\nu}\_h\right) + \left(\nabla \cdot \left(\frac{\mathbf{u}\_h^{n+1} - \mathbf{u}\_h^{n-1}}{2\Delta t}, \boldsymbol{\nu}\_h\right), \nabla \cdot \boldsymbol{\nu}\_h\right) + a\_f\left(\frac{\mathbf{u}\_h^{n+1} + \mathbf{u}\_h^{n-1}}{2}, \boldsymbol{\nu}\_h\right) \\ \quad - \left(\frac{p\_h^{n+1} - p\_h^{n-1}}{2}, \nabla \cdot \boldsymbol{\nu}\_h\right) + g\left(\boldsymbol{\nu}\_h, \boldsymbol{\delta}\_h^n\right)\_{\Gamma} = \left(\boldsymbol{f}\_f^n, \boldsymbol{\nu}\_h\right), \\ \left(q\_h, \nabla \cdot \frac{\mathbf{u}\_h^{n+1} - \mathbf{u}\_h^{n-1}}{2}\right) = 0, \end{cases} \tag{52}$$

$$\mathrm{gS}\_0\left(\frac{\phi\_h^{n+1} - \phi\_h^{n-1}}{2\Delta t}, \phi\_h\right) + a\_p\left(\frac{\phi\_h^{n+1} + \phi\_h^{n-1}}{2}, \boldsymbol{\psi}\_h\right)\_h - \mathbf{g}\left(\mathbf{u}\_h^n, \boldsymbol{\psi}\_h\right)\_\Gamma = \mathbf{g}\left(\boldsymbol{f}\_p^n, \boldsymbol{\psi}\_h\right). \tag{53}$$

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 *κ 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.

In general, the structure model reads as:

$$\begin{cases} \rho\_s \left( \frac{\partial \mathbf{u}\_s}{\partial t} + (\mathbf{u}\_s \cdot \nabla) \mathbf{u}\_s \right) - \operatorname{div} \mathbf{T}(\mathbf{u}\_s, p\_s) = \mathbf{f}\_s, & \forall \mathbf{x} \in \Omega\_t, \\\ \operatorname{div} \mathbf{u}\_t = \mathbf{0}, \quad \forall \mathbf{x} \in \Omega\_t. \end{cases} \tag{54}$$

Here, *ρ<sup>s</sup>* is the density of the structure, *u<sup>s</sup>* is the structure velocity, *ps* is the structure pressure. In structure mechanics, the displacement *d* is usually used as a primary variable ( \_ *d* ¼ *us*), 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 Γ, which may be described as

$$\begin{cases} \mathfrak{u}\_f = \mathfrak{u}\_s, & \text{on } \Gamma, \\ \displaystyle T\left(\mathfrak{u}\_f, p\_f\right)\mathfrak{n} = T(\mathfrak{u}\_s, p\_s)\mathfrak{n}, \text{on } \Gamma. \end{cases} \tag{55}$$

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

Here, *n* denotes one of *n<sup>f</sup>* and *ns*. 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:

$$\begin{cases} \rho\_s \varepsilon \partial\_t \dot{\mathbf{d}} + \mathbf{L}\_s \left( \mathbf{d}, \dot{\mathbf{d}} \right) = -T \left( \mu\_f, p\_f \right) \mathbf{n}, \quad \text{on} \quad \Gamma, \\\ \mathbf{u}\_s = \dot{\mathbf{d}} \quad \text{on} \quad \Gamma. \end{cases} \tag{56}$$

Here, *ε* denotes the structure thickness, *L<sup>s</sup> d*, \_ *d* � � 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 *β*-scheme, which is actually decoupled in the sense that computations are realized subdomain by subdomain [13, 21]. The derivation of the *β*scheme is based on operator splitting.

Partitioned Robin-Neumann scheme:

1.Given the initial solution *u*<sup>0</sup> *<sup>f</sup>* ,*p*<sup>0</sup> *<sup>f</sup>* and *<sup>d</sup>*0.

2.For *m* ¼ 0, 1, 2, 3, … , *N* � 1,

• **Fluid step**: find *u<sup>m</sup>*þ<sup>1</sup> *<sup>f</sup>* and *<sup>p</sup><sup>m</sup>*þ<sup>1</sup> *<sup>f</sup>* such that

$$\begin{cases} \frac{\rho\_f}{\Delta t} \left( \mathbf{u}\_f^{m+1} - \mathbf{u}\_f^m \right) - \operatorname{div} T \left( \mathbf{u}\_f^{m+1}, p\_f^{m+1} \right) = \mathbf{0} & \text{in} \quad \Omega\_f, \\\ \operatorname{div} \mathbf{u}\_f^{m+1} = \mathbf{0} & \text{in} \quad \Omega\_f, \\\ T \left( \mathbf{u}\_f^{m+1}, p\_f^{m+1} \right) \mathbf{n} + \frac{\rho\_i \varepsilon}{\Delta t} \mathbf{u}\_f^{m+1} = \frac{\rho\_i \varepsilon}{\Delta t} \dot{\mathbf{d}}^m + T \left( \mathbf{u}\_f^m, p^m \right) \mathbf{n} & \text{on} \quad \Gamma. \end{cases} (57)$$

• **Structure step**: find **d***<sup>m</sup>*þ<sup>1</sup> *<sup>s</sup>* such that *Decoupling Techniques for Coupled PDE Models in Fluid Dynamics DOI: http://dx.doi.org/10.5772/intechopen.105997*

$$\begin{cases} \frac{\rho\_{\boldsymbol{\beta}}\boldsymbol{\varepsilon}}{\Delta t} \left( \dot{\boldsymbol{d}}^{m+1} - \dot{\boldsymbol{d}}^{m} \right) + \mathbf{L}\_{\boldsymbol{s}} \left( \mathbf{d}^{m+1}, \dot{\mathbf{d}}^{m+1} \right) = -\mathbf{T} \left( \mathbf{u}\_{f}^{m+1}, p\_{f}^{m+1} \right) \mathbf{n} & \text{on} \quad \Gamma, \\\dot{\mathbf{d}}^{m+1} = \frac{1}{\Delta t} \left( \mathbf{d}^{m+1} - \mathbf{d}^{m} \right) & \text{on} \quad \Gamma, \end{cases} \tag{58}$$

3.End

Here, *<sup>ρ</sup>s<sup>ε</sup>* <sup>Δ</sup>*<sup>t</sup>* 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 *β*-**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 *u*<sup>0</sup> *<sup>f</sup>* ,*p*<sup>0</sup> *<sup>f</sup>* and *<sup>d</sup>*0.

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

• **Structure step:** find *u*~*<sup>m</sup>*þ<sup>1</sup> *<sup>s</sup>* such that

$$\begin{cases} \rho\_s \dot{\boldsymbol{u}}\_s^{m+1} - \boldsymbol{u}\_s^m + \mathbf{L}\_s \left( \mathbf{d}^{m+1}, \dot{\mathbf{d}}^{m+1} \right) = -\beta \sigma\_f \left( \mathbf{u}\_f^m, p\_f^m \right) \mathbf{n} & \text{on } \Gamma, \\\ \dot{\mathbf{d}}^{m+1} = \ddot{\mathbf{u}}\_s^{m+1}, \mathbf{d}^{m+1} = \mathbf{d}^m + \Delta t \ddot{\mathbf{u}}\_s^{m+1} & \text{on } \Gamma. \end{cases} \tag{59}$$

$$\bullet \text{ Fluid step: find } \mathfrak{u}\_f^{m+1}, p\_f^{m+1} \text{ and } \mathfrak{u}\_s^{m+1} \text{ such that}$$

$$\begin{cases} \frac{\rho\_f}{\Delta t} \left( \mathbf{u}\_f^{m+1} - \mathbf{u}\_f^m \right) - \operatorname{div} \sigma\_f \left( \mathbf{u}\_f^{m+1}, p\_f^{m+1} \right) = \mathbf{0} & \text{in } \Omega\_f, \\\\ \operatorname{div} \mathbf{u}\_f^{m+1} = \mathbf{0} & \text{in } \Omega\_f, \\\\ \rho\_i e \frac{\mathbf{u}\_i^{m+1} - \tilde{\mathbf{u}}\_i^{m+1}}{\Delta t} = -\sigma\_f \left( \mathbf{u}\_f^{m+1}, p\_f^{m+1} \right) \mathbf{n} + \beta \sigma\_f \left( \mathbf{u}\_f^m, p\_f^m \right) \mathbf{n} & \text{on } \Gamma, \\\\ \mathbf{u}\_f^{m+1} = \mathbf{u}\_i^{m+1} & \text{on } \Gamma. \end{cases} (60)$$

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 *β*scheme, since similar performance is observed for both the Robin-Neumann scheme and the *β*-scheme in numerical experiments. Furthermore, the decoupled multirate *β*-scheme can be extended to more general FSI problems involving nonlinearity, irregular domains, and large structural deformations.

**Figure 3.** *An illustration of a multirate time stepping technique.*

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 *β*-scheme is nothing else but the *β*-scheme itself when the time-step ratio *r* ¼ 1.

Multirate time-stepping *β*-scheme:

1.Given the initial solution *u*<sup>0</sup> *<sup>f</sup>* ,*p*<sup>0</sup> *<sup>f</sup>* and *<sup>d</sup>*0.

	- **Structure steps:** for *m* ¼ *mk*, *mk* þ 1, *mk* þ 2, … , *mk*þ<sup>1</sup> � 1,

$$\begin{cases} \rho\_s \mathbf{\dot{u}}\_s^{m+1} - \mathbf{u}\_s^m + \mathbf{L}\_s \left( \mathbf{d}^{m+1}, \mathbf{\dot{d}}^{m+1} \right) = -\beta \sigma\_f \left( \mathbf{u}\_f^{m\_k}, p\_f^{m\_k} \right) \mathbf{n} & \text{on } \Gamma, \\\ \mathbf{\dot{d}}^{m+1} = \mathbf{\ddot{u}}\_s^{m+1}, \mathbf{d}^{m+1} = \mathbf{d}^m + \Delta t\_s \mathbf{\ddot{u}}\_s^{m+1} & \text{on } \Gamma. \end{cases} \tag{61}$$

$$\begin{cases} \begin{aligned} &\text{\* } \text{Fluid step:}\\ &\frac{\rho\_{f}}{\Delta t\_{f}} \left( \mathbf{u}\_{f}^{m\_{k+1}} - \mathbf{u}\_{f}^{m\_{k}} \right) - \text{div}\,\sigma\_{f} \left( \mathbf{u}\_{f}^{m\_{k+1}}, p\_{f}^{m\_{k+1}} \right) = \mathbf{0} & \text{in } \Omega\_{f},\\ &\text{div}\,\mathbf{u}\_{f}^{m\_{k+1}} = \mathbf{0} & \text{in } \Omega\_{f},\\ &\rho\_{i}\varepsilon \frac{\mathbf{u}\_{f}^{m\_{k+1}} - \tilde{\mathbf{u}}\_{i}^{m\_{k+1}}}{\Delta t\_{f}} = -\sigma\_{f} \left( \mathbf{u}\_{f}^{m\_{k+1}}, p\_{f}^{m\_{k+1}} \right) \mathbf{n} + \beta \sigma\_{f} \left( \mathbf{u}\_{f}^{m\_{k}}, p\_{f}^{m\_{k}} \right) \mathbf{n} & \text{on } \Gamma,\\ &\mathbf{u}\_{f}^{m\_{k+1}} = \mathbf{u}\_{s}^{m\_{k+1}} & \text{on } \Gamma. \end{aligned} \tag{62}$$

3.End.

#### **3.3 Numerical experiment**

In this subsection, we present numerical experiments to demonstrate the convergence and stability performance of the multirate *β*-scheme (60)–(61) for coupling a Stokes flow with a thin-walled structure by using a benchmark model. The model

#### *Decoupling Techniques for Coupled PDE Models in Fluid Dynamics DOI: http://dx.doi.org/10.5772/intechopen.105997*

consists of a *2-D* rectangular fluid domain Ω*<sup>f</sup>* ¼ ½ �� 0, *L* ½ � 0, *R* with *L* ¼ 6 and *R* ¼ 0*:*5 and a *1-D* structure domain Γ ¼ ½ �� 0, *L 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-gramsecond (CGS) system of units.

The physical parameters are set as follows: *ρ<sup>f</sup>* ¼ 1*:*0, *ρ<sup>s</sup>* ¼ 1*:*1, *μ* ¼ 0*:*035; *L<sup>s</sup> d*, \_ *d* <sup>¼</sup> *<sup>c</sup>*1*∂*<sup>2</sup> *<sup>x</sup><sup>d</sup>* <sup>þ</sup> *<sup>c</sup>*0*d*, where *<sup>c</sup>*<sup>1</sup> <sup>¼</sup> *<sup>E</sup><sup>ε</sup>* 2 1ð Þ <sup>þ</sup>*<sup>ν</sup>* , *<sup>c</sup>*<sup>0</sup> <sup>¼</sup> *<sup>E</sup><sup>ε</sup> <sup>R</sup>*<sup>2</sup> <sup>1</sup>�*ν*<sup>2</sup> ð Þ with *<sup>ε</sup>* <sup>¼</sup> <sup>0</sup>*:*1, the Poisson ratio *<sup>ν</sup>* <sup>¼</sup> <sup>0</sup>*:*5 and the Young's modulus *<sup>E</sup>* <sup>¼</sup> <sup>0</sup>*:*<sup>75</sup> � <sup>10</sup>6. A pressure-wave

$$P(t) = P\_{\max} \left( 1 - \cos \left( 2t\pi / T^\* \right) \right) / 2 \quad \text{with} \quad P\_{\max} = 2 \cdot 10^4,\tag{63}$$

is prescribed on the fluid inlet boundary for *<sup>T</sup>*<sup>∗</sup> <sup>¼</sup> <sup>5</sup> � <sup>10</sup>�<sup>3</sup> (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 *d* ¼ **0** at *x* ¼ 0 and *x* ¼ 6.

The first experiment is set up to compare the Robin-Neumann scheme with the *β*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 *β*-scheme, together with the coupled implicit scheme for reference, where the

**Figure 4.** *Geometrical configuration.*

#### **Figure 5.**

*Comparisons of the numerical results obtained by the coupled implicit scheme, the RN scheme, and the β-scheme under the setting: h* <sup>¼</sup> <sup>0</sup>*:*<sup>05</sup> *and* <sup>Δ</sup>*ts* <sup>¼</sup> <sup>10</sup>�<sup>4</sup>*.*

mesh size and the time step size are *<sup>h</sup>* <sup>¼</sup> <sup>0</sup>*:*05 and <sup>Δ</sup>*<sup>t</sup>* <sup>¼</sup> <sup>10</sup>�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 *β*-scheme for investigating the multirate time-step technique.

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 *<sup>h</sup>* <sup>¼</sup> <sup>0</sup>*:*1 and <sup>Δ</sup>*<sup>t</sup>* <sup>¼</sup> <sup>10</sup>�<sup>5</sup> . In addition, we observe that, when <sup>Δ</sup>*ts* Δ*tf* is further increased to be <sup>Δ</sup>*ts* <sup>Δ</sup>*tf* <sup>¼</sup> 5 or <sup>Δ</sup>*ts* <sup>Δ</sup>*tf* ¼ 10, there are substantial numerical instability. This screens out the possibility of using a larger time step size in the structure part.

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 Δ*ts* and *h* while varying the time-step ratio *r* ¼ 1, 5, 10, 20, 50. **Figure 7** displays the computed displacements at *t* ¼ 0*:*015 with the structure time step size <sup>Δ</sup>*ts* <sup>¼</sup> <sup>10</sup>�<sup>5</sup> , 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 *β*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 *β*-scheme is still stable even the time-step ratio is reasonably large.

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 *β*-scheme, and the solution by the multirate *β*-scheme with *r* ¼ 10. By comparing the results, we see that the multirate *β*-scheme provides a very good approximation.

In order to examine the order of convergence, we start with *h* ¼ 0*:*1 and Δ*ts* ¼ 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:

#### **Figure 6.**

*Comparison of the β-scheme and the multirate β-schemes with two different time-step ratios (h* ¼ 0*:*1 *and* <sup>Δ</sup>*<sup>t</sup>* <sup>¼</sup> <sup>10</sup>�<sup>5</sup> *are fixed).*

*Decoupling Techniques for Coupled PDE Models in Fluid Dynamics DOI: http://dx.doi.org/10.5772/intechopen.105997*

**Figure 8.**

*Fluid pressure distribution at t* ¼ 0*:*005, 0*:*010, 0*:*015 *obtained by the coupled implicit scheme (top), the multirate β-scheme with r* ¼ 1 *(middle) and r* ¼ 10 *(bottom) with h* ¼ 0*:*01 *and* Δ*ts* ¼ 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.*

$$\left\{ h, \Delta t\_i \right\}^i = \left\{ \mathbf{0}.\mathbf{1} \cdot \left( \mathbf{0}.5 \right)^i, \mathbf{0}.\mathbf{0}\mathbf{0}\mathbf{0}\mathbf{1} \cdot \left( \mathbf{0}.25 \right)^i \right\}, \ i = \mathbf{0}, \mathbf{1}, \mathbf{2}, \mathbf{3}, \mathbf{4}. \tag{64}$$

We compare the numerical solutions of the multirate *β*-scheme with the reference solution. The reference solution is computed by using the coupled implicit scheme with a high space–time grid resolution *<sup>h</sup>* <sup>¼</sup> <sup>3</sup>*:*<sup>125</sup> � <sup>10</sup>�<sup>3</sup> , <sup>Δ</sup>*<sup>t</sup>* <sup>¼</sup> <sup>10</sup>�<sup>6</sup> � � as that in [14]. In the multirate scheme, *r* ¼ 1 and *r* ¼ 10. The relative errors of the primary variables (*u<sup>f</sup>* , *pf* and *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

#### **Figure 9.**

*Relative errors of primary variables with the spacing h and time step size* Δ*ts in (64). (a) Relative error of uf, (b) Relative error of pf, (c) Relative error of d.*

and the time step size are refined once. Therefore, the multirate *β*-scheme is of a second order in *h* and a first order in *t*.

Finally, in order to demonstrate the advantage of the multirate *β*-scheme, we compare in **Table 3** the CPU times of the concerned numerical algorithms under various settings. We fixed <sup>Δ</sup>*ts* <sup>¼</sup> <sup>10</sup>�<sup>5</sup> and varying the mesh sizes as *<sup>h</sup>* <sup>¼</sup> <sup>1</sup> <sup>10</sup> , <sup>1</sup> <sup>20</sup> , <sup>1</sup> <sup>40</sup> , <sup>1</sup> <sup>80</sup> , <sup>1</sup> 160. From the table, it is observed that the multirate *β*-scheme takes much less computational cost than that of the coupled implicit scheme, particularly when *r* is large.


#### **Table 3.**

*CPU times (in seconds) for the coupled implicit scheme and the multirate β-scheme (with r* ¼ 1 *or* 10*) under different settings of mesh sizes (*Δ*ts* <sup>¼</sup> <sup>10</sup>�<sup>5</sup> *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.

## **Author details**

Mingchao Cai<sup>1</sup> , Mo Mu<sup>2</sup> \* and Lian Zhang<sup>3</sup>

1 Department of Mathematics, Morgan State University, Baltimore, MD, USA

2 Department of Mathematics, The Hong Kong University of Science and Technology, Hong Kong

3 In-Chao Institute Ltd, Shenzhen, China

\*Address all correspondence to: mamu@ust.hk

© 2022 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

## **References**

[1] Beavers G, Joseph D. Boundary conditions at a naturally permeable wall. Journal of Fluid Mechanics. 1967;**30**(1): 197-207

[2] Saffman PG. On the boundary condition at the surface of a porous medium. Studies in Applied Mathematics. 1971;**50**(2):93-101

[3] Jäger W, Mikelić A. On the interface boundary condition of Beavers, Joseph, and Saffman. SIAM Journal on Applied Mathematics. 2000;**60**:1111-1127

[4] Cai M, Mu M, Xu J. Numerical solution to a mixed Navier-Stokes/Darcy model by the two-grid approach. SIAM Journal on Numerical Analysis. 2009;**47** (5):3325-3338

[5] Cesmelioglu A, Girault V, Rivière B. Time-dependent coupling of Navier-Stokes and Darcy flows. ESAIM: Mathematical Modelling and Numerical Analysis. 2013;**47**:539-554

[6] Girault V, Riviére B. DG approximation of coupled Navier-Stokes and Darcy equations by Beaver-Joseph-Saffman interface condition. SIAM Journal on Numerical Analysis. 2009;**47**: 2052-2089

[7] Layton WJ, Schieweck F, Yotov I. Coupling fluid flow with porous media flow. SIAM Journal on Numerical Analysis. 2003;**40**(6):2195-2218

[8] Layton W, Lenferink H. A multilevel mesh independence principle for the Navier-Stokes equations. SIAM Journal on Numerical Analysis. 1996;**33**(1):17-30

[9] Cai M. Modeling and Numerical Simulation for the Coupling of Surface Flow with Subsurface Flow, [PhD thesis]. Hong Kong University of Science and Technology; 2008

[10] Li Z. An augmented cartesian grid method for Stokes¨CDarcy fluid¨Cstructure interactions. International Journal for Numerical Methods in Engineering. 2016;**106**(7): 556-575

[11] Mu M, Xu J. A two-grid method of a mixed stokes-Darcy model for coupling fluid flow with porous media flow. SIAM Journal on Numerical Analysis. 2007;**45** (5):1801-1813

[12] Mu M, Zhu X. Decoupled schemes for a non-stationary mixed Stokes-Darcy model. Mathematics of Computation. 2010;**79**(270):707-731

[13] Bukac M, Muha B. Stability and convergence analysis of the extensions of the kinematically coupled scheme for the fluid-structure interaction. SIAM Journal on Numerical Analysis. 2016;**54**(5): 3032-3061

[14] Fernández M, Mullaert J, Vidrascu M. Explicit Robin–Neumann schemes for the coupling of incompressible fluids with thin-walled structures. Computer Methods in Applied Mechanics and Engineering. 2013;**267**:566-593

[15] Fernández M, Mullaert J, Vidrascu M. Generalized Robin–Neumann explicit coupling schemes for incompressible fluid-structure interaction: Stability analysis and numerics. International Journal for Numerical Methods in Engineering. 2015;**101**(3):199-229

[16] Hou Y. Optimal error estimates of a decoupled scheme based on two-grid finite element for mixed StokesCDarcy model. Applied Mathematics Letters. 2016;**57**:90-96

[17] Huang P, Cai M, Wang F. A Newton type linearization based two grid method for coupling fluid flow with porous media flow. Applied Numerical Mathematics. 2016;**106**:182-198

[18] Zuo L, Hou Y. A decoupling twogrid algorithm for the mixed Stokes-Darcy model with the Beavers-Joseph interface condition. Numerical Methods for Partial Differential Equations. 2014; **30**(3):1066-1082

[19] Zuo L, Hou Y. Numerical analysis for the mixed Navier–Stokes and Darcy problem with the Beavers–Joseph interface condition. Numerical Methods for Partial Differential Equations. 2015; **31**(4):1009-1030

[20] Zhang T, Yuan J. Two novel decoupling algorithms for the steady Stokes-Darcy model based on two-grid discretizations. Discrete and Continuous Dynamical Systems—Series B. 2014;**19** (3):849-865

[21] Čanić S, Muha B, Bukač M. Stability of the kinematically coupled β-scheme for fluid-structure interaction problems in hemodynamics. International Journal of Numerical Analysis and Modeling. 2015;**12**(1):54-80

[22] Rybak I, Magiera J, Helmig R, Rohde C. Multirate time integration for coupled saturated/unsaturated porous medium and free flow systems. Computational Geosciences. 2015;**19**(2): 299-309

[23] Zhang L, Cai M, Mu M. A multirate approach for fluid-structure interaction computation with decoupled methods. Communications in Computational Physics. 2020;**27**(4):1014-1031

[24] Badia S, Codina R. Unified stabilized finite element formulations for the stokes and the Darcy problems. SIAM Journal on Numerical Analysis. 2009; **47**(3):1971-2000

[25] Burman E, Hansbo P. A unified stabilized method for Stokes' and Darcy's equations. Journal of Computational and Applied Mathematics. 2007;**198**(1): 35-51

[26] Cai M, Mu M. A multilevel decoupled method for a mixed Stokes/ Darcy model. Journal of Computational and Applied Mathematics. 2012;**236**(9): 2452-2465

[27] Discacciati M, Miglio E, Quarteroni A. Mathematical and numerical models for coupling surface and groundwater flows. Applied Numerical Mathematics. 2002;**43**(1):57-74

[28] Discacciati M, Quarteroni A. Convergence analysis of a subdomain iterative method for the finite element approximation of the coupling of Stokes and Darcy equations. Computing and Visualization in Science. 2004;**6**(2–3): 93-103

[29] Discacciati M, Quarteroni A, Valli A. Robin-Robin domain decomposition methods for the Stokes-Darcy coupling. SIAM Journal on Numerical Analysis. 2007;**45**(3):1246-1268

[30] Rivière B, Yotov I. Locally conservative coupling of Stokes and Darcy flows. SIAM Journal on Numerical Analysis. 2005;**42**(5):1959-1977

[31] Girault V, Raviart PA. Finite Element Methods for Navier-Stokes Equations, Theory and Algorithms, Springer Series in Computational Mathematics. Vol. Vol. 5. Berlin: Springer; 1986

[32] Elman H, Silvester D, Wathen A. Finite Elements and Fast Iterative Solvers: With Applications in Incompressible Fluid Dynamics. Oxford University Press; 2014

*Decoupling Techniques for Coupled PDE Models in Fluid Dynamics DOI: http://dx.doi.org/10.5772/intechopen.105997*

[33] Nield D, Bejan A. Convection in Porous Media. Vol. Vol. 3. Springer; 2006

[34] Badea L, Discacciati M, Quarteroni A. Numerical analysis of the Navier-Stokes/Darcy coupling. Numerische Mathematik. 2010;**115**(2):195-227

[35] Quarteroni A, Valli A. Domain Decomposition Methods for Partial Differential Equations. Oxford University Press; 1999

[36] Layton W, Tobiska L. A two-level method with backtracking for the Navier-Stokes equations. SIAM Journal on Numerical Analysis. 1998;**35**(5):2035- 2054

[37] Brezzi F, Fortin M. Mixed and Hybrid Finite Element Methods. New York: Springer–Verlag; 1991

[38] Taylor S, Hood P. A numerical solution of the Navier-Stokes equations using the finite element technique. Computers & Fluids. 1973;**1**:73-100

[39] Cai M, Mu M, J. Xu preconditioning techniques for a mixed Stokes/Darcy model in porous media applications. Journal of Computational and Applied Mathematics. 2009;**233**(2):346-355

[40] Kay D, Loghin D, Wathen A. A preconditioner for the steady-state Navier-Stokes equations. SIAM Journal on Scientific Computing. 2002;**24**(1): 237-256

[41] Layton W, Lee H, Peterson J. Numerical solution of the stationary Navier-Stokes equations using a multilevel finite element method. SIAM Journal on Scientific Computing. 1998; **20**:1-12

[42] Dai X, Cheng X. A two-grid method based on Newton iteration for the

Navier-Stokes equations. Journal of Computational and Applied Mathematics. 2008;**220**(1):566-573

[43] Cai M, Huang P, Mu M. Some multilevel decoupled algorithms for a mixed Navier-Stokes/darcy model. Advances in Computational Mathematics. 2017:1-31

[44] Adams RA. Sobolev Spaces. New York: Academic Press; 1975

[45] Layton W, Tran H, Trenchea C. Analysis of long time stability and errors of two partitioned methods for uncoupling evolutionary groundwater– surface water flows. SIAM Journal on Numerical Analysis. 2013;**51**(1):248-272

[46] Layton W, Tran H, Xiong X. Longtime stability of four methods for splitting the evolutionary Stokes-Darcy problem into Stokes and Darcy subproblems. Journal of Computational and Applied Mathematics. 2012;**236**:3198- 3217

[47] Kubacki M. Uncoupling evolutionary groundwater–surface water flows using the Crank-Nicolson Leap-Frog method. Numerical Methods for Partial Differential Equations. 2013;**29**: 1192-1216

[48] Kubacki M, Moraiti M. Analysis of a second-order, unconditionally stable, partitioned method for the evolutionary Stokes–Darcy model. International Journal of Numerical Analysis and Modeling. 2015;**12**:704-730

[49] Chen W, Gunzburger M, Sun D, Wang X. Efficient and long-time accurate second-order methods for Stokes–Darcy system. SIAM Journal on Numerical Analysis. 2013;**51**:2563-2584

[50] Causin P, Gerbeau J, Nobile F. Added-mass effect in the design of partitioned algorithms for fluid– structure problems. Computer Methods in Applied Mechanics and Engineering. 2005;**194**(42):4506-4527

[51] Xu J, Yang K. Well-posedness and robust preconditioners for discretized fluid–structure interaction systems. Computer Methods in Applied Mechanics and Engineering. 2015;**292**: 69-91

[52] Wu Y, Cai X. A fully implicit domain decomposition based ale framework for three-dimensional fluid–structure interaction with application in blood flow computation. Journal of Computational Physics. 2014;**258**: 524-537

[53] Badia S, Nobile F, Vergara C. Fluid– structure partitioned procedures based on Robin transmission conditions. Journal of Computational Physics. 2008; **227**:7027-7051

## *Edited by Bruno Carpentieri and Aamir Shahzad*

Fusion power may offer a long-term energy supply with an uninterrupted power delivery, a high power-generation density, and no greenhouse gas emissions, contributing to preventing the worst effects of climate change and making an enduring contribution to future energy supply. However, the intense conditions inside a fusion power plant (extreme temperatures and high magnetic fields necessary for nuclear fusion) call for addressing several potential problems. These include the development of new materials with extremely high heat tolerances and low enough vapor pressure and the design of mechanical structures that can withstand the electromagnetic force generated as well as feedback controllers to measure and counteract the unstable modes of evolution of the plasma, to name a few. The future of nuclear fusion as an efficient alternative energy source depends largely on techniques that enable us to control these instabilities. Mathematical modelling and physical experiments attempt to overcome some of the hindrances posed by these complexities. This book provides a comprehensive overview of the current state of the art in this fascinating and critically important field of pure and applied physics, mathematics, and engineering, presenting some of the most recent developments in theory, modelling, algorithms, experiments, and applications.

Published in London, UK © 2022 IntechOpen © Petr Ciz / Dollarphotoclub

Advances in Fusion Energy Research - From Theory to Models, Algorithms,

Advances in

Fusion Energy Research

From Theory to Models, Algorithms,

and Applications

*Edited by Bruno Carpentieri* 

*and Aamir Shahzad*

and Applications