Multiscale dimension reduction for flow and transport problems in thin domain with reactive boundaries

Mathematical models in thin domains occur in many real-world applications, scientific and engineering problems. Fluid flow and transport in thin tube structures are widely used in biological applications, for example, to simulate blood flow in vessels [1, 2, 3]. In engineering problems, flow simulation is used to study fluid flow in complex pipe structures, for example in pipewise industrial installations, wells in oil and gas industry, heat exchangers, etc. In reservoir simulations, thin domains are related to the fractures that usually have complex geometries with very small thickness compared to typical reservoir sizes [4, 5, 6]. Such problems are often considered with complex interaction processes with surrounding media or walls. In many applications, these problems are transformed to reduced (e.g., one) dimensional problems via some type of apriori postulated models. Our goal is to present an alternative approach to analytical model reduction by using multiscale basis functions.

For the applicability of the convenient numerical methods for simulations of such problems, a very fine grid should be constructed to resolve the geometry’s complex structure on the grid level. Moreover, a very small domain thickness provides an additional complexity in the grid construction for thin and long domains. For a fast and accurate solution of the presented problem, a homogenization (upscaling) technique or multiscale models are used, which are based on constructing the lower dimensional coarse grid models. The asymptotic behavior of the solutions in thin domains is intensively studied. In [7], the authors consider the problem of complete asymptotic expansion for the time-dependent Poiseuille flow in a thin tube. In [8], the method of asymptotic partial decomposition of a domain (MAPDD) was presented to reduce the computational complexity of the numerical solution of such problems. This method combines the three-dimensional description in some neighborhoods of bifurcations and the one-dimensional description out of these small subdomains, and it prescribes some special junction conditions at the interface between these 3D and 1D submodels. Numerical results were presented in [2] for the method of asymptotic partial domain decomposition for thin tube structures with Newtonian and non-Newtonian flows in large systems of vessels. Our goal is to provide an alternative systematic approach to model reduction for thin domain problems that can add complexity via additional multiscale basis functions.

Model reduction techniques usually depend on a coarse grid approximation, which can be obtained by discretizing the problem on a coarse grid and choosing a suitable coarse-grid formulation of the problem. In the literature, several approaches have been developed to obtain the coarse-grid formulation for the problems in heterogeneous domains, including multiscale finite element method [9, 10, 11], mixed multiscale finite element method [12, 13], generalized multiscale finite element method [14, 15, 16, 17, 18], multiscale mortar mixed finite element method [19], multiscale finite volume method [20, 21], variational multiscale methods [22, 23, 24], and heterogeneous multiscale methods [25] and etc. The non - conforming multiscale method is considered for the solution of the Stokes flow problem in a heterogeneous domain in [26, 27]. In [28, 29], we presented the Generalized Multiscale Finite Element Method (GMsFEM) for the solution of the flow problems in perforated domains with continuous multiscale basis functions. GMsFEM shows a good accuracy for solving the nonlinear (non - Newtonian) fluid flow problems [30]. In [31], we presented the Discontinuous Galerkin Generalized Multiscale Finite Element Method (DG-GMsFEM) for the solving the two - dimensional problems in perforated domains.

