一种隐式-显式求解方案以及适用于裂隙多孔介质中广义Forchheimer流动的自适应多尺度近似方法

《Computers & Mathematics with Applications》：An implicit-explicit scheme and adaptive multiscale approximation of generalized Forchheimer flow in fractured porous media

　　玛丽亚·瓦西里耶娃 | 娜娜·阿乔亚·姆布罗 | 德尼斯·斯皮里多诺夫 | 塔蒂亚娜·季莫费耶娃
德克萨斯A&M大学-科珀斯克里斯蒂分校数学与统计系，美国德克萨斯州科珀斯克里斯蒂

摘要
在本文中，我们研究了一种适用于裂隙多孔介质的广义Forchheimer流动模型。我们基于有限体积方法构建了一个半离散系统，并在其中嵌入了裂隙模型。对于非线性半离散问题的时间近似，我们考虑了两种方案：（i）线性化隐式方案（LI）和（ii）显式-隐式-空方案（EIN）。在LI方案中，我们使用前一时间级的线性化结果在每个时间级生成一个线性系统，这需要每个时间级更新问题算子。在EIN方案中，我们利用算子的加性表示来分离线性部分，并对非线性部分进行显式近似，从而为每个时间层生成线性系统。我们将这两种方案以通用（标准）形式表达，并给出了稳定性估计。我们证明了，通过仔细选择算子的加性表示，EIN方案对于非线性流动模型满足稳定性所需的性质。接下来，我们通过在粗网格上构建精确的多尺度近似来减小线性离散系统的规模，在此过程中EIN方案为多尺度基的构建提供了有效的算子选择。多尺度空间的构建基于解决局部线性谱问题，通过引入多个宏观尺度变量来准确捕捉高对比度渗透系数。我们展示了基于主要局部流动模式的自适应方法能够准确地近似裂隙多孔介质中的流动过程。为了减少与非线性流动模式相关的误差，我们基于局部残差对多尺度空间进行了自适应局部多尺度增强。研究表明，几次自适应迭代可以显著减少误差，并为广义Forchheimer流动模型生成高度精确的解。

引言
裂隙多孔介质中的流体流动是一个研究热点，其在水文地质学、石油工程和环境科学等领域具有重要意义。因此，理解流体通过这些复杂地质结构的动态至关重要。通常，描述压力梯度与流体速度之间线性关系的达西定律是研究多孔介质流动的主要框架。然而，这种线性关系在较高速度下常常失效，此时非达西效应变得显著[1]。这些偏差使得在异质系统中（如裂隙储层）准确建模流体流动变得困难。因此，研究者们转向使用非线性方程来描述多孔介质中的流动行为[2]，[3]，[4]。特别是Forchheimer方程组在捕捉线性达西近似未考虑的物理现象方面受到了广泛关注[5]，[6]，[7]，[8]，[9]，[10]。

在本文中，我们关注广义Forchheimer方程[5]，[11]。该模型通过引入速度场与压力梯度之间的非线性关系来扩展了动量方程。它结合了经验或理论定律，例如两项式、幂律和三项式Forchheimer方程，这些方程用与速度大小及其幂次成比例的项修改了达西定律。这些额外项反映了依赖于流体性质、渗透率和孔隙度的经验常数或张量，使模型能够适用于广泛的流动范围。所得关系可以通过实验数据、均质化方法或理论论证推导出来，这使得该公式非常适合捕捉多孔介质中的惯性和高速度效应[5]。然而，将这些模型应用于裂隙多孔介质会因内在复杂性而面临计算挑战。为了解决这些问题，提出了双连续体或多连续体模型[12]，[13]，[14]，[15]。这些模型通过将裂隙和多孔基质概念化为两个重叠的连续体来简化裂隙多孔介质的表示。通过区分每种连续体的独特性质及其相互作用，双连续体模型为模拟异质系统中的流体流动提供了一个稳健的框架。此外，低维建模已被证明是捕捉复杂系统关键特征的有效方法。例如，研究文章[16]，[17]，[18]展示了这些技术如何在保持关键物理行为的同时简化裂隙多孔介质的计算建模。

广义Forchheimer模型的复杂性，尤其是在裂隙多孔介质中，需要使用数值方法进行实际应用。传统的多孔介质流体流动问题数值方案，如有限差分法和有限元法，由于其稳健性和简单性而被广泛使用。例如，混合有限元法已应用于单相流体在多孔介质中的流动[19]，而块中心有限差分方案则专为非线性达西-Forchheimer方程设计[20]。在这两种情况下，都已严格证明了解的存在性和唯一性，并为相应的离散化方案提供了误差估计。这些方法将多孔介质离散为细网格，并在每个网格点求解控制方程。尽管对于小尺度问题有效，但由于需要高分辨率来捕捉细尺度异质性，这些传统方法在大规模模拟时往往计算成本高昂且不实用。为了解决这些限制，人们进行了大量研究以开发模型阶数降低技术，以大幅减少计算域的同时保留问题的关键特征。其中最常用的技术是上采样和多尺度方法。上采样技术通过采用通量平均程序来降低计算复杂性，计算有效参数以构建粗网格近似[21]，[22]，[23]，[24]，[25]，[26]。已经提出了多种方法从详细的微观描述中推导出有效的宏观性质，以确保准确和高效的大规模模拟[11]，[28]。例如，已经在[29]中研究了低雷诺数流动的多项式过滤定律的上采样。这些方法能够在保持关键物理行为的同时高效模拟复杂系统。然而，多尺度方法基于传统的有限元或有限体积方法，通过结合局部基函数来捕捉每个元素内的细尺度异质性[30]，[31]，[32]，[33]，[34]，[35]，[36]。

