pkuthss/doc/example/chap/chap1.tex

% Copyright (c) 2014,2016,2018 Casper Ti. Vector
% Public domain.

\chapter{绪论}

\section{研究背景}
各向异性扩散问题在众多科学领域均有广泛应用, 例如辐射流体动力学(Radiation Hydrodynamics, RHD), 磁流体动力学(Magneto-Hydrodynamics, MHD), 油藏模拟以及等离子体物理学.
在这些应用领域中, 扩散方程通常与其他复杂的物理模型耦合求解\parencite{Yuan2009}, 这使得设计扩散方程的数值方法颇具难度.



例如, 在用于油藏模拟的多孔介质模型中\parencite{Bear2013}, 扩散系数代表介质的渗透率, 通常是一个各向异性、非均匀不连续的张量.
多孔介质中的两相流模型在油藏工程中起着至关重要的作用, 例如, 在二次采油中, 注入水以从油藏中驱油.
在某些简化假设下, 描述石油储层中流体流动的数学模型包括压力的扩散方程和饱和度的非线性双曲方程.
其中, 压力的扩散方程是一个含有各向异性间断系数的二阶椭圆方程\parencite{Chen2006}.
在这样的应用背景下, 压力扩散方程与其他物理模型耦合形成多物理过程.


另一个例子是在多介质拉格朗日辐射流体动力学的数值模拟中,
多介质辐射流体力学中的扩散方程描述辐射在空间的传输过程\parencite{Lindl1995},
辐射扩散系数是温度和密度的强非线性函数.
在激波两侧, 密度是不连续的, 因此扩散系数具有很强的间断.
此外, 网格随流体移动, 会发生变形. 扩散过程与流体过程耦合, 直接使用变形网格求解扩散方程.
在这种情况下, 数值求解扩散方程需要仔细考虑离散格式可能引入的误差(参见参数锁定现象的相关文章\parencite{Babuska1992,Manzini2007}).

离散格式的设计应当基于实际应用和物理背景的需求.
理想的格式应尽可能多地具备优良特性,
比如经典的稳定性、局部守恒性、仅包含单元中心未知量的局部模板、简单性、鲁棒性、单调性、对称正定线性系统、允许任意扩散张量和网格等等.
然而, 据我们所知, 不存在一种格式能够满足上述所有特性.
在我们看来, 精度、鲁棒性以及允许任意扩散张量和网格是其中最为基础的特性.


在上述背景下, 为近似求解偏微分方程, 我们选择有限体积(Finite Volume, FV)方法.
长期以来, 有限体积格式在辐射流体动力学的数值模拟程序中得到了广泛应用\parencite{Ding2008,TRHD2015}.
特别是对于守恒律问题, 守恒型有限体积法更容易获得熵解\parencite{Lax-Wendroff1960,Hou1994}.




在多物理耦合问题中, 这类各向异性扩散方程的离散格式通常面临三大挑战: 
\begin{itemize}
  \item 多种多样的扩散系数. 在油藏模拟领域的多孔介质模型中, 扩散系数即介质的渗透率, 通常是各向异性、非均匀的、间断的张量. 介质渗透率的各向异性强度用最大特征值与最小特征值的比值定义, 在倾斜的岩层中各向异性现象非常明显.
  \item 大变形大长宽比的非结构网格. 在实际应用中计算区域通常十分复杂, 难以用单层单块结构网格描述, 需要使用多块结构网格拼接或者使用非结构网格. 小尺度的裂缝周围通常采用大长宽比的网格, 并且使用非结构拓扑关系过渡到其他区域.
  \item 三维情形计算量大. 扩散方程通常采用隐式求解, 需要求解大规模稀疏线性代数方程组. 例如对于三维结构网格, 如果每个方向剖分1000份, 那么总的网格单元数将达到10亿. 在三维情形, 计算效率和性能成为了一个重要的问题\parencite{FVCA6}.
\end{itemize}

从以上三方面挑战可以看出, 各向异性扩散方程的离散的主要问题是精度、健壮性和计算效率.
我们聚焦于辐射流体力学中的辐射扩散方程数值求解, 精度、健壮性和计算效率是我们研究的出发点.



在这样的应用背景下, 有限体积方法是一种常用的数值方法.
有限体积法不但能适用于结构和非结构网格, 而且能保持数值通量的局部守恒性, 即越过每一条网格边的法向流是连续的. 当所求的问题中通量起的作用十分重要时, 有限体积法更能显示出它的优越性.