In this paper, we consider flow and transport processes in thin structures with reactive boundaries. The mathematical model is described using Stokes equations and the unsteady convection-diffusion equation with non - homogeneous boundary conditions. Non-homogeneous boundary conditions occur in many applications. For example, the pore-scale modeling and simulation of reactive flow in porous media have many applications in many branches of science such as biology, physics, chemistry, geomechanics, and geology [32, 33, 34]. In order to handle the complex geometry of the walls and non-homogeneous boundary conditions on them, we present additional spectral multiscale basis functions. In this work, we continue developing the multiscale model reduction techniques for problems with multiscale features and developing the generalization of the techniques for DG-GMsFEM. In our previous work [35, 36], we considered elliptic problems in perforated domains and constructed additional basis functions to capture non-homogeneous boundary conditions on perforations. In [35], we proposed a non-local multi-continua (NLMC) method for Laplace, elasticity, and parabolic equations with non-homogeneous boundary conditions on perforations. In [36], we considered the Continuous Galerkin Generalized Multiscale Finite Element Method (CG-GMsFEM) for problems in perforated domains with non-homogeneous boundary conditions, where we constructed one additional basis function for local domains with perforations. Recently, we extended this technique for solving unsaturated flow problems in heterogeneous domains with rough boundary [37]. In this work, we consider DG-GMsFEM and present additional spectral basis functions for rough non-homogeneous boundaries in transport problems. The concept is based on the separation of the snapshots for each feature and shares a lot of similarities with multiscale methods for fractured media presented in our previous works [38, 39]. Our work is also motivated by a recently developed Edge GMsFEM, where multiscale basis functions are constructed for each coarse grid interface [40, 41].

In this work, we use the DG-GMsFEM for constructing multiscale basis functions for problems in thin domains. The goals of this paper are in constructing the general approach for problems in complex thin geometries with an accurate approximation of the velocity space and transport processes. We construct local multiscale basis functions to generate a lower-dimensional model on a coarse grid. In DG-GMsFEM for flow problem, we start with constructing the snapshot space that captures possible flows between coarse cell interfaces. After constructing the snapshot space, we perform a dimension reduction by a solution of the local spectral problems. For the pressure approximation, we use piecewise constant functions on the coarse grid. For the transport problem, we construct multiscale basis functions for each interface between coarse grid cells and present additional basis functions to capture non-homogeneous boundary conditions on walls. The presented snapshot spaces can accurately capture processes on the rough wall boundaries with non-homogeneous boundary conditions on them. We present the results of the numerical simulations for three test geometries for two and three-dimensional problems.

As we mentioned earlier, our goal is to investigate the use of multiscale and generalized multiscale methods for dimension reduction for problems in this domains. Many existing approaches propose analytical or semi-analytical reduced-dimensional problems, where the dimensions are determined apriori. Our idea is to use generalized multiscale method and identify the dimension. The proposed approach combined with aposteriori error estimate can further identify the dimension across thin layers and allow obtaining more accurate solution.

The paper is organized as follows. In Section 2, we describe the problem formulation and the fine-scale approximation. In Section 3, we present the multiscale method for flow and transport processes in thin structures with rough reactive boundaries. In Section 4, we present numerical results. The paper ends with a conclusion.

Let $\Omega$ be a thin domain with multiscale features, where the thickness is small compared to the domain length (See Figure 1 for an illustration). We will consider the following flow and transport equations in the thin domain $\Omega$

where $\mu$ is the fluid viscosity, $\rho$ is the fluid density, $D$ is the diffusion coefficient, $c$ is the concentration, $u$ and $p$ are the velocity and pressure.

The system (1) is equipped with following initial conditions

Moreover, We consider the following boundary conditions for the flow problem

and for the transport problem

where $n$ is the unit outward normal vector on $\partial\Omega$, $\mathcal{I}$ is the identity matrix, $\Gamma_{in}$ is the inflow boundary, $\Gamma_{out}$ is the outflow boundary, $\Gamma_{w}$ is the reactive boundary of the thin domain, $\Gamma_{w}\cup\Gamma_{in}\cup\Gamma_{out}=\partial\Omega$ (see Figure 1).

In order to solve the problem (1), we generate an unstructured grid (fine grid) and use a finite element approximation for the spatial discretization. Let $\mathcal{T}^{h}$ be a fine-grid partition of the domain $\Omega$ given by

