饱和多孔介质渗流的强-弱不连续统一非局部模型

摘要

非均匀饱和多孔介质流体传输模拟关键在于需要考虑由孔隙、裂隙以及不同渗透性介质组成的强–弱不连续面，但这与经典连续介质力学局部流体模型所定义的偏微分方程是不相容的，从而导致饱和多孔介质渗流模拟的困难。基于统一变分近场动力学的基本思想，提出了一种向量态函数用于描述饱和多孔介质流体的非局部传输作用，并在此基础上建立饱和多孔介质流体输运的非局部空间积分–时间微分型控制方程。该非局部流体输运模型统一描述了流体在强–弱不连续面的力学行为，避免了传统局部力学模型在不连续界面处的导数无定义性，且从理论上证明，当非局部流体模型中的非局部作用半径趋于零时，非局部模型可退化为局部模型。为了消除非局部模型内在的零能模式问题，基于罚函数方法提出一种全隐式的数值求解格式以保证数值计算的精确性。最后通过算例分析，模拟了三维非均匀不连续多孔介质中流体的传输过程，揭示流体穿越物质界面、孔隙和裂隙的力学行为。为研究饱和多孔介质中流体在不连续界面传输的力学行为提供了一种分析模型。

关键词

非均匀饱和多孔介质; 非局部模型; 界面传输; 裂隙; 孔隙

岩土体是含有孔隙和裂隙的非均匀多孔介质，这类多孔介质中流体传输过程需要穿越由介质中孔隙、裂隙和物质界面组成的强–弱不连续面，这与经典连续介质力学模型中描述局部流体模型的偏微分方程是不相容的，从而导致需要额外的数值技术来描述流体在这些不连续界面的力学行为^［1-3］。

目前兴起的近场动力学（Peridynamics）理论为该问题分析提供了一种新的思路。近场动力学是由Silling等^［4-5］提出的一种新的非局部理论，其控制方程是积分方程，考虑连续体内物质点之间非局部相互作用，避免了经典连续介质模型中偏微分方程带来的空间求导奇异性^［6-7］。目前近场动力学已经被广泛应用于包括脆性材料断裂^［8-9］、黏弹性断裂^［10-11］和多物理场耦合断裂^［12］等各类非连续问题。将原本描述断裂等不连续问题的近场动力学理论应用于饱和多孔介质流体传输问题，其关键在于如何建立多孔介质中流体传输的非局部形式守恒律方程和本构律方程。

近期已有一些学者通过修改近场动力学方程以建立孔隙介质中流体传输的控制方程，如Ni等^［13］和Bobaru等^［14］首先将键型近场动力学模型中能量对应的思想扩展到热传导和流体传输问题，建立了热–流耦合的非局部模型；同样地，Oterkus等^［15］也建立了一种热–流耦合的非局部键型近场动力学模型；随后Katiyar等^［16-17］基于常规态型近场动力学模型推导出多孔介质中流体传输的非局部形式控制方程，并研究了流体穿越介质中不连续界面的力学行为。然而，上述研究主要针对键型和常规态型近场动力学模型，由于物质点之间键的相互作用是独立的，从而导致其推导的非局部扩散模型无法描述流体传输的非线性力学行为。此外，目前所发展的非局部传输模型主要基于键型模型中能量对应的思想，与经典连续介质力学模型之间的对应关系不够清晰。

基于统一变分近场动力学^［18］的基本思想，本文提出一种表征饱和多孔介质中流体非局部传输作用的流量向量态函数，统一描述流体在强–弱不连续面传输过程中的力学行为，避免局部力学模型在不连续界面处出现的导数无定义性，并建立相应非局部形式的质量守恒方程、线动量守恒方程和本构方程。从理论证明了当非局部流体模型中的非局部作用半径趋于零时，该非局部模型将退化为经典连续介质力学中的Darcy流体模型。另外，提出一种基于罚函数方法的全隐式数值求解格式，旨在消除非局部模型内在的零能模式难题，保证数值结果的精确性。

1 渗流传输的非局部控制方程

1.1　非局部质量守恒方程

在传统局部模型中，水压力梯度是驱动多孔介质中流体传输的主要因素，故质量守恒方程以偏微分方程形式呈现。基于非常规态型近场动力学中的力态公式^［6-7］，本文提出了流量态向量的概念，并给出了非均匀岩土介质中流体传输的非局部形式质量守恒方程。

图1为非均匀饱和多孔介质中渗流输运的非局部模型示意。

图1  非均匀饱和多孔介质中渗流输运的非局部模型示意

Fig. 1  Schematic of the nonlocal fluid transport model in the heterogeneous saturated porous media

假设多孔介质中流体处于恒温条件下的弱可压缩状态，流体密度和压力之间满足关系如式（1）：

$$\rho (\mathit{x},t)={\rho }^{\mathrm{o}}(\mathit{x})[\mathrm{1}+c_{\mathrm{f}}(p(\mathit{x},t)-p^{\mathrm{o}}(\mathit{x}))]$$

