\documentclass[11pt]{article}
\usepackage{fullpage}
\usepackage{setspace}
\usepackage{parskip}
\usepackage{titlesec}
\usepackage[section]{placeins}
\usepackage{xcolor}
\usepackage{breakcites}
\usepackage{lineno}
\usepackage{hyphenat}
\usepackage{times}
\PassOptionsToPackage{hyphens}{url}
\usepackage[colorlinks = true,
linkcolor = blue,
urlcolor = blue,
citecolor = blue,
anchorcolor = blue]{hyperref}
\usepackage{etoolbox}
\makeatletter
\patchcmd\@combinedblfloats{\box\@outputbox}{\unvbox\@outputbox}{}{%
\errmessage{\noexpand\@combinedblfloats could not be patched}%
}%
\makeatother
\renewenvironment{abstract}
{{\bfseries\noindent{\abstractname}\par\nobreak}\footnotesize}
{\bigskip}
\titlespacing{\section}{0pt}{*3}{*1}
\titlespacing{\subsection}{0pt}{*2}{*0.5}
\titlespacing{\subsubsection}{0pt}{*1.5}{0pt}
\usepackage{authblk}
\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\providecommand{\tightlist}{\setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}%
\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}
\usepackage[english]{babel}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\usepackage[utf8]{inputenc}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{mathtools}
\usepackage{textcomp}
\usepackage{array}
\usepackage{graphicx}
\usepackage{xcolor}
\usepackage{xspace}
%\usepackage{ulem}
%\usepackage[numbers,sort&compress]{natbib}
\usepackage{hyperref}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\bibliographystyle{rsc}
%\linespread{2.0}
\makeatletter
\newcommand*{\overtabline}{%
\noalign{%
% normal "baselineskip" in tabular is height + depth of \@arstrutbox
\vskip-.5\dimexpr\ht\@arstrutbox+\dp\@arstrutbox\relax
% default line thickness is 0.4pt
\vskip-.2pt\relax
{\color{red}\hrule}
\vskip-.2pt\relax
\vskip+.5\dimexpr\ht\@arstrutbox+\dp\@arstrutbox\relax
}%
}
\makeatother
\newcommand{\todo}[1]{\textcolor{blue}{#1}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\renewcommand{\vec}[1]{{\ensuremath{\mathchoice
{\mbox{$\displaystyle\mathbf{#1}$}}
{\mbox{$\textstyle\mathbf{#1}$}}
{\mbox{$\scriptstyle\mathbf{#1}$}}
{\mbox{$\scriptscriptstyle\mathbf{#1}$}}}}}%
\newcommand{\tens}[1]{{\ensuremath{\mathchoice
{\mbox{$\displaystyle\mathsf{#1}$}}
{\mbox{$\textstyle\mathsf{#1}$}}
{\mbox{$\scriptstyle\mathsf{#1}$}}
{\mbox{$\scriptscriptstyle\mathsf{#1}$}}}}}%
\newcommand{\cvec}{\ensuremath{\vec{c}}}
\newcommand{\evec}{\ensuremath{\vec{e}}}
\newcommand{\fvec}{\ensuremath{\vec{f}}}
\newcommand{\jvec}{\ensuremath{\vec{j}}}
\newcommand{\nvec}{\ensuremath{\vec{n}}}
\newcommand{\rvec}{\ensuremath{\vec{r}}}
\newcommand{\uvec}{\ensuremath{\vec{u}}}
\newcommand{\vvec}{\ensuremath{\vec{v}}}
\newcommand{\xvec}{\ensuremath{\vec{x}}}
\newcommand{\zetavec}{\ensuremath{\pmb{\zeta}}}
\newcommand{\Fvec}{\ensuremath{\vec{F}}}
\newcommand{\e}{\epsilon}
\newcommand{\1}{\ensuremath{\tens{I}}}
\newcommand{\Cop}{\ensuremath{\mathcal{C}}}
\newcommand{\Iop}{\ensuremath{\mathcal{I}}}
\newcommand{\Lop}{\ensuremath{\mathcal{L}}}
\newcommand{\Pitens}{\ensuremath{\tens{\Pi}}}
\newcommand{\Omegaop}{\ensuremath{\mathsf{\Omega}}}
\newcommand{\Lambdaop}{\ensuremath{\mathsf{\Lambda}}}
\newcommand{\eq}{\text{eq}}
\newcommand{\nq}{\text{neq}}
\renewcommand{\d}{\ensuremath{\partial}}
\renewcommand{\O}{\ensuremath{\mathcal{O}}}
\newcommand{\therefore}{\item[$\rightarrow$]}
\newcommand{\ie}{\emph{i.e.}\xspace}
\newcommand{\eg}{\emph{e.g.}\xspace}
\newcommand{\cf}{\emph{cf.}\xspace}
\newbox\tmpbox
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{document}
\title{Multiscale Simulations of Porous Media: From Toy Models to Materials
Models}
\author[1]{Fang Wang}%
\author[2]{Ulf D. Schiller}%
\affil[1]{Affiliation not available}%
\affil[2]{Clemson University}%
\vspace{-1em}
\date{\today}
\begingroup
\let\center\flushleft
\let\endcenter\endflushleft
\maketitle
\endgroup
\selectlanguage{english}
\begin{abstract}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Multiscale modeling and simulation techniques are transforming the way
we can address questions regarding design, characterization, and
optimization of novel materials. Advanced computational models that
incorporate realistic geometries of porous media can be used as tools
to predict flow and transport phenomena. Recent developments in
mesoscopic and pore-scale modeling include workflows that combine
experimental information and direct modeling into an integrated
multiscale approach. This paper surveys the progress, challenges and
future directions in predictive modeling and simulation of
multiphysics phenomena in porous media.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
\end{abstract}%
\sloppy
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Introduction}
\label{sec:intro}
Porous media are materials that contain pores, \ie, void spaces that
are embedded in the condensed phase of the material. Examples include
rocks, clay, ceramics, membranes, and biological tissue. Porous media
continue to attract research interest as novel micro- and nanoscale
engineering applications are emerging, \eg, in the biomedical
\cite{Khanafer2006} and energy sectors \cite{Tarascon2001}. These
applications often pose multiple physics problems. For example, the
utilization of nonwoven fibrous membranes as filters or detectors
depends on morphology, wettability, and interactions of dissolved
components. A central research problem is the determination of
effective properties that describe the behavior of a porous medium on
macroscopic scales, \ie, the properties that are relevant and
observable in laboratory or field experiments.
The challenge in describing the macroscopic dynamics effectively lies
in the intricate structure of porous media. Three examples of
different types of porous media are shown in
Fig.~\ref{fig:porous-media}. The pore sizes typically span a range of
length scales, giving rise to complex morphology and topology that
often exhibit significant heterogeneity and anisotropy. Consequently,
the dynamics on the microscale are highly non-trivial and strongly
influenced by the geometry of the pore space. To predict macroscopic
transport properties thus requires both an accurate representation of
the geometry and sufficiently sophisticated models of the physical
interactions
Considerable efforts have been spent on developing advanced imaging
methods that can capture the complex structure of the pore
space. Experimental techniques such as X-ray computed tomography and
scanning electron microscopy are now capable of providing
three-dimensional images that resolve pores on the scale of microns
and smaller \cite{Blunt2001,Xiong2016}. Such images provide
information on the local properties of pore space at very high level
of detail \cite{Gueven2017}. This information can be used to construct
computational models that can predict the flow and transport through
the experimentally recorded pore space.
There are essentially two types of modeling approaches, namely,
pore-scale \cite{Blunt2013} and pore network models
\cite{Xiong2016}. While pore-scale models address the smallest scale
and aim at predicting effective parameters that capture the specific
interactions and surface effects for individual pores, pore network
models are mesoscale models that adopt effective models for single
pores and predict the flow and transport through larger networks of
pores. Pore-network models bridge between the pore-scale and the
macroscale and rely on models that discard irrelevant features of the
microscopic dynamics. In contrast, pore-scale models resolve a greater
amount of detail and are therefore more accurate but also
computationally more expensive. With increasing performance of
available computing resources, however, the computational demands of
pore-scale models have become less of a concern. On the other hand,
direct bridging between the pore-scale and the macroscale implicates a
need for improved pore-scale models that capture essential details of
physico-chemical processes with sufficient accuracy and fidelity.
A number of computational approaches can be applied to porous media,
including particle based methods and density functional techniques
\cite{Meakin2009}. In recent years, the lattice Boltzmann method has
gained in popularity as a versatile mesoscopic method that can be used
to model a range of physical processes beyond single-phase fluid
dynamics. Multiphase models \cite{Shan1993,Shan1994} and
electrokinetics solvers \cite{Capuani2004}, among others, have been
successfully implemented in the lattice Boltzmann framework. Moreover,
the underlying grid enables a straightforward representation of
complex geometries by dividing the grid into solid matrix and pore
space. The lattice Boltzmann method has been used to simulate a range
of transport phenomena in a variety of different porous media,
see~\eg~\cite{Nabovati2009,Boek2010,Ghassemi2011,Gueven2017}.
%
Such multiscale modeling and simulation can be a valuable aide in
designing systems using functionalized porous media and tailored
geometries. However, careful validation and optimization of computer
models are imperative and require close feedback with
experimentation. The path forward thus leads to an integrated approach
that connects theory and experiments through multiscale modeling and
data analysis.\selectlanguage{english}
\begin{figure}
%\begin{minipage}{.49\textwidth}
% \includegraphics[width=\textwidth]{sandstone}\par
% %{\scriptsize [R. Hilfer, \emph{Adv. Chem. Phys.} \textbf{92}, 299-424 (1996)]}
%\end{minipage}\\
\begin{minipage}{.33\textwidth}
\includegraphics[height=.16\textheight]{pores-two-phase}
%{\scriptsize [X. H. Yang et al., \emph{J. Phys. D} \textbf{46} 125305 (2013)]}
\end{minipage}%
\begin{minipage}{.33\textwidth}
\includegraphics[height=.16\textheight]{10971_2017_4387_Fig4_HTMLb}
%{\mbox{}\hfill\scriptsize Image courtesy of L. Jacobsohn}
\end{minipage}%
\begin{minipage}{.33\textwidth}
\includegraphics[height=.16\textheight]{jeecs_013_03_031005_f003}
\end{minipage}
\caption{{Examples of different types of porous media. (left) Aluminum
foam (porosity $\phi=0.76$) as a typical close-celled medium with
non-spherical pores. Reproduced from \cite{Yang2013}. Permission
pending. (middle) SEM micrograph of an electrospun YAG fiber membranes. Reproduced from \cite{Chen2017}. Permission pending.
(right) Mesostructure of lithium cobalt oxide reconstructed from
from a three-dimensional FIB/SEM image stack. The segmented active
particles are shown in different colors. Reproduced from
\cite{Roberts2016}. Permission pending.}}
\label{fig:porous-media}
\end{figure}
In this review, we highlight the basic building blocks of an
integrated approach for multiscale modeling of porous media. The
remainder of the paper is organized as follows. In section
\ref{sec:background}, we highlight a few selected applications of
porous media that are of timely interest. In section
\ref{sec:mesoscale}, we discuss the mesoscale structure of porous media
and how it can be recorded using advanced imaging techniques. Section
\ref{sec:theory} introduces the basic mathematical formulations that
are used to describe transport phenomena in porous media, and section
\ref{sec:modeling} reviews multiscale modeling approaches for
numerical simulation of multiphysics properties across the scales. We
summarize in section \ref{sec:conclusions} and close with a discussion
of challenges and future directions for multiscale modeling of porous
media.
%For example, petroleum recovery from reservoir rocks relies on
%multiphase flows with a complicated interplay of capillary pressures,
%permeability, and wettability of the porous formation.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Complex flows in porous media}
\label{sec:background}
Porous media are a class of materials that contain pores, \ie, void
spaces that are embedded in the material phase. Porous media
encompass, \emph{inter alia}, rocks, clay, ceramics, membranes, and
biological tissue. Three examples of porous media are shown in
Fig.~\ref{fig:porous-media}. Since the void space is typically filled
with fluids (gases or liquids), flow and transport through porous
media are of particular interest. However, the pore sizes typically
span a range of length scales and the resulting complex geometry leads
to highly non-trivial flow patterns. In addition, the dynamic behavior
is often governed by an interplay of different physical mechanisms
such as viscous friction, interfacial tension, and charge
interactions. The transport coefficients emerge from an intricate
competition of hydrodynamic, thermodynamic, and electrokinetic effects
in the bulk and at the surface of the complex geometry. General models
to obtain the transport coefficients as a function of the bulk
properties of the medium and the fluid are not available. The
prediction of flow phenomena in porous media thus remains a formidable
task, and is a genuine multiscale and multiphysics problem. For
engineering applications, on the other hand, the rich phenomenology of
porous media offers ample opportunities to design systems using
functionalized materials and tailored geometries. Here, we exemplarily
highlight a few applications of porous media that are of timely
interest.
A classical application of porous media is filtration. Aerosol
filtration is relevant in a range of sectors including civil and
automotive engineering, and health. Similarly, water purification is
among the most pressing needs for global health
\cite{Shannon2008,Amin2014}. Filtration is based on the removal of
contaminants from gaseous or liquid fluids streaming through a porous
medium, with typical pore sizes on the micron scale. The efficiency of
filtration media depends primarily on the pore size and the flow
resistance, and there is typically a trade off between the collection
efficiency and the required pressure drop \cite{Theron2017}. However,
the detailed filtration characteristics depend on a number of factors
such as the morphology of the medium, pore size distribution, and
surface functionalization. One important class of filtration media are
fibrous membrane fabrics which can be further divided into woven and
nonwoven fabrics. Woven fabrics consist of a regular pore structure,
nevertheless the texture varies depending on the type of weave. This
can significantly affect the fluid flow, and the flow pattern in the
basic fabric pore have an influence on the overall flow resistance
\cite{Lu1996}. Theoretical or numerical predictions of the
permeability or flow resistance thus have to take the details of the
tortuous geometry into account. A numerically computed flow field
through an array of fibers representing a model fibrous membrane is
shown in Fig.~\ref{fig:fiber-flow}. Nonwoven fabrics are a low-cost
alternative thanks to inexpensive manufacturing methods such as melt
blowing and spun bonding. However, their characterization is more
complicated due to local variations in fiber orientation, pore size
distribution, etc. Nonwoven fabrics thus exhibit a more heterogeneous
and anisotropic morphology which can be both advantageous and limiting
for filtration \cite{Gregor2003}. The irregular structure of the void
space makes it more difficult to determine the tortuosity which
influences the flow through the interstices and the filtration
performance \cite{Vallabh2010}. Characterization of the porous
geometry of nonwoven fibrous media \cite{Theron2017} is thus an
important aspect in gaining a better understanding of the macroscopic
flow properties, which may ultimately guide improvement of the
filtration capabilities.\selectlanguage{english}
\begin{figure}
%\includegraphics[width=.49\textwidth]{ordered}
\includegraphics[width=.45\textwidth]{regular-geometry}
%\includegraphics[width=.45\textwidth]{L_200_R_13}
%\includegraphics[width=.49\textwidth]{L_200_R_13}
\includegraphics[width=.45\textwidth]{velocity}
%\includegraphics[width=.3\textwidth]{porous}
\caption{{Flow through an array of fibers representing a model
fibrous membrane. The fibers are arranged in layers where each
layer is rotated by $90^\circ$ with respect to the previous
one. Simulations were performed using the lattice Boltzmann method
implemented in the LB3D software package
\cite{Schmieschek2017}. The flow field is visualized with arrow
glyphs where the color indicates the velocity magnitude. The flux
can be measured for varying pressure gradient in order to
determine the permeability in Darcy's law.}}
\label{fig:fiber-flow}
\end{figure}
Another active area of research is hydrocarbon recovery from
reservoirs. Here the scales range from the millimeter grain size to
meters on the reservoir scale and to kilometers on the field
scale. The oil recovery process is governed by viscous friction,
surface wetting, and multiphase interactions between hydrocarbons and
injected fluids \cite{Zitha2011}. In order to extract hydrocarbons
from the reservoir, the viscous drag in the porous medium has to be
overcome by the pressure gradient. This is basically the simplest
realization of Darcy's law (see section \ref{sec:theory}), where the
permeability is a property of the porous geometry. The natural
pressure is typically only sufficient to recover about 20\% of
hydrocarbons from the reservoir \cite{Zitha2011}. While pressure
depletion can be delayed by expansion of aquifers or gas caps, the
most common method to sustain pressure is waterflooding. Besides
maintaining the pressure, the injected water turns the fluid into a
multiphase mixture, where the oil displacement depends on factors such
as interfacial tension and surface wettability. The multiphase
behavior opens up various routes for enhanced oil recovery that can
broadly be divided into thermal and nonthermal recovery methods. While
thermal methods mainly aim at decreasing the oil viscosity through
heating, nonthermal methods modify the mixture composition in order to
tune the multiphase interactions so to increase displacement
efficiency. For example, surfactant flooding is used to tune the
interfacial tension and wettability, which can lead to a reduction in
capillary pressure allowing more oil to be extracted from the
pores. Surfactants can also be mixed with gas to create foams that
improve the mobility of the residual oil. Modifying the rock/fluid
interactions can also increase the displacement efficiency. Apart from
targeting the surface chemistry, one can introduce ionic interactions
to activate electrokinetic transport in response to electric fields
\cite{Ghazanfari2012}. For instance, a change in ionic composition has
been linked to improved oil displacement for low-salinity
waterflooding \cite{Seccombe2008,Webb2008}, however, the underlying
mechanism is not fully understood.
Electrokinetic transport in porous media involves additional physical
mechanisms that are also highly relevant for batteries and fuel
cells. In Li-ion batteries, a voltage is present due to free energy
difference of lithium intercalation in the porous cathode and anode,
respectively. Electric current results from the transport of electrons
in the circuit while the Li-ions undergo transport from the anode to
the cathode. The electrodes are composite porous media consisting of
active particles, binder, and conductive. The pore space is filled
with electrolyte, and during discharge the lithium ions migrate to the
surface of the active particles and through the electrolyte. The
capacity and conductivity depend strongly on the properties of the
porous electrodes, primarily the cathode material, which is often
lithium cobalt oxide. Over the last one or two decades, considerable
efforts to improve performance and safety have focused on the
electrochemistry \cite{Tarascon2001}. One of the performance
degradation issues is capacity fade due to loss of lithium to side
reactions. This can occur in particular in the solid-electrolyte
interphase and while it can be addressed by modifying the surface
chemistry, the effect of the pore structure on capacity fade has
received little attention. Another concern is the swelling of the
active particles upon Lithium intercalation which can degrade the
performance and ultimately lead to failure \cite{Garrick2014,Abada2016}.
In order to better understand the non-equilibrium dynamics during
charge and discharge, the coupling between electrochemistry,
thermodynamics, and mechanics needs to be
considered \cite{Roberts2016,Smith2017}. To this end, multiphysics
models can be used to simulate the behavior of realistic porous
geometries. Roberts et al. \cite{Roberts2016} have performed
simulations of experimentally reconstructed mesostructures and found
that both the particle shape and the structure of the porous network
can significantly affect the capacity and mechanical performance of
the battery. The detailed functional dependence of macroscopic
properties on the mesoscale structure, \eg, pore size distribution,
anisotropy, etc., is the subject of ongoing research using multiscale
models \cite{Lee2016,Kashkooli2016,Turner2015}.
The above examples illustrate how the multiscale and multiphysics
nature of porous media gives rise to complex flow behavior. Porous
media are indeed multiscale systems in a twofold way: First, the
morphology of the medium and the pore size distribution introduce a
range of mesoscopic length scales inherent in the geometry. Second,
the physical interactions may span time- and length-scales of
different order of magnitude, \eg, the diffusion in the active
particles occurs on shorter times than the diffusion through the
electrolyte. Bridging the microscale and the macroscale thus does not
only require multiscale models to capture the relevant physical
mechanisms, but also experiment specific representations of the porous
geometry that accurately reflect the morphology of the porous medium
across scales.
%\begin{itemize}
%\item General aspects of flow in porous media
%\item Structure-property relations
%\item Important applications
%\begin{itemize}
%\item Filters/Detectors
%\item Oil recovery (hydrocarbon extraction)
%\item Fracture dissolution
%\item Batteries (electrokinetics)
%\end{itemize}
%\item Opportunities for materials science
%\end{itemize}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Some setences are duplicate, which can be deleted. The rest can be put forward.
\textcolor{green}{The investigation associated with hydraulic properties in fracture porous systems
is one of the most challenging problems that engineers have to face. The unique features
displayed by transport phenomena and flow in fractures include the many orders of magnitude
change in permeability and flow rate, compared with granular porous media. Research
into microscopic fracture geometry and detection methods has mainly been motivated by
the issues related to water quality, underground storage reservoirs and repositories
of nuclear waste. Laboratory measurements of water flow in single, saturated fractures
is complicated due to the anisotropy of natural fractures. Quantitative analysis of
solute transport is more difficult to elucidate even with detailed geological
characterization in fractured rock systems \cite{Hadermann,1996}. There is a clear
need for reliable model from simulation studies to couple the geometrical and hydraulic
properties in large scale. The transport behavior of single phase flow can be further
studied in filter media. The importance of designing efficient water filtration systems
in the area of environmental engineering is well know and served as a motivation for
research on hydraulic properties of fibrous filters. Pores of filter membrane range
in size over several orders of magnitude due to the multi-scale of particle size in
water treatment. The conventional filter typically are used for the removal of particles
larger than 1-10 $\mu m$. \cite{Chang,2006} conducted experimental studies to demonstrate
that non-woven membranes bioreactor system can be used in wastewater treatment by proper
selection of its pore size. Nanofiber filter prepared by electrospinning technique provides
capability to resist penetration of chemical and biological pollution and accelerate the
progress of research of aerosol filtration. 3D structure model of manufactured polyurethane
filter has been proposed by \cite{Sambaer,2011} to predict the filtration efficiency as a
function of particle size.}\
%The following paragraphs are duplicate with previous part.
\textcolor{blue}{A new direction that have attracted much attention in energy engineering is the optimization
of performance in fuel cells. Many attempts have been made to develop a model for PEMFCs
that account for 3D multiphase flow originating from electrochemical reactions
\cite{Kone,2017}\cite{Tamayol,2011}. The multiphase transport is composed of the motion of
reactants, hydrogen and oxygen flow, and the removal of excess produces water. The great
potential of PEMFCs as a clean power source has increased the needs of research on
permeability of gas diffusion layers. The electrochemical performance of lithium-ion
batteries is very similar to fluid flow behavior, in ways that mesoscale network of
active material particles is in complex geometry with small particle size. The result
of coupled electrochemical, mechanical and thermal simulations on reconstructed mesosctures
is capable of explaining electronic transport mechanisms\cite{Roberts,2016}. Biomedical
engineering is another major area that demands the mimic of electrokinetic transport.
\cite{Faraji,2011} quantitatively represents the contribution of $\zeta$-potential
to solute transport in brain tissue by using organotypic hippocampal slice cultures as
a representative for porous brain tissue.}\
\textcolor{blue}{Another active area of research is hydrocarbon recovery from reservoirs. Here the scales
range from the mm-size of sandstone grains to several meters for the reservoir scale and
to kilometers on the field scale. The oil recovery process is governed by viscous friction,
surface interaction, and multiphase interactions between hydrocarbons and injected fluids \cite{Zitha}.}
\textcolor{blue}{Modeling the complex flow behavior in porous media presents a challenge because the
description of dependence on pore structure of macroscopic properties requires high
computational effort. The problem of viscous flow past a sold matrix is part of the
study of the relative motion of a mixtures of a fluid and solid bodies. In considering
a realistic porous medium which exhibits a serpentine pore sequence with complex
interconnections, the flow path is also tortuous. Unconnected pores do not contribute
to fluid flow. A macroscopic homogeneous media have uniform properties with the change
of representation elementary volume(REV), which means the pore structure parameters,
such as permeability or porosity will not change in a series of samples of increasing
size. A second reason for the complexity relates to the three dimensional geometry in
which a statistically meaningful representation of flow field is obtained by solving
the full three-dimensional flow equation. While a simplified flow scheme in a saturated
single phase flow can provide quantitative simulation results, most engineering applications
involving fuel cells and oil extractions generally contains multiphase and multi-components,
water-air glow, water-oil flow. A fundamental feature of simulation in immiscible fluids is
the introduction of interfacial energy and equilibrium distribution of multiple phases. For
precise description of pore-scale transport phenomena in porous medium, both the fluid properties
in terms of its rheology and pore structure parameters should be taken into account.}\\
\subsection*{Characteristics of Complex Flow}
The physical properties of flow though porous media have a great influence on the hydraulic conductivity of porous system. Some recent research on single phase flow is receiving much attention due to its application in real industrial problems. \citealt{tamayol_single_nodate} studied the through-plane gas permeability of carbon paper in gas diffusion layers of proton exchange membrane fuel cells (PEMFCs). They found a reverse relation between permeability and the PTFE content, which is useful in design of PEMFCs. Transverse air flow has also been used in the experiment work of measuring permeability of nonwoven fabrics to investigate the effect of fabric thickness\citep{vallabh_new_2010} . Most of the research in the last decades have been concentrated on multiphase and multicomponent flows due to the necessity of complicated models for a wide range of transport phenomena in porous media. A major improvement of lattice Boltzmann method in simulating multiple phases and components fluid is achieved by incorporating interpaticle potential and surface tension. Work by \citealt{boek_lattice-boltzmann_2010} presents simulation results of relative permeability of binary immiscible fluids in the Bentheimer rock sample. The phase seperation and wettability situations are modelled with lattice Botzmann method. Another fluid properties affects the computed permeability by LBM is viscosity. Viscosity effect has been studied by invoking different LBE scheme\citep{cho_permeability_2013}\citep{pan_evaluation_2006}
\subsection*{Pore Geometry}
It is difficult to characterize the void space structure due to its complex nature. The distinct pores shape is intrinsically a result of available choices in particle shapes, arrangement and surface roughness. Many efforts have been made to describe the influence by pore shape. The narrow channel in dense fibrous systems can be viewed as a bunch of interconnected capillary tubes from lubrication theory \citep{yazdchi_validity_nodate}. They considered the orientation of narrow channels at different porosities and an isotropic distributions of orientation is observed. \citet{azhdari_heravi_investigation_2015} constructed two dimensional porous media packing with random rectangular particles, which give rise to straight rectangular channels. The measurement of tortuous flow path is simple in such two-dimensional media. Flow configuration by specifying spherical particles in face-centered cube(FCC) and body-centered cube(BCC) has been studied by \cite{Eshghinejadfard_2016}, which provides a highly symmetric flow field. The steady-state fluid velocity in a smooth channel can be easily observed. However, most realistic porous media like soil and metal foam is composed of coarse and random particles. The rough surfaces due to sediments of rocks and clays may have a significant impact on the flow behavior. \citet{chukwudozie_pore_2013} investigated the effect of grain roughness and found that higher roughness height results in more tortuous flow paths.
\subsection*{Applications}
The object of current work is to present a prospective study of complex flow in porous media, in several respects. This fundamental model-based approach has close connection to other disciplines and provides theoretical basis for materials designing. The most concentrated research in this field is concerning multiphase flow which has important technical application, most notably in oil extraction from petroleum reservoirs. Understanding the multiphase flow in the reservoirs at the microscopic scale helps researchers to predict the relative permeability as a function of material properties. Another new direction that have attracted much attention in energy engineering is the performance of fuel cells. Many attempts have been made to develop a model for PEMFCs that account for 3D multiphase flow\citep{kone_three-dimensional_2017}\citep{tamayol_-plane_2011}. The electrochemical performance of lithium-ion batteries is very similar to fluid flow behavior, in ways that mesoscale network of active material particles is in complex geometry with small particle size. The result of coupled electrochemical, mechanical and thermal simulations on reconstructed mesosctures is capable of explaining electronic transport mechanisms\citep{roberts_insights_2016}. Biomedical engineering is another major area that demands the mimic of electrokinetic transport. Faraji, 2011 quantitatively represents the contribution of $\zeta$-potential to solute transport in brain tissue by using organotypic hippocampal slice cultures as a representative for porous brain tissue.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\section*{Mathematical models for porous media}
\section*{Mathematical modeling of transport phenomena in porous media}
\label{sec:theory}
Most transport phenomena require to solve a boundary value problem on
the complex geometry of pore space. For instance, incompressible fluid
flow is governed by the Stokes equations
%
\begin{align}
\eta \nabla^2 \uvec(\rvec) &= \nabla P(\rvec) , \\
\nabla\cdot\uvec(\rvec) &= 0 ,
\end{align}
%
where $\uvec(\rvec)$ denotes the flow velocity, $P(\rvec)$ the
pressure, and $\eta$ the dynamic viscosity of the fluid. The boundary
condition on the surface of pore space is typically the no slip
boundary condition
%
\begin{equation}
\uvec(\rvec_s) = 0 .
\end{equation}
%
Due to the complicated geometry of the surface, obtaining analytical
solutions for porous media is often not feasible. Therefore,
simplified descriptions based on empirical relationships are commonly
used. The most widely used formulation is \emph{Darcy's law}
\cite{Darcy1856,Dullien1992}
%
\begin{equation}
\frac{1}{A} \vec{q} = - \frac{1}{\eta}\,\mathsf{K}\cdot\nabla P
\end{equation}
%
which describes the discharge rate $\vec{q}$ through a cross-section
$A$ of the porous medium subject to a hydraulic gradient. The linear
Darcy's law is valid for stationary low Reynolds number flow of an
incompressible Newtonian fluid with dynamic viscosity $\eta$. The
\emph{permeability tensor} $\mathsf{K}$ is an effective physical
property (with units of area) that depends on the geometry and
topology of pore space. Darcy's law can be derived from the Stokes
equations by means of an asymptotic expansion in a small parameter
$\epsilon$ that represents the ratio of microscopic to macroscopic
length scales, \cf \cite{Hilfer1996}. In an isotropic porous medium,
the permeability tensor $\mathsf{K}$ reduces to a scalar permeability
$k$.
The determination of the permeability $k$ as a function of the
geometry and topology of the porous medium has been a central
objective in many investigations and remains an active research
topic. A widely referenced approach is the Kozeny-Carman model
\cite{Kozeny1927,Carman1937,Carman1956,Dullien1992}, which is based on the
capillary tube model. The discharge rate of $N$ capillary tubes of
radius $a_i$ and lengths $L_i$ is $Q=\sum_{i=1}^{N} \pi a_i^4 \Delta
P/(8\eta L_i)$, and the average permeability of a cubic sample of size
$L$ is thus
%
\begin{equation}
\langle k \rangle = \frac{\pi}{8} \frac{N}{L} \left\langle \frac{a^4}{L} \right\rangle .
\end{equation}
%
The Kozeny-Carman model assumes that a mean-field approximation is
valid in which case the average permeability can be written as
%
\begin{equation}
\langle k \rangle
%= \frac{\pi}{8} \frac{N}{L} \frac{\langle a \rangle^4}{\langle L \rangle}
= \frac{\langle \phi \rangle \langle a \rangle^2}{8 \langle T \rangle^2}
= \frac{\langle \phi \rangle^3}{2\langle T \rangle^2 \langle S \rangle^2} .
\end{equation}
%
In the mean-field approximation, the ratio $\langle\phi\rangle/\langle
S\rangle$ becomes equal to the average hydrodynamic radius
$R_h=\langle a\rangle/2$ \cite{Hilfer1996}.
The semi-empirical Kozeny-Carman permeability $k$ for monodisperse
granular media is commonly written in the
form \cite{Pan2001,Ahmadi2011,Duda2011}
%
\begin{equation}
k = \frac{1}{C \, T^2 S^2} \frac{\phi^3}{(1-\phi)^2} %= \frac{d_p^2}{36C T^2} \frac{\phi^3}{(1-\phi)^2} ,
\end{equation}
%
where $C$ is the Kozeny-Carman (KC) constant. The KC constant depends
on the specific structure of the porous medium and a considerable body
of work has been devoted to determining permeability relations for
various types of porous media \cite{Ahmadi2011,Xu2008}. However, in
most cases, only the value of $CT^2$ can be extracted and the
Kozeny-Carman constant and the tortuosity $T$ remain undetermined. The
tortuosity is used in various contexts, and the difference between the
geometric and hydraulic, electric, diffusional, etc. tortuosities has
been recognized
\cite{Matyka2008,Duda2011,Valdes-Parada2011,Ghanbarian2013,Allen2017}.
Accordingly, no unified theory for tortuosity is available
\cite{Matyka2008}, and numerical and experimental correlations carry
some uncertainty regarding their applicability to different types of
porous media and porosity range. An analytical expression for
tortuosity of sphere arrays can for example be derived based on a
representative elementary volume (REV) model \cite{Ahmadi2011} or by
means of fractal geometry \cite{Xu2008}. In practice, the
tortuosity-porosity relation is commonly cast in the power law form
$T-1=p(1-\phi)^\gamma$, however, the values of $\gamma$ extracted from
empirical trends remain non-universal \cite{Allen2017}.
%
Table \ref{tab:permeability} and \ref{tab:tortuosity} list various
correlations for permeability and tortuosity as a function of porosity
for different porous media. We should note that the table only lists a
selection of the most common correlations used in the literature.
%
Most studies have focused on the scalar Kozeny-Carman permeability
without attempting to clarify the tensorial nature of the permeability
in Darcy's law. This is just one aspects where computational models
can provide additional insight by revealing tensorial
permeability-porosity relationships.\selectlanguage{english}
\begin{table}\footnotesize
\begin{tabular}{cc>{$}l<{$}ll}
\hline
Porous structure & Method & \text{Permeability-porosity relation} & Remarks & References \\
\hline
Granular media & empirical & k = \frac{d_p^2}{36C T^2} \frac{\phi^3}{(1-\phi)^2} & Kozeny-Carman, $d_p$: average grain diameter & \cite{Carman1956,Dias2006,Ahmadi2011} \\
Porous media & empirical & k = C d_p^2 \phi^\gamma & Rumpf-Gupte, $d_p$: average grain diameter & \cite{Rumpf1971,Bourbie1987} \\
%Porous media && k = C d^2 \phi^n && \cite{Bourbie} \\
%Monodisperse sphere packing & experimental & k = \frac{d^2_p}{180}\frac{\phi^3}{(1-\phi)^2} && \cite{Ahmadi2011} \\
Square particles & numerical & k = \frac{\phi^3}{CT^2S^2} && \cite{Koponen1997} \\
%Mixed granular beds & experimental & k = \frac{d^2_{av}\phi^3}{36K_0T^2(1-\phi)^2} && \cite{Dias2006} \\
Textile assembly & empirical & k = \frac{d_f^2}{16C} \frac{\phi^3}{(1-\phi)^2} & $d_f$: filament diameter & \cite{McGregor1965} \\
Glass and fiber & experimental & k = \frac{\phi^{n+1}}{C(1-\phi)^n} && \cite{Rodriguez2004}\footnotemark[1] \\
Fiber mats & theoretical & k = C \frac{\phi^n}{1-\phi} && \cite{Costa2006}\footnotemark[1] \\
%Fibrous media & numerical & k = C \beta \delta_\text{min}^2 \tau && \cite{Tamayol2011} \\
Random fibers & numerical & k = C r_f^2 \left(\sqrt{\frac{1-\phi_c}{1-\phi}}-1\right)^\gamma & $r_f$: fiber radius & \cite{Gebart1992,Clague2000,Nabovati2009} \\
\hline
\end{tabular}
\footnotetext[1]{As cited in \cite{Xu2008}.}
\caption{{Permeability-porosity relations proposed in the literature
for different porous media. Here $C$ is a constant permeability
factor (sometimes called Kozeny-Carman constant), $S$ is the
specific internal surface, and $T$ is the tortuosity of the porous
medium. $\gamma$ is an (empirical) scaling exponent.}}
\label{tab:permeability}
\end{table}\selectlanguage{english}
\begin{table}\footnotesize
\begin{tabular}{cc>{$}c<{$}ll}
\hline
Porous structure & Type & \text{Tortuosity-porosity relation} & Remarks & References \\
\hline
Granular beds & electric & T \propto \phi^{-\gamma} & Archie & \cite{Archie1942,Dullien1992,Dias2006} \\
Isotropic representative unit cell & diffusive & T = \frac{\phi}{1-(1-\phi)^\gamma} & 2D: $\gamma=\frac{1}{2}$, 3D: $\gamma=\frac{2}{3}$ & \cite{Plessis1988,Plessis1991,Plessis2010} \\
Sediments & diffusive & T-1 = p (1-\phi) & sand: $p=2$, clay: $p=3$ & \cite{Iversen1993,Matyka2008,Duda2011} \\
Sandy sediments & diffusive & T = \sqrt{1+p(1-\phi)} & & \cite{Iversen1993} \\
Packed beds & hydraulic & T = 1 - p \ln\phi & spheres: $p=\frac{1}{2}$, cylinders: $p=\frac{2}{3}$ & \cite{Weissberg1963,Comiti1989,Barrande2007,Matyka2008,Mauret1997,Duda2011,Koponen1996} \\
Freely overlapping obstacles & hydraulic & T-1=p(1-\phi)^\gamma & & \cite{Iversen1993,Matyka2008,Duda2011,Koponen1996} \\
Freely overlapping obstacles & hydraulic & T-1\propto R\frac{S}{\phi} & $R$: hydraulic radius, $S$: specific surface & \cite{Matyka2008} \\
Granular media & hydraulic & T = \sqrt{\frac{2}{3}\frac{\phi}{1-p(1-\phi)^{2/3}}+\frac{1}{3}} & & \cite{Ahmadi2011} \\
\hline
\end{tabular}
\caption{{Tortuosity-porosity relations proposed in the literature for
different porous media. Here $\gamma$ is an (empirical) scaling
exponent and $p$ is an (empirical) shape factor.}}
\label{tab:tortuosity}
\end{table}
% Rumpf and Gupte
% Ergun
% limited accuracy
% breakdown for low and high porosities
% local porosity theory
\subsubsection*{Electrokinetic transport}
Electrokinetic transport in a porous medium filled with an electrolyte
is described by the continuity equation and the Nernst-Planck equation
%
\begin{align}\label{eq:nernst-planck}
\frac{\partial \rho_i}{\partial t} + \nabla \cdot \vec{j}_i &= 0 \\
\vec{j}_i &= \rho_i \uvec - D_i^* \nabla \rho_i - \mu_i^* z_i \rho_i \nabla \Phi
\end{align}%\todo{Needs checking.}
%
where $\rho_i$ and $\vec{j}_i$ are the density and flux of species $i$
in the fluid phase, respectively, $\uvec$ is the Darcy velocity of the
electrolyte, $D_i^*$ and $\mu_i^*$ are the effective diffusion
coefficient and mobility, and $\Phi$ is the electrostatic potential.
%
In a porous medium, the effective transport coefficients are obtained
by renormalizing the microscopic values with the tortuosity, \ie,
$D_i^*=D_i/T^2$ and $\mu_i^*=\mu_i/T^2$.
%and the flux $\vec{j}^*=\phi \vec{j}$ is adjusted to account for the porosity.
The electrostatic potential is described by the Poisson equation
$\nabla\cdot(\epsilon\nabla\Phi)=-\sum_i z_i e \rho_i$, where $z_i$ is
the valency of species $i$, $e$ is the electron charge, and
$\epsilon$ the permittivity. The Poisson equation needs to be solved
subject to boundary conditions on the surface of pore space which can
prescribe a surface charge. The surface charges will lead to formation
of a Debye layer which affects the potential of mean force and the
charge transport. The investigation of the structure of the Debye
layer in porous media is a direction for further research that may
reveal how microscopic details affect the macroscopic transport of
ions.
%
It should be noted that the Nernst-Planck equation holds for
(semi-)dilute electrolytes. For concentrated electrolytes
eq. \eqref{eq:nernst-planck} needs to be replaced with Stefan-Maxwell
coupled fluxes \cite{Smith2017}.
\subsubsection*{Immiscible displacement}
Multiphase transport of species in a binary mixture is governed by the
interplay of viscous, capillary, and volumetric (pressure, gravity)
forces. One the pore scale, the flow of the immiscible components is
described by the incompressible Navier-Stokes equations. At the
fluid-fluid interface, the velocities of the components are equal and
the stress-balance $\sigma_i = \sigma_j + 2\gamma_{ij} \kappa$, where
$\sigma_i$ is the fluid stress of species $i$, $\gamma_{ij}$ the
surface tension between $i$ and $j$, and $\kappa$ is the curvature of
the interface.
%
On the macroscopic scale, multiphase flow in a porous medium is described by a generalization of Darcy's law
%
\begin{align}
\frac{1}{A} \vec{q}_i &= - \sum_j \mathsf{K}^\text{rel}_{ij} \frac{1}{\eta_j} \mathsf{K}\,\nabla P_j \\
\phi \frac{\partial S_i}{\partial t} &= \nabla\cdot\uvec_i ,
\end{align}
%
where $\eta_i$ and $\nabla P_i$ are the dynamic viscosity and partial
pressure of species $i$, respectively, and
$\mathsf{K}^\text{rel}_{ij}$ is the \emph{relative permeability} of
species $i$ with respect to species $j$, and the saturations $S_i$
satisfy $\sum_i S_i=1$. The importance of the microscopic wetting
behavior on the multiphase transport was highlighted by Hilfer
\cite{Hilfer1996} along with the dimensional analysis of the governing
equations \cite{Hilfer1996b}. By denoting the breakthrough
pressure from the capillary pressure, those authors elucidated the
process dependence of the desaturation as a function of the capillary
pressure and explained observed differences for oil recovery between
laboratory and field experiments \cite{Morrow1990}. The calculation of
accurate relative permeabilities and capillary pressures is an active
area of research where multiscale simulations can provide further
insights.
%\begin{itemize}
%\item analytical method in simple geometry (spherical and square)
%\item theoretical tortuosity-porosity relations from previous work, Formation factor?)
%\item
%\end{itemize}
Permeability is an essential parameter in probing the fluid flow behavior of porous media due to its direct effect on bulk properties. Analytical and numerical study of hydraulic permeability have steadily received much attention to determine the effect of porous structure on hydraulic properties. A simple empirical relation in macroscopic terms, Darcy's law, states that the fluid flux is proportional to the pressure gradient. In the creeping flow regime, Newtonian fluid flow is governed by Darcy's law which can be written as:\
\begin{equation}
\label{eqn:drag}
__=\frac{-\kappa}{\mu}\nabla p
\end{equation}
where $____$ is the volume averaged fluid velocity in the entire flow domain and $\mu$ is the dynamic viscosity of the fluid. The proportionality constant, $\kappa$, denotes the permeability of the medium. Darcy's law is valid in flows where fluid velocity is low so that the intertial force is much smaller than viscous forces. In this situation, Reynold number is very small $(Re<1)$ and forchheimer drag is negligible. Permeability is independent of fluid properties including viscosity, saturation and density. Only geometrical features of the solid matrix contribute to the intrinsic permeability.
A broadly used semiempirical relationship between permeability and porous medium properties, Kozeny- Carman(KC) formula(Carman, 1937), read:
\begin{equation}
\label{eqn:drag}
k=\frac{\phi^3}{\beta T^2 S^2}
\end{equation}
where $\beta$ is the Kozeny constant that depend on the channel shape. The other three geometrical parameters, $\phi$, $T$ and $S$, denotes porosity, hydraulic tortuosity and specific surface area per unit volume of particles, respectively. While KC equation is designed particularly for granular random packings, the validity of KC equation for flow through unidirectional fibrous media is reported in the work of Gutowski, 1987. Kozeny constant is an experimental based value and cannot be derived straightforwardly.\
The most important structure parameter affecting the permeability is porosity, a measure of relative volume of void space. Fluid can only transport through interconnected pores which define the concept of effective porosity. Dead-end pores or isolated pores barely contribute to flows and fluid in them is stagnant in fact. The advantage of using porosity to characterize permeability is that porosity is experimental measurable and can be determined by voxel resolution in simulation with less computational work. There are many methods to experimentally measure porosity, including determining bulk volume of a sample directly or measuring the electrical properties of pore space. The relationship between permeability and porosity has been extensively studied since the pioneer work of Kozeny, 1927. In recent years, considerable progress has been made to elucidate the effect of porosity based on computational method. At low porosity, capillary bundle approximation is a proper approach for modeling porous materials, in which the flow path is assumed to be aligned and tortuous capillary tubes. KC equation is based on this dense packing where the stokes equation leads to the Poiseuille flow through pipes. However, at higher porosity, cell method is appropriate by dividing the void space into independent cells without the influence of fiber shape. Yazdchi, 2011 proposed a modified KC coefficient for intermediate porosity range and showed that permeability of random and regular porous media is the same at high porosity.\
More specific parameters that affect permeability includes particle size, aspect ratio, specific surface, representive element volume ( REV), pore throat, pore size distribution.
When performing the flow simulation, the impact of size of the numerical domain cannot be ignored because it corresponds to different pore network's randomness. For a smaller domain size, the distribution of permeability is broader with respect to porosity. Specific surface of monosized sphere packing can be estimated from the pore size distribution. The dependence on specific surface of permeability is revealed by particles' diameter. Pore size distribution can be determined by the mercury porosimetry technique and by X-ray micro-tomography. Many simulation models incorporate fiber radius or particle diameter in predicting permeability. Nabovati,2009 confirmed that the impact of fiber aspect ratio is significant in fibrous media when it is smaller than 6 by LBM simulation. Heravi, 2015 explored the accuracy at different aspect ratio. In addition, pore throat mean diameter are shown to be an essential parameter in that higher mean pore throat diameter corresponds to higher permeability by Gueven, 2017.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Mesoscale structure of porous media}
\label{sec:mesoscale}
Virtually all materials contain internal surfaces and phase boundaries
and thus are porous media at the microstructure level. The length
scales where this occurs can vary, \eg, from nanometers for porous
silicon to meters for rocks and soils. Moreover, the shape, size, and
connectivity of the pores gives rise to a variety of complex
microstructures. In most cases, a multi-level hierarchy of pore
structures exists, and local variations give rise to heterogeneous
morphology on larger scales. The length scale where the size and shape
of pores and connecting throats is observed is commonly referred to as
the \emph{mesoscale}. The complexity of the mesoscale architecture
influences the macroscopic properties of porous media, and the
development of predictive models thus relies on the accurate
characterization of the mesoscale structure.
The geometry of a porous medium is given by the partition of the
sample into the solid matrix $\text{M}$ and the pore space
$\text{P}$. However, on macroscopic scales it is impractical to
specify the geometry in full detail, and general characterizations of
the geometry are sought. Commonly used characterizations of porous
media include the porosity $\phi$, hydraulic radius $r_h$, and
tortuosity $T$. The porosity is defined as volume fraction of pore
space, \ie, the ratio $\phi=V_\text{P}/V_\text{bulk}$. In a porous
medium filled with a single phase fluid, the porosity can also be
expressed in terms of the densities of the phases
$\phi=(\rho_\text{M}-\rho_\text{bulk})/(\rho_\text{M}-\rho_\text{P})$. Porosity
can be measured experimentally using, \eg, gas expansion or mercury
intrusion porosimetry \cite{Hilfer1996,Xiong2016}. The latter can
measure the open porosity only, \ie, the pore space that is connected
to the surface of the sample.
%
Another geometric quantity of interest is the specific internal
surface $S=A_\text{P}/V_\text{bulk}$. The specific internal surface
yields an inverse length and quantifies the surface to volume ratio of
the porous medium. It can be measured using similar experimental
techniques as for porosity, \eg, gas adsorption isotherms. %This sentence is duplicated with the latter one which is also marked red
\textcolor{red}{The interpretation of the isotherms by means of the Brunauer-Emmett-Teller
(BET) or Barret-Joiner-Halenda (BJH) methods involve certain
assumptions about the pore geometry and thus only appropriate for
specific types of porous materials \cite{Xiong2016}.} A global
characteristic radius of the pore space can be defined as the volume
of pore space divided by its area $R_c=\phi
V_\text{bulk}/A_\text{P}$. The porosity, specific internal surface,
and the characteristic radius are global properties of porous
media. They do not capture local variations in the pore space which
are characterized by local porosity and pore size. In order to
quantify the local pore characteristics, the pore space is divided
into individual pores by defining pore throats at the local minima of
the hydraulic radius ratio $R_h=S/P$, \ie, the ratio of the
cross-sectional area $S$ and the perimeter $P$ of the cross section
\cite{Kwiecien1990,Dullien1992}. The partition of pore space allows to
define a pore size distribution function $\Pi(r)$ as the probability
of finding a pore with size smaller than $r$. An alternative to pore
size distributions are local porosity distributions which are based on
partitioning the porous medium into localized measurement cells
\cite{Hilfer1996}. The definition of specific internal surface carries
over naturally to local specific internal surface distributions. Since
the local porosity distribution depends on the partitioning into
measurement cells, the size of the measurement cells can be varied in
order to determine a correlation length up to which local variations
are important. Beyond the correlation length the porous medium can be
treated as homogeneous, however, if the porous medium is heterogeneous
on all scales, a finite correlation length does not exist.
While the porosity and pore size distributions describe local
morphology in terms of volume and area ratios, they do not capture the
connectivity of the pore space. The latter is important for transport
phenomena that are primarily dependent on the connectedness of the
pores by pore throats. The topological properties can be characterized
by pore coordination number and the pore coordination spectrum, \ie,
the distribution of pores connected by a certain number of pore
throats. If two opposite surfaces of the porous medium (or the local
measurement cell) are connected by a path through the pore space, the
porous medium is referred to as (locally) percolating. The percolation
threshold depends on the specific partitioning and percolation model
\cite{Hilfer1996}.
%
A simplified measure of the topology of porous media is the
\emph{tortuosity} that characterizes the length of the pathways
through pore space. It is defined as the ratio $T=L_i/L$ of the length
$L_i$ of a complex path through a porous medium over a straight-line
distance $L$. The tortuosity is strongly dependent on the properties
of pore space, including porosity, specific surface, connectivity and
their respective distributions, in a non-trivial way
\cite{Panda1995,Garrouch2001}. Hence in practice a number of empirical
models have been suggested that are often specific to certain types of
porous media and experimental conditions.
\subsection*{Simplified geometry models}
The interpretation of experimental measurements is often based on
simplified models of porous media such as the capillary tube model. In
the capillary tube model, the pore space is modeled as an array of
cylindrical tubes that do not intersect each other. The porosity of
the capillary tube model is $\phi=\pi/L^3 \sum_{i=1}^{N} L_i a_i^2$,
where N is the number of tubes in a porous medium of dimension $L$,
and $L_i$ and $a_i$ are the length and the radius of the $i$-th tube,
respectively. Similarly, the specific internal surface is given by
$S=2\pi/L^3\sum_{i=1}^{N} L_i a_i$. If the lengths and radii are
randomly distributed, the average porosity and average specific
internal surface can be written as
%
\begin{align}
\langle \phi \rangle &= \pi \frac{N}{L^2} \langle T \rangle \left\langle a^2 \right\rangle \\
\langle S \rangle &= \pi^2 \frac{N}{L^2} \langle T \rangle \left\langle a \right\rangle ,
\end{align}
%
where $\langle T \rangle$ is the average tortuosity and $N/L^2$ is the
number of capillaries per unit cross-section.
Another simple model is the spherical grain model which models either
the matrix or the pore space by randomly placed spheres. For a matrix
of non-overlapping spherical grains, the porosity and specific
internal surface are given by
%
\begin{align}
\langle \phi \rangle &= 1 - \frac{4}{3} \pi R^3 \rho \\
\langle S \rangle &= 4\pi R^2 \rho ,
\end{align}
%
where $R$ and $\rho$ are the radius and number density of the spheres,
respectively. Instead of randomly placing the spheres, regular
arrangements on simple cubic, bcc, and fcc lattices of spheres can
also be used.
\subsection*{Experimental characterization of porous media}
There are various ways to characterize porous media
experimentally. Two common approaches are porosimetry and
imaging. Mercury porosimetry is arguably the most common method to
quantify the geometry of pore space \cite{Hilfer1996,Xiong2016}. It is
based on the injection of mercury into the evacuated pore space as a
function of external pressure. Since mercury has a high surface
tension and is highly nonwetting, larger pores are filled first and
the amount of injected mercury is an indicator of the fraction of
pores down to a certain size that are filled. The pore size
distribution can be extracted by relating the size of the pores to the
capillary pressure of mercury through the Lucas-Washburn relation
\cite{Lucas1918,Washburn1921}. This interpretation of the measured
injected volume is based on the capillary tube model and thus assumes
that all pores are equally accessible. Hence mercury porosimetry can
only be used to measure the open pore space and does not provide any
connectivity information. More complex capillary network models can be
used to provide more accurate pore size distributions for porous media
with complex microstructure, \cf \cite{Xiong2016}.
The surface area of pore space can be probed by gas
adsorption \cite{Sing2001}. Like mercury porosimetry it is based on
surface tension and pressure, however, in gas adsorption smaller pores
are probed first. Depending on the targeted pore size and ambient
conditions, adsorptives such as nitrogen, argon or carbon dioxide are
used to determine the physical adsorption isotherms. \textcolor{red}{A number of
theories are available to interpret the measured isotherms, \eg, the
Brunauer-Emmett-Teller (BET) \cite{Brunauer1938} or Barret-Joiner-Halenda
(BJH) \cite{Barrett1951} method.} The interpretation methods rely on specific
assumptions on the shape of the pores and the nature of the adsorbate,
therefore the accuracy of the extracted pore size distribution can be
limited for porous media whose microstructure differs from the
assumptions \cite{Xiong2016}. For complex and hierarchical
microstructures, newer interpretation models have been developed that
use simulations methods such as molecular dynamics and density
functional theory \cite{Meakin2009}. These methods promise to improve the
accuracy for a wider range of pore sizes and geometries \cite{Sing2001}.
\subsubsection*{Imaging techniques}
Modern imaging methods are capable of capturing three-dimensional
representations of the full geometry and topology of porous media. For
instance, X-rays allow to scan the internal structure of a
three-dimensional sample by rotating the sample and recording the
absorption in all directions. The standard technique for this type of
non-destructive three-dimensional microscopy is X-ray computed
micro-tomography (micro-CT) \cite{Ketcham2001,Blunt2013}. The image
resolution of micro-CT scanners depends on the beam set-up, the
detectors, and the sample size. Laboratory systems have a typical
resolution of up to 50 \textmu{}m, and some systems are capable of
providing submicron resolution in the range of 50 nm
\cite{Xiong2016}. Synchroton sources can also be used for CT imaging
and provide resolutions around 1 \textmu{}m. The advantage of
synchroton based CT is that monochromatic beams are available and the
beam is collimated, thus reducing aberration artifacts. After
acquisition the recorded absorption data is processed and a stack of
two-dimensional images is generated. This stack can be segmented to
obtain a three-dimensional representation of the solid matrix or
pore-space. Three examples of pore-space representations obtained from
micro-CT images of porous carbonate samples are shown in
Fig.~\ref{fig:pore-network}. The main limitation of micro-CT is the
available resolution, and the reconstructed three-dimensional
representation may miss the smallest details of the pore space
\cite{Blunt2013}. In addition, insufficient contrast in the
absorption data can lead to inaccurate representation of features in
the reconstructed images. Some image artifacts such as noise can be
accounted for during processing by applying filters and thresholds
that have been calibrated for the specific CT setup. However, micro-CT
can not provide any information beyond the sample size such that
heterogeneities on larger scales will not be captured.\selectlanguage{english}
\begin{figure}
\includegraphics[width=.49\textwidth]{{{1-s2.0-S0309170812000528-gr2}}}
\includegraphics[width=.49\textwidth]{{{1-s2.0-S0309170812000528-gr3}}}
\caption{{Visualization of pore space and corresponding pore-network
representations for three porous carbonates (Estaillades, Ketton,
Mount Gambier). The pore space images were reconstructed from
three-dimensional micro-CT images of the samples. The networks
show the pores as spheres and the pore throats as cylinders, where
the sizes are derived from the pore-space images. Reproduced from
\cite{Blunt2013}. Permission pending.}}
\label{fig:pore-network}
\end{figure}
Higher resolution in the 10 nm range can be achieved using scanning
electron microscopy (SEM). Most porous media can sufficiently be
resolved at this resolution, however, only two-dimensional images can
be captured and no information about the three-dimensional topology
can be extracted. This can be mitigated by combining SEM with focused
ion beams (FIB). The FIB is used to cut the sample into fine slices
which can be imaged with SEM to provide small scale pore features. The
information provided by the combined FIB/SEM can augment data obtained
from micro-CT to generate a three-dimensional registration of pore
space. FIB/SEM is a destructive imaging technique, but it can provide
the best resolution of pore space covering the geometry and topology
in the whole range from cm to nm \cite{Xiong2016}. The combination of
different imaging techniques paves the way to more realistic geometry
representations that can be integrated with multiscale modeling and
simulation.
%\begin{itemize}
%\item{Granular materials}
%\item{previous work on granular media(electrokinetic)}
%\item{Fibrous materials}
%\item{previous work on fibrous media}
%\item{Packed beds}
%\item{Rocks etc.}
%\end{itemize}
%A measure of the geometric complexity of a porous medium. Tortuosity is a ratio that characterizes the convoluted pathways of fluid diffusion and electrical conduction through porous media. In the fluid mechanics of porous media, tortuosity is the ratio of the length of a streamline—a flow line or path—between two points to the straight-line distance between those points.
%Tortuosity is thus related to the ratio of a fluid's diffusion coefficient when it is not confined by a porous medium to its effective diffusion coefficient when confined in a porous medium.
%Tortuosity is also related to the formation factor, which is the ratio of electrical resistivity of a conductive fluid in a porous medium to the electrical resistivity of the fluid itself.
%In some literature, tortuosity denotes the square of the ratio defined above, whereas in other literature, the term tortuosity factor is used for the square of the ratio.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Multiscale modeling of porous media}
\label{sec:modeling}
The central objective of modeling porous media is to be able to
predict the phenomenological and quantitative behavior of a porous
medium on macroscopic scales, \ie, the scales that are observed in
laboratory or field experiments. Due to the complex geometry of pore
space and the hierarchy of length scales, a fully resolved model of
the porous medium is usually not feasible. Therefore, macroscopic
models of porous media use effective properties such as the
permeability in Darcy's law. The equations governing the macroscopic
variables arise from averaging over the microscopic features of pore
space. They are thus a coarse-grained representation of the
microscopic dynamics, much like coarse-grained models of complex
molecules arise from averaging over fast degrees of freedom. The
difference is that in the latter case the average is taken in time,
while for porous media the average is performed over a spatial domain.
The averaging process is referred to as \emph{homogenization} and is
based on the notion of a \emph{representative elementary volume}
(REV). The REV is a subdomain of the porous medium which contains
sufficient information such that the average yields a representative
effective property. The average of a quantity $\theta$ over
phase $i$ is written as
%
\begin{equation}
\langle \theta_i(\rvec) \rangle = \int_\text{REV} \theta(\rvec') \chi_i(\rvec') \omega(\rvec'-\rvec) d\rvec' ,
\end{equation}
%
where $\chi_i(\rvec)$ is the indicator function for phase $i$, and
$\omega(\rvec)$ is a suitable weight function. Standard choices for
the weight function are a constant or Gaussian function. The size of
the REV is a mesoscopic length scale beyond which the variance of the
effective properties decays and the average becomes a constant. For
instance, averaging Darcy's law over REV with the local pressure
gradient $\nabla P(\rvec)= s(\rvec) \nabla P$ yields
%
\begin{equation}
\langle \vvec \rangle = - \frac{1}{\eta} \left\langle s k \right\rangle \nabla P ,
\end{equation}
%
where $K=\langle s k \rangle$ is the effective permeability. Note that
in order to determine the effective permeability, the statistical
distribution of the microscopic permeability $k(\rvec)$ and the
scaling factor $s(\rvec)$ for the local pressure gradient need to be
known. The process of transitioning from the microscale description to
the coarse-grained description by determining the effective
permeability is called \emph{upscaling}. The microscopic
heterogeneities are absorbed into the effective properties and the
porous medium becomes homogeneous on the coarse-grained scale. Various
upscaling methods exist \cite{Chen2006,Renard1997}, however, in
many situations they face difficulties regarding accuracy and
robustness. In recent years, upscaling has increasingly become
augmented by numerical methods such as the heterogeneous multiscale
methods (HMM) \cite{E2003}. The development of novel multiscale
methods has led to robust, accurate, and efficient procedures that
outperform the classical upscaling techniques in many practical
applications. Today, the determination of effective properties and
macroscopic flow predictions thus relies more and more on numerical
multiscale modeling.
\subsection*{Pore network models}
A simplified representation of pore space consists of a network of
idealized pores connected by pore throats. The pores and throats
constitute the nodes and bonds of a pore network and are assumed to
have simple geometric shapes for which the permeability (and other
properties) is known. In this sense, pore network models are a
coarse-grained representation that retains the essential features on a
mesoscale. Instead of averaging Darcy's equation, the flow through the
network elements is described by a discretized Laplace equation
%
\begin{equation}\label{eq:pore-network}
\sum_j \frac{k_{ij}}{\eta} (P_i - P_j) = 0 ,
\end{equation}
%
where $P_i$ is the pressure $P(\rvec_i)$ at element $i$ of the
network, and $k_{ij}=k((\rvec_i+\rvec_j)/2)$ is the permeability of
the throat connecting $i$ and $j$. The discrete equations
\eqref{eq:pore-network} can be solved using, \eg, successive
overrelaxation or conjugate gradient methods. The boundary conditions
give rise to a source term on the right hand side if $i$ is a boundary
node. Pore network models date back to the seminal work of Fatt
\cite{Fatt1956} who used a random resistor network to study two-phase
capillary displacement. Since then, significant advances have been
made and a variety of more sophisticated network models have been
developed. On the one hand, these aim at including more details of the
pore scale physics such as wettability and relative permeability, or
improved physical models for non-Newtonian displacement and non-linear
extensions of Darcy's law \cite{Blunt2013}. There are also ongoing
efforts to develop pore-network models for electrokinetic flows
\cite{Obliger2013,Obliger2014}. On the other hand, better
representations of the geometry and topology of pore space have been
sought. These include more accurate modeling of the morphology of the
throats in order to capture the influence of surface wetting in
multiphase flows \cite{Blunt2002,Blunt2013,Xiong2016}. The main aim,
however, is to accurately represent the spatial correlations and
connectivity of pore space which are not captured by statistically
mapped networks. For this reason, randomly generated representations
of porous media are often not predictive \cite{Mehmani2015} and are of
limited value in multiscale modeling.
\subsubsection*{Network reconstruction}
Constructing the pore network for a particular representation of a
porous medium is a key element in pore-network modeling. The first
realistic models for granular media were based on simulated grain
packing and diagenesis \cite{Bryant1992,Bryant1993a,Bryant1993b}. Pore
centers are identified in the void space and the topology is created
by connecting four or less adjacent pores. However, grain based
network reconstruction is limited to porous media that consist of
spherical grains.
%
Alternative approaches are based on reconstruction of porous media
from experimentally acquired images of pore space. The maximal ball
method places spheres on every voxel and grows them until they touch
the solid matrix. The largest inscribed spheres are identified as
pores while smaller spheres represent pore throats. While this method
can identify the pores, it tends to produce high coordination numbers
and the pore throats may include cascades of small
spheres \cite{Xiong2016}.
%
In contrast, the erosion-dilation method and the medial axis method
use certain transformations to reduce the pore space to its
topological skeleton. The pores are identified based on a tesselation
of the sample space. While these methods preserve the topology of pore
space, the identification of pores remains ambiguous and is sensitive
to resolution and image quality
\cite{Xiong2016}. Fig.~\ref{fig:pore-network} shows three examples of
pore-network representations consisting of pores and pore throats
whose sizes were determined from the corresponding micro-CT images.
%
In general, the network reconstruction methods suffer from the lack of
a rigorous distinction of pores and pore throats. The reconstructed
morphology may depend on this definition and on the specific
discretization and resolution. Therefore, different pore-network
models can be constructed for the same image data, and it can only be
validated a posteriori whether a representation yields accurate
physical predictions. Hence the comparison to results from direct
models and experimental data is essential to determine whether the
pore-network predictions of the relevant properties are sensitive to
the particular reconstruction strategy \cite{Xiong2016}.
\subsection*{Direct models}
Direct models have the advantage that they do not seek a simplified
representation of the porous medium but use the voxelized
three-dimensional geometry of pore space to solve the governing
equations on a discrete grid. Since the only limitation of this
approach is the resolution and quality of the image, direct models are
ideally suited to investigate the physics of porous media down to the
pore scale \cite{Biswal2009}. In principle, any numerical method can
be applied to solve the governing equations and classical approaches
include finite-difference, finite-element methods, and density
functional methods \cite{Meakin2009}. However, most of these methods
are plagued by the difficulty of incorporating the complex boundary
conditions of pore space. In addition, the treatment of interfaces in
multiphase systems requires additional effort in order to accurately
locate and track the interface. The lattice Boltzmann method (LBM), on
the other hand, allows incorporation of arbitrary boundary shapes with
relative ease and does not require any interface tracking in
multiphase systems. Therefore, the LBM is arguably the most popular
method for direct modeling of porous media.
\subsubsection*{Lattice Boltzmann methods}
The lattice Boltzmann method is a lattice kinetic method that
reproduces the incompressible Navier-Stokes equation on hydrodynamic
length and time scales. Unlike particle-based methods, it represents
the fluid by a set of discrete (mass) distribution functions on
regular cubic lattice. These distribution functions are updated
according to the lattice Boltzmann equation (LBE)
%
\begin{equation}\label{eq:lbe}
f_i(\mathbf{x}+\mathbf{c}_i\, h, t+h)
%=f_i^\star(\mathbf{x},t)
=f_i(\mathbf{x}, t) - \sum_j \Lambdaop_{ij} \left[ f_i(\xvec,t) - f_i^\eq(\xvec,t) \right] + \sum_j \left( \delta_{ij} - \frac{1}{2} \Lambdaop_{ij} \right) g_j(\xvec,t) ,
\end{equation}
%
where $\cvec_i$ are velocity vectors connecting the discrete lattice
points $\xvec$, and $h$ is a discrete time step. The collision matrix
$\Lambdaop_{ij}$ relaxes the distributions $f_i$ towards their
equilibrium $f_i^\eq$. The equilibrium $f_i^\eq$ and the force term
$g_i$ are given in terms of low-velocity expansions
%
\begin{align}\label{eq:lbm-eq}
f_i^\text{eq}(\rho, \mathbf{u}) &= w_i\,\rho\left(1 + \frac{\mathbf{u}\cdot\mathbf{c}_i}{c_\text{s}^2} + \frac{\vec{u}\vec{u}:(\vec{c}_i\vec{c}_i-c_\text{s}^2\mathsf{1})}{2c_\text{s}^4} \right) ,\\
g_i &= w_i \left(\frac{\vec{F}\cdot\vec{c}_i}{c_\text{s}^2} + \frac{\vec{u}\vec{F}:(\vec{c}_i \vec{c}_i - c_\text{s}^2 \mathsf{1})}{2c_\text{s}^4}\right) ,
\end{align}
%
where $c_s$ is a (pseudo-)speed of sound that depends on the specific
lattice model. The density $\rho$, momentum density
$\rho \mathbf{u}$, and momentum flux $\mathsf{\Pi}$ of the fluid are the
zeroth, first, and second moments~\cite{Duenweg2009}
%
\begin{align}
\rho &= \sum_i \, f_i, & \rho \mathbf{u} &= \sum_i \, \mathbf{c}_i f_i + \frac{1}{2} \vec{F}, & \mathsf{\Pi} &= \frac{1}{2} \sum_i \left( f_i+f_i^\star \right)\vec{c}_i\vec{c}_i ,
\end{align}
%
where the momentum flux $\mathsf{\Pi}=\mathsf{\Sigma}+\mathsf{\sigma}$ contains
the Euler stress $\mathsf{\Sigma} = \rho \mathbf{u} \mathbf{u} +
p \mathbf{I}$ and the deviatoric stress tensor $\mathsf{\sigma} = \rho \nu
[\nabla \mathbf{u} + (\nabla \mathbf{u})^\top]$. A popular
choice for the collision matrix is the BGK approximation
$\Lambdaop_{ij}=\tau^{-1} \delta_{ij}$ in which case the viscosity is
$\nu=c_s^2(\tau-h/2)$.
%The link between the LBE (Eq.~\ref{lbe}) and the macroscopic
%Navier-Stokes equation can be established through a Chapman-Enskog
%expansion \cite{Succi:2001, Duenweg:2009, kruger2016lattice}.
The geometry of pore space can be modeled on the grid by marking the
nodes that are located within the solid matrix. At the boundary
between fluid nodes and solid nodes, the LBE \eqref{eq:lbe} is
augmented by the bounce-back boundary conditions which reflects a
distribution when it impinges the boundary, \ie,
$f_{\bar{\imath}}(\mathbf{x}, t + h) = f^\star_i(\mathbf{x}, t)$, where
$\mathbf{c}_{\bar{\imath}} = - \mathbf{c}_i$ denotes the reversal of
direction. More involved boundary conditions have been developed to
improve accuracy \cite{Bouzidi2001,Ginzburg2003}, however, the bounce-back
rule remains arguably the most widely applied boundary condition.
A multiphase mixture can be modeled in the LBM by including a second
set of distributions that evolves on the same lattice. A short-range
interaction force is applied between component $\sigma$ and $\sigma'$,
which is commonly written in terms of a Shan-Chen pseudo-potential
\cite{Shan1993,Shan1994}
%
\begin{equation}
\mathbf{F}^\sigma(\mathbf{x}) = -\psi^\sigma(\mathbf{x}) \sum_{\sigma'} g_{\sigma \sigma'} \sum_{i} w_i \psi^{\sigma'}(\mathbf{x} + \mathbf{c}_i h) \mathbf{c}_i h
\end{equation}
%
where $g_{\sigma \sigma'} = g_{\sigma' \sigma}$ is the interaction
strength between components $\sigma$ and $\sigma'$, and the
pseudo-potential $\psi^\sigma(\mathbf{x})$ is a typically linear or
exponential function of the component density
$\rho^\sigma(\mathbf{x})$. The short-range interactions induce a
surface tension at the interface between the components, and the
interaction strength $g_{\sigma\sigma'}$ controls the phase separation
of the components.
%
Thanks to the kinetic basis and the mesoscopic scale of the LBM, it
can be coupled to a range of other physical dynamics if they are
amenable to a lattice representation on a similar scale. For instance,
Capuani et al. \cite{Capuani2004} have developed a link-flux algorithm
to solve the electrokinetic equations on the same lattice and couple
the ionic fluxes to the hydrodynamic flow field. This approach can be
applied to simulate pore scale electrokinetic flow in porous media
\cite{Rotenberg2008,Rotenberg2013}.
%
For more comprehensive reviews of various lattice Boltzmann methods,
we refer the reader
to \cite{Duenweg2009,Schmieschek2017}.
\subsubsection*{Geometry generation and reconstruction}
There are essentially two ways to obtain the porous geometry for
direct modeling. The first is based on computer-generated structures
that mimic the morphology and topology of a certain type of porous
medium. The seminal example are grain-based models for diagenesis and
compaction of granular media and sandstones
\cite{Bryant1992,Bryant1993a,Bryant1993b}. They model the natural
formation process by growing a random packing of spheres until they
overlap and percolate. The packing is compacted by moving the centers
of the spheres closer together in the direction of compression. Grain
size distributions can be modeled by using spheres of different
radius. While grain-based models where originally used to construct
pore-network models, they can be used straightforwardly in direct
models. Instead of spheres, other basic shapes such as ellipsoids,
cubes, parallelepipeds, and fibers can be used.
Random arrangements of fibers can be used as models for paper,
textiles, and non-woven materials. For instance, Koponen et
al. \cite{Koponen1998} have used random fiber sheets to study the
permeability of paper and non-woven fabrics. Ordered and disordered
fiber arrangements were compared by Claque et al. \cite{Clague2000} in
wall-bounded and unbounded flows. Nabovati et al. have studied the
influence of fiber radius, curvature, and aspect ratio on the
permeability of non-overlapping and overlapping random fiber
arrangements. In our research, we employ a workflow using the GeoDict
software \cite{GeoDict2015} to generate models of non-woven fibrous
media. The fibers are placed randomly but only a limited range of
directions is allowed, giving rise to an anisotropic arrangement that
mimics non-woven mats produced by electrospinning
\cite{Chen2017}. Both \textcolor{orange}{overlapping} and non-overlapping arrangements
can be created to cover a range of porosities. Two examples of
anisotropic, random fiber arrangements are shown in
Fig.~\ref{fig:random-fibers}. An animated visualization of the
three-dimensional structures is available in the supplementary
materials. In contrast to stochastic models for porous media,
process-based models for generation of porous media yield a better
representation of the porous medium beyond porosity. While random
generation can recover the local porosity, it typically fails to give
an accurate representation of the connectivity which can significantly
deteriorate the predictive capabilities of stochastic models. Modeling
of the formation process gives a more realistic representation of pore
space, however, the process-based simulation can be difficult to
realize for specific types of porous media.\selectlanguage{english}
\begin{figure}
\includegraphics[width=.49\textwidth]{{{r18_phi_0.796_new}}}
\includegraphics[width=.49\textwidth]{{{r12_phi_0.695}}}
%\movie[externalviewer]{%
% \includegraphics[height=6\baselineskip]{L_200_R_12}
%}{avi/r12_svp30.mpg}
%\movie[externalviewer]{%
% \includegraphics[height=6\baselineskip]{L_200_R_12}
%}{avi/r18_svp20.mpg}
%\includemedia{\includegraphics{L_200_R_12}}{avi/r18_svp20.swf}
%\includemedia{\includegraphics{L_200_R_12}}{avi/r12_svp30.swf}
%\includemedia[
% width=0.4\linewidth,
% height=0.3\linewidth,
% activate=pageopen,
% addresource=img/r18_svp20.mp4,
% flashvars={source=img/r18_svp20.mp4}
%]{}{VPlayer.swf}
%{\flushright{\footnotesize Generated with GeoDict.}}
\caption{{Fibrous membranes with varying porosity generated with the
software GeoDict \cite{GeoDict2015}. (left) Non-overlapping fibers
with a porosity of $\phi=0.796$. (right) Overlapping fibers with a
porosity of $\phi=0.695$.}}
\label{fig:random-fibers}
\end{figure}
The second approach is based on experimental imaging techniques such
as micro-CT/XRCT and FIB/SEM. Image reconstruction can deliver the
most realistic representations of porous media by digitizing the full
three-dimensional pore space. Heterogeneities are naturally included
and the reconstructed geometry yields accurate pore spectra and
connectivity information. However, the extraction the necessary
information from the experimentally recorded images is not trivial and
requires careful filtering and segmentation techniques. Gueven et
al. \cite{Gueven2017} have digitized CT images of porous sintered
glass bead packings. They applied several filters to the raw data
before segmenting the geometry. Filtering is necessary to adjust the
gray-scale distribution such that the glass beads and the pore space
can reliably be identified. The filters included the removal of peaks
(bright spots), phase contrast improvement to enhance the edges of the
glass beads, and noise reduction and thresholding to generate a
binarized image. The resulting binary map was used to analyze the
particle size distribution, pore throat distribution, porosity, and
tortuosity. The workflow for this image analysis is shown in
Fig.~\ref{fig:segmentation}. The voxelized geometry was used to perform
lattice Boltzmann simulations to determine the permeability of the
sample and in different REVs. The resolution of the CT images was 16
\textmu{}m and domains of up to $1024\times1024\times2048$ voxels were
considered. The results demonstrate that direct modeling can predict
the local variations in permeability with good accuracy as long as a
representative sample size is used.\selectlanguage{english}
\begin{figure}
\includegraphics[width=.99\textwidth]{10035_2017_705_Fig4_HTML}
\caption{{Flowchart of the image analysis workflow employed by Gueven
et al. \cite{Gueven2017} based on XRCT images of sintered glass bead
systems. Several filters are applied to the raw images in order to
obtain a binary gray scale distribution. The binarized image was
segmented to analyze the bead distribution. The pore space was
analyzed by inverting the binarized image. Reproduced from
\cite{Gueven2017}. Permission pending.}}
\label{fig:segmentation}
\end{figure}
The use of FIB/SEM image reconstruction for direct modeling was
demonstrated by Roberts et al. \cite{Roberts2016}. They binarized
images of lithium-cobalt-oxide and segmented the active particle,
binder, and electrolyte phases. The identification of individual
particles was facilitated by a watershed algorithm. The aim of the
work was to perform finite-element simulations of the coupled
electrochemical-mechanical-thermal equations governing the
intercalation of lithium in the porous material. Therefore, instead of
using the voxel data directly, a tetrahedral mesh was created and
decomposed into a phase conforming mesh according to the segmentation
of the original image. The smallest mesh size used by Roberts et
al. was 0.0625 \textmu{}m resulting in more than $100\times10^6$
elements, and second-order accuracy of the results was verified. An
analysis of the required REV of the reconstructed images is also
possible \cite{Kanit2003}.
These examples serve to show that the progress in imaging techniques
makes it possible to reconstruct high-fidelity direct models of porous
media. Limitations in image quality can partly be overcome using
digital image processing techniques, and direct modeling can reach the
experimental resolution provided the numerical techniques are
sufficiently efficient and computationally feasible.
%\begin{itemize}
%\item Hybrid models
%\end{itemize}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Conclusions and future directions}
\label{sec:conclusions}
In this contribution we have reviewed the development of modeling
approaches for transport phenomena in porous media. Advances in this
research area are driven by the need to better understand the physical
mechanisms that are at play across the scales in porous media. These
mechanism impact the performance of important technological
applications such as water filtration, hydrocarbon recovery, and
energy storage. The heterogeneous structure of porous media is of
particular importance and leads to a complex interplay of multiscale
and multiphysics effects.
Modern imaging methods such as micro-CT/XRCT and FIB/SEM are now
capable of probing the geometry of pore space with sub-micron
resolution. At the same time, advances in high performance computing
have made simulations with millions and billions of voxels
feasible. The combination of these developments enables direct
simulation of transport in porous media with unprecedented accuracy,
and direct modeling of sample sizes on the scale of laboratory
experiments seem almost possible. Fig.~\ref{fig:workflow} illustrates
an integrated workflow that combines multiscale modeling and
simulation with the reconstruction of realistic geometries from
experimental images.\selectlanguage{english}
\begin{figure}
\includegraphics[width=.99\textwidth]{workflow}
\caption{{Illustration of an integrated workflow for multiscale
modeling of porous media. Experimental images are used to
characterize the porous medium and to generate or reconstruct
realistic geometries. Pore-scale simulations are used to determine
the local permeability and other effective properties. These
inform theoretical analysis and upscaling approaches to model the
porous medium on larger scales. The whole process is iterative and
two-way communication between smaller and larger scales is used in
hybrid methods. High performance computing enables multiscale
multiphysics simulations of porous media with improved predictive
capabilities.}}
\label{fig:workflow}
\end{figure}
However, direct modeling is not without challenges. One problem that
remains despite the increased resolution of experimental imaging is
the limited ability to obtain information on the surface
characteristics at the pore scale. For instance, the wettability of
the porous medium may be heterogeneous and the contact angle may vary
in different pores. Similarly, surface charges may vary throughout the
porous medium and give rise to a heterogeneous thickness of the Debye
layer. Inferring the wettability or surface charges on the pore scale
remains difficult and hence simulations typically use rather generic
model assumptions.
%
Another concern is the validity of the physical models across the
scales. Most governing equations contain certain model parameters that
have to be determined from constitutive relations or by some other
means. Such relations are typically based on simplifying assumptions
and the models are virtually always imperfect. For instance, Darcy's
law breaks down at very low or very high porosities, and the linear
dependence of the flux on the pressure gradient is only valid under
creeping flow conditions. Direct modeling offers the opportunity to
validate simulation models by comparing numerical predictions and
experimental measurements. Multiscale modeling is particularly suited
in this context since the same observable quantity can be tackled with
different methods and at different scales. This allow to determine
effective parameters and devise coarse-grained models that can be used
at larger scales. However, this typically requires statistical
analysis of considerable amounts of data from different scales, and
effective parameters usually carry some uncertainty. Approaches from
data science and machine learning can prove valuable to incorporate
all available information. This will also require a better integration
of modeling approaches and data analysis, since a purely data driven
approach is prone to overfitting and limited in the physical insight
it can provide.
A particular challenge arises when the separation of scales breaks
down, \ie, when the macroscopic flow dynamics feeds back into the
interactions at the pore scale and the system becomes dynamically
coupled. This occurs, \eg, in reactive flows where the molecular
kinetics is limited by mixing on the pore scale which depends on flow
rates. Overlap of scales can also be promoted by the presence of
microscopic heterogeneities \cite{Lu2005}. If no scale separation
exists, a hierarchical approach that samples the pore scale dynamics
to extract effective properties for the macro scale is not
appropriate. This limitation is addressed by \emph{hybrid multiscale
methods} which simulate a given system on different scales
simultaneously and implement a two-way coupling between the two
instances to exchange information between the micro- and the
macroscale. An example of a hybrid approach is the heterogeneous
multiscale method (HMM), which uses a small scale model to generate
data from which constitutive relations for a macroscale model are
extracted on the fly.
Novel hybrid and multiphysics models are required to address systems
with a dynamically evolving pore structure. A change in the geometry
may be due to chemical reactions leading to erosion or deposition at
the surface, or due to mechanical deformation or fracture of the solid
matrix. In both cases, the evolution of pore space needs to be
included in the model which leads to multiphysics models that
potentially broaden the spectrum of relevant length and time
scales. Implementing consistent couplings between the deformation
mechanics of the solid matrix and reactive flows in porous media
relies on accurate incorporation of reaction rates and kinetics on the
pore scale, which typically increases the computational costs of the
overall numerical algorithm.
Whereas hybrid multiscale methods are computationally challenging due
to the two-way communication between models that operate at disparate
scales, they will likely gain in popularity on emerging exascale
computing systems. It is anticipated that exascale systems will be
massively parallel and highly heterogeneous, such that novel computing
patterns are required which may prove particularly fruitful for hybrid
algorithms \cite{Alowayyed2017}.
%
Overall, modeling of porous media continues towards more sophisticated
models as more complex applications emerge, \eg, in biomedical
engineering \cite{Khanafer2006}. The combination of advanced imaging
techniques with multiscale modeling facilitates numerical simulation
of realistic geometries which can help turn simulation models of
porous media into predictive tools of wide applicability.
%\begin{itemize}
%\item limitations
%\item identification of physics from images
%\item models
%\item data
%\item multiscale
%\item Hybrid models
%\item Challenges from Blunt
%\end{itemize}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Acknowledgments}
This work has been supported by a faculty start-up grant from the College of Engineering, Computing and Applied Sciences at Clemson University. Clemson University is acknowledged for generous allotment of compute time on Palmetto cluster.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\selectlanguage{english}
\FloatBarrier
\bibliographystyle{unsrt}
\bibliography{bibliography/converted_to_latex.bib%
}
\end{document}
__