where $N_{cell}^{h}$ is the number of fine grid cells. We use $\mathcal{E}^{h}$ to denote the set of facets in $\mathcal{T}^{h}$ with $\mathcal{E}^{h}=\mathcal{E}^{h}_{o}\cup\mathcal{E}^{h}_{b}$, where $\mathcal{E}^{h}_{o}$ is the set of interior facets and $\mathcal{E}^{h}_{b}$ is the set of boundary facets with $\mathcal{E}^{h}_{b}=\mathcal{E}^{h}_{b,in}\cup\mathcal{E}^{h}_{b,w}\cup\mathcal{E}^{h}_{b,out}$ (see Figure 2). We use the notations $K$ and $E$ to denote a generic cell and facet in the fine grid $\mathcal{T}^{h}$ (see Figure 2). We define the jump $[u]$ and the average $\{u\}$ of a function $u$ on interior facet by

where $u_{+}=u|_{K^{+}}$, $u_{-}=u|_{K^{-}}$, $K^{+}$ and $K^{-}$ are the two cells sharing the facet $E$. Note that, we define $[u]=u|_{E}$ and $\{u\}=u|_{E}$ for $E\in\mathcal{E}^{h}_{b}$.

We define the fine-scale velocity space

which contains functions that are piecewise linear in each fine-grid element $K$ and are discontinuous across coarse grid edges. For the pressure, we use the space of piecewise constant functions

The fine-scale space for concentration is the following

Using these spaces, we have the following IPDG variational formulation with implicit time approximation for the approximation of (1):

• •

  Flow problem: find $(u_{h},p_{h})\in V_{h}\times Q_{h}$ such that

  $$\begin{split}\frac{1}{\tau}m^{u}(u_{h}-\check{u}_{h},v)+a^{u}_{\text{DG}}(u_{h},v)+b_{\text{DG}}(p_{h},v)&=l^{u}(v),\quad\forall v\in V_{h},\\ b_{\text{DG}}(u_{h},q)&=l^{p}(q),\quad\forall q\in Q_{h},\end{split}$$

  (2)

  where

  $$\begin{split}a^{u}_{\text{DG}}(u,v)&=\sum_{K\in\mathcal{T}^{h}}\int_{K}\mu\nabla u\cdot\nabla v\,dx\\ &-\sum_{E\in\mathcal{E}^{h}/\mathcal{E}^{h}_{b,out}}\int_{E}\Big{(}\{\mu\nabla u\cdot n\}\cdot[v]+\{\mu\nabla v\cdot n\}\cdot[u]-\frac{\gamma_{u}}{h}\{\mu\}[u]\cdot[v]\Big{)}\,ds,\\ m^{u}(u,v)&=\sum_{K\in\mathcal{T}^{h}}\int_{K}\rho\,u\,v\,dx,\\ b_{\text{DG}}(u,p)&=-\sum_{K\in\mathcal{T}^{h}}\int_{K}p\,\nabla u\,dx+\sum_{E\in\mathcal{E}^{h}/\mathcal{E}^{h}_{b,out}}\int_{E}p\,[u]\cdot n\,ds,\\ l^{u}(v)&=\sum_{E\in\mathcal{E}^{h}_{b,in}}\int_{E}\left(\frac{\gamma_{u}}{h}\,\mu v-(\mu\nabla v\cdot n)\right)\cdot g\,ds,\\ l^{p}(q)&=\sum_{E\in\mathcal{E}^{h}_{b,in}}\int_{E}(g\cdot n)\,q\,dx,\end{split}$$

  and $\gamma_{u}$ is the penalty perm.

