@@ -8,7 +8,7 @@ This extends the updated performance study from \cite{klinkovsky:2022meshes} wit
The section is organized as follows.
First, we introduce the general formulation of the PDE system solved by \emph{NumDwarf} in \cref{sec:mhfem:formulation} and summarize the numerical scheme in \cref{sec:mhfem:numerical scheme}.
In\cref{sec:mhfem:mcwhdd}, we present the results of a benchmark problem that was used to verify the numerical scheme and implementation by means of comparison against known semi-analytical solutions.
\Cref{sec:mhfem:two-phase flow} summarizes the mathematical model of incompressible two-phase flow in porous media and\cref{sec:mhfem:mcwhdd} presents the results of a benchmark problem that was used to verify the numerical scheme and implementation by means of comparison against known semi-analytical solutions.
The same benchmark problem was used to also evaluate the performance of the solver using various computational parameters.
The benchmarking methodology is described in \cref{sec:mhfem:benchmarking methodology} and the benchmark results are presented in \cref{sec:mhfem:computational results}.
@@ -60,21 +60,230 @@ The matrix of each linear system is generally non-symmetric due to the upwind st
In this work, we use conforming unstructured homogeneous meshes denoted as $\mathcal K_h$ for the discretization of the domain and the lowest-order Raviart--Thomas--Nédélec function space $\mathbf{RTN}_0(\mathcal K_h)$.
The set of faces of the mesh $\mathcal K_h$ will be denoted as $\mathcal E_h$ and the set of faces of an individual element $K \in\mathcal K_h$ will be denoted as $\mathcal E_K$.
\subsection{Two-phase flow in porous media}
\label{sec:mhfem:two-phase flow}
This section summarizes the mathematical model of incompressible two-phase flow in porous media.
See \cite{fucik:thesis,bastian:thesis} for details on the theory behind the model.
The two-phase incompressible flow in porous media can be described macroscopically by the system of partial differential equations
where $\Phi$~[-] is the material porosity, $K$~[\si{\metre\squared}] is the material permeability, and $\vec g$~[\si{\metre\per\second\squared}] is the gravitational acceleration vector.
Symbols related to the phases in the system are denoted by subscripts $\alpha\in\{w,n\}$.
In particular, the quantities $\rho_\alpha$~[\si{\kilogram\per\metre\cubed}], $S_\alpha$~[-], $\lambda_\alpha$~[\si{\metre\second\per\kilogram}], and $p_\alpha$~[\si{\pascal}] represent the $\alpha$-phase density, saturation, mobility, and pressure, respectively.
The densities $\rho_\alpha$ are assumed to be constant, the mobilities are defined as $\lambda_\alpha= k_{r,\alpha}(S_\alpha)/\mu_\alpha$, where $k_{r,\alpha}(S_\alpha)$~[-] is the relative $\alpha$-phase permeability given by a specific empirical model and $\mu_\alpha$~[\si{\kilogram\per\metre\per\second}] is the $\alpha$-phase dynamic viscosity, and the quantities $S_\alpha$, $p_\alpha$ for $\alpha\in\{w,n\}$ are the four unknown variables in \cref{eq:mhfem:two-phase flow}.
The model is closed by algebraic constraints for the phase pressures and saturations
\begin{subequations}
\label{eq:mhfem:algebraic constraints}
\begin{align}
p_n - p_w &= p_c , \\
S_w + S_n &= 1 ,
\end{align}
\end{subequations}
where $p_c$~[\si{\pascal}] is the capillary pressure that is related to the wetting phase saturation $S_w$ according to $p_c = p_c(S_w)$, which is a formula supplied by a specific empirical model.
The function $p_c(S_w)$ must be invertible to allow the calculation of saturation from given capillary pressure.
Hence, multiple formulations of the problem are possible by selecting different pairs of primary unknown variables, such as $\{p_w,p_n\}$, $\{p_c,p_n\}$, or $\{p_w,S_w\}$.
The empirical model for capillary pressure used in this work is the Brooks--Corey model \cite{brooks:1964}
where $p_\ast$~[\si{\pascal}] and $\lambda_\ast$~[-] are model parameters for given porous material, and $S_{w,e}$ is the effective saturation of the wetting phase. The effective saturation of the $\alpha$-phase is defined by
The \cref{eq:mhfem:two-phase flow} governing the two-phase incompressible flow in porous media can be transformed to the form of \cref{eq:NumDwarf} by taking $n =2$ and setting the general coefficients appropriately for the selected pair of primary unknown variables.
Note that formulations using $S_w$ or $S_n$ as a primary unknown are not admissible for the numerical scheme, because phase saturations may become discontinuous across material interfaces in heterogeneous porous media.
The $\{p_w,p_n\}$ formulation is used in this work, i.e., $\vec Z =[p_w, p_n]^T$, and the coefficients in \cref{eq:NumDwarf} are set as follows:
In this section, we present the results of the generalized McWhorter--Sunada problem, a special case of two-phase flow in porous media with known semi-analytical solution \cite{mcwhorter1990exact,fucik2016multidimensional} which was used in \cite{fucik:2019NumDwarf} to verify the numerical scheme.
In this section, we present the results of the generalized McWhorter--Sunada problem, a special case of incompressible two-phase flow in porous media described by \cref{eq:mhfem:two-phase flow,eq:mhfem:algebraic constraints} with known semi-analytical solution \cite{mcwhorter1990exact,fucik2016multidimensional}.
The geometrical configuration of the problem and boundary conditions are shown in \cref{fig:mcwhdd:domain}.
The computational domain $\Omega\subset\mathbb R^D$ is a unit square in the two-dimensional case and a unit cube in the three-dimensional case.
A homogeneous porous medium is assumed in the domain and the effect of gravity is neglected.
Initially, at time $t =0$, the entire domain is filled with water at the reference pressure $p_w = p_{\mathrm{ref}}$ and $p_n$ is calculated from $p_w$ and the constant initial water saturation $S_{w,\mathrm{ini}}$.
The non-aqueous phase liquid (NAPL) is being injected into the domain since $t =0$ until the final time $t_{\max}$ through a point source located in the origin $\vec x =\vec0$.
The volumetric rate of injection, denoted as $Q_0$~[\si{\metre\tothe{D}\per\second}], is prescribed by
\begin{equation}
Q_0 = Q_0(t) = A_D t^{\frac{D-2}{2}},
\end{equation}
where $D \in\{2,3\}$ denotes the spatial dimension and $A_D$~[\si{\metre\tothe{D}\second\tothe{-D/2}}] is an injection rate parameter.
\inline{problem description}
As noted in \cite{fucik:2019NumDwarf,klinkovsky:2017thesis}, the boundary condition involving the point source injection cannot be treated exactly in the numerical scheme and it must be approximated using boundary faces of the numerical mesh adjacent to the origin $\vec x =\vec0$.
The part of the domain boundary $\Gamma=\partial\Omega$ that is used for the approximation of the injection is illustrated in \cref{fig:mcwhdd:Gamma_in} and denoted as $\Gamma_{\mathrm{in}}$.
The injection rate $Q_0$ is enforced on $\Gamma_{\mathrm{in}}$ for all $t \in[0, t_{\max}]$ as
\begin{equation}
% \int\limits_{\Gamma_{\mathrm{in}}} \left( - \lambda_n \tensor K (\nabla p_n - \rho_n \vec g) \right) \cdot \vec n = - Q_0(t),
\int\limits_{\Gamma_{\mathrm{in}}}\lambda_n K \nabla p_n \cdot\vec n = Q_0(t),
\end{equation}
where $\vec n$ is the outward unit normal vector on $\Gamma_{\mathrm{in}}$.
The boundary condition for the wetting phase on $\Gamma_{\mathrm{in}}$ remains the same, i.e., $\nabla p_w \cdot\vec n =0$.
Injection rate parameter -- 2D case &$A_2$& e-5 &\si{\metre\squared\per\second}\\
\phantom{Injection rate parameter} -- 3D case &$A_3$& e-7 &\si{\metre\cubed\second\tothe{-3/2}}\\
Final time &$t_{\max}$& 2e4 &\si{\second}\\\bottomrule
\end{tabular}
\caption{Parameters of the mathematical model used in the generalized McWhorter--Sunada problem.}
\label{tab:mcwhdd:parameters}
\end{table}
\subsubsection{Verification results}
For the purpose of this thesis, the computational domain, boundary conditions, and physical parameters were chosen identically to \cite{fucik:2019NumDwarf}.
The simulations were computed on the meshes listed in \cref{tab:meshes} which are the same as the triangular and tetrahedral meshes used in \cite{fucik:2019NumDwarf}.
All simulations were computed in \ic{double} precision.
For each mesh $\mathcal K_h$, we define the error $E_{h,S_n}= S_{n,\mathrm{ex}}(t_{\max})- S_{n,h}(t_{\max})$ in terms of the injected phase saturation $S_n$ as the difference between the semi-analytical solution $S_{n,\mathrm{ex}}$ and numerical solution $S_{n,h}$ on the mesh $\mathcal K_h$ at time $t_{\max}$.
The $L_p$ norm in $\Omega\subset\mathbb R^D$ of the error $\norm{E_{h,S_n}}_p$ is evaluated for $p \in\{1,2\}$ as
\begin{equation}
\label{eq:mhfem:norm additivity}
\norm{E_{h,S_n}}_p = \left( \int\limits_\Omega\abs{E_{h,S_n}(\vec x)}^p \dif\vec x \right)^\frac1p
The integrals over mesh cells $K \in\mathcal K_h$ are first transformed to integrals over the reference cell, and then the Fubini theorem is applied to obtain one-dimensional integrals that are computed numerically using the Gauss--Kronrod quadrature (G7--K15) with adaptive interval subdivision, which is implemented in the \CPPeleven library \ic{qdt}~\cite{munoz:2014,qdt}.
\inline{define $E_{h,S_n}$ and $eoc_{S_n,p}$}
\inline{describe the EOC tables}
Given two approximate solutions $S_{n,h_1}$ and $S_{n,h_2}$ on meshes $\mathcal K_{h_1}$ and $\mathcal K_{h_2}$, respectively, the order of convergence of the numerical scheme with respect to $S_n$ in the $L_p$ norm can be approximated as
where the acronym $EOC$ stands for experimental order of convergence.
In order to investigate the convergence and accuracy of the numerical scheme using the experimental order of convergence, several simulations were computed on a series of triangular and tetrahedral meshes.
The meshes are the same as those used in \cite{fucik:2019NumDwarf} and their properties are listed in \cref{tab:meshes}.
All simulations were computed in \ic{double} precision.
\Cref{tab:mhfem:EOC 2D,tab:mhfem:EOC 3D} show the resulting $L_1$ and $L_2$ norms of the error, the corresponding experimental orders of convergence, as well as the mesh size parameters $h$ and constant time steps $\tau$ that were used in the simulations.
The results of the numerical analysis indicate that the numerical scheme converges with the first order of accuracy in both dimensions and both $L_1$ and $L_2$ norms.
@@ -144,7 +144,7 @@ The diffusive flux is given by the Fick's law \cite{bird:2002} as
\begin{equation}\label{eq:transport:Fick law}
\vec J_\alpha = - \rho D \nabla\omega_\alpha,
\end{equation}
where $D$ [\si{\metre\squared\per\second}] is the molecular diffusivity coefficient and $\omega_\alpha=\rho_\alpha/\rho$ [-] is the mass fraction of the component $\alpha$ in the mixture.
where $D$\todo{$D$ was used to denote the spatial dimension in previous chapters} [\si{\metre\squared\per\second}] is the molecular diffusivity coefficient and $\omega_\alpha=\rho_\alpha/\rho$ [-] is the mass fraction of the component $\alpha$ in the mixture.
\inline{to deal with LBM compressibility ($\nabla\cdot\vec v \ne0$), we would likely have to solve the general equation (\cref{eq:transport:rho_alpha_2} below) or do different manipulations with it}
\inline{the continuity eq is $\nabla\cdot\vec v =0$ here, but LBM solves $\frac{\partial\rho}{\partial t}+\nabla\cdot(\rho\vec v)=0$ (according to LBMAT)}
@@ -356,7 +356,7 @@ However, water vapor is still released diffusively from the sides where the no-s
\inline{Connect to the sections in \cref{chapter:numerical methods}}
The transport equation \cref{eq:transport} has been incorporated into the mathematical framework of the \emph{NumDwarf} solver described in \cref{sec:MHFEM}.
Both \cref{eq:transport:conservative} and \cref{eq:transport:non-conservative} can be converted to the form \eqref{eq:NumDwarf} by taking $n=1$, $\vec Z =[\phi]$, $\matrix N =[1]$, $\vec m =[1]$, $\vec D =[D]$, $\vec w =\vec0$, $\vec r =\vec0$, $\vec f =[0]$, %[f_\phi]$,
Both \cref{eq:transport:conservative} and \cref{eq:transport:non-conservative} can be converted to the form of \cref{eq:NumDwarf} by taking $n=1$, $\vec Z =[\phi]$, $\matrix N =[1]$, $\vec m =[1]$, $\vec D =[D]$, $\vec w =\vec0$, $\vec r =\vec0$, $\vec f =[0]$, %[f_\phi]$,
and depending on the equation:
\begin{itemize}
\item$\vec u =\vec0$, $\vec a =[\vec v]$ for \cref{eq:transport:conservative},
