孫宇麒，禹海濤

（1.同濟(jì)大學(xué)土木工程學(xué)院，上海 200092；2.成都理工大學(xué)地質(zhì)災(zāi)害防治與地質(zhì)環(huán)境保護(hù)國家重點(diǎn)實(shí)驗(yàn)室，四川成都 610059；3.同濟(jì)大學(xué)土木工程防災(zāi)國家重點(diǎn)實(shí)驗(yàn)室，上海 200092）

巖土體是含有孔隙和裂隙的非均勻多孔介質(zhì)，這類多孔介質(zhì)中流體傳輸過程需要穿越由介質(zhì)中孔隙、裂隙和物質(zhì)界面組成的強(qiáng)-弱不連續(xù)面，這與經(jīng)典連續(xù)介質(zhì)力學(xué)模型中描述局部流體模型的偏微分方程是不相容的，從而導(dǎo)致需要額外的數(shù)值技術(shù)來描述流體在這些不連續(xù)界面的力學(xué)行為[1-3］。

目前興起的近場動(dòng)力學(xué)（Peridynamics）理論為該問題分析提供了一種新的思路。近場動(dòng)力學(xué)是由Silling等[4-5］提出的一種新的非局部理論，其控制方程是積分方程，考慮連續(xù)體內(nèi)物質(zhì)點(diǎn)之間非局部相互作用，避免了經(jīng)典連續(xù)介質(zhì)模型中偏微分方程帶來的空間求導(dǎo)奇異性[6-7］。目前近場動(dòng)力學(xué)已經(jīng)被廣泛應(yīng)用于包括脆性材料斷裂[8-9］、黏彈性斷裂[10-11］和多物理場耦合斷裂[12］等各類非連續(xù)問題。將原本描述斷裂等不連續(xù)問題的近場動(dòng)力學(xué)理論應(yīng)用于飽和多孔介質(zhì)流體傳輸問題，其關(guān)鍵在于如何建立多孔介質(zhì)中流體傳輸?shù)姆蔷植啃问绞睾懵煞匠毯捅緲?gòu)律方程。

近期已有一些學(xué)者通過修改近場動(dòng)力學(xué)方程以建立孔隙介質(zhì)中流體傳輸?shù)目刂品匠?，如Ni等[13］和Bobaru等[14］首先將鍵型近場動(dòng)力學(xué)模型中能量對應(yīng)的思想擴(kuò)展到熱傳導(dǎo)和流體傳輸問題，建立了熱-流耦合的非局部模型；同樣地，Oterkus等[15］也建立了一種熱-流耦合的非局部鍵型近場動(dòng)力學(xué)模型；隨后Katiyar等[16-17］基于常規(guī)態(tài)型近場動(dòng)力學(xué)模型推導(dǎo)出多孔介質(zhì)中流體傳輸?shù)姆蔷植啃问娇刂品匠?，并研究了流體穿越介質(zhì)中不連續(xù)界面的力學(xué)行為。然而，上述研究主要針對鍵型和常規(guī)態(tài)型近場動(dòng)力學(xué)模型，由于物質(zhì)點(diǎn)之間鍵的相互作用是獨(dú)立的，從而導(dǎo)致其推導(dǎo)的非局部擴(kuò)散模型無法描述流體傳輸?shù)姆蔷€性力學(xué)行為。此外，目前所發(fā)展的非局部傳輸模型主要基于鍵型模型中能量對應(yīng)的思想，與經(jīng)典連續(xù)介質(zhì)力學(xué)模型之間的對應(yīng)關(guān)系不夠清晰。

基于統(tǒng)一變分近場動(dòng)力學(xué)[18］的基本思想，本文提出一種表征飽和多孔介質(zhì)中流體非局部傳輸作用的流量向量態(tài)函數(shù)，統(tǒng)一描述流體在強(qiáng)-弱不連續(xù)面?zhèn)鬏斶^程中的力學(xué)行為，避免局部力學(xué)模型在不連續(xù)界面處出現(xiàn)的導(dǎo)數(shù)無定義性，并建立相應(yīng)非局部形式的質(zhì)量守恒方程、線動(dòng)量守恒方程和本構(gòu)方程。從理論證明了當(dāng)非局部流體模型中的非局部作用半徑趨于零時(shí)，該非局部模型將退化為經(jīng)典連續(xù)介質(zhì)力學(xué)中的Darcy流體模型。另外，提出一種基于罰函數(shù)方法的全隱式數(shù)值求解格式，旨在消除非局部模型內(nèi)在的零能模式難題，保證數(shù)值結(jié)果的精確性。

1 滲流傳輸?shù)姆蔷植靠刂品匠?/h2>

1.1 非局部質(zhì)量守恒方程

在傳統(tǒng)局部模型中，水壓力梯度是驅(qū)動(dòng)多孔介質(zhì)中流體傳輸?shù)闹饕蛩?，故質(zhì)量守恒方程以偏微分方程形式呈現(xiàn)?；诜浅Ｒ?guī)態(tài)型近場動(dòng)力學(xué)中的力態(tài)公式[6-7］，本文提出了流量態(tài)向量的概念，并給出了非均勻巖土介質(zhì)中流體傳輸?shù)姆蔷植啃问劫|(zhì)量守恒方程。

圖1為非均勻飽和多孔介質(zhì)中滲流輸運(yùn)的非局部模型示意。

圖1 非均勻飽和多孔介質(zhì)中滲流輸運(yùn)的非局部模型示意Fig.1 Schematic of the nonlocal fluid transport model in the heterogeneous saturated porous media

假設(shè)多孔介質(zhì)中流體處于恒溫條件下的弱可壓縮狀態(tài)，流體密度和壓力之間滿足關(guān)系如式（1）：

式中：ρ(x)為流體密度，kg·m-3；x為流體物質(zhì)點(diǎn)的位置矢量，m；t為時(shí)間變量，s；ρo(x)為參考流體密度，kg·m-3；cf為流體壓縮系數(shù)，Pa-1；p(x，t)為流體壓力，Pa；po(x)為參考流體壓力，Pa。流體傳輸?shù)姆蔷植抠|(zhì)量守恒方程可表示為

式中：-U[x′，t]x′-x和-U[x，t]x-x′分別表示在流體物質(zhì)點(diǎn)x′和x處的流量態(tài)矢量，kg·m-6·s-1；Hx表示流體物質(zhì)點(diǎn)x的非局部作用區(qū)域；q(x，t)為所考慮多孔介質(zhì)區(qū)域內(nèi)流體源的流量函數(shù)，kg·m-3·s-1.需要說明的是，流量態(tài)矢量中x′-x和x-x′表征了流體x′和x之間的傳輸作用。

流體物質(zhì)點(diǎn)之間非局部傳輸作用的流量態(tài)矢量定義如下：

式中：u(x，t)為流體流速，m·s-1；‖x-x‖′為x′和x之間的歐氏距離；ω(‖x′-x‖)為核函數(shù)，其定義如式（4）：

式中：x′的積分區(qū)域?yàn)閤的非局部作用區(qū)域Hx。