• •

  Transport problem: find $c_{h}\in P_{h}$ such that

  $$\begin{split}\frac{1}{\tau}m^{c}(c_{h}-\check{c}_{h},r)+c^{c}_{\text{DG}}(c_{h},r)+a^{c}_{\text{DG}}(c_{h},r)&=l^{c}(r),\quad\forall r\in P_{h},\end{split}$$

  (3)

  where

  $$\begin{split}a^{c}_{DG}(c,r)&=\sum_{K\in\mathcal{T}^{h}}\int_{K}D\nabla c\cdot\nabla r\,dx\\ &-\sum_{E\in\mathcal{E}^{h}_{0}\cup\mathcal{E}^{h}_{b,in}}\int_{E}\Big{(}\{D\nabla c\cdot n\}[r]+\{D\nabla r\cdot n\}[c]-\frac{\gamma_{c}}{h}\{D\}[c][r]\Big{)}\,ds\\ &+\sum_{E\in\mathcal{E}^{h}_{b,w}}\int_{E}\alpha\,c\,r\,ds,\\ m^{c}(c,r)&=\sum_{K\in\mathcal{T}^{h}}\int_{K}c\,r\,dx,\\ c^{c}_{DG}(c,r)&=\sum_{K\in\mathcal{T}^{h}}\int_{K}(u_{h}\,c)\cdot\nabla r\,dx+\sum_{E\in\mathcal{E}^{h}_{0}\cup\mathcal{E}^{h}_{b,out}}\int_{E}(\tilde{u}_{+}c_{+}-\tilde{u}_{-}c_{-})\,[r]\,ds,\\ l^{c}(r)&=\sum_{E\in\mathcal{E}^{h}_{b,in}}\int_{E}\left(\frac{\gamma_{c}}{h}\,D\,r-\{D\nabla r\cdot n\}\right)\,c_{in}\,dx+\sum_{E\in\mathcal{E}^{h}_{w}}\int_{E}\alpha\,c_{w}\,r\,ds,\end{split}$$

  and $\gamma_{c}$ is the penalty perm and $\tilde{u}=(u_{h}\cdot n+|u_{h}\cdot n|)/2$.

In the above systems (2) and (3), $\check{u}_{h}$ and $\check{c}_{h}$ are the solutions from the previous time step and $\tau$ is the given time step size.

We can write the above discrete systems (2) and (3) as follows.

• •

  Flow problem:

  $$\begin{split}\frac{1}{\tau}\begin{pmatrix}M^{u}_{h}&0\\ 0&0\end{pmatrix}\begin{pmatrix}u_{h}-\check{u}_{h}\\ p_{h}-\check{p}_{h}\end{pmatrix}+\begin{pmatrix}A^{u}_{h}&B_{h}^{T}\\ B_{h}&0\end{pmatrix}\begin{pmatrix}u_{h}\\ p_{h}\end{pmatrix}=\begin{pmatrix}F_{h}^{u}\\ F_{h}^{p}\end{pmatrix}.\end{split}$$

  (4)

• •

  Transport problem:

  $$\frac{1}{\tau}M^{c}_{h}(c_{h}-\check{c}_{h})+(A^{c}_{h}+C^{c}_{h}(u_{h}))c_{h}=F^{c}_{h}.$$

  (5)

  Here, the matrix $C^{c}_{h}(u_{h})$ depends on the function $u_{h}$.

In the next section, we will present the proposed multiscale method for the solution of the flow and transport problems that used to reduce the system size. In the multiscale method, we solve problems in local domains with various boundary conditions to form a snapshot space and use a spectral problem in the snapshot space to perform the required dimension reduction.

Let $\mathcal{T}^{H}$ be a coarse-grid partition of the domain $\Omega$ with mesh size $H$ (see Figure 3)

where $N_{cell}^{H}$ is the number of coarse grid cells (local domains). We use $\mathcal{E}^{H}$ to denote the set of facets in $\mathcal{T}^{H}$ with $\mathcal{E}^{H}=\mathcal{E}^{H}_{o}\cup\mathcal{E}^{H}_{b}$. For the construction of the coarse grid approximation, we use the Discontinuous Galerkin Generalized Multiscale Finite Element Method (DG-GMsFEM). In DG-GMsFEM, the multiscale basis functions are supported in each coarse cell $K_{i}$[28, 29, 31]. We define $V_{H}$ as the multiscale velocity space and $P_{H}$ as the multiscale space for concentration. For the pressure approximation, we use the piecewise constant function space $Q_{H}$ over the coarse cells.

We denote the multiscale spaces for concentration and velocity as