出于理论和应用需求的驱动, 针对扩散问题的有限体积格式已开展了大量研究工作. 关于近期的发展动态,
可参阅文献\parencite{Handbook2000,Droniou2014,Camier2016,FVCA5,FVCA6}及其所引用的参考文献.
近年来, 数值模拟中的三维效应备受关注, 人们已提出了若干三维格式. 在此, 我们将重点关注三维单元中心有限体积格式,
因其每个单元仅有一个自由度, 所以在数值模拟中得到了广泛应用.


有限体积格式可以保持法向流的局部守恒,
在油藏模拟问题的压力扩散方程中, 扩散系数和压力梯度的乘积在网格面的法向投影称为法向流, 与流体速度在网格面的法向投影成正比, 而速度是两相流数值模拟中的重要物理量, 因此有限体积方法是我们关注的离散方法. 使用有限体积方法求解扩散方程已广泛应用于许多领域. 例如用有限体积方法求解辐射流体中的扩散方程, 需要求解流体、辐射、粒子输运等多个物理过程, 有限体积方法广泛地应用于该领域的数值模拟程序, 例如MARED\parencite{Ding2008}、TRHD\parencite{TRHD2015}以及CHIC\parencite{Breil2011}. 然而在油藏模拟商业软件中, 通常使用基于结构网格的差分格式, 例如CMG\parencite{CMG2019}以及ECLIPSE\parencite{GeoQuest2014}.


通常, 在油藏模拟器中, 使用简单的两点通量近似(Two-point flux approximation, TPFA)方法来求解压力扩散方程, 并且使用简单的一阶迎风方法来求解饱和度的双曲对流方程\parencite{Ewing1983, Arnold2002, Peaceman1977}. 尽管这种组合策略易于实现并且计算高效, 但它有某些缺点. 一方面, TPFA方法不能适当地模拟复杂几何形状, 这些复杂几何形状是从储层的断层和倾斜的地质建模得到的. 另一方面, 一阶迎风方法引入的数值耗散不仅使得饱和度前沿过度弥散, 也容易产生所谓的网格方向效应(Grid Orientation Effect, GOE), 也就是数值解对网格线方向的强烈依赖\parencite{Eymard2013, Lamine2010}. 为了解决这些问题, 文献中提出了许多策略, 涉及不同的数值方法组合, 例如袁益让和他的学术团队提出了一些稳定有效的数值方法, 并在胜利油田的实际生产中得到了应用\parencite{Yuan2014, Yuan2019}. 因此, 对各向异性多孔介质两相流模型应用更加鲁棒、更加高分辨的有限体积格式是一个很有意义的课题.



我们主要研究扩散方程有限体积格式, 因此, 精度、鲁棒性、计算效率是我们首要满足的三要素. 基于这三要素, 我们将聚焦于线性精确的单元中心型有限体积格式. 文\parencite{Wu2010,wjm2005}首次提出线性精确(Linearity-preserving)思想, 并用该思想重新推导了广泛应用于辐射流体数值模拟的九点格式\parencite{Li1980}. 线性精确思想作为一种在扭曲网格上构造扩散方程高精度有限体积格式的启示性途径, 要求当问题的精确解是分片线性函数并且扩散系数是分片常数时, 格式中的所有离散步骤精确成立. 许多研究工作证实\parencite{Shashkov1996,Gao2011,Contreras2016}, 线性精确格式对任意各向异性扩散系数以及任意非结构网格均具有二阶精度.

在过去的几年里, 线性精确准则已被用作在二维网格上推导一些单元中心格式的重要工具\parencite{Wu2010, Gao2011}.
这些工作的共同特点是, 定义在单元顶点处的辅助未知量的插值算法在格式的构造中起着重要作用,
这需要考虑顶点周围的几何拓扑结构, 因此使得将这些格式推广到三维情况成为一项困难的工作.
2012年, 文\parencite{Wu2012}中提出了一种全新的线性精确格式,
该格式利用了在单元边上的所谓调和平均点\parencite{Agelas2009}.
尽管这种格式很容易推广到三维网格, 并且其保极值的对应格式也很容易构造\parencite{Gao2013},
但新格式在某些单元边上可能不存在调和平均点, 这在一定程度上使这类格式的发展受到了影响.
文献\parencite{Zhang2020}已证明, 调和平均点可能位于单元面之外, 这会导致间断系数问题和各向异性扩散问题的精度降低. 因此, 在线性精确的单元中心格式中为辅助未知量寻求简单且稳健的插值算法是一个有趣的问题.


