# Monolithic approach for coupling Stokes-Darcy in severe regimes

Lara Abouorm, Julien Bruchon, Nicolas Moulin, Sylvain Drapier

► To cite this version:

Lara Abouorm, Julien Bruchon, Nicolas Moulin, Sylvain Drapier. Monolithic approach for coupling Stokes-Darcy in severe regimes. 11e colloque national en calcul des structures, CSMA, May 2013, Giens, France. hal-01722059

HAL Id: hal-01722059

https://hal.science/hal-01722059

Submitted on 2 Mar 2018

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# Monolithic approach for coupling Stokes-Darcy in severe regimes

Lara ABOUORM*, Julien BRUCHON, Nicolas MOULIN, Sylvain DRAPIER

LGF UMR CNRS 5307, Centre SMS, École Nationale Supérieure des Mines de Saint-Étienne
{abouorm,bruchon,nmoulin,drapier}@emse.fr

**Résumé** — This paper presents a numerical strategy for coupling Stokes and Darcy equations in severe regimes (low permeability and curved geometries). P1/P1 elements, stabilized with ASGS method (Algebraic SubgridScale) are used both in Stoke sand Darcy's domains. Two level-set functions are used to describe the interface between Stokes and Darcy, as well as the moving flow front. The accuracy of the presented strategy is shown through various two- and three-dimensional numerical simulations.

**Mots clés** — Stokes-Darcy, ASGS method, resin infusion

The Stokes - Darcy coupling problem has been studied by many researchers in many field of engineering. The main strategies found in the literature to solve the coupling problem are either “decoupled” strategies or “unified” strategies. A decoupled approach consists in using two different discrete spaces to solve Stokes and Darcy’s equations [1]. Conversely, a unified or monolithic approach consists in using the same finite element spaces [2], [12], [8]. The approach proposed in this paper is a unified strategy for coupling Stokes and Darcy’s equations.

Different interpolations orders are known to satisfy the stability for the Stokes and Darcy systems separately. For example, we have the MINI-element (also known as P1 + /P1), Taylor-hood ( [10]) element, penalized stabilized methods [3] and Variational MultiScale Methods for Stokes equations; Least squares methods [4], Raviart-Thomas (RT) and Brezzi-Douglas-Marini elements [5], penalized stabilized methods and Variational MultiScale method [6] for Darcy’s equations.

It is argued that finite element discretizations based on the same finite element spaces for both region will have some advantages in terms of code implementation, but also on the user point of view (only one unstructured mesh to construct). However, the difficulty is that the stable Stokes elements will not be stable in the Darcy region and that most of stable Darcy elements will not be convergent when coupling to Stokes elements. Different unified strategies are used in literature, as P1/P0 with pressure stabilization. This coupling is introduced in [2]. Taylor-Hood elements or Mini- elements for Stokes with BDM or RT elements in Darcy. This coupling is developed in [7]. Mini-elements for Stokes with HVM in Darcy domain, this coupling is developed in [9].

All these methods demonstrate accuracy and convergence in simple cases (permeability equal to 1 m²), but present also difficulties : P1/P0 used in [2] requires special interpolation orders and constant pressure which is not a general case. Taylor-Hood elements with Raviart-Thomas elements used in [7] is expensive in memory cost and CPU time. In addition, the method developed in [9] is the only monolithic approach, used in severe regimes : low permeability, down to 10⁻¹⁵ m² and curved geometries. However, this method shows spurious oscillations and deterioration of convergence in these severe cases. For this reason, we have to introduce a new stabilized finite element methods for Stokes and Darcy coupled problems, able to deal with severe cases.

## 1 mathematical model