respectively, where $N_{u}=\text{dim}(V_{H})$ is the number of basis functions for velocity, $N_{c}=\text{dim}(P_{H})$ is the number of basis functions for concentration. For the pressure, we use the space of piecewise constant functions over the coarse grid cells

where $N_{p}=\text{dim}(Q_{H})$ is equal to the number of coarse grid cells ($N_{p}=N_{cell}^{H}$).

For the coarse grid approximation, we have the following variational formulation:

• •

  Flow problem: find $(u_{H},p_{H})\in V_{H}\times Q_{H}$ such that

  $$\begin{split}\frac{1}{\tau}m^{u}(u_{H}-\check{u}_{H},v)+a^{u}_{\text{DG}}(u_{H},v)+b_{\text{DG}}(p_{H},v)&=l^{u}(v)\quad\forall v\in V_{H},\\ b_{\text{DG}}(u_{H},q)&=l^{p}(q),\quad\forall q\in Q_{H}.\end{split}$$

  (6)

• •

  Transport problem: find $c_{H}\in P_{H}$ such that

  $$\begin{split}\frac{1}{\tau}m^{c}(c_{H}-\check{c}_{H},r)+a^{c}_{\text{DG}}(c_{H},r)+c^{c}_{\text{DG}}(c_{H},r)&=l^{c}(r),\quad\forall r\in P_{H}.\end{split}$$

  (7)

Next, we consider the construction of the multiscale basis functions for velocity and concentration.

In this section, we will present some numerical results. We will use the following three computational domains (Figure 13) to demonstrate the performance of our method:

• Geometry 1 with fine grid that contains 17350 cells. Coarse grid contains 10 local domains.

• Geometry 2 with fine grid that contains 18021 cells. Coarse grid contains 20 local domains.

• Geometry 3 with fine grid that contains 15094 cells. Coarse grid contains 10 local domains.

In order to construct structured coarse grids, we explicitly add lines between local domains (coarse grid cells) in geometry construction. For unstructured coarse grid, local domains have a rough interface between local domains. Computational domains with fine and coarse grids are presented in Figure 13 for Geometry 1, 2, and 3. We use Gmsh to construct computational geometries and unstructured fine grids [46]. The numerical implementation is based on the FEniCS library [47].

To investigate the presented multiscale method for solving problems in thin domains with non-homogeneous boundary conditions, we consider the following test:

• •

  Test 1 for different geometries of computational domain. We consider Geometries 1, 2 and 3 (see Figure 13). We consider problems with non-homogeneous Robin boundary conditions for concentration.

• •

  Test 2 for different boundary conditions and diffusion coefficients ($D=0.01$, $0.1$ and $1$). We consider Dirichlet and Neumann non-homogeneous boundary conditions on the wall boundary $\Gamma_{w}$. In this test, we also investigate different types of multiscale basis functions in detail.

• •

  Test 3 for structured and unstructured coarse grids. We investigate the performance of the multiscale method with a rough interface between local domains (unstructured coarse grid). We also investigate the influence of the multiscale velocity accuracy to the concentration errors.

For flow problem, we set $\mu=1$, $\rho=1$ and $u_{0}=0$. We set $g=(\tilde{g},0)$ for Geometry 1,2 and $g=(0,0,\tilde{g})$ for Geometry 3 as inflow boundary condition $\Gamma_{in}$, where $\tilde{g}(x)=u_{in}n^{-1}(n+2)(1-(r/r_{max})^{n})$ with $n=2$ and $u_{in}=1$ [48, 3]. Here $r$ is the distance to center point $x_{0}$ and $r_{max}$ is the radius of left boundary $\Gamma_{in}$, where $x_{0}=(0,0.05)$, $r_{max}=0.05$ for Geometry 1, $x_{0}=(0,1)$, $r_{max}=0.025$ for Geometry 2, and $x_{0}=(0,0,1)$, $r_{max}=0.035$ for Geometry 3.

We calculate relative errors in $L^{2}$ norm in percentage