Multiscale finite element method for Stokes-Darcy model

Let us consider the following mixed model for coupling a fluid flow and a porous media flow in a bounded domain $\Omega\subset\mathbf{R}^{d},d=2,3$. Here $\Omega=\Omega_{f}\cup\Gamma\cup\Omega_{p}$, where $\Omega_{f}$ and $\Omega_{p}$ are two disjoint, connected and bounded domains occupied by fluid flow and porous media flow and $\Gamma=\bar{\Omega}_{f}\cap\bar{\Omega}_{p}$ is the interface. For simplicity, we assume $\partial\Omega_{p}$ and $\partial\Omega_{f}$ are smooth enough in the rest of this paper. We denote $\Gamma_{f}=\partial\Omega_{f}\cap\partial\Omega,\Gamma_{p}=\partial\Omega_{p}\cap\partial\Omega$ and we also denote by $\mathbf{n}_{p}$ and $\mathbf{n}_{f}$ the unit outward normal vectors on $\partial\Omega_{p}$ and $\partial\Omega_{f}$, respectively. See Fig. 1 for a sketch.

The fluid motion in the fluid region $\Omega_{f}$ is governed by the Stokes equations

where

are the stress tensor and the deformation rate tensor, $v>0$ is the kinetic viscosity and $\mathbf{g}_{f}$ is the external force.

We assume that the porous region possesses multiple scales. The fluid motion in the porous medium region $\Omega_{p}$ is governed by

where $\mathbb{K}_{\epsilon}$ denotes the hydraulic conductivity in $\Omega_{p}$, which is a positive symmetric tensor and includs multiscale information where $\epsilon$ is a small parameter and $g_{p}$ is a source term. The first equation is the saturated flow model and the second equation is the Darcy’s law. Here $\phi_{p}=z+\frac{p_{p}}{\rho g}$ is the piezometric (hydraulic) head, where $p_{p}$ represents the dynamic pressure, $z$ the height from a reference level, $\rho$ the density and $g$ the gravitational constant, and $\mathbf{u}_{d}$ is the flow velocity in the porous medium which is proportional to the gradient of $\phi_{p}$, namely, the Darcy’s law.

Combining the two equations in (2.2), we get the equation for the piezometric head, which we will refer to simply as the Darcy equation:

Eqs. (2.1) and (2.3) are completed and coupled together by the following boundary conditions:

and the interface conditions on $\Gamma$ :

where $\bm{\tau}_{i},i=1,\cdots,d-1$ is an orthonormal basis of the tangential space on $\Gamma$ and $g$ the gravitational acceleration and we assumed that $g=1$ in the following. $\alpha$ is an experimentally determined parameter and $\Pi$ represents the permeability, which has the following relation with the hydraulic conductivity, $\mathbb{K}_{\epsilon}=\frac{\Pi g}{v}$. The first interface condition (2.4) describes the mass conservation and the second equation (2.5) represents the balance of momentum. The third interface condtion (2.6) is called Beavers-Joseph-Saffman condition, which means the tangential components of the normal stress force is proportional to the tangential components of the fluid velocity[5].

Furthermore, we assume $\mathbb{K}_{\epsilon}$ is symmetric, perodic, and uniformly elliptic. There exist two constants $\lambda_{\max}>0,\lambda_{\min}>0$ such that

Also, we assume

For simplicity, we consider 2D problems through the full text. We reserve $\Omega$ for a domain (bounded and open set) with Lipschitz boundary and $d$ for spacial dimension $(d=2)$. The Einstein summation convention is adopted, means summing repeated indexes from 1 to $d$. The Sobolev spaces $W^{k,p}$ and $H^{k}$ are defined as usual (see e.g., [6]) and we abbreviate the norm Sobolev space $H^{k}(D)$ as $\|\cdot\|_{k,D}$.

Later on we need to introduce some Hilbert spaces

The space $L^{2}(D)$, where $D=\Omega_{\mathrm{f}}$ or $\Omega_{\mathrm{p}}$, is equipped with the usual $L^{2}$-scalar product $(\cdot,\cdot)$ and $L^{2}$-norm $\|\cdot\|_{L^{2}(D)}$. The spaces $H_{\mathrm{f}}$ and $H_{\mathrm{p}}$ are equipped with the following norms:

Let us denote