Let us consider a bounded domain Ω, made up of two non overlapping subdomains Ωs and Ωd separated by a surface Γ = ∂Ωs ∩ ∂Ωd. Index s is used to denote everything that concerns the purely fluid part (Stokes domain) and index d to denote everything that concerns the porous medium (Darcy's domain). The Stokes equations are then expressed as :

<footnote>*Auteur correspondant</footnote>
Find the velocity $\mathbf{v}_s$ and the pressure $p_s$ fields defined on $\Omega_s$, solution of

$$ 
\begin{aligned}
-div(2\mu\dot{\varepsilon}(\mathbf{v}_s)) + \nabla p_s &= \mathbf{f}_s \quad \text{in} \quad \Omega_s \\
-div \mathbf{v}_s &= h_s \quad \text{in} \quad \Omega_s \\
\mathbf{v}_s &= \mathbf{v}^1 \quad \text{on} \quad \Gamma_{s,D} \\
\boldsymbol{\sigma} \cdot \mathbf{n}_s &= \mathbf{t} \quad \text{on} \quad \Gamma_{s,N}
\end{aligned}
\qquad (1) $$

where $\mathbf{n}_s$ is the unit vector normal to the boundary of $\Omega_s$, $\mathbf{t}$ is the stress vector to be prescribed on $\Gamma_{s,N}$, $\mathbf{v}^1$ the velocity to be prescribed on $\Gamma_{s,D}$, $\mu$ denotes the viscosity, $\dot{\varepsilon}$ is the Eulerian strain rate tensor, and $\boldsymbol{\sigma}_s$ is the Cauchy stress tensor. If the fluid is incompressible then $h_s = 0$.

The Darcy equations are expressed as :

Find the velocity $\mathbf{v}_d$ and the pressure $p_d$ fields defined on $\Omega_d$, solution of

$$ 
\begin{aligned}
\frac{\mu}{k} \mathbf{v}_d + \nabla p_d &= \mathbf{f}_d \quad \text{in} \quad \Omega_d \\
-div \mathbf{v}_d &= h_d \quad \text{in} \quad \Omega_d \\
\mathbf{v}_d \cdot \mathbf{n}_d &= 0 \quad \text{on} \quad \Gamma_{d,D} \\
p_d &= p_{\text{ext}} \quad \text{on} \quad \Gamma_{d,N}
\end{aligned}
\qquad (2) $$

where $k$ is the permeability tensor reduced to a scalar in the isotropic case considered here, $p_{\text{ext}}$ is a pressure to be prescribed on $\Gamma_{d,N}$ and $\mathbf{n}_d$ is the outward unit vector normal to the boundary of $\Omega_d$.

In Stokes-Darcy coupled problem, boundary conditions must be considered on the interface $\Gamma$ of normal $\mathbf{n} = \mathbf{n}_s$. These conditions are :

**Continuity of normal velocity** :

$$ \mathbf{v}_s \cdot \mathbf{n}_s + \mathbf{v}_d \cdot \mathbf{n}_d = 0 \quad \text{on} \quad \Gamma \qquad (3) $$

**Continuity of the fluid normal stress** :

$$ \mathbf{n} \cdot \boldsymbol{\sigma}_s \cdot \mathbf{n} = \mathbf{n} \cdot \boldsymbol{\sigma}_d \cdot \mathbf{n} \quad \text{on} \quad \Gamma \qquad (4) $$

**Beaver-Joseph-Saffman condition** :

$$ 2\mathbf{n} \cdot \dot{\mathbf{\epsilon}}(\mathbf{v}_s) \cdot \tau_j = \frac{\alpha}{\sqrt{k}} (\mathbf{v}_s \cdot \tau_j) \qquad (5) $$

where $\alpha$ is a dimensionless parameter, so-called slip coefficient and $\tau_j$ are the unit vectors tangential to the interface.

## 2 weak formulation

In order to solve the Stokes-Darcy coupled problem by a finite element method, the weak formulation has to be established. We present the weak formulation of Stokes and Darcy separately. The weak formulation of the coupled problem is obtained by summing up the weak formulation of Stokes and Darcy taking into consideration interface conditions. For the sake of simplicity, we choose to write the $L^2$ inner product in $\Omega_{d,s}$ as $\langle, \rangle$.

The spaces of velocity, pressures and test functions are defined by :

$$ L^2(\Omega_i) = \{q, \int_{\Omega_i} q^2 d\Omega < \infty\} $$
$$H^1(\Omega_i)^m = \{q \in L^2(\Omega_i)^m, \nabla q \in L^2(\Omega_i)^{m \times m}\}$$

$$H_{\Gamma_{i,D}}^{1}(\Omega_i)^m = \{q \in H^1(\Omega_i)^m \mid q = 0 \text{ on } \Gamma_{i,D}\}$$

with $i = s$ or $i = d$ and $m$ is the dimension equal to 2 or 3.

The variational formulation of Stokes problem consists in finding a velocity - pressure pair $[v, p]$ such that:

$$B_s([v, p], [w, q]) = L_s([w, q]) \qquad (6)$$

where $w$ and $q$ are weighting functions defined in $H_{\Gamma_{s,D}}^1(\Omega_s)^m$ and $L^2(\Omega_s)$ respectively. The bilinear form $B_s$ and the linear form $L_s$ are defined in Stokes region by

$$\begin{aligned} B_s([v, p][w, q]) &= 2\mu <\dot{\epsilon}(v) : \dot{\epsilon}(w)> - <p, \nabla \cdot w> - <q, \nabla \cdot v> \\ L_s([w, q]) &= <f_s, w> + <h_s, q> + <\sigma_n, w> \end{aligned} \qquad (7)$$

where $\sigma_n$ is the stress vector.

The variational formulation of Darcy's problem consists in finding a velocity-pressure pair $[v, p]$ such that:

$$B_d([v, p], [w, q]) = L_d([w, q]) \qquad (8)$$

where $w$ and $q$ are weighting functions defined in $H_{\Gamma_{d,D}}^1(\Omega_d)^m$ and $L^2(\Omega_d)$.

The bilinear form $B_d$ and the linear form $L_d$ are defined in Darcy by:

$$\begin{aligned} B_d([v, p][w, q]) &= \frac{\mu}{k} <v, w> - <p, \nabla \cdot w> - <q, \nabla \cdot v> \\ L_d([w, q]) &= <f_d, w> + <h_d, q> - <\sigma_n, w> \end{aligned} \qquad (9)$$

The mixed formulation of the Stokes-Darcy problem is established by considering a velocity field $v$ defined on $\Omega$ and a pressure field $p$ defined on $\Omega$ such that $v/\Omega_i = v_i$ and $p/\Omega_i = p_i$ with $i = s, d$. The mixed weak formulation of Stokes-Darcy is obtained by summing up Equations 7 and 9 and taking into consideration the conditions imposed on the Stokes-Darcy interface. Hence, the variational formulation of the Stokes-Darcy coupled problem consist in finding $[v, p] \in V \times Q$ such that:

$$B_c([v, p], [w, q]) = L_c([w, q])$$

where $[w, q]$ are weighting functions defined in $V \times Q$.

$$V = H_{\Gamma_{s,D}}^1 \times H_{\Gamma_{d,D}}^1(\mathrm{div}, \Omega_d)$$

$$Q = L^2(\Omega_s) \times L^2(\Omega_d)$$

The bilinear form $B_c$ and the linear form $L_c$ are defined by:

$$\begin{aligned} B_c([v, q], [w, q]) &= \int_{\Omega} 2\mu \dot{\epsilon}(v) : \dot{\epsilon}(w) H_s d\Omega \\ &+ \int_{\Omega} \frac{\mu}{k} v \cdot w H_d d\Omega \\ &- \int_{\Omega} p \operatorname{div} w d\Omega - \int_{\Omega} q \operatorname{div} v d\Omega + \int_{\Gamma} \frac{\alpha \mu}{\sqrt{k}} v \cdot w H_d d\Gamma \\ L_c([v, p], [w, q]) &= <f, w> + <h, q> \end{aligned} \qquad (10)$$

where $(f, h)$ are defined by $(f_d, h_d)$ in Darcy and $(f_s, h_s)$ in Stokes. $H_i$ is an Heaviside function equal to 1 in domain $i$ and vanishing elsewhere.
### 3 finite element method for Stokes-Darcy problem

The whole computational domain $\Omega \subset \mathbb{R}^m$ is discretized with one single unstructured mesh. This mesh is made up of triangles when $m = 2$ and of tetrahedrons for $m = 3$. Let $V_h$ and $Q_h$ be the finite element spaces of the continuous piecewise linear functions $\mathbf{v}_h$ and $p_h$. The Galerkin approximation of both the Stokes and the Darcy problems requires the use of velocity-pressure interpolation that satisfies the adequate inf-sup condition. Different interpolation pairs are known to satisfy this condition for each problem independently, but the key issue is to find interpolations that satisfy both at same time. In this paper, we choose the use of stabilized finite element methods. The philosophy of the stabilized methods is to strengthen classical variational formulations so that discrete approximation which would otherwise be unstable becomes stable and convergent. One of the stabilized finite element method for Stokes-Darcy problem is VMS (Variational Multiscale Method) [6], [11]. The basic idea of this method is to approximate the effect of the component of the continuous solution which cannot be captured by the finite element solution. It consists in splitting the continuous solution for velocity and pressure into two components, one coarse corresponding to the finite element scale $[\mathbf{v}_h, p_h]$, and a finer component corresponding to lower scale $[\mathbf{v}', p'] : \mathbf{v} = \mathbf{v}_h + \mathbf{v}'$ and $p = p_h + p'$, $V \times Q = (V_h \times Q_h) \oplus (V' \times Q')$. Invoking this decomposition in the continuous problem for both the solution and test functions, we get the two scales systems :

$$B_c([\mathbf{v}_h, p_h]), [\mathbf{w}_h, q_h]) + B_c([\mathbf{w}', p'], [\mathbf{w}_h, q_h]) = L_c([\mathbf{w}_h, q_h]) \quad (11)$$