（1）

式中：$\rho (\mathit{x})$为流体密度，$\mathrm{k}\mathrm{g}\cdot {\mathrm{m}}^{-\mathrm{3}}$；$\mathit{x}$为流体物质点的位置矢量，$\mathrm{m}$；t为时间变量，$\mathrm{s}$；${\rho }^{\mathrm{o}}(\mathit{x})$为参考流体密度，$\mathrm{k}\mathrm{g}\cdot {\mathrm{m}}^{-\mathrm{3}}$；$c_{\mathrm{f}}$为流体压缩系数，$\mathrm{P}{\mathrm{a}}^{-\mathrm{1}}$；$p(\mathit{x},t)$为流体压力，$\mathrm{P}\mathrm{a}$；$p^{\mathrm{o}}(\mathit{x})$为参考流体压力，$\mathrm{P}\mathrm{a}$。流体传输的非局部质量守恒方程可表示为

（2）

式中：和分别表示在流体物质点${\mathit{x}}^{\text{'}}$和x处的流量态矢量，$\mathrm{k}\mathrm{g}\cdot {\mathrm{m}}^{-\mathrm{6}}\cdot {\mathrm{s}}^{-\mathrm{1}}$；$H_{\mathit{x}}$表示流体物质点x的非局部作用区域；$q(\mathit{x},t)$为所考虑多孔介质区域内流体源的流量函数，kg·m^-3·s^-1.需要说明的是，流量态矢量中和表征了流体${\mathit{x}}^{\text{'}}$和x之间的传输作用。

流体物质点之间非局部传输作用的流量态矢量定义如下：

（3）