1.2 非局部線動(dòng)量守恒方程與本構(gòu)方程

為了得到完整描述非均勻多孔介質(zhì)中流體傳輸?shù)姆蔷植磕Ｐ?，還需定義流體傳輸?shù)姆蔷植烤€動(dòng)量守恒方程和本構(gòu)方程，其非局部動(dòng)量守恒方程定義為

式中：μ為流體黏度；-P[x]為x處的非局部有效滲透率張量，由式（7）非局部本構(gòu)方程確定。

式中：Pab為飽和巖土介質(zhì)的絕對滲透率張量；kr(x)和kr(x′)分別為x′和x的相對滲透率系數(shù)。?Nx?-[x，t]為流體勢函數(shù)的非局部梯度，如式（8）：

式中：-?[x′，t]x′-x和-?[x，t]x-x′為x′和x處的非局部流體勢標(biāo)量態(tài)函數(shù)。

式中：?(x，t)=p(x，t)-ρ(x)gz(x)為物質(zhì)點(diǎn)位置矢量的局部勢流函數(shù)，本文忽略由于重力產(chǎn)生的流體勢函數(shù)ρ(x)gz(x)。

式（1）、式（5）和式（6）構(gòu)成了非均勻飽和多孔介質(zhì)中流體傳輸?shù)娜靠刂品匠?，為了求解該方程，還需引入非局部壓力邊界Γpg和流量邊界Γph，如圖1所示，采用文獻(xiàn)[16］所提出的邊界施加方法。

1.3 與經(jīng)典Darcy流模型的一致性證明

本文的非局部滲流模型引入非局部作用半徑來表征流體物質(zhì)點(diǎn)之間的傳輸作用。接著理論證明，當(dāng)物質(zhì)點(diǎn)之間的作用半徑趨于零時(shí)，非局部模型退化為局部的Darcy流體模型。

基于Taylor公式，Hx中物質(zhì)點(diǎn)處x′的局部流量

式中：?x(ρ(x，t)u(x，t))表示x處的流量梯度；r(‖x′-x‖2)表示Taylor公式中的高階余項(xiàng)。

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

式中：Δ(ρu)(x′，x，t)=ρ(x′，t)u(x′，t)-ρ(x，t)u(x，t)為x′和x處流量的差值；?x(ρ(x)u(x，t))為流量函數(shù)的局部梯度。從式（11）可知，當(dāng)物質(zhì)點(diǎn)的非局部作用區(qū)域Hx的半徑δx趨于零時(shí)，非局部質(zhì)量守恒方程將退化為局部質(zhì)量守恒方程。由此定義流量的非局部梯度算子為

同理，對Hx中x′的勢流函數(shù)?(x′，t)進(jìn) 行Taylor展開，可以得到

將式（12）代入式（7），可得滲流的局部線動(dòng)量守恒方程和非局部質(zhì)量守恒方程如式（14）：

式中：Δ?(x′，x，t)=?(x′，t)-?(x，t)為x′和x處流體勢函數(shù)的差值；?x?(x，t)為流體勢函數(shù)在x處的具體梯度。注意到x處的非局部滲透率張量和局部滲透率張量有如式（15）的關(guān)系：

結(jié)合式（13）和式（14）可知，當(dāng)物質(zhì)點(diǎn)的非局部作用半徑δx趨于零時(shí)，非局部的線動(dòng)量守恒方程、本構(gòu)方程與局部線動(dòng)量守恒方程、本構(gòu)方程可通過如式（16）等式建立對應(yīng)關(guān)系。

2 非局部變分框架與全隱式數(shù)值求解格式

2.1 非局部變分公式與物質(zhì)點(diǎn)離散

不連續(xù)伽遼金方法[19-20］和物質(zhì)點(diǎn)離散方法是高效求解非局部模型的有效數(shù)值離散方法。將待求解多孔介質(zhì)區(qū)域離散為Np個(gè)物質(zhì)點(diǎn){xI}Np I=1，每個(gè)物質(zhì)點(diǎn)的體積為ΔVI，則待求解的多孔介質(zhì)區(qū)域Ω可用物質(zhì)點(diǎn)表征為

xI的非局部作用區(qū)域HxI內(nèi)，其物質(zhì)點(diǎn)可用局部編號表示為

與有限元方法求解偏微分方程的離散弱形式相同，定義離散試驗(yàn)函數(shù)空間和離散測試函數(shù)空間如式（19）：

則非局部傳輸模型的離散弱形式可表示為

式中：δp=δp(x，t)為壓力的變分函數(shù)；gp(xI，t)為非局部壓力邊界Γpg上的壓力函數(shù)；hp=hp(x，t)為

式中：下標(biāo)dis表示歐氏空間的離散內(nèi)積；ΔVxI為xI的體積；N?p為非局部流量邊界Γph內(nèi)的物質(zhì)點(diǎn)數(shù)；

2.2 零能模態(tài)控制

由于非局部梯度（11）內(nèi)在的零能模式將導(dǎo)致數(shù)值解答產(chǎn)生嚴(yán)重的數(shù)值震蕩問題[21］，通過引入罰函數(shù)控制零能模式的出現(xiàn)，定義xI處的罰函數(shù)如式（23）：

式中：λi(xI)(i=1，...，nd)為物質(zhì)點(diǎn)xI的形狀張量K(xI)的特征值；α為罰函數(shù)因子，取1×10-3～1×103。zp(xJ，xI，t)為xI和xJ之間非一致部分壓力變化，定義如式（25）：

對離散罰函數(shù)（22）分別求一階和二階變分導(dǎo)數(shù)，可以得到罰函數(shù)荷載向量和罰函數(shù)剛度矩陣

2.3 隱式數(shù)值求解格式

為能精確模擬流體在介質(zhì)中界面、裂隙和孔隙處傳輸行為，在推導(dǎo)罰函數(shù)荷載向量和剛度矩陣的基礎(chǔ)上，采用Newton-Raphson隱式方法進(jìn)行求解[22-23］。結(jié)合式（19）和式（27），可得最終求解的非線性方程如式（29）：

式中：Δtn為時(shí)間步長；pn=[p1(tn)，p2(tn)，...，pNp(tn)]T。矩陣B和C(p)、向量f(p)的定義為

其中BI、Gp，I和DI定義為

3 模型驗(yàn)證

為驗(yàn)證本文提出的非局部流體模型在模擬不連續(xù)面上流體力學(xué)行為的精確性，采用二維高滲透差問題作為驗(yàn)證模型。幾何模型和邊界如圖2所示，模型的長、寬各為0.02m，流體的黏度系數(shù)取μ=1kg·m-1·s-1。模型的滲透性分布由式（33）控制：

圖2 高滲透差問題模型示意Fig.2 2D porous media with a high-contrast permeability

假設(shè)整個(gè)模型中存在單一流體源q(r)=

在數(shù)值計(jì)算中，物質(zhì)點(diǎn)間距取ΔX=ΔY=2.5×10-5m，物 質(zhì) 點(diǎn) 作 用 半 徑 為罰函數(shù)因子α=5κ。圖3a為非局部模型計(jì)算結(jié)果，圖3b為數(shù)值計(jì)算結(jié)果與解析解誤差。由圖3可知，本文提出的非局部流體模型能精確捕捉流體在穿越物質(zhì)界面上的不連續(xù)力學(xué)行為，同時(shí)也證明了非局部模型解空間中內(nèi)蘊(yùn)的不連續(xù)性，及其在模擬含不連續(xù)界面多孔介質(zhì)流體傳輸問題上的有效性。

圖3 非局部模型數(shù)值計(jì)算結(jié)果與解析解答對比Fig.3 Comparison of the numerical results of the nonlocal model with analytical solutions

4 三維含孔隙非均勻多孔介質(zhì)中的流體傳輸模擬

由于多孔介質(zhì)的非均勻性與孔隙分布特征，三維情況下流體穿越不連續(xù)界面的力學(xué)行為更加復(fù)雜，因此研究三維條件多孔介質(zhì)流體傳輸問題有重要意義。本算例考慮長L=1m、寬W=1m、高H=1m的多孔介質(zhì)，在多孔介質(zhì)兩側(cè)各有一個(gè)半圓形空腔，幾何模型和邊界如圖4所示。

圖4 含半圓形空腔三維多孔介質(zhì)、多孔幾何模型示意[23］Fig.4 Schematic of the geometric model for 3D porous media containing cavities[23]

流體黏度μ=2.5kg·m-1·s-1，模擬區(qū)域的孔隙度?=0.55，基質(zhì)滲透性為κp=1m2，并隨機(jī)分布滲透性為κn=1×10-4m2的顆粒物，定義流體源分布函數(shù)為

5 結(jié)語

圖5 非均勻飽和多孔介質(zhì)中非局部流量和流體速度矢量場分布Fig.5 Nonlocal fluid flux and velocity field distribution in the saturated heterogeneous porous media containing two cavities at a time

基于非常規(guī)態(tài)型近場動(dòng)力學(xué)模型及其變分原理創(chuàng)建了非均勻多孔介質(zhì)流體傳輸模擬的強(qiáng)-弱不連續(xù)統(tǒng)一非局部模型，能有效模擬多孔介質(zhì)中流體在穿越孔隙、裂隙和物質(zhì)界面時(shí)的不連續(xù)力學(xué)行為，通過與解析解對比驗(yàn)證了模型的精確性，并基于三維復(fù)雜多孔介質(zhì)中流體傳輸?shù)哪M分析了多孔介質(zhì)中流體傳輸時(shí)的不連續(xù)力學(xué)行為。得到的結(jié)論如下：

（1）基于近場動(dòng)力學(xué)力態(tài)函數(shù)提出了流量態(tài)函數(shù)和相應(yīng)的流體非局部控制方程，可以描述多孔介質(zhì)中流體傳輸?shù)姆蔷植啃?yīng)，并有效捕捉了多孔介質(zhì)中強(qiáng)-弱不連續(xù)處流體的輸運(yùn)，且當(dāng)物質(zhì)點(diǎn)的非局部作用趨于零時(shí)，流體的非局部控制方程將退化為經(jīng)典連續(xù)介質(zhì)力學(xué)中的Darcy流體模型，從而建立了多孔介質(zhì)流體傳輸模擬的局部模型與非局部模型之間的聯(lián)系。

（2）建立的非局部流體輸運(yùn)方程具有與經(jīng)典連續(xù)介質(zhì)力學(xué)統(tǒng)一的變分原理，不僅適用于無網(wǎng)格離散，也同樣適用于基于網(wǎng)格的數(shù)值方法，為耦合局部模型與非局部模型提供了新的耦合途徑。

（3）針對非局部模型存在的零能模態(tài)問題，提出了一種罰函數(shù)方法并推導(dǎo)出全隱式數(shù)值方法，能有效消除非局部模型數(shù)值計(jì)算存在的數(shù)值震蕩問題，精確模擬了多孔介質(zhì)中流體傳輸?shù)姆沁B續(xù)力學(xué)行為，為研究流體在多孔介質(zhì)中不連續(xù)界面?zhèn)鬏數(shù)牧W(xué)行為提供了理論依據(jù)。

作者貢獻(xiàn)聲明：

孫宇麒：模型構(gòu)建、數(shù)據(jù)分析呈現(xiàn)及論文撰寫。

禹海濤：項(xiàng)目負(fù)責(zé)人、指導(dǎo)模型構(gòu)建及數(shù)據(jù)分析，論文修改。