$$B_c([\mathbf{v}_h, p_h], [\mathbf{w}', p']) + B_c[\mathbf{v}', p'], [\mathbf{w}', q']) = L_c([\mathbf{w}', q']) \quad (12)$$

for all $[\mathbf{w}_h, q_h] \in V_h \times Q_h$ and $[\mathbf{w}', q'] \in V' \times Q'$. After approximating (12) with an algebraic formulation, by introducing the operator of projection $P'$ onto $V'$, the approximated fields $[\mathbf{v}', p']$ are taken into account in the finite element problem (11). Hence, we get the stabilized forms of the bilinear and linear forms in Stokes, Darcy, and Stokes-Darcy coupled problem. The stabilized problem in Stokes can be written as follow :

Find $[\mathbf{v}_h, p_h] \in V_h \times Q_h$ such as

$$B_{s,stable}([\mathbf{v}_h, p_h], [\mathbf{w}_h, q_h]) = L_{s,stable}([\mathbf{w}_h, q_h]) \quad (13)$$

The bilinear stabilized form $B_{s,stable}$ and the linear stabilized form $L_{s,stable}$ are defined by

$$\begin{align}
B_{s,stable}([\mathbf{v}_h, p_h], [\mathbf{w}_h, q_h]) &= B_s([\mathbf{v}_h, p_h], [\mathbf{w}_h, q_h]) \tag{14} \\
&\quad + \tau_q \sum_K \langle P'_{h,p}(\nabla \cdot \mathbf{v}_h), \nabla \cdot \mathbf{w}_h \rangle + \tau_v \sum_K \langle P'_{h,u}(-\mu \Delta \mathbf{v}_h + \nabla p_h), \mu \Delta \mathbf{w}_h - \nabla q_h \rangle \nonumber \\
L_{s,stable}([\mathbf{w}_h, q_h]) &= L_s([\mathbf{w}_h, q_h]) \nonumber \\
&\quad + \tau_q \sum_K \langle P'_{h,p}(h_s), \nabla \cdot \mathbf{w}_h \rangle + \tau_v \sum_K \langle P'_{h,p}(f_s), \mu \Delta \mathbf{w}_h - \nabla q_h \rangle \nonumber
\end{align}$$

