1. Coupling Stokes and Darcy equations

1.1. Problem setting

Definition

Let \(\Omega\subset\mathbb R^d\) be a domain split into two subdomains \(\Omega_f\) and \(\Omega_p\) such that :

\[ \begin{array}{l} \Omega=\overline{\Omega_f}\cup\overline{\Omega_p}\\ \Omega_f\cap\Omega_p=\emptyset\\ \overline{\Omega_f}\cap\overline{\Omega_p}=\Gamma \end{array}\]

We call \(\Gamma\) the interface.

Definition

We denote by \(\underline n_f\text{ and }\underline n_p\) the respective unit outward normal vector fields along the boundary of \(\Omega_f\text{ and }\Omega_p\) respectively.

Remark

We have : \(\underline n_f = -\underline n_p\text{ on }\Gamma\).

Definition

For each \(i\in\{1,...,d-1\}\), we define an unit vector field \(\tau_i\) tangent to \(\Gamma\) such that :

\[ \forall\underline x\in\Gamma,\quad \left(\underline n_f,\underline\tau_1,...,\underline\tau_{d-1}\right)\text{ in an orthonormal basis for }\mathbb R^d\]

We aim to solve the Stokes equation in domain \(\Omega_f\) and the Darcy equation in domain \(\Omega_p\), subject to boundary conditions on \(\partial\Omega=(\partial\Omega_f\backslash\Gamma)\cup(\partial\Omega_p\backslash\Gamma)\) together with coupling relations along \(\Gamma\). We denote by \(\Gamma_f\text{ and }\Gamma_p\) the outer boundary of \(\Omega_f\text{ and }\Omega_p\) respectively.

We use the following notations :

  • \(\underline u_f\text{ and }p_f\) are the fluid velocity and fluid pressure respectively within the Stokes domain \(\Omega_f\),

  • \(\underline u_p\text{ and }\varphi\) are the fluid velocity and piezometric head respectively within the Darcy domain \(\Omega_p\),

  • \(\underline{\underline D}(\underline u)\) is the symmetrized gradient defined by :

\[ \underline{\underline D}(\underline u)=\frac12(\nabla\underline u+\nabla\underline u^T)\]

  • \(\underline{\underline T}(\underline u,p)\) is the Cauchy stress tensor defined by :

\[ \underline{\underline T}(\underline u,p)=2\nu\underline{\underline D}(\underline u)-p\underline{\underline I}_d\]

1.2. Equations

Writing Stokes and Darcy problems on their respective domains yields the system :