裂隙多孔介质的多尺度和异质特性给流体流动建模和多尺度方法构建带来了独特挑战。裂隙可以增强或抑制流体流动，显著影响渗透性和传输性质。这些挑战通过在各局部域引入多个多尺度基函数得到解决[37]，[38]，[39]。最近的研究还将多尺度方法与非局部多连续体理论结合，以提高精度并提供物理意义上的模型。例如，在[40]，[41]，[42]中，作者使用局部基函数构建了裂隙介质的非局部传递性质，从而得到了能够准确模拟裂隙与基质之间相互作用的低阶模型。进一步的进展包括采用隐式-显式时间近似方案来提高多连续体方法的效率[43]，[44]。在[45]中，提出了一种混合广义多尺度方法来解决达西-Forchheimer方程。作者使用局部细网格解决方案构建速度空间，同时在粗网格上使用简单的分段常数函数表示压力空间。在[46]中，为裂隙多孔介质中的基本Forchheimer流动引入了非局部上采样框架，该方法通过局部基函数识别多连续体参数，并在过采样区域构建非局部传递和有效性质。

数值不稳定性是求解非线性、高阶或多尺度偏微分方程（PDEs）时的常见挑战。为此，开发了各种稳定技术。特别是[47]提出了一种边界积分方法，将高阶、非线性和非局部项重新表述为线性、常数系数，从而在傅里叶空间中有效地消除了典型的稳定性约束。他们的方法在捕捉Kelvin-Helmholtz不稳定性下的复杂动态时显示出了稳健性和准确性。类似的成功也适用于不稳定的分层Boussinesq流体，突显了该方法在多种界面问题中的有效性。[48]提出了一种更通用的方法，通过在PDE的右侧添加和减去精心选择的项来引入稳定性：一项作为阻尼部分并隐式处理，另一项则显式处理。这种方法已成功应用于稳定平均曲率流动、Kuramoto-Sivashinsky和Rayleigh-Taylor不稳定性等方程，证明了其广泛应用性。在此基础上，EIN方法及其变体（如Variable-Coefficient EIN (VC-EIN) 方法被提出用于高阶扩散和扩散方程[49]。这些方法通过隐式-显式离散化添加和减去线性最高阶导数项来提高稳定性。VC-EIN方法通过在计算域内变化系数进一步改进了这种方法，提高了高变异性区域的准确性。结合局部离散Galerkin (LDG) 空间离散化，这些方法实现了高阶精度和稳定性。[50]的最新工作提出了一种高效方案，用于使用EIN方法解决多尺度、非线性、高对比度扩散问题。非线性项被分解为线性和阻尼部分，分别由隐式和显式时间推进方案处理。为了捕捉线性部分的多尺度特性，他们采用了部分显式分裂方案，并构建了合适的多尺度子空间以加快计算。在这种设置中，高对比度子空间被隐式离散化，而其余子空间则显式处理，确保了稳定性和收敛性。EIN方法在管理非线性方面以及使用部分显式分裂方案处理线性多尺度分量方面的新颖应用是关键贡献。EIN方案将非线性PDE算子分为线性和阻尼部分，并基于一般的隐式-显式 (IMEX) 方案和算子加性表示[51]。IMEX为在具有挑战性的PDE中平衡稳定性和计算成本提供了一种方法[52]，[53]。通常，IMEX方法将PDE算子分为扩散部分（隐式）和对流或反应部分（显式），允许更大的时间步长并降低刚性[54]。同样，加性算子差分方案将PDE分解为多个子算子，通常按空间变量或物理过程进行，并分别解决每个子问题[51]。尽管有这些进展，将这些稳定技术应用于裂隙多孔介质中的广义Forchheimer流动模型仍然有限。裂隙多孔系统具有独特的挑战，包括强烈的非线性、多尺度相互作用和耦合的裂隙-基质动力学。

本文提出了一种用于裂隙多孔介质中广义非线性流动模型的高效计算算法。该非线性模型采用混合维度公式表示，并使用带有嵌入裂隙模型的传统有限体积方法进行近似。所提出的方法基于算子的加性表示，仔细分离线性和非线性部分，确保了隐式-显式方案的稳定性。我们比较了传统LI方案和提出的EIN方案的迭代和非迭代方法。这两种方案都在通用框架中进行了展示，并提供了稳定性估计，证明了它们在解决裂隙多孔介质中的广义Forchheimer流动模型时的稳健性和适用性。接下来，我们通过在粗网格上构建精确的多尺度近似来减小线性离散系统的规模。EIN方案为多尺度基的构建提供了有效的算子选择。多尺度空间的构建基于解决局部线性谱问题，通过引入多个宏观尺度变量来准确捕捉高对比度渗透系数。为了减少与非线性流动模式相关的误差，我们基于局部残差对多尺度空间进行了自适应局部多尺度增强。本文为离线和在线多尺度方法提出了误差分析。数值结果证明了多尺度隐式-显式方案在二维和三维公式中三个测试案例中的适用性。