where $B_s([\mathbf{v}_h, p_h], [\mathbf{w}_h, q_h])$ and $L_s([\mathbf{w}_h, q_h])$ are defined in (7). $\tau_q$ and $\tau_v$ are the stabilization parameters (obtained by Fourier analysis) that we compute as

$$\begin{align}
\tau_q &= c_1 \mu \nonumber \\
\tau_v &= c_1 \mu h_k^2 \tag{15}
\end{align}$$

where $c_1$ is an algorithmic constant and $h_k$ is the size of mesh. $P'_{h,u}$ is the broken $L^2$ projection onto $V'$ and $P'_{h,p}$ is the broken $L^2$ projection onto $Q'$. The simplest approach is to take $P'_h$ as the identity operator when acting on the FE residual. Assuming this, we obtain the stabilized method called Algebraic Subgrid
Scale (ASGS) method [6, 11]. Invoking this expression of $P'_h$, we get the following stabilized finite form $B_s$ and $L_s$ in Stokes

$$
\begin{align}
B_{s,\text{stable}}([\mathbf{v}_h, p_h], [\mathbf{w}_h, q_h]) &= B_s([\mathbf{v}_h, p_h], [\mathbf{w}_h, q_h]) \tag{16} \\
&\quad + \tau_q \sum_K \langle \nabla \cdot \mathbf{v}_h, \nabla \cdot \mathbf{w}_h \rangle + \tau_v \sum_K \langle -\mu \Delta \mathbf{v}_h + \nabla p_h, \mu \Delta \mathbf{w}_h - \nabla q_h \rangle \nonumber \\
L_{s,\text{stable}}([\mathbf{w}_h, q_h]) &= L_s([\mathbf{w}_h, q_h]) \nonumber \\
&\quad + \tau_q \sum_K \langle h, \nabla \cdot \mathbf{w}_h \rangle + \tau_v \sum_K \langle f, \mu \Delta \mathbf{w}_h - \nabla q_h \rangle \nonumber
\end{align}
$$