\[ \left\{\begin{array}{rll} -\nabla\cdot\underline{\underline T}(\underline u_f,p_f) & =\underline f_f & \text{ in }\Omega_f\\ \nabla\cdot\underline u_f & = g_f & \text{ in }\Omega_f\\ -\nabla\cdot\underline{\underline\kappa}\nabla\varphi & =f_p & \text{ in }\Omega_p \end{array}\right.\]

1.3. Boundary conditions

The same essential/natural splitting as in the sole Darcy problem occurs here. One thus has to split the outer boundary \(\partial\Omega\) in four relatively open disjoint pieces. The subscripts \(f,p\) denote objects relative to the Stokes and The Darcy problem respectively, and the subscripts \(n,e\) denote respectively natural and essential boundary condition.

Definition
\[ \left.\begin{array}{l} \overline{\Gamma_{f,n}}\cup\overline{\Gamma_{f,e}}=\partial\Omega_f\backslash\Gamma\\ \Gamma_{f,n}\cap\Gamma_{f,e}=\emptyset\\ \overline{\Gamma_{p,n}}\cup\overline{\Gamma_{p,e}}=\partial\Omega_p\backslash\Gamma\\ \Gamma_{p,n}\cap\Gamma_{p,e}=\emptyset\\ \end{array}\right.\]

One can impose :

\[ \left\{\begin{array}{rll} \underline u_f & = \underline u_{f,e} & \text{ on }\Gamma_{f,e} \\ \underline{\underline T}(\underline u_f,p_f)\underline n_f & = \underline T_{f,n} & \text{ on }\Gamma_{f,n} \\ \varphi & =\varphi_{e} & \text{ on }\Gamma_{p,e} \\ (-\underline{\underline\kappa}\nabla\varphi )\cdot\underline n_p & = u_{p,n} & \text{ on }\Gamma_{p,n} \end{array}\right.\]

1.4. Coupling conditions

We impose the following coupling conditions along the interface \(\Gamma\).

  • Continuity of the normal velocity, i.e. \(\underline u_f\cdot\underline n_f + \underline u_p\cdot\underline n_p = 0\), or :

\[ \underline u_f\cdot\underline n_f = -(\underline{\underline\kappa}\nabla\varphi)\cdot\underline n_f\quad\text{ on }\Gamma\]

  • Continuity of the normal component of the Cauchy stress :

\[ -\underline n_f\cdot(\underline{\underline T}(\underline u_f,p_f)\underline n_f)=g\varphi\quad\text{ on }\Gamma\]

  • Beavers-Joseph condition :

\[ (\underline u_f - \underline u_p)\cdot\underline\tau_i+\alpha\underline\tau_i\cdot(\underline{\underline T}(\underline u_f,p_f)\underline n_f)=0\quad\text{ on }\Gamma\quad i\in\{1,...,d-1\}\]

where \(\alpha\) is proportional to the square root of the permeability of \(\Omega_p\) and inversely proportional to the square root of the kinematic viscosity \(\nu\) and depends on the characteristic length of the pores in \(\Omega_p\).

It is common to neglect the tangential components of \(\underline u_p\), thus the Beavers-Joseph condition reduces to the Beavers-Joseph-Saffman condition :

\[ \underline u_f\cdot\underline\tau_i+\alpha\underline\tau_i\cdot(\underline{\underline T}(\underline u_f,p_f)\underline n_f)=0\quad\text{ on }\Gamma\quad i\in\{1,...,d-1\}\]

See D. A. Nield, The Beavers-Joseph boundary condition and related matters : a historical and critical note, Transp. Porous Med. (2009) 78:537-540 and M. Hoffmann, The Navier-Stokes-Darcy problem, Master’s Thesis (2013). This identity comes from empirical observations.

1.5. Weak formulation

We basically just write independently the weak formulation for Darcy and Stokes problem, then make the coupling condition appear by substitution.

Remark

In the following, function restrictions on a boundary are to be understood in terms of trace operator.

Definition

In order to write the weak formulation, we build the following function spaces from the boundary conditions :

\[ \begin{array}{l} H_f=\left\{\underline v=(v_i)_{i=1,...,d}\in(H^1(\Omega_f))^d\quad|\quad\left.\underline v\right|_{\Gamma_{f,e}}=0\right\}\\ H_p=\left\{\psi\in H^1(\Omega_p)\quad|\quad\left.\psi\right|_{\Gamma_{p,e}}=0\right\} \end{array}\]

in which we will pick our test functions. We endow them with the usual norms :

\[ \begin{array}{l} \|\underline v\|_{1,\Omega_f}=\left(\sum_{i=1}^d\|v_i\|_2^2+\int_{\Omega_f}(\nabla\underline v:\nabla\underline v)\right)^{\frac12}\\ \|\psi\|_{1,\Omega_p}=\left(\|\psi\|_2^2+\int_{\Omega_p}\nabla\psi\cdot\nabla\psi\right)^{\frac12} \end{array}\]

inherited from the scalar products :

\[\begin{array}{l} \langle\underline v,\underline w\rangle_{\Omega_f}=\left(\sum_{i=1}^d\left(\int_{\Omega_f}|v_i||w_i|+\int_{\Omega_f}(\nabla\underline v:\nabla\underline w)\right)\right)^{\frac12}\\ \langle\varphi,\psi\rangle_{\Omega_p}=\left(\int_{\Omega_p}|\varphi||\psi|+\int_{\Omega_p}\nabla\varphi\cdot\nabla\psi\right)^{\frac12}\\ \end{array}\]
Definition

Moreover, we endow \(X=H_f\times H_p\) with the norm :

\[\|\overline{\underline w}\|_X=\left(\|\underline w\|_{1,\Omega_f}^2+\|\psi\|_{1,\Omega_p}^2\right)^\frac12\]

From Stokes' equation, we have to take the divergence of the stress tensor, in particular of the deformation tensor, times a test function \(\underline v\in H_f\). We use the well-known identity :

\[\nabla\cdot\left(D(\underline u)\underline v\right)=(\nabla\cdot D(\underline u))\cdot\underline v + D(\underline u):\nabla\underline v\]

where \(A:B=Tr(AB^T)\) denotes the generalized scalar product. Since \(D(\underline u)\) is symmetric, we can replace \(\nabla\underline v\) by \(D(\underline v)\), according to the following relations : \(Tr(AB^T)=Tr(BA^T)=Tr(BA)=Tr(AB)\) with \(A\) symmetric and any \(B\).

We define the following bilinear forms :

\[a_f(\underline u,\underline v) = 2\nu\int_{\Omega_f}\underline{\underline D}(\underline u):\underline{\underline D}(\underline v)\quad;\quad b_f(\underline v,p) = -\int_{\Omega_f}p\nabla\cdot\underline v\quad;\quad a_p(\varphi,\psi)=\int_{\Omega_p}\nabla\psi\cdot\underline{\underline\kappa}\nabla\varphi\]

1.5.1. Stokes

Let \(\underline v\in H_f\). Multiplying and integrating gives :

\[\int_{\Omega_f}(-\nabla\cdot\underline{\underline T}(\underline u_f,p_f))\cdot\underline v=\int_{\Omega_f}\underline f\cdot\underline v\]

We now focus on the left hand side term.

\[ \begin{array}{rl} \displaystyle\int_{\Omega_f}(-\nabla\cdot\underline{\underline T}(\underline u_f,p_f))\cdot\underline v =&\displaystyle\int_{\Omega_f}(\nabla\cdot(p_f\underline{\underline I}_d-2\nu \underline{\underline D}(\underline u_f)))\cdot\underline v\\ =&\displaystyle\int_{\Omega_f}\nabla p_f\cdot\underline v-2\nu(\nabla\cdot\underline{\underline D}(\underline u_f))\cdot\underline v\\ =&\displaystyle\int_{\Omega_f}\nabla\cdot(p_f\underline v)-\int_{\Omega_f}p_f\nabla\cdot\underline v+2\nu\int_{\Omega_f}\underline{\underline D}(\underline u_f):D(\underline v)-2\nu\int_{\Omega_f}\nabla\cdot(\underline{\underline D}(\underline u_f)\cdot\underline v)\\ =&a_f(\underline u_f,\underline v)+b_f(\underline v,p_f)-\displaystyle\int_{\Omega_f}\nabla\cdot\left(\underline{\underline T}(\underline u_f,p_f)\underline n_f\cdot\underline v\right)\\ =&a_f(\underline u_f,\underline v)+b_f(\underline v,p_f)-\displaystyle\int_{\Gamma_{f,e}\cup\Gamma_{f,n}\cup\Gamma}\underline v\cdot(\underline{\underline T}(\underline u_f,p_f)\underline n_f) \end{array}\]

Notice that the \(\Gamma_{f,e}\) part vanishes since \(\left.\underline v\right|_{\Gamma_{f,e}}=0\). Writing \(\displaystyle\underline v=(\underline v\cdot\underline n_f)\underline n_f+\sum_{i=1}^{d-1}(\underline v\cdot\underline\tau_i)\underline\tau_i\) on the boundary gives :

\[\int_{\Gamma}-\underline v\cdot(\underline{\underline T}(\underline u_f,p_f)\underline n_f)=-\int_\Gamma(\underline v\cdot\underline n_f)(\underline{\underline T}(\underline u_f,p_f)\underline n_f)\cdot\underline n_f-\sum_{i=1}^{d-1}\int_\Gamma(\underline v\cdot\underline\tau_i)(\underline{\underline T}(\underline u_f,p_f)\underline\tau_i)\cdot\underline n_f\]

allowing us to apply the coupling conditions :

\[\int_\Gamma-\underline v\cdot(\underline{\underline T}(\underline u_f,p_f)\underline n_f)=\int_\Gamma(\underline v\cdot\underline n_f)g\varphi + \frac1\alpha\sum_{i=1}^{d-1}\int_\Gamma(\underline v\cdot\underline\tau_i)(\underline u_f-\underline u_p)\cdot\underline\tau_i\]

so that the total weak Stokes equation is :

\[a_f(\underline u_f,\underline v)+b_f(\underline v,p_f)+\int_{\Gamma}\left((\underline v\cdot\underline n_f)g\varphi + \frac1\alpha\sum_{i=1}^{d-1}(\underline v\cdot\underline\tau_i)(\underline u_f-\underline u_p)\cdot\underline\tau_i\right) = \int_{\Omega_f}\underline f\cdot\underline v - \int_{\Gamma_{f,n}}\underline v\cdot\underline T_{f,n}\]

We may often use Saffmann’s approximation and a zero Neumann condition on the outer boundary, yielding :

\[a_f(\underline u_f,\underline v)+b_f(\underline v,p_f)+\int_{\Gamma}\left((\underline v\cdot\underline n_f)g\varphi+\frac1\alpha\sum_{i=1}^{d-1}(\underline v\cdot\underline\tau_i)(\underline u_f\cdot\underline\tau_i)\right) = \int_{\Omega_f}\underline f\cdot\underline v-\int_{\Gamma_{f,n}}\underline v\cdot\underline T_{f,n}\]

1.5.2. Darcy

Since \(\underline u_p=\underline{\underline\kappa}\nabla\varphi\) and \(\nabla\cdot\underline u_p=0\), we rewrite Darcy’s equation :

\[-\nabla\cdot(\underline{\underline\kappa}\nabla\varphi)=0\text{ on }\Omega_p\]

Let us write its weak formulation :

\[ \begin{array}{rl} \displaystyle\int_{\Omega_p}-\nabla\cdot(\underline{\underline\kappa}\nabla\varphi)\psi=&\displaystyle\int_{\Omega_p}\underline{\underline\kappa}\nabla\varphi\cdot\nabla\psi-\int_{\Omega_p}\nabla\cdot(\underline{\underline\kappa}\psi\nabla\varphi)\\ =&\displaystyle a_p(\varphi,\psi)-\int_{\Gamma_{p,e}\cup\Gamma_{p,n}\cup\Gamma}\underline{\underline\kappa}\psi\nabla\varphi\cdot\underline n_p\\ =&\displaystyle a_p(\varphi,\psi)+\int_{\Gamma}\underline{\underline\kappa}\psi\nabla\varphi\cdot\underline n_f + \int_{\Gamma_{p,n}}\underline{\underline\kappa}\psi\nabla\varphi\cdot\underline n_f\quad\text{ thanks to b.c.}\\ =&\displaystyle a_p(\varphi,\psi)-\int_{\Gamma}\psi(\underline u_f\cdot\underline n_f) - \int_{\Gamma_{p,n}}\psi u_{p,n}\quad\text{ thanks to c.c.} \end{array}\]

so the weak Darcy is :

\[\displaystyle a_p(\varphi,\psi)-\int_{\Gamma}\psi(\underline u_f\cdot\underline n_f) = \int_{\Gamma_{p,n}}\psi u_{p,n}\]

1.5.3. Sum up

\[ \left\{\begin{array}{rll} a_f(\underline u_f,\underline v)+b_f(\underline v,p_f)-\displaystyle\int_{\Gamma_{f,e}\cup\Gamma_{f,n}\cup\Gamma}\underline v\cdot(\underline{\underline T}(\underline u_f,p_f)\underline n_f)=&\displaystyle\int_{\Omega_f}\underline f\cdot\underline v&\text{ in }\Omega_f\\ \displaystyle a_p(\varphi,\psi)-\int_{\Gamma_{p,e}\cup\Gamma_{p,n}\cup\Gamma}\underline{\underline\kappa}\psi\nabla\varphi\cdot\underline n_p=&0&\text{ in }\Omega_p\\ \underline u_f=&\underline u_{f,e}&\text{ on }\Gamma_{f,e}\\ \underline u_p\cdot\underline n_p=&u_{p,n}&\text{ on }\Gamma_{p,n}\\ \varphi=&\varphi_e&\text{ on }\Gamma_{p,e}\\ (-\underline{\underline\kappa}\nabla\varphi)\cdot\underline n_p=&u_{p,n}&\text{ on }\Gamma_{p,n}\\ \underline u_p\cdot\underline n_f=&\underline u_f\cdot\underline n_f&\text{ on }\Gamma\\ \underline\tau_i\cdot(\underline{\underline T}(\underline u_f,p_f)\underline n_f)=&\displaystyle\frac1\alpha(\underline u_f-\underline u_p)\cdot\underline\tau_i&\text{ on }\Gamma\text{ for }i=1,..,d-1\\ -\underline n_f\cdot(\underline{\underline T}(\underline u_f,p_f)\underline n_f)=&g\varphi&\text{ on }\Gamma \end{array}\right.\]

Adding up these equations yields the weak formulation of this coupled Stokes-Darcy problem :

\[ \left\{\begin{array}{l} \text{find }\underline u_f\in H_f\text{, }p_f\in L^2(\Omega_f)\text{ and }\varphi\in H_p\text{ s.t. for all }\underline v\in H_f\text{, }q\in L^2(\Omega_f)\text{ and }\psi\in H_p :\\ \displaystyle a_f(\underline u_f,\underline v)+b_f(\underline v,p_f)+ga_p(\varphi,\psi)+\int_\Gamma g\varphi(\underline v\cdot\underline n_f)-\int_\Gamma g\psi(\underline u_f\cdot\underline n_f)+\int_\Gamma\sum_{i=1}^{d-1}\frac1\alpha(\underline u_f\cdot\underline\tau_i)(\underline v\cdot\underline\tau_i)=\int_{\Omega_f}\underline f\cdot\underline v - \int_{\Gamma_{f,n}}\underline v\cdot\underline T_{f,n} - \int_{\Gamma_{p,n}}\psi u_{p,n}\\ b_f(\underline u_f,q)=0 \end{array}\right.\]