本文的结构如下：第2节介绍了裂隙多孔介质中广义Forchheimer流动的数学模型。第3节介绍了基于有限体积法和两种时间近似方案的空间半离散化方法。接下来，我们以通用（典型）形式描述所提出的方案，对两种方法的稳定性进行了估计，并证明了所考虑的非线性流模型在适当选择算子表示方式的情况下满足稳定性所需的属性。第4节通过引入自适应离线多尺度空间来阐述粗网格近似方法。第5节提出了针对非线性流问题的自适应局部在线多尺度空间增强方法。我们对离线和在线多尺度粗网格方法都建立了误差分析。第6节展示了在二维和三维 fractured porous media 模型中三个测试案例的数值结果。最后，第7节给出了结论。

**数学模型**
对于可压缩流体，在域 Ω 内有质量守恒方程：
$$
\partial(\varphi\rho)\partial t + \nabla \cdot (\rho u) = f, \quad x \in \Omega,
$$
其中 $\varphi$ 是孔隙率，$\rho$ 是流体密度，$u$ 是流体速度，$p$ 是压力，$f$ 是给定的源项。例如，当 $\varphi = \varphi(x, p)$ 且 $\rho$ 为常数（适用于可压缩岩石和不可压缩流体）时，可得到：
$$
\partial\varphi\rho\partial t = \rho\nabla\varphi\partial t = \rho\nabla\varphi\frac{\partial p}{\partial \rho}\frac{\partial p}{\partial t} = \rho c \frac{\partial p}{\partial t},
$$
其中 $c = \varphi_0 c^r$ 且 $c_r = \frac{1}{\varphi_0\nabla\varphi\frac{\partial p}{\partial \rho}$ 是岩石的压缩性。

现在我们有：
$$
c \frac{\partial p}{\partial t} + \nabla \cdot u = f, \quad x \in \Omega.
$$

**离散问题**
本节使用有限体积方法（FVM）介绍了多连续介质流动问题的细网格近似。对于非线性流动问题的时间近似，我们构建了 EIN [47]、[48]、[50]、[56] 方法。

**多尺度空间上的粗网格近似**
在本节中，我们描述了在粗网格上构建多尺度近似的方法。我们构造了一个粗网格 $TH$，定义为：
$$
TH = \bigcup_i K_i,
$$
其中 $K_i$ 是粗网格中的单元格。接下来，通过合并粗网格顶点周围的单元格来创建一个局部粗邻域域 $\omega_i$：
$$
\omega_i = \{ \bigcup_j(K_j : x_i \in K_j \in TH) \},
$$
其中 $i = 1, \ldots, N_c$，$N_c$ 是粗网格顶点的数量。

在每个 $\omega_i$ 内，我们解决一个局部耦合的广义特征值问题，并使用 dominant 特征函数来构建局部解。

**多尺度空间的在线增强**
局部谱基函数是为算子的线性部分构建的，可以离线预先计算。为了处理问题的非线性和动态特性，我们进行了多尺度空间的在线增强 [39]、[63]、[64]。

根据 (26)，我们在粗网格上每个时间步解决以下线性方程组：
$$
BH(pmsn)pH_{n+1} = bH_n(pmsn),
$$
其中：
$$
\begin{align*}
&\text{LI-MS: } & BH(pmsn) = MH + \tau AH(pmsn) \quad \text{和} \quad bH_n(pmsn) = bH_n = \tau FH + MHpH_n, \\
&\text{EIN-MS: } & BH(pmsn) = MH + \tau AH(1) \quad \text{和} \quad bH_n(pmsn) = \tau F.
\end{align*}
$$

**数值结果**
我们考虑了二维和三维测试问题，分别在域 $\Omega = [0,1] \times [0,1]$ 和 $\Omega = [0,1] \times [0,1] \times [0,0.1]$ 中进行了研究，其中包含 14 条裂缝（图 1）。如前几节所述，对于细网格近似，我们在域 $\Omega$ 内的结构化网格上使用有限体积法对空间变量进行近似。具体来说，对于二维问题，我们使用了一个 400 × 400 的网格，其中每个单元格为正方形；裂缝网格包含 1226 个单元格和 1248 个面。细尺度问题的自由度为 $DOF_h$。

**结论**
我们考虑了在 fractured porous media 中的广义 Forchheimer 流动问题，其数学模型以混合维度形式给出。为了数值求解该问题，我们使用有限体积法（FVM）在细网格上对空间进行了半离散化，并嵌入了裂缝模型（EFM）。为了减小半离散系统的规模，我们在粗网格上构建了多尺度近似。我们引入了自适应离线多尺度空间和自适应在线增强机制。