我们对于多孔介质两相流数值模拟, 采用经典的隐式压力显式饱和度法(IMplicit Pressure Explicit Saturation, IMPES)求解该方程组. IMPES方法最初由Sheldon等\parencite{Sheldon1959}和Stone等\parencite{Stone1961}提出. 它已广泛应用于求解多孔介质中两相流动的非线性耦合系统\parencite{Chen2004, Chen2019, Kou2010, Silva2016}. 在标准的IMPES方法中, 将饱和度归一化条件和达西定律代入各相的两个质量守恒定律并且将所得的方程求和, 就得到压力扩散方程. 对压力方程中除了压力之外的其他变量显式处理, 以消除非线性耦合, 然后隐式求解压力方程. 只要得到压力, 就可以显式更新达西速度和两相的饱和度. 重复这个过程直到模拟结束. 在IMPES算法框架下, 每个时间步的压力扩散方程和饱和度对流方程依次求解, 这就允许扩散格式和对流格式进行任意组合. 我们将使用线性精确的单元中心型有限体积格式求解各向异性多孔介质两相流问题中的压力扩散方程, 同时结合双曲守恒律中广泛使用的二阶单调MUSCL格式\parencite{Lamine2010, Xie2017, Contreras2016,Blazek2015}求解饱和度方程, 进而形成格式可扩展的两相流模拟器. 我们将为这类问题的数值模拟提供一些新的离散格式. 同时, 我们的研究成果可推广应用于例如金属铸造等两相驱替的工程问题.

综上所述, 在多介质辐射流体以及各向异性多孔介质两相流问题背景下, 对含有任意各向异性扩散系数的扩散方程, 发展适用于任意非结构网格的线性精确有限体积格式, 具有十分重要的理论意义和应用价值. 这也是开展我们研究的基本立意所在.


\section{国内外研究现状}


在理论和应用需求驱动下, 关于各向异性扩散方程的有限体积格式已有大量的研究工作. 2000年出版的专著\parencite{Handbook2000}系统地介绍了扩散方程有限体积方法, 文\parencite{Droniou2014}综述了目前常用的扩散方程有限体积格式. 连续两届“复杂应用有限体积方法国际研讨会”分别公布了二、三维扩散问题的典型算例和典型网格\parencite{FVCA6, FVCA5}.



有限体积格式按未知量定义的位置可以分为以下三类: 混合-杂交型, 节点型以及单元中心型. 我们主要关注单元中心型格式. 以下将综述国内外关于各向异性扩散方程单元中心型有限体积格式的进展.


\subsection{混合-杂交型格式}
混合-杂交型格式是指最终的离散系统包含单元中心未知量、节点未知量、边心未知量、面心未知量、面法向通量等未知变量的两种或两种以上的有限体积格式.

F. Hermeline在2000年提出一种对偶有限体积法(Discrete Duality FiniteVolume, DDFV)\parencite{Hermeline2000}, 该方法的基本思想是在原网格和对偶网格分别建立格式, 耦合求解. 在三维情形, DDFV格式分为两种类型, 一种是CeVe-DDFV\parencite{Hermeline2007,Hermeline2009}, 包含单元中心和顶点两种未知量, 另一种是CeVeFE-DDFV\parencite{Coudiere2011}, 包含单元中心、顶点、边心和面心四种未知量. 这类格式的优点是守恒、无条件满足强制性, 特别地, 如果网格与扩散系数满足正交性条件, 所得的线性代数方程组是对称正定的. 这类格式的缺点是未知量多, 计算量大, 并且在扩散系数间断情形精度较低, 需要用更为精确的方法重构梯度\parencite{Boyer2008}.

单元泛函极小化格式CFM\parencite{Yin2010}可视为文\parencite{Li1995}第3章第7节的全局泛函极小化格式的改进版本.
基于CFM的区域分解算法及分析可参见\parencite{Yin2014}. 文\parencite{Yin2015}提出一簇CFM格式,
满足线性精确性质并且所得线性方程组是对称正定的, 在弱星形区域假设下给出误差估计.
目前为止, 这类格式还没有三维版本.

支撑算子方法(Support Operators Method, SOM)在\parencite{Shashkov1995}中有详细的描述, 现在称为Mimetic有限差分方法\parencite{Lipnikov2014b}. SOM方法具有二阶精度, 但是形成稠密的对称正定代数方程组, 模板非局部, 且没有显式的通量表达式. 1998年, Morel等提出局部支撑算子方法LSOM\parencite{Morel1998}, 导出一个稀疏的对称正定线性方程组, 而且保持了二阶精度, 但其代价是引入网格边中点量, 计算量较大. 2001年Morel等发展了六面体网格的局部支撑算子方法\parencite{Morel2001}.