By using ASGS method, the stabilized problem in Darcy can be written as follow :

Find $[\mathbf{v}_h, p_h] \in V_h \times Q_h$ such as :

$$B_{d,\text{stable}}([\mathbf{v}_h, p_h][\mathbf{w}_h, q_h]) = L_{d,\text{stable}}([\mathbf{w}_h, q_h]) \quad (17)$$

$$
\begin{align}
B_{d,\text{stable}}([\mathbf{v}_h, p_h], [\mathbf{w}_h, q_h]) &= B_d([\mathbf{v}_h, p_h], [\mathbf{w}_h, q_h]) \nonumber \\
&\quad + \tau_p \sum_K \langle \nabla \cdot \mathbf{v}_h, \nabla \cdot \mathbf{w}_h \rangle + \tau_u \sum_K \langle \frac{\mu}{k} \mathbf{v}_h + \nabla p_h, -\frac{\mu}{k} \mathbf{w}_h - \nabla q_h \rangle \nonumber \\
L_{d,\text{stable}}([\mathbf{w}_h, q_h]) &= L_d([\mathbf{w}_h, q_h]) \nonumber \\
&\quad + \tau_p \sum_K \langle h_d, \nabla \cdot \mathbf{w}_h \rangle + \tau_u \sum_K \langle \mathbf{f}_d, -\frac{\mu}{k} \mathbf{w}_h - \nabla q_h \rangle \tag{18}
\end{align}
$$

where $B_d([\mathbf{v}_h, p_h], [\mathbf{w}_h, q_h])$ and $L_d([\mathbf{w}_h, q_h])$ are defined in 9. $\tau_p$ and $\tau_u$ are the stabilization parameters, that we compute as

$$
\begin{align}
\tau_p &= c_p \frac{\mu}{k} l_p^2 \nonumber \\
\tau_u &= (c_u \frac{\mu}{k} l_u)^{-1} h_k^2 \tag{19}
\end{align}
$$

with $c_p$ and $c_u$ algorithmic constants. $l_u$ and $l_p$ are length scales, chosen as $(L_0 h_k)^2$, where $L_0$ is a characteristic length of domain and $h_k$ is the size of mesh.

For Stokes and Darcy flow coupled through the interfaces, the problem stabilized with ASGS can be written as

Find $[\mathbf{v}_h, p_h] \in V_h \times Q_h$ such as

$$B_{c,\text{stable}}([\mathbf{v}_h, p_h][\mathbf{w}_h, q_h]) = L_{c,\text{stable}}([\mathbf{w}_h, q_h]) \quad (20)$$