Hence, we use the notational convention that $\underline{\mathbf{u}}=(\mathbf{u}_{f},\phi_{p})$ and $\underline{\mathbf{v}}=(\mathbf{v}_{f},\psi_{p})$. They all belong to $\mathbf{U}$: for $\mathbf{g_{f}}\in\mathbf{L}^{2}(\Omega_{f})$ and $g_{p}\in L^{2}(\Omega_{p})$, find $(\underline{\mathbf{u}},p_{f})\in\mathbf{U}\times Q_{f}$ such that $\forall(\underline{\mathbf{v}},q_{f})\in\mathbf{U}\times Q_{f}$ and $\mathbf{U}^{\prime}$ is the dual space of $\mathbf{U},P_{\tau}(\cdot)$ is the projection onto the local tangential plane that can be explicitly expressed as $P_{\tau}\left(\mathbf{v}_{f}\right)=\mathbf{v}_{f}-\left(\mathbf{v}_{f}\cdot\mathbf{n}_{f}\right)\mathbf{n}_{f}$

where

We know that there exists a positive constant $\beta>0$ such that the following Ladyzhenskaya-Babuš kaBrezzi (LBB) condition holds:

We give the finite element approximation of this model. For any given small parameter $h>0$, we construct the regular triangulations $\mathcal{T}_{h},\mathcal{T}_{fh}$ and $\mathcal{T}_{ph}$ of $\Omega,\Omega_{f}$ and $\Omega_{p}$.

Let $\textbf{V}_{fh}\subset\textbf{V}_{f}$ ,$X_{ph}\subset X_{p}$,and $Q_{fh}\subset Q_{f}$ be finite element spaces such that the space pair $(\textbf{V}_{fh},Q_{fh})$ satisfies the discrete LBB condition : there exists a constant $\beta>0$,independent of mesh ,such that

Here we choose MINI finite element pair for $(\textbf{V}_{fh},Q_{fh})$ and multiscale finite element for $V_{ph}$. Define

We define Lagrange element space $\mathcal{P}:=\left\{u\in C(\bar{\Omega}),\forall T\in\mathcal{T}_{h},\left.u\right|_{T}\in\mathcal{P}_{1}(T)\right\}$ and consider a triangulation $\mathcal{T}_{ph}$ and base functions $\left(\psi^{i}\right)_{1\leq i\leq Nb_{\text{vertex }}}\in\mathcal{P}$ which consists of the globally continuous on $\mathcal{T}_{ph}$ and affine on each triangle $K\subset\mathcal{T}_{ph}$ functions which satisfy: $\psi^{i}\left(x_{j}\right)=\delta_{ij}$,where $\left(x_{j}\right)_{1\leq j\leq Nb_{\text{vertex }}}$ are the vertices of the coarse mesh ($Nb_{\text{vertex }}$ the set of interior vertices of the coarse mesh) and $\delta_{ij}$ is the Kronecker symbol. The multiscale finite element basis functions $\left(\eta_{i}^{MsFEM}\right)_{1\leq i\leq Nb_{\text{vertex }}}$ which take into account the multiscale of the coefficients. More precisely we compute $\eta_{i,K}^{MsFEM}$ , such that for any triangle $K\subset\mathcal{T}_{ph}$,

In practice, we do not have access to $\eta_{i,K}^{MsFEM}$. We build $\eta_{i,K}^{MsFEM,h_{finer}}$, an approximation of $\eta_{i,K}^{MsFEM}$ on a finer embedded grid of mesh size $h_{finer}$, with $\mathcal{P}_{1}$ FE. See Fig.2.

Therefore, we can get the multiscale finite element basis function $\{\eta_{i}^{MsFEM}=\cup_{K\in\mathcal{K}^{h}}\eta_{i,K}^{MsFEM}\}$ and a simple observation tells that $\eta_{i}^{MsFEM}$ is locally support.

Then we can introduce

Therefore, the coupled finite element Galerkin approximation of reads:

Coupled Finite Element Scheme: find $\underline{\mathbf{u}}_{h}=\left(\mathbf{u}_{fh},\phi_{ph}\right)\in\mathbf{U}_{h},p_{fh}\in Q_{fh}$ such that for any $\underline{\mathbf{v}}_{h}=\left(\mathbf{v}_{fh},\phi_{ph}\right)\in\mathbf{U}_{h}$ and $q_{fh}\in Q_{fh}$

The Darcy region exhibits multiscale characteristics. In order to better elucidate our model, we introduce the model from [22]

For simplicity, we assume the source term $f$ belongs to $L^{2}(\Omega)$, and the boundary term $g\in L^{2}\left(\Gamma_{N}\right)$.

Then, we introduce the multiscale expansion technique in deriving homogenized equations. We look for $\phi_{\epsilon}(x)$ in the form of asymptotic expansion

where the functions $\phi_{j}(x,y)$ are periodic in $y$ with period 1 .

Here, we referenced some definitions and theorems from [22].

We need an interpolation operator. If the triangulation $\mathcal{T}_{h}$ is regular, then there exist an interpolation operator $\mathcal{I}:H^{1}(\Omega)\mapsto\mathcal{P}$ , $\phi_{\epsilon,\mathcal{I}}\in X_{ph}$.