混合-杂交格式通常未知量多, 计算量大, 在三维情形, 必须使用并行计算提升计算效率.


\subsection{节点型格式}
我们提及一些成熟的格式, 例如Bank等提出的Box格式\parencite{Bank1987}, 李荣华等提出的有限体积元格式(广义差分方法)\parencite{Zhang2015,Li1999}, 这类格式不需要辅助未知量插值, 在结构网格上可以达到任意高阶. 邬吉明等提出线性精确节点格式VLPS\parencite{Wu2016}, 满足局部守恒、对称正定、线性精确以及无条件稳定. 董倩楠等\parencite{Dong2020}证明了VLPS格式与最低阶虚拟元方法在星形域上总刚度矩阵相同,
并为这两种方法提出统一的保正、守恒的后处理方法. 苏帅等将VLPS格式推广至三维并且做了保正修正\parencite{Su2019}.

目前, 节点型格式处于蓬勃发展阶段, 距离实用化仍然有一定距离. 对于多介质辐射流体力学问题, 常用的计算方法是将未知量定义在单元中心. 在大规模并行计算情形, 对于单元中心量的通信算法相对于节点量的通信算法更为简单\parencite{Liu2019}.


\subsection{单元中心型格式}


最简单且应用最广泛的有限体积格式是线性两点通量近似(Two-point flux approximation, TPFA, 见文献\parencite{Handbook2000}第3章).
TPFA方法相当稳健且易于实现, 但在变形网格和全扩散张量的情况下,
其精度有显著损失.





为了在具有任意扩散张量的一般网格上获得相容的离散,
二十世纪九十年代中后期, 挪威学者Aavatsmark等\parencite{Aavatsmark1998}和英国学者Edwards等\parencite{Edwards1998}提出二维多边形网格上的线性多点流近似方法(Multipoint Flux approximation, MPFA-O), 并广泛应用于油藏模拟中.
该方法最初用于二维问题, 并且很容易扩展到三维\parencite{Aavatsmark2006}.
该格式的优点是保持流连续, 只含单元中心未知量, 且有局部的模板, 导出稀疏线性方程组, 可以直接与流体计算耦合, 易于推广至三维情形\parencite{Aavatsmark2002}. 早期的MPFA-O格式的缺点是, 对于强各向异性扩散问题,
它容易产生非物理振荡, 导致计算精度较低.

为了改进MPFA-O方法, 研究人员随着时间的推移开发了许多MPFA格式, 例如MPFA-L\parencite{Aavatsmark2008}和EMPFA\parencite{EMPFA2008}. 邬吉明等\parencite{Wu2010}将MPFA格式的边中点改为边上的动点, 类似MPFA-O方法构造一簇线性精确的广义“多点流”格式GMPFA, 并且分析了局部线性系统的可解性. 最近, 经典的MPFA方法与高阶迎风格式相结合, 用于模拟多孔介质中的两相流动\parencite{Edwards2010,Lamine2010}.



还有一些满足单调性准则的以单元为中心的有限体积格式值得在此介绍.
文献\parencite{Potier2005}提出了一种在三角形网格上的非线性两点通量近似格式,
文献\parencite{Kapyrin2007}对该方法进行了三维扩展并进行了分析.
文献\parencite{Danilov2009}的作者提出了一种在多面体网格上的单调的单元中心有限体积格式,
当扩散张量光滑且边界条件为Dirichlet边界条件时, 该格式无需插值.
文献\parencite{Nikitin2010}将此方法推广到了对流扩散方程. 单调格式通常是非线性的, 这是为此需要付出的代价.


许多三维单元中心有限体积格式\parencite{Eymard2012,Sun2013,Gao2013,Xie2018,Wang2017}是基于调和平均点\parencite{Agelas2009}推导出的,
其调和平均点上的辅助未知量可通过面两侧的单元中心未知量进行插值.
在\parencite{Zhang2020,Gao2013}中表明, 调和平均点可能会偏离相应的单元面, 这可能会导致各向异性和间断系数扩散问题的精度损失.
此外, 由于调和平均点还取决于扩散系数, 在求解非线性扩散问题时, 其位置在每次非线性迭代中都会发生变化,
从而导致相关格式的动态模板.



本文所关注的是菱形格式, 1980年, 李德元等基于积分插值法提出了任意多边形网格上求解扩散方程的守恒格式\parencite{Li1980}. 将每条网格边上的法向流表示为相邻两个单元中心量以及这条边上的两个节点量的线性组合. 在结构四边形网格上, 该格式包含五个单元中心量和四个节点量, 而节点量用周围四个单元中心量插值消去, 最终得到的格式具有九点模板,
该格式在许多研究中得到探讨,
人们称之为九点格式\parencite{Sheng2008,Wu2010,Gao2011}.
一些欧洲学者\parencite{Coudiere1999,Manzini2007,Bertolazzi2004}也提出类似的格式,
称之为菱形格式(Diamond格式).
菱形格式的名称源于通量近似所采用的菱形模板.
最近, 一些巴西学者将相关格式称为MPFA-D格式\parencite{Contreras2016, Galindez2020a, Ricardo2021},
并将该格式应用于油藏模拟.








九点格式的优点是简洁、容易实现并且计算开销较小, 长期以来人们将九点格式应用于大型辐射(磁)流体力学程序, 例如MARED程序\parencite{Ding2008}等.
这类方法通常将面的法向通量表示为该面两侧单元中心未知量以及定义在该面顶点处的一些辅助未知量的特定的组合.
为了获得纯单元中心格式, 应通过单元中心未知量消去顶点未知量.
顶点插值算法被认为是实现菱形格式最佳精度的关键.
当网格扭曲程度比较大或者扩散系数出现间断时, 早期的九点格式精度达不到二阶, 原因在于其中的节点量采用简单的算术平均加权, 精度较低.


为了克服这个问题, 有两种解决途径.

一种途径是构造避免节点插值的单元中心型格式. 例如骆龙山等\parencite{Luo2017a}使用网格边上的两个点作为辅助未知量, 按照九点格式的推导过程推出一个法向流的表达式, 再结合MPFA-O方法推导网格边上的辅助点插值公式. 为了完全避免辅助量的插值, 美国Los Alamos国家实验室的研究人员\parencite{Lipnikov2009,Lipnikov2010}提出无需插值的单元中心型有限体积格式, 并且将该格式应用于扩散方程和对流扩散方程.

另一种途径是构造二阶精度的节点插值算法或者称为加权算法. 对于光滑扩散系数问题, 法国学者Y. Coudi\`{e}re等提出的最小二乘插值算法\parencite{Coudiere1999}, 该算法十分简洁并且具有二阶精度, 可直接推广至三维情形\parencite{Coudiere2011}, 但是对于间断系数问题, 该方法精度较低. 原因在于最小二乘插值算法没有考虑扩散系数的影响, 它总是假设解的梯度是连续的, 而事实上扩散系数间断处, 解的梯度也是间断的.

近年来, 一些二阶精度的节点插值算法被提出. 在二维情形, 一些常用的节点插值算法如下. 盛志强等\parencite{Sheng2008}基于泰勒展开构造二阶精度插值算法. 高志明等\parencite{Gao2011}提出两种线性精确的显式节点加权算法, 允许任意形式的扩散张量, 不仅非间断依赖而且非网格拓扑依赖. 其中的第二种算法称为LPEW2, 已被许多学者使用\parencite{TRHD2015,Contreras2016, Yao2012}, 有较高的国际影响力.

在三维情形, 引入调和平均点的插值算法是许多学者的替代选择, 例如文\parencite{Xie2018,Wang2017,Eymard2012}. 然而文\parencite{Zhang2020}指出, 调和平均点可能落在单元面外部, 导致精度下降, 而且对于非线性系数的扩散问题, 每次迭代都要重新计算调和平均点的位置, 导致格式的模板不固定.

在三维情形, 构造不依赖调和平均点的插值算法十分困难. 将鲁棒的LPEW2插值算法\parencite{Gao2011}扩展到三维非常困难, 直到最近, 一些巴西学者\parencite{Ricardo2021}使用与LPEW2类似的技术, 推导出四面体网格上具有显式权重的插值算法.
2017年, 来翔等\parencite{Lai2017}提出四面体网格上的非线性节点加权算法. 2023年, 骆龙山和董成\parencite{Luo2023}提出三维非结构网格中的扩展最小二乘插值算法, 克服了经典最小二乘插值算法在间断系数情形精度下降的问题, 该插值算法对于多种多样的扩散系数都能保持二阶精度.

综上, 各向异性扩散方程的单元中心型有限体积格式已有大量研究工作, 基于线性精确思想, 将二维节点插值算法推广至三维, 是研究扩散方程高精度、高置信、高效数值离散的有效手段, 这也是我们研究的出发点.


国外已有许多学者将这类格式应用于油藏模拟, 国内较少使用这类格式求解各向异性多孔介质两相流问题.

应用非结构网格上的线性精确单元中心型有限体积格式, 求解各向异性多孔介质两相流压力方程, 进而模拟两相流问题, 是一种十分有意义的集成创新.
% vim:ts=4:sw=4