$$
\begin{align}
B_{c,\text{stable}}([\mathbf{v}_h, p_h], [\mathbf{w}_h, q_h]) &= 2\mu \int_{\Omega} H_s \dot{\varepsilon}(\mathbf{v}_h) : \dot{\varepsilon}(\mathbf{w}_h) + \frac{\mu}{k} \int_{\Omega} H_d \mathbf{w}_h \mathbf{v}_h \nonumber \\
&\quad - \int_{\Omega} \nabla \cdot \mathbf{w}_h p_h - \int_{\Omega} \nabla \cdot \mathbf{v}_h q_h + \tau_{p,c} \int_{\Omega} \nabla \mathbf{v}_h \nabla \mathbf{w}_h \nonumber \\
&\quad + \int_{\Gamma} \alpha \frac{\mu}{\sqrt{k}} (\mathbf{v}_h \cdot \mathbf{\tau})(\mathbf{w}_h \cdot \mathbf{\tau}) + \tau_{u,c} \int_{\Omega} \langle -\Delta \mathbf{v}_h + \frac{\mu}{k} \mathbf{v}_h + \nabla p_h, \mu \Delta \mathbf{w}_h - \frac{\mu}{k} \mathbf{w}_h - \nabla q_h \rangle \nonumber \\
L_{c,\text{stable}}([\mathbf{w}_h, q_h]) &= \langle f, \mathbf{w}_h \rangle + \langle h, q_h \rangle \nonumber \\
&\quad + \tau_{p,c} \langle h_c, \nabla \cdot \mathbf{w}_h \rangle + \tau_{u,c} \langle \mathbf{f}_c, \mu \Delta \mathbf{w}_h - \frac{\mu}{k} \mathbf{w}_h - \nabla q_h \rangle \tag{21}
\end{align}
$$

where $\tau_{p,c}$, $\tau_{u,c}$ are the stabilization parameters, that we compute as
$$ 
\begin{aligned}
\tau_{p,c} &= c_p \frac{\mu}{k} l_p^2 + c_1 \mu \\
\tau_{u,c} &= (c_1 \mu + c_u \frac{\mu}{k} l_u)^{-1} h_k^2
\end{aligned} 
\qquad (22) $$

where $H_i$ is a Heaviside function equal to 1 in domain $i$ and vanishing elsewhere. The surface integral $\int_{\Gamma} \alpha_k^\mu(\mathbf{v}_h \cdot \mathbf{\tau})(\mathbf{w}_h \cdot \mathbf{\tau})$ is turned into a volume integral for calculation simplification.

## 4 Numerical Validation

We have investigated the convergence of ASGS method for a coupled Stokes - Darcy flow by using the manufactured solution method. The theoretical convergence orders are well recovered, even they are not presented in this paper. In this section, the accuracy and robustness of the proposed method are evaluated in two situations. First, the case of a flow perpendicular to the Stokes - Darcy interface, as sketched in Figure 1 a), is investigated. We compare our results with another monolithic approach using MINI-element (P1 + /P1) in Stokes and HVM (Hughes Variational Multi-scale) for the stabilization in Darcy. Figure 1 b) compares the normal velocity $v_y$ obtained with ASGS to the one obtained by P1+/P1-HVM method. Theoretically, normal velocity should be equal to $10^{-9}$ m/s for a permeability $k = 10^{-14}$m$^2$. Figure 1 b) shows that the normal velocity computed with ASGS method is continuous through the Stokes - Darcy interface, and does not exhibit spurious oscillations contrary to the results obtained with P1+/P1-HVM method.

The robustness of ASGS method is demonstrated on more realistic geometries of interest for real life, especially in the context of the manufacturing processes of composite materials. Two- and three-dimensional simulations are carried out with a very small height of the Stokes domain with respect to the Darcy's domain height, and a very small permeability ($k < 10^{-14}$ m$^2$). Figure 2 a) and 2 b) show the isovvalues of pressure obtained for the 2D and 3D simulations obtained with ASGS method. Despite the curved interface, velocity and pressure fields computed with the monolithic approach do not exhibit any spurious oscillations.