In the offline phase, we solve the multiscale basis functions in parallel. This approach is adopted to enhance the efficiency of basis function generation.

• •

  The grid is partitioned to obtain the required coarse mesh $K\in\mathcal{T}_{h}$.

• •

  For every $K\in\mathcal{T}_{h}$, we build $\mathcal{T}h_{finer}^{K}$ and parallelly solve $-\nabla\cdot\left(\mathbb{K}^{\epsilon}(x/\epsilon)\right)\nabla\eta_{i}=0$ in $K,\eta_{i}=\varphi_{i}$ on $\partial K$ with $\varphi_{i}$ the standard $\mathrm{P}1$ basis function and store $\eta_{i}$. See Fig.4.

  

• •

  Compute and store $A_{i,j}^{\text{loc }}=\int_{K}\left(\mathbb{K}^{\epsilon}\nabla\eta_{i}\right)\left(\nabla\eta_{j}\right)$ and $B_{i}=B_{i}+\int_{K}\eta_{i}\times f_{p}$ with $A^{loc}$ stifness matrix and $B$ the generic RHS.

• •

  Output: Assemble and store $A$ the stiffness matrix associted with the new basis $\eta_{H}.$

In the online phase, we employ the Robin-Robin Stokes-Darcy algorithm. In this process, we utilize the multiscale basis functions generated in the offline stage. This algorithm is a decoupled algorithm, and here are the main steps:

• •

  Given the initial condtion $\mathbf{u}_{f}^{0},p_{f}^{0},\varphi_{p}^{0}$, $\epsilon_{iter}$,$\epsilon_{error}$,$\gamma_{f},\gamma_{p}$.

• •

  For $k=1,2,\ldots$, independently solve the Stokes and Darcy systems with Robin boundary conditions. More precisely, find $(\mathbf{u}_{fh}^{k+1},p_{fh}^{k+1})$ such that for all $\mathbf{v}_{fh}\in\mathbf{V}_{fh}$ and $q_{fh}\in Q_{fh}$ it is

  $$a_{f}\left(\mathbf{u}_{fh}^{k},\mathbf{v}_{fh}\right)+b_{f}\left(\mathbf{v}_{fh},p_{fh}^{k}\right)+\gamma_{f}\left(\mathbf{u}_{fh}^{k}\cdot\mathbf{n}_{f},\mathbf{v}_{fh}\cdot\mathbf{n}_{f}\right)+\alpha\left(P_{\tau}\mathbf{u}_{fh}^{k},P_{\tau}\mathbf{v}_{fh}\right)=\left(\eta_{f}^{k},\mathbf{v}_{fh}\cdot\mathbf{n}_{f}\right)+\left(\mathbf{g}_{f},\mathbf{v}_{fh}\right),$$

  $$b_{f}\left(\mathbf{u}_{fh}^{k},q_{fh}\right)=0.$$

• •

  Find $\varphi_{ph}^{k+1}$ such that for all $\psi_{ph}\in Q_{ph}$ it is

  $$\gamma_{p}a_{p}\left(\phi_{ph}^{k},\psi_{ph}\right)+\left(g\phi_{ph}^{k},\psi_{ph}\right)=\left(\eta_{p}^{k},\psi_{ph}\right)+\left(g_{p},\psi_{ph}\right).$$

• •

  Update $\eta_{p}^{k+1}$ and $\eta_{f}^{k+1}$

  $$\eta_{f}^{k+1}=a\eta_{p}^{k}+bg\phi_{ph}^{k},\text{where}\quad a=\frac{\gamma_{f}}{\gamma_{p}},b=-1-a$$

  $$\eta_{p}^{k+1}=c\eta_{f}^{k}+d\mathbf{u}_{fh}^{k}\cdot\mathbf{n}_{f},\text{where}\quad c=-1,d=\gamma_{f}+\gamma_{p}$$

• •

  Compute $\epsilon_{iter}=\|\mathbf{u}_{fh}^{k+1}-\mathbf{u}_{fh}^{k}\|^{2}+\|p_{fh}^{k+1}-p_{fh}^{k}\|^{2}+\|\varphi_{ph}^{k+1}-\varphi_{ph}^{k}\|^{2}$

• •

  If $\epsilon_{iter}>\epsilon_{error}$, then $k:=k+1$.

• •

  Output: get the soultion $\mathbf{u}_{fh}^{N},p_{fh}^{N},\varphi_{ph}^{N}$.

The offline phase involves a parallel method for generating multiscale basis functions. In the online phase, a Robin-Robin online Stokes-Darcy method is employed to obtain the solution. The algorithm offers advantages such as parallel construction of basis functions, which enhances reusability, and ensures high computational accuracy and efficiency. Additionally, the method is easily implementable for generalization and maintenance.