式中：$\mathit{u}(\mathit{x},t)$为流体流速，$\mathrm{m}\cdot {\mathrm{s}}^{-\mathrm{1}}$；$\Vert \mathit{x}-{\mathit{x}}^{\text{'}}\Vert$为${\mathit{x}}^{\text{'}}$和$\mathit{x}$之间的欧氏距离；$\omega (‖{\mathit{x}}^{\text{'}}-\mathit{x}‖)$为核函数，其定义如式（4）：

$$\omega (‖{\mathit{x}}^{\text{'}}-\mathit{x}\Vert )=\left\{\begin{array}{cc}\frac{\mathrm{1}}{\mathrm{1}-{\mathrm{e}}^{-\mathrm{9}}}({\mathrm{e}}^{-\frac{\Vert {\mathit{x}}^{\text{'}}-\mathit{x}{\Vert }^{\mathrm{2}}}{{\delta }^{\mathrm{2}}}}-{\mathrm{e}}^{-\mathrm{9}}),& {\mathit{x}}^{\text{'}}\in H_{\mathit{x}}\\ \mathrm{0},& {\mathit{x}}^{\text{'}}\notin H_{\mathit{x}}\end{array}\right.$$

（4）

$\mathit{K}(\mathit{x})$为x的形状张量，定义为

$$\mathit{K}(\mathit{x})=\underset{{\mathit{H}}_{\mathit{x}}}{\int }\mathit{\omega }(‖{\mathit{x}}^{\text{'}}-\mathit{x}\Vert )({\mathit{x}}^{\text{'}}-\mathit{x})\otimes ({\mathit{x}}^{\text{'}}-\mathit{x})\mathrm{d}{V}_{{\mathit{x}}^{\text{'}}}$$

（5）

式中：${\mathit{x}}^{\text{'}}$的积分区域为x的非局部作用区域$H_{\mathit{x}}$。

1.2　非局部线动量守恒方程与本构方程

为了得到完整描述非均匀多孔介质中流体传输的非局部模型，还需定义流体传输的非局部线动量守恒方程和本构方程，其非局部动量守恒方程定义为

$$\underline{\mathit{u}}[\mathit{x},t]=-\frac{\mathrm{1}}{\mu }\underline{\mathit{P}}[\mathit{x}]{\nabla }_x^N\underline{\varphi }[\mathit{x},t]$$

（6）

式中：μ为流体黏度；$\underline{\mathit{P}}[\mathit{x}]$为$\mathit{x}$处的非局部有效渗透率张量，由式（7）非局部本构方程确定。

$$\underline{\mathit{P}}[\mathit{x}]=\mathrm{2}\frac{k_{\mathrm{r}}(\mathit{x})k_{\mathrm{r}}(\mathit{x}\mathrm{\text{'}})}{k_{\mathrm{r}}(\mathit{x})+k_{\mathrm{r}}(\mathit{x}\mathrm{\text{'}})}{\mathit{P}}_{\mathrm{a}\mathrm{b}}$$

（7）

式中：${\mathit{P}}_{\mathrm{a}\mathrm{b}}$为饱和岩土介质的绝对渗透率张量；$k_{\mathrm{r}}(\mathit{x})$和$k_{\mathrm{r}}(\mathit{x}\mathrm{\text{'}})$分别为${\mathit{x}}^{\text{'}}$和$\mathit{x}$的相对渗透率系数。${\nabla }_{\mathit{x}}^N\underline{\varphi }[\mathit{x},t]$为流体势函数的非局部梯度，如式（8）：

（8）

式中：和为${\mathit{x}}^{\text{'}}$和x处的非局部流体势标量态函数。

（9）

式中：$\varphi (\mathit{x},t)=p(\mathit{x},t)-\rho (\mathit{x})gz(\mathit{x})$为物质点位置矢量的局部势流函数，本文忽略由于重力产生的流体势函数$\rho (\mathit{x})gz(\mathit{x})$。

式（1）、式（5）和式（6）构成了非均匀饱和多孔介质中流体传输的全部控制方程，为了求解该方程，还需引入非局部压力边界${\Gamma }_{\mathrm{g}}^{\mathrm{p}}$和流量边界${\Gamma }_{\mathrm{h}}^{\mathrm{p}}$，如图1所示，采用文献［16］所提出的边界施加方法。

1.3　与经典Darcy流模型的一致性证明

本文的非局部渗流模型引入非局部作用半径来表征流体物质点之间的传输作用。接着理论证明，当物质点之间的作用半径趋于零时，非局部模型退化为局部的Darcy流体模型。

基于Taylor公式，$H_{\mathit{x}}$中物质点处${\mathit{x}}^{\text{'}}$的局部流量$\rho ({\mathit{x}}^{\text{'}})\mathit{u}({\mathit{x}}^{\text{'}},t)$可表示为

$$\rho ({\mathit{x}}^{\text{'}},t)\mathit{u}({\mathit{x}}^{\text{'}},t)=\rho (\mathit{x},t)\mathit{u}(\mathit{x},t)+{\nabla }_x(\rho (\mathit{x})\mathit{u}(\mathit{x},t))\cdot ({\mathit{x}}^{\text{'}}-\mathit{x})+r(‖{\mathit{x}}^{\text{'}}-\mathit{x}{\Vert }^{\mathrm{2}})$$

（10）

式中：${\nabla }_x(\rho (\mathit{x},t)\mathit{u}(\mathit{x},t))$表示$\mathit{x}$处的流量梯度；$r(‖{\mathit{x}}^{\text{'}}-\mathit{x}{\Vert }^{\mathrm{2}})$表示Taylor公式中的高阶余项。

将 式 （9） 代 入 式 （1） ， 得 到

（11）

式中：$\mathrm{\Delta }(\rho \mathit{u})({\mathit{x}}^{\text{'}},\mathit{x},t)=\rho ({\mathit{x}}^{\text{'}},t)\mathit{u}({\mathit{x}}^{\text{'}},t)-\rho (\mathit{x},t)\mathit{u}(\mathit{x},t)$为${\mathit{x}}^{\text{'}}$和$\mathit{x}$处流量的差值；${\nabla }_{\mathit{x}}(\rho (\mathit{x})\mathit{u}(\mathit{x},t))$为流量函数的局部梯度。从式（11）可知，当物质点的非局部作用区域$H_{\mathit{x}}$的半径${\delta }_x$趋于零时，非局部质量守恒方程将退化为局部质量守恒方程。由此定义流量的非局部梯度算子为

（12）

同理，对$H_{\mathit{x}}$中${\mathit{x}}^{\text{'}}$的势流函数$\varphi ({\mathit{x}}^{\text{'}},t)$进行Taylor展开，可以得到

$$\begin{array}{l}{}_{}\varphi ({\mathit{x}}^{\text{'}},t)=\varphi (\mathit{x},t)+{\nabla }_x\varphi (\mathit{x},t)\cdot ({\mathit{x}}^{\text{'}}-\\ \mathit{x})+r(‖{\mathit{x}}^{\text{'}}-\mathit{x}{\Vert }^{\mathrm{2}})\end{array}$$

（13）

将式（12）代入式（7），可得渗流的局部线动量守恒方程和非局部质量守恒方程如式（14）：

（14）

式中：$\mathrm{\Delta }\varphi ({\mathit{x}}^{\text{'}},\mathit{x},t)=\varphi ({\mathit{x}}^{\text{'}},t)-\varphi (\mathit{x},t)$为${\mathit{x}}^{\text{'}}$和$\mathit{x}$处流体势函数的差值；${\nabla }_x\varphi (\mathit{x},t)$为流体势函数在$\mathit{x}$处的具体梯度。注意到$\mathit{x}$处的非局部渗透率张量和局部渗透率张量有如式（15）的关系：

$$\begin{array}{c}\underline{\mathit{P}}[\mathit{x}]=\mathrm{2}\frac{k_{\mathrm{r}}(\mathit{x})k_{\mathrm{r}}(\mathit{x}\mathit{\text{'}})}{k_{\mathrm{r}}(\mathit{x})+k_{\mathrm{r}}(\mathit{x}\mathrm{\text{'}})}{\mathit{P}}_{\mathrm{a}\mathrm{b}}\to \\ \mathit{P}(\mathit{x})=k_{\mathrm{r}}(\mathit{x}){\mathit{P}}_{\mathrm{a}\mathrm{b}},{\delta }_x\to \mathrm{0}\end{array}$$

（15）

结合式（13）和式（14）可知，当物质点的非局部作用半径${\delta }_x$趋于零时，非局部的线动量守恒方程、本构方程与局部线动量守恒方程、本构方程可通过如式（16）等式建立对应关系。

$$\begin{array}{l}\mathit{u}(\mathit{x},t)=-\frac{\mathrm{1}}{\mu }\underline{\mathit{P}}[\mathit{x}]{\nabla }_{\mathit{x}}^N\underline{\varphi }[\mathit{x},t]=\\ -\frac{\mathrm{1}}{\mu }\underline{\mathit{P}}[\mathit{x}]{\nabla }_{\mathit{x}}\varphi (\mathit{x},t)+r(\Vert {\mathit{x}}^{\text{'}}-\mathit{x}\Vert )\to \\ -\frac{\mathrm{1}}{\mu }\mathit{P}(\mathit{x}){\nabla }_{\mathit{x}}\varphi (\mathit{x},t)\end{array}$$

（16）

2 非局部变分框架与全隐式数值求解格式

2.1　非局部变分公式与物质点离散

不连续伽辽金方法^［19-20］和物质点离散方法是高效求解非局部模型的有效数值离散方法。将待求解多孔介质区域离散为$N_p$个物质点$\{{\mathit{x}}_I{\}}_{I=\mathrm{1}}^{N_p}$，每个物质点的体积为$\mathrm{\Delta }{V}_I$，则待求解的多孔介质区域$\Omega$可用物质点表征为

$$\Omega ={\displaystyle \underset{I=\mathrm{1}}{\overset{N_p}{\cup }}}\Delta V_I$$

（17）

${\mathit{x}}_I$的非局部作用区域$H_{{\mathit{x}}_I}$内，其物质点可用局部编号表示为

$$N_I=\left\{\mathrm{1},...,k,...,n_I\right\}$$

（18）

与有限元方法求解偏微分方程的离散弱形式相同，定义离散试验函数空间和离散测试函数空间如式（19）：

$$\begin{array}{l}{\mathit{S}}_t^{\mathrm{d}\mathrm{i}\mathrm{s}}=\{p(\mathit{x},t)\in L_{\mathrm{2}}(\Omega \times [\mathrm{0},T])|\\ p(\mathit{x},t)=g^{\mathrm{p}}(\mathit{x},t)\mathrm{i}\mathrm{n}{\Gamma }_{\mathrm{g}}^{\mathrm{p}}\}\\ {\mathit{W}}_t^{\mathrm{d}\mathrm{i}\mathrm{s}}=\{\delta p(\mathit{x},t)\in L_{\mathrm{2}}(\Omega \times [\mathrm{0},T])|\\ \delta p(\mathit{x},t)=\mathrm{0}\mathrm{i}\mathrm{n}{\Gamma }_{\mathrm{g}}^{\mathrm{p}}\}\end{array}$$

（19）

则非局部传输模型的离散弱形式可表示为

（20）

式中：$\mathrm{\delta }p=\mathrm{\delta }p(\mathit{x},t)$为压力的变分函数；$g^{\mathrm{p}}({\mathit{x}}_{\mathrm{I}},\mathrm{t})$为非局部压力边界${\Gamma }_g^p$上的压力函数；$h^{\mathrm{p}}=h^{\mathrm{p}}(\mathit{x},t)$为非局部流量边界${\Gamma }_h^p$的流体流量函数。，${}_{\Omega ,\mathrm{d}\mathrm{i}\mathrm{s}}$和${}_{{\Gamma }_{\mathrm{h}}^{\mathrm{p}},\mathrm{d}\mathrm{i}\mathrm{s}}$定义如式（21）：

（21）

式中：下标dis表示欧氏空间的离散内积；$\mathrm{\Delta }{V}_{{\mathit{x}}_{\mathit{I}}}$为${\mathit{x}}_I$的体积；${\tilde{N}}_p$为非局部流量边界${\Gamma }_{\mathrm{h}}^{\mathrm{p}}$内的物质点数；${\nabla }_{x_I}^N\mathrm{\delta }\underline{p}={\nabla }_{x_I}^N\mathrm{\delta }\underline{p}[\mathit{x},t]$为压力变分函数在${\mathit{x}}_I$处的非局部梯度，定义如式（22）：

$${\nabla }_{x_I}^N\underline{p}={\displaystyle \sum _{J=\mathrm{1}}^{n_I}}\omega (‖{\mathit{x}}_I-{\mathit{x}}_J‖)\omega (‖{\mathit{x}}^{\text{'}}-\mathit{x}\Vert )(p({\mathit{x}}_J,t)-p({\mathit{x}}_I,t)){\mathit{K}}^{-\mathrm{1}}({\mathit{x}}_J)\cdot ({\mathit{x}}_J-{\mathit{x}}_I)\mathrm{\Delta }{V}_{{\mathit{x}}_J}$$

（22）

式中：$\mathrm{\Delta }{V}_{x_J}$为物质点${\mathit{x}}_J$的体积。

2.2　零能模态控制

由于非局部梯度（11）内在的零能模式将导致数值解答产生严重的数值震荡问题^［21］，通过引入罚函数控制零能模式的出现，定义${\mathit{x}}_I$处的罚函数如式（23）：

$$W_p^{\mathrm{d}\mathrm{i}\mathrm{s}}({\mathit{x}}_I,t)=\frac{\beta ({\mathit{x}}_I)\mathrm{\Delta }{V}_{{\mathit{x}}_I}}{\mathrm{2}}{\displaystyle \sum _{J=\mathrm{1}}^{n_I}}{\omega }_{{\mathit{x}}_I{\mathit{x}}_J}{z}_p^{\mathrm{2}}({\mathit{x}}_J,{\mathit{x}}_I,t)\mathrm{\Delta }{V}_{{\mathit{x}}_J}$$

（23）

式中：${\omega }_{{\mathit{x}}_I{\mathit{x}}_J}=\omega (‖{\mathit{x}}_J-{\mathit{x}}_I‖)$为${\mathit{x}}_I$和${\mathit{x}}_J$之间的影响函数；$\beta ({\mathit{x}}_I)$为正则化罚函数因子，定义为

$$\beta ({\mathit{x}}_I)=\frac{\alpha }{\sqrt[]{{\displaystyle \sum _{i=\mathrm{1}}^{n_{\mathrm{d}}}}{\lambda }_i^{\mathrm{2}}({\mathit{x}}_I)}}$$

（24）

式中：${\lambda }_{\mathrm{i}}({\mathit{x}}_I)(i=\mathrm{1},...,n_{\mathrm{d}})$为物质点${\mathit{x}}_I$的形状张量$\mathit{K}({\mathit{x}}_I)$的特征值；$\alpha$为罚函数因子，取$\mathrm{1}\times {\mathrm{10}}^{-\mathrm{3}}\sim \mathrm{1}\times {\mathrm{10}}^{\mathrm{3}}$。$z_p({\mathit{x}}_J,{\mathit{x}}_I,t)$为${\mathit{x}}_I$和${\mathit{x}}_J$之间非一致部分压力变化，定义如式（25）：

$$z_p({\mathit{x}}_I,{\mathit{x}}_J,t)={\nabla }_{\mathit{x}}^N\underline{p}[{\mathit{x}}_I,t]\cdot \mathrm{\Delta }{\mathit{x}}_{IJ}-\mathrm{\Delta }p({\mathit{x}}_I,{\mathit{x}}_J,t)$$

（25）

式中：$\mathrm{\Delta }{\mathit{x}}_{IJ}={\mathit{x}}_J-{\mathit{x}}_I$；$\mathrm{\Delta }p({\mathit{x}}_I,{\mathit{x}}_J,t)=p({\mathit{x}}_J,t)-p({\mathit{x}}_I,t)$。

对离散罚函数（22）分别求一阶和二阶变分导数，可以得到罚函数荷载向量和罚函数刚度矩阵

$$\begin{array}{l}{\mathit{R}}_p({\mathit{x}}_I,t)=\beta ({\mathit{x}}_I)\mathrm{\Delta }{V}_{{\mathit{x}}_I}\\ {\left[-{\displaystyle \sum _{J=\mathrm{1}}^{n_I}}{r}_{I\mathrm{1}}^p({\mathit{x}}_I,t),\cdots ,r_{Ik}^p({\mathit{x}}_I,t),\cdots ,r_{I{n}_I}^p({\mathit{x}}_I,t)\right]}^{\mathrm{T}}\\ {\mathit{K}}_p({\mathit{x}}_I,t)=\beta ({\mathit{x}}_I)\mathrm{\Delta }{V}_{{\mathit{x}}_I}{\mathbf{S}}_I-\\ \beta ({\mathit{x}}_I)\mathrm{\Delta }{V}_{{\mathit{x}}_I}{\mathit{G}}_I^{\mathrm{T}}\mathit{K}({\mathit{x}}_I){\mathit{G}}_I\end{array}$$

（26）

式中：$r_{IJ}^p({\mathit{x}}_I,t)={\omega }_{{\mathit{x}}_I{\mathit{x}}_J}(\mathrm{\Delta }p({\mathit{x}}_I,{\mathit{x}}_J,t)-{\nabla }_{x_I}^N\underline{p}\mathrm{\Delta }{\mathit{x}}_{IJ})\mathrm{\Delta }{V}_{{\mathit{x}}_J}$；${\mathbf{S}}_I$和${\mathit{G}}_I$的定义为

$$\begin{array}{l}{\mathit{S}}_I=\left[\begin{array}{cccccc}{\displaystyle \sum _{J=\mathrm{1}}^{n_I}}{S}_J& -S_{\mathrm{1}}& \dots & -S_k& \dots & -S_{n_I}\\ -S_{\mathrm{1}}& S_{\mathrm{1}}& \dots & \mathrm{0}& \dots & \mathrm{0}\\ ⋮& ⋮& & ⋮& & ⋮\\ -S_k& \mathrm{0}& \dots & S_k& \dots & \mathrm{0}\\ ⋮& ⋮& & & & \\ -S_{n_I}& \mathrm{0}& \dots & \mathrm{0}& \dots & S_{n_I}\end{array}\right]\\ {\mathit{G}}_I=\left[\begin{array}{cccccc}-{\displaystyle \sum _{k=\mathrm{1}}^{n_I}}{b}_k^x& b_{\mathrm{1}}^x& \cdots & b_k^x& \cdots & b_{n_I}^x\\ -{\displaystyle \sum _{k=\mathrm{1}}^{n_I}}{b}_k^y& b_{\mathrm{1}}^y& \cdots & b_k^y& \cdots & b_{n_I}^y\\ -{\displaystyle \sum _{k=\mathrm{1}}^{n_I}}{b}_k^z& b_{\mathrm{1}}^z& \cdots & b_k^z& \cdots & b_{n_I}^z\end{array}\right]\end{array}$$

（27）

式中：$S_k={\omega }_{{\mathit{x}}_I{\mathit{x}}_k}\mathrm{\Delta }{V}_{{\mathit{x}}_k}$；${\mathit{b}}_k={\left[b_k^x,b_k^y,b_k^z\right]}^{\mathrm{T}}=\omega (‖{\mathit{x}}_k-{\mathit{x}}_I‖)({\mathit{x}}_k-{\mathit{x}}_I{)}^{\mathrm{T}}{\mathit{K}}^{-\mathrm{1}}({\mathit{x}}_I)\mathrm{\Delta }{V}_{{\mathbf{x}}_k}$$(k=\mathrm{1},...,n_I)$。

因此，总体罚函数荷载向量和刚度矩阵为

$${\mathit{F}}_p^g(t)=\stackrel{N_p}{\underset{I=\mathrm{1}}{\mathit{A}}}{\mathit{R}}_p({\mathit{x}}_I,t)$$

$${\mathit{K}}_p^g(t)=\stackrel{N_p}{\underset{I=\mathrm{1}}{\mathit{A}}}{\mathit{K}}_p({\mathit{x}}_I,t)$$

（28）

2.3　隐式数值求解格式

为能精确模拟流体在介质中界面、裂隙和孔隙处传输行为，在推导罚函数荷载向量和刚度矩阵的基础上，采用Newton-Raphson隐式方法进行求解^［22-23］。结合式（19）和式（27），可得最终求解的非线性方程如式（29）：

$$\mathit{B}\frac{{\mathit{p}}^{n+\mathrm{1}}-{\mathit{p}}^n}{\mathrm{\Delta }{t}^n}+\mathit{C}({\mathit{p}}^{n+\mathrm{1}}){\mathit{p}}^{n+\mathrm{1}}=\mathit{f}({\mathit{p}}^{n+\mathrm{1}})$$

（29）

式中：$\mathrm{\Delta }{t}^n$为时间步长；${\mathit{p}}^n=[p_{\mathrm{1}}(t_n),p_{\mathrm{2}}(t_n),...,p_{N_p}(t_n{)]}^{\mathrm{T}}$。矩阵$\mathit{B}$和$\mathit{C}(\mathit{p})$、向量$\mathit{f}(\mathit{p})$的定义为

$$\begin{array}{l}\mathit{B}(\mathit{p})=\stackrel{N_p}{\underset{I=\mathrm{1}}{\mathit{A}}}{\mathit{B}}_I\\ \mathit{C}(\mathit{p})=\stackrel{N_p}{\underset{I=\mathrm{1}}{\mathit{A}}}{\mathit{G}}_{p,I}^T{\mathit{D}}_I{\mathit{G}}_{\varphi ,I}\\ \mathit{f}(\mathit{p})={\mathit{f}}_{\mathrm{e}\mathrm{q}}(\mathit{p})-{\mathit{R}}_p^g(t)\end{array}$$

（30）

其中${\mathit{B}}_{\mathrm{I}}$、${\mathit{G}}_{p,I}$和${\mathit{D}}_I$定义为

$$\begin{array}{l}{\mathit{B}}_I=\left[\begin{array}{ccc}{B}_I^{\mathrm{1}}& & \\ & \ddots & \\ & & B_I^{n_I}\end{array}\right]\\ {\mathit{G}}_{p,I}=\left[-{\displaystyle \sum _{J=\mathrm{1}}^{n_I}}{\mathit{g}}_p^J,{\mathit{g}}_p^{\mathrm{1}},...,{\mathit{g}}_p^{n_I}\right]\\ {\mathit{D}}_I=\left[\begin{array}{ccc}{D}_I^{\mathrm{1}}& & \\ & \ddots & \\ & & D_I^{n_I}\end{array}\right]\end{array}$$

（31）

式中：$B_I^J=c_{\mathrm{f}}\rho ({\mathit{x}}_J)\mathrm{\Delta }{V}_{x_J}$；${\mathit{g}}_p^J={\mathit{K}}^{-\mathrm{1}}({\mathit{x}}_J)({\mathit{x}}_J-{\mathit{x}}_I)\mathrm{\Delta }{V}_{x_J}$；$D_I=-\frac{\mathrm{1}}{\mu }\underline{\mathit{P}}[{\mathit{x}}_J]\mathrm{\Delta }{V}_{{\mathit{x}}_J}$，$J=\mathrm{1,2},...,n_I$。${\mathit{f}}_{\mathrm{e}\mathrm{q}}(\mathit{p})$可定义为

$${\mathit{f}}_{\mathrm{e}\mathrm{q}}(\mathit{p})=$$$$[\mathrm{\Delta }{V}_{x_{\mathrm{1}}}h({\mathit{x}}_{\mathrm{1}},t)$$

$$\mathrm{\Delta }{V}_{x_{\mathrm{2}}}h({\mathit{x}}_{\mathrm{2}},t),...,\mathrm{\Delta }{V}_{x_{N_p}}h({\mathit{x}}_{N_p}{,t)]}^{\mathrm{T}}$$

（32）

式中：$h({\mathit{x}}_I,t)=\mathrm{0},{\mathit{x}}_I\in {\Gamma }_h^p$。

3 模型验证

为验证本文提出的非局部流体模型在模拟不连续面上流体力学行为的精确性，采用二维高渗透差问题作为验证模型。几何模型和边界如图2所示，模型的长、宽各为0.02m，流体的黏度系数取$\mu =\mathrm{1}\mathrm{k}\mathrm{g}\cdot {\mathrm{m}}^{-\mathrm{1}}\cdot {\mathrm{s}}^{-\mathrm{1}}$。模型的渗透性分布由式（33）控制：

$$\kappa =\left\{\begin{array}{cc}{\kappa }_{\mathrm{1}},& r\le r_{\mathrm{o}}\\ {\kappa }_{\mathrm{o}},& r>r_{\mathrm{o}}\end{array}\right.$$

（33）

式中：$r=\sqrt[]{x^{\mathrm{2}}+y^{\mathrm{2}}}$；$r_{\mathrm{o}}=\mathrm{0.005}\mathrm{m}$；${\kappa }_{\mathrm{1}}=\mathrm{0.001}{\mathrm{m}}^{\mathrm{2}}$；${\kappa }_{\mathrm{o}}=\mathrm{1}{\mathrm{m}}^{\mathrm{2}}$。

图2  高渗透差问题模型示意

Fig. 2  2D porous media with a high-contrast permeability

假设整个模型中存在单一流体源$q(r)=-\mathrm{9}r\mathrm{k}\mathrm{g}\cdot {\mathrm{m}}^{-\mathrm{1}}\cdot {\mathrm{s}}^{-\mathrm{1}}$，压力边界条件为$p(r)=\frac{r^{\mathrm{3}}}{{\kappa }_{\mathrm{o}}}+(\frac{\mathrm{1}}{{\kappa }_{\mathrm{1}}}-{\frac{\mathrm{1}}{\kappa }}_{\mathrm{o}})r_{\mathrm{o}}^{\mathrm{3}}\mathrm{P}\mathrm{a}$。在上述条件下，局部流体模型的解析解定义为

$$p(r)=\left\{\begin{array}{cc}\frac{r^{\mathrm{3}}}{{\kappa }_{\mathrm{1}}},& r\le r_{\mathrm{o}}\\ \frac{r^{\mathrm{3}}}{{\kappa }_{\mathrm{o}}}+(\frac{\mathrm{1}}{{\kappa }_{\mathrm{1}}}-{\frac{\mathrm{1}}{\kappa }}_{\mathrm{o}})r_{\mathrm{o}}^{\mathrm{3}},& r>r_{\mathrm{o}}\end{array}\right.$$

（34）

在数值计算中，物质点间距取$\mathrm{\Delta }X=\mathrm{\Delta }Y=\mathrm{2.5}\times {\mathrm{10}}^{-\mathrm{5}}\mathrm{m}$，物质点作用半径为$\delta =\mathrm{2}\sqrt[]{\mathrm{\Delta }{X}^{\mathrm{2}}+\mathrm{\Delta }{Y}^{\mathrm{2}}}$，罚函数因子$\alpha =\mathrm{5}\kappa$。图3a为非局部模型计算结果，图3b为数值计算结果与解析解误差。由图3可知，本文提出的非局部流体模型能精确捕捉流体在穿越物质界面上的不连续力学行为，同时也证明了非局部模型解空间中内蕴的不连续性，及其在模拟含不连续界面多孔介质流体传输问题上的有效性。

图3  非局部模型数值计算结果与解析解答对比

Fig. 3  Comparison of the numerical results of the nonlocal model with analytical solutions

4 三维含孔隙非均匀多孔介质中的流体传输模拟

由于多孔介质的非均匀性与孔隙分布特征，三维情况下流体穿越不连续界面的力学行为更加复杂，因此研究三维条件多孔介质流体传输问题有重要意义。本算例考虑长$L=\mathrm{1}\mathrm{m}$、宽$W=\mathrm{1}\mathrm{m}$、高$H=\mathrm{1}\mathrm{m}$的多孔介质，在多孔介质两侧各有一个半圆形空腔，几何模型和边界如图4所示。

图4  含半圆形空腔三维多孔介质、多孔几何模型示意^［23］

Fig. 4  Schematic of the geometric model for 3D porous media containing cavities^[23]

流体黏度$\mu =\mathrm{2.5}\mathrm{k}\mathrm{g}\cdot {\mathrm{m}}^{-\mathrm{1}}\cdot {\mathrm{s}}^{-\mathrm{1}}$，模拟区域的孔隙度$\varphi =\mathrm{0.55}$，基质渗透性为${\kappa }_{\mathrm{p}}=\mathrm{1}{\mathrm{m}}^{\mathrm{2}}$，并随机分布渗透性为${\kappa }_{\mathrm{n}}=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{4}}{\mathrm{m}}^{\mathrm{2}}$的颗粒物，定义流体源分布函数为

$$q(x,y,t)={\mathrm{e}}^{-\mathrm{9}\frac{\sqrt[]{x^{\mathrm{2}}+y^{\mathrm{2}}+z^{\mathrm{2}}}}{\sqrt[]{L^{\mathrm{2}}+W^{\mathrm{2}}+H^{\mathrm{2}}}}}\mathrm{c}\mathrm{o}\mathrm{s}(\mathrm{9}t)$$

（35）

对模拟区域施加压力边界条件$p(x,y,z,t)=\mathrm{0}$，物质点间距取$\mathrm{\Delta }X=\mathrm{\Delta }Z=\mathrm{\Delta }Y=\mathrm{0.01}\mathrm{m}$，物质点非局部作用半径取$\delta =\mathrm{1.25}\sqrt[]{\mathrm{\Delta }{X}^{\mathrm{2}}+\mathrm{\Delta }{Y}^{\mathrm{2}}+\mathrm{\Delta }{Z}^{\mathrm{2}}}$，罚函数因子设为$\alpha =\frac{\mathrm{1.0}\times {\mathrm{10}}^{-\mathrm{3}}{\kappa }_{\mathrm{p}}}{\mu }$。图5为t为1.25s、2.50s、3.75s、5.00s时非均匀饱和多孔介质中非局部流量和流体速度矢量场分布情况。可以看到，由于低渗透性颗粒物质的存在，流体将沿着这些颗粒形成复杂“通道”进行流动，从而数值模拟结果进一步说明了所提出的非局部模型在描述流体传输不连续力学行为方面的有效性。

图5  非均匀饱和多孔介质中非局部流量和流体速度矢量场分布

Fig.5  Nonlocal fluid flux and velocity field distribution in the saturated heterogeneous porous media containing two cavities at a time

5 结语

基于非常规态型近场动力学模型及其变分原理创建了非均匀多孔介质流体传输模拟的强–弱不连续统一非局部模型，能有效模拟多孔介质中流体在穿越孔隙、裂隙和物质界面时的不连续力学行为，通过与解析解对比验证了模型的精确性，并基于三维复杂多孔介质中流体传输的模拟分析了多孔介质中流体传输时的不连续力学行为。得到的结论如下：

（1） 基于近场动力学力态函数提出了流量态函数和相应的流体非局部控制方程，可以描述多孔介质中流体传输的非局部效应，并有效捕捉了多孔介质中强–弱不连续处流体的输运，且当物质点的非局部作用趋于零时，流体的非局部控制方程将退化为经典连续介质力学中的Darcy流体模型，从而建立了多孔介质流体传输模拟的局部模型与非局部模型之间的联系。

（2） 建立的非局部流体输运方程具有与经典连续介质力学统一的变分原理，不仅适用于无网格离散，也同样适用于基于网格的数值方法，为耦合局部模型与非局部模型提供了新的耦合途径。

（3） 针对非局部模型存在的零能模态问题，提出了一种罚函数方法并推导出全隐式数值方法，能有效消除非局部模型数值计算存在的数值震荡问题，精确模拟了多孔介质中流体传输的非连续力学行为，为研究流体在多孔介质中不连续界面传输的力学行为提供了理论依据。

作者贡献声明

孙宇麒：模型构建、数据分析呈现及论文撰写。

禹海涛：项目负责人、指导模型构建及数据分析，论文修改。