Finally, the mechanical problem with ASGS method is coupled with the level-set problem used to capture the resin flow front to simulate the transient fluid flow in Stokes-Darcy domains. A 3D simulation corresponding to a more complex part, essentially to illustrate the 3D resin flow is conducted. The results of this simulation are shown in figure 3. Moreover, the simulation is relatively robust since the flow front is not disturbed when it passes through the interface between Stokes and Darcy.

<figure><img src="image_1.png" /><figcaption>Fig. 1 – (a)Boundary conditions of a flow perpendicular to Stokes - Darcy interface. (b) Comparison of normal velocity between P1+/P1-HVM and ASGS methods for $x = 2.5$ ($k = 10^{-14}$m$^2$, $\mu = 1$Pa.s, $\alpha = 1$).</figcaption></figure>
<figure><img src="image_2.png" /><figcaption>Fig. 2 – (a) Isovalues of pressure for 2D simulation of a flow with a curved interface. (b) Isovalues of pressure for 3D simulation of a flow with a curved interface.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 3 – 3D simulation of a transient flow into a Stokes-Darcy domain.</figcaption></figure>

# 5 Conclusions

The monolithic approach has been developed to solve Stokes-Darcy coupled problem. A stabilized finite element method has been proposed to stabilize Stokes, Darcy and Stokes-Darcy coupled problem in the case where the Brezzi-Babuska stability condition is not fulfilled. This stabilization method is based on a variational multiscale technique called ASGS method. Convergence of this method was validated and compared with another monolithic approach [9], showing the accuracy of the results. Also, it was shown that the method could be used for more complex geometries and with a wide range of permeability, which is essential for the industrial applications.

# Références

[1] B. Riviere and I. Yotov, *Locally conservative coupling of Stokes and Darcy flows*, SIAM on Numerical Analysis, 42, 5, pp.1959-1977, 2005.

[2] E. Burmana and P. Hansbob, *A unified stabilized method for Stokes and Darcy equations*, Journal of Computational and Applied Mathematics, 198, pp.35 - 51, 2007

[3] J. Douglas and J. wang, *An absolutely Stabilized Finite Element Method for the Stokes problem*, Mathematics of computation, 52, 186, pp.495 - 508, April 1989.

[4] P. Bochev and M. Gunzburger, *Least-squares finite element methods for elliptic equations*, SIAM J. Numer. Anal., pp.789 - 837, 404, 1998.

[5] P.A.Raviart and J.M.Thomas, *A mixed-finite element method for second order elliptic problems*, Springer, New York, Mathematical aspects of the finite element method, Lecture Notes in Mathematics, 1977.

[6] S. Badia and R. Codina, *Stabilized Continuous and discontinuous Galerkin techniques for Darcy flow*, Computer Methods in Applied Mechanics and Engineering, 199, 25-28, pp.654-667, 2010.
[7] W.J.Layton and F.schieweck and I.Yotov, *coupling fluid flow with porous media flow*, SIAM J. Numer. Anal, pp.2195-2218,40,6,2003.

[8] T. Arbogast and D.S. Brunson *A computational method for approximation a Darcy-Stokes system governing a vuggy porous medium*, Computational Geoscience, 11,3, pp.207-218, 2007

[9] G. Pacquaut, and J. Bruchon, and N. Moulin, and S. Drapier, *Combining a level set method and mixed stabilized P1/P1 formulation for coupling stokes-Darcy flows*, International Journal for Numerical Methods in Fluids, 69, pp.459-480,2012.

[10] C.Taylor and P.Hood, *A numerical solution of the Navier-Stokes equations using the finite element*, Computers and Fluids, 1,1, pp 73-100, 1973.

[11] S.Badia and R.Codina, *Unified stabilized finite element formulation for the stokes and the Darcy problems*, SIAM journal on Numerical Analysis, 47 :1971â2000, 2008.

[12] T.Karper and K.A.Mardal and R.Winther, *Unified finite element discretizations of coupled Darcy-Stokes flow*. *Numerical Methods for Partial Differential Equations*,SIAM on Numerical Analysis,25,2,pp.311-326,2009.