彈塑性多孔介質(zhì)流固耦合新理論:混合耦合理論

徐麗陽, 王 鍇, 丁 智, 徐日慶, 陳曉輝

(1.浙江大學(xué) 濱海和城市巖土工程研究中心,杭州 310058;2.浙江省城市地下空間開發(fā)工程技術(shù)研究中心,杭州 310058;3.浙大城市學(xué)院 土木工程系,杭州 310015;4.浙江省城市盾構(gòu)隧道安全建造與智能養(yǎng)護(hù)重點(diǎn)實(shí)驗(yàn)室,杭州 310015;5.北京師范大學(xué) 水科學(xué)研究院,北京 100875;6.利茲大學(xué) 土木工程系,英國利茲 LS2 9JT)

1 引 言

在全球氣候變化和雙碳政策的大背景下,亟需發(fā)展新能源和減碳減排技術(shù),如開發(fā)利用地?zé)帷㈨搸r氣、天然氣、水合物、核電和地質(zhì)封存二氧化碳[1]。使得多孔介質(zhì)中固體的變形和流體的輸運(yùn)問題受到越來越多研究者的廣泛關(guān)注[2,3]。然而,在多孔介質(zhì)中建立多場耦合模型是復(fù)雜的,面臨的主要挑戰(zhàn)是,多孔介質(zhì)內(nèi)部孔隙的尺度范圍一般在10-9m～10-2m,而固體顆粒的尺度范圍一般在宏觀尺度,如黏土顆粒直徑通常小于0.002 mm (2 μm),有時(shí)會(huì)小到納米級(jí)別,具有高度的塑性和黏性。因此,多孔介質(zhì)中的流固耦合現(xiàn)象跨越了宏觀尺度到納米尺度,即不僅要考慮宏觀尺度(>10-7m)的流固耦合行為[4],如泥石流、軟土的彈性變形、塑性變形和大孔隙中水的滲流;還需要考慮納米級(jí)(10-9m～10-7m)的流固耦合作用[5],如軟土微觀孔隙結(jié)構(gòu)的變化、吸水膨脹、毛細(xì)作用和吸附等。如果忽略這些微細(xì)觀作用,對(duì)力場和流場的預(yù)測會(huì)產(chǎn)生較大誤差,難以解釋一些現(xiàn)象,如白堊巖遇水弱化[6](水滲透到白堊巖中時(shí),水的滲透力可能會(huì)在顆粒之間產(chǎn)生沖擊力,導(dǎo)致顆粒分離和破碎,從而降低了巖石的強(qiáng)度和穩(wěn)定性)、北海Ekofisk油田的海底沉降[7](高孔隙率白堊的塑性變形使得Ekofisk油藏壓實(shí),造成了上覆海底的沉降)。

開展多場耦合分析有許多種思路,可歸納為力學(xué)方法和混合理論方法兩類主要方法。力學(xué)方法是基于Terzaghi有效應(yīng)力原理[8,9]和Biot固結(jié)理論[10]進(jìn)行建模,其本質(zhì)是宏觀唯象關(guān)系。但對(duì)納米級(jí)的多場耦合行為,力學(xué)方法只能基于經(jīng)驗(yàn)和試驗(yàn)結(jié)果進(jìn)行大量直覺性假設(shè),因此在數(shù)學(xué)和物理上都不夠嚴(yán)謹(jǐn)[11,12]?；旌侠碚摲椒ㄓ蒚ruesdell等[13,14]提出,該方法基于熱力學(xué)建模,為多孔介質(zhì)中的流體輸運(yùn)行為提供了更精確的理論方程。由于該方法假設(shè)固相和流相為獨(dú)立體,需要對(duì)各相建立連續(xù)方程和動(dòng)量守恒方程,并確定相界面上各相間的相互作用。然而納米級(jí)的各相由于尺度太小,相間的作用力是很難獲得的,只能依賴主觀經(jīng)驗(yàn)和特定的假設(shè),限制了混合理論方法進(jìn)一步考慮納米尺度的作用[15]。

為了避免力學(xué)方法和混合理論方法的局限性,在多場耦合分析中實(shí)現(xiàn)考慮納米尺度的作用,Chen等[16,17]提出了混合耦合理論(Mixture-Coupling Theory)。該理論源自Heidug等[18]對(duì)具有吸水膨脹效應(yīng)的飽和多孔介質(zhì)的研究,基于多孔介質(zhì)力學(xué)和非平衡熱力學(xué),利用Helmholtz自由能和耗散熵產(chǎn)推導(dǎo)出固體變形能與耗散過程之間的關(guān)系[19]?；旌像詈侠碚摰膬?yōu)點(diǎn)是具有熱力學(xué)一致性且推導(dǎo)方法嚴(yán)謹(jǐn),可以建立不同尺度的平滑聯(lián)系,適合于構(gòu)建多尺度多場復(fù)雜耦合模型[20-23]。特別地,混合耦合理論還可以考慮大變形問題,建立彈性、塑性和固體流(泥石流)的統(tǒng)一理論。

然而,目前基于混合耦合理論建立的多場耦合模型,均假設(shè)固體的變形是彈性的,尚未考慮固體骨架塑性變形的影響。當(dāng)固體骨架發(fā)生塑性變形時(shí),固體顆粒會(huì)重新排列(固體顆粒的尺度范圍為<0.002 m,少量在納米尺度),儲(chǔ)存在固體顆粒之間的自由能會(huì)變化[24,25];同時(shí),固體顆粒重新排列,會(huì)改變流體流動(dòng)的孔隙通道(10-9m～10-2m),造成孔壓的變化。這些變化一方面引起宏細(xì)觀的固體骨架有效應(yīng)力變化和微觀的孔隙結(jié)構(gòu)的變化;另一方面又反過來影響微觀的孔隙流體流動(dòng)和壓力的分布[26]。

本文首次將混合耦合理論從彈性變形假設(shè)推廣到彈塑性變形,推導(dǎo)了飽和多孔介質(zhì)的彈塑性流固耦合控制方程(即HM模型),模型實(shí)現(xiàn)了考慮從宏觀尺度到納米尺度的流固耦合行為,也為后續(xù)在彈塑性多孔介質(zhì)中耦合更多納米尺度的物理化學(xué)作用提供了參考。

2 平衡方程

混合耦合理論是基于非平衡熱力學(xué)的方法,假設(shè)滿足局域平衡條件,即飽和多孔介質(zhì)整體雖然是變化的,但是局部可以近似視為平衡,并滿足平衡體系中的熱力學(xué)關(guān)系。在飽和多孔介質(zhì)中任意選擇代表單元體(REV)(圖1),要求REV在宏觀尺度上足夠小,則其內(nèi)部近似均勻;REV在微觀尺度上要足夠大,包含有足夠多的分子,以進(jìn)行統(tǒng)計(jì)分析。因此,飽和多孔介質(zhì)系統(tǒng)整體的描述方程可采用局部形式的微分方程表示。假設(shè)REV具有體積V(單位:m3),由固體骨架、孔隙流體和邊界S(單位:m2)組成,只允許孔隙流體的流入/流出。

圖1 代表單元體(REV)

2.1 質(zhì)量平衡方程

飽和多孔介質(zhì)中固相和液相的質(zhì)量平衡方程可以統(tǒng)一表示為

(1)

(2)

式中 ?t為時(shí)間導(dǎo)數(shù),為梯度算子。因此,由式(1)可以得到固相質(zhì)量平衡方程的局部形式為

(3)

液相質(zhì)量平衡方程的局部形式為

(4)

式中Iw為液體單位面積的質(zhì)量通量(單位:kg/(m2·s)-1)。

(5)

2.2 Helmholtz自由能平衡方程

2.2.1 內(nèi)能平衡

飽和多孔介質(zhì)的內(nèi)能平衡方程為

(6)

則式(6)的局部形式為

(7)

2.2.2 熵平衡

飽和多孔介質(zhì)的熵平衡方程為

(8)

式中Iη為REV與周圍環(huán)境交換的熵通量(單位:J/(K·m2·s)),γ為單位體積單位時(shí)間內(nèi)的熵產(chǎn)[28](單位:J/(K·m3·s))。方程(8)的局部形式就是REV的熵密度平衡方程,

(9)

2.2.3 自由能平衡方程

(10)

3 本構(gòu)關(guān)系

3.1 熵產(chǎn)

下一步需要確定REV的耗散過程,即確定方程(9)的γ。本文假設(shè)固體經(jīng)歷了不可逆的彈塑性變形,因此,REV中有兩個(gè)耗散機(jī)制,即固體骨架的塑性變形,以及當(dāng)流體通過多孔骨架時(shí)在固體/流體相邊界產(chǎn)生的摩擦。其中,土壤在受力變形時(shí),骨架顆粒重新排列產(chǎn)生功,一部分功是可恢復(fù)的彈性功(當(dāng)土壤受的外力解除后,彈性功會(huì)完全釋放),另一部分功是塑性功。假設(shè)所有的塑性功都會(huì)耗散,參考文獻(xiàn)[29-31]的研究得到塑性變形耗散熵產(chǎn)的表達(dá)式

(11)

0≤Tγ2=-Iw·μw=-u·pw

(12)

因此,耗散過程的總熵產(chǎn)為

(13)

(14)

(15)

(16)

(17)

綜上,得到耦合矩陣如下,

(18)

(19)

方程(19)的σd表達(dá)式與經(jīng)典的耗散應(yīng)力與塑性應(yīng)變率的關(guān)系[31,33]一致,即較高的應(yīng)變率通常會(huì)導(dǎo)致材料更快地發(fā)生變形,則固體顆粒更容易發(fā)生滑移,從而產(chǎn)生更多的耗散力;方程(19)的u表達(dá)式與經(jīng)典的達(dá)西定律完全一致。以上驗(yàn)證了方程(14,15)的正確性。

3.2 連續(xù)介質(zhì)力學(xué)方程

依據(jù)經(jīng)典的連續(xù)介質(zhì)力學(xué)理論,可以定義固體骨架的變形狀態(tài)。一些基本的關(guān)系式為

T=JF-1σF-T,J=dV/dV0

(20)

式中F為變形梯度,X為參考組構(gòu)REV所在的位置,且在t時(shí)刻的位置(當(dāng)前位置)是x,E為Green應(yīng)變,T為第二Piola-Kirchhoff應(yīng)力(PK2),I為單位張量,J為函數(shù)的雅可比(Jacobian)矩陣,是函數(shù)輸出相對(duì)于其輸入的變化率,V0為參考組構(gòu)REV的體積,V為當(dāng)前組構(gòu)REV的體積。

假設(shè)試驗(yàn)土樣一直保持著力學(xué)平衡,因此,在不考慮外部體力的情況下,平衡方程為·σ=0。將方程(13)代入Helmholtz自由能密度平衡方程(10)得

(21)

方程(21)描述了當(dāng)前組構(gòu)REV的自由能平衡方程,在式(21)的兩邊乘以J,可推導(dǎo)得到參考組構(gòu)的自由能平衡方程為

(22)

3.3 孔隙水的Helmholtz自由能

根據(jù)經(jīng)典熱力學(xué)可推導(dǎo)得到飽和多孔介質(zhì)中,孔隙水的 Helmholtz自由能密度Ψpore為[15]

(23)

式中p為孔隙水壓力(單位:Pa)。根據(jù)Gibbs-Duhem方程,不考慮溫度變化,對(duì)方程進(jìn)行時(shí)間微分,推導(dǎo)得

(24)

將方程(23)兩邊對(duì)時(shí)間求偏導(dǎo),再將方程(24)代入,可得

(25)

3.4 固體骨架的Helmholtz自由能

(26)

3.5 本構(gòu)方程

方程(26)采用v作為變量,但在巖土工程中,使用p作為變量比較方便,因此采用Legendre變換為

(27)

方程(27)兩邊對(duì)時(shí)間求偏導(dǎo),再將方程(26)代入,可得方程為

(28)

W可理解為REV的總變形能,方程(28)表明,W為E,εp和p的函數(shù)。因此,T,σd和ν的表達(dá)式為

(29)

(30)

因此,應(yīng)力T、耗散應(yīng)力σd和孔隙體積分?jǐn)?shù)v的本構(gòu)方程為

(31)

(32)

(33)

式中 下標(biāo)i,j,k,l=1,2,3。參數(shù)Lijkl,Sijkl,Mij,Rijkl,Bij和Q分別計(jì)算為

(34)

4 彈塑性多孔介質(zhì)HM模型控制方程

方程(31～33)給出了Tij,σdij和v的一般形式耦合方程。這些方程適用于非線性、大變形和各向異性條件下的一般情況。本文主要用到Tij和v,為了簡化討論和關(guān)注土骨架塑性變形的影響,只討論方程(31,33),并作了如下假設(shè)。

(35)

(36)

ζ=1-(K/Ks)

(37)

此外,Q為孔隙的壓縮性能,即

(38)

式中 當(dāng)假設(shè)固體顆粒為不可壓縮時(shí)(Ks=∞),Q=0。

4.1 流固耦合控制方程

4.1.1 固相的控制方程

根據(jù)上述假設(shè),方程(31)可簡化為

(39)

由于REV保持力學(xué)平衡,不考慮體力,因此應(yīng)力滿足平衡條件?σij/?x=0,推導(dǎo)得到固相的控制方程為

(40)

式中D為彈塑性剛度系數(shù),由于本文研究的流固耦合行為需要考慮土骨架彈塑性變形,采用修正Cam-Clay模型描述彈塑性本構(gòu)關(guān)系,因此,D可以由式(41)計(jì)算[34](詳細(xì)推導(dǎo)過程見附錄)。

(41)

(42)

此時(shí),力學(xué)平衡方程包含有塑性的影響和流體的影響

4.1.2 液相的控制方程

(43)

(44)

(45)

(46)

4.1.3 創(chuàng)新分析

本文提出的HM本構(gòu)框架由方程(18,40,46)組成。這是一個(gè)統(tǒng)一的模型,將彈性變形、塑性變形和液體滲透都納入了混合耦合理論框架,實(shí)現(xiàn)了對(duì)從宏觀尺度到微觀尺度的流固耦合行為的表征與分析。可以用于不同的場景實(shí)現(xiàn)更準(zhǔn)確的預(yù)測,如模擬基坑降水、石油和天然氣在地下儲(chǔ)層中的流動(dòng)等,特別是固體顆粒尺度大量處于納米級(jí)的多孔介質(zhì)的彈塑性流固耦合情況。此外,方程(18)精確描述了固體塑性變形對(duì)塑性耗散應(yīng)力和液體輸運(yùn)方程的影響,方程(40,46)精確描述了固體塑性變形對(duì)力學(xué)平衡方程和流體均衡方程的影響。建立的HM模型具有準(zhǔn)確的推導(dǎo)過程和多種形式的比較驗(yàn)證,在數(shù)學(xué)和物理上都很嚴(yán)謹(jǐn)。

本文提出的HM-塑性模型與Chen等[15]采用非平衡熱力學(xué)方法提出的模型相比,力學(xué)變形方程和流體傳輸方程都有相似的形式,但是Chen等[15]的模型假設(shè)土體只發(fā)生彈性變形,沒有考慮土的塑性變形。同時(shí),如果不考慮塑性變形,可將方程(40,46)退化成為力學(xué)方法推導(dǎo)得到的經(jīng)典Biot方程[10]。此外,考慮塑性變形的HM方程與文獻(xiàn)[26,31,35]提出的方程相比,具有類似的形式。最后,本文的數(shù)值模擬部分對(duì)HM-塑性模型也進(jìn)行了驗(yàn)證。

5 數(shù)值模擬

5.1 數(shù)值模型簡介

基于Heidug等[18]設(shè)計(jì)的實(shí)驗(yàn)裝置,建立一個(gè)長度為0.2 m、高度為0.1 m的二維軸對(duì)稱模型(圖2)。模型中的試樣是飽和軟土(0.1 m×0.2 m),用上述建立的彈塑性流固耦合模型表征,試樣固定在兩個(gè)堅(jiān)硬且無摩擦的上下兩個(gè)邊界之間,因此試樣只能產(chǎn)生水平方向(x方向)的位移。選擇試樣的中心點(diǎn)作為觀測點(diǎn)C(圖2)?；贑OMSOL Multiphysics平臺(tái)采用傳統(tǒng)有限元方法計(jì)算。將彈塑性流固耦合模型的計(jì)算結(jié)果與李培超等[26]的彈塑性HM模型的計(jì)算結(jié)果進(jìn)行對(duì)比分析。同時(shí),選用劍橋模型表征軟土的力學(xué)行為,模型參數(shù)列入表1。

圖2 幾何模型和邊界條件

模型的邊界和初始條件如下。

(1) 邊界條件。邊界A是自由的,而邊界B是固定的(位移=0)。兩個(gè)邊界都是可透水的,可以讓水在飽和軟土試樣中流進(jìn)和流出。上下邊界是不可透水的固定邊界,不允許水的流動(dòng),也不能產(chǎn)生豎向位移。

(3) 模擬的開始。在t=0時(shí)刻,自由邊界A上的水壓力從5 MPa突然增加到10 MPa,并且在后續(xù)的分析時(shí)間段里(即t>0)都保持這個(gè)水壓力數(shù)值;固定邊界B的水壓力一直保持在5 MPa。

表1 考慮塑性的流固耦合模型參數(shù)

5.2 數(shù)值模擬結(jié)果

圖3是觀測點(diǎn)C的孔隙率、孔隙水壓力p、有效應(yīng)力和水平位移隨時(shí)間的發(fā)展趨勢(shì)。本文模型和李培超等[26]的模型計(jì)算結(jié)果是相似的。本文模型計(jì)算的孔隙率比文獻(xiàn)模型的略大(圖3(a)),體現(xiàn)了飽和軟土發(fā)生塑性變形時(shí)土顆粒重排,使得孔隙結(jié)構(gòu)變化。

該塑性變形會(huì)影響孔隙水壓力。由于李培超等[26]的模型中,忽略重力項(xiàng)的水輸運(yùn)方程為

(47)

外部荷載相同的情況下,根據(jù)有效應(yīng)力原理可知,孔隙水壓力降低,相應(yīng)的有效應(yīng)力將增加。如圖3(c)所示,本文模型的有效應(yīng)力比文獻(xiàn)模型略大;同時(shí),本文模型的水平位移也比文獻(xiàn)模型略大,如圖3(d)所示。

圖3 觀測點(diǎn)C的性質(zhì)參數(shù)隨時(shí)間的變化曲線

通過上述對(duì)比分析可知,如果忽略塑性變形的影響,對(duì)力場和流場的預(yù)測會(huì)產(chǎn)生較大誤差。這是由于多孔介質(zhì)中的流固耦合現(xiàn)象跨越了宏觀尺度到納米尺度,本文模型通過混合耦合理論推導(dǎo)得到的HM-塑性模型,可以考慮納米尺度的流固耦合作用。

6 結(jié) 論

本文提出了一個(gè)彈塑性多孔介質(zhì)流固耦合模型,在同一個(gè)理論框架內(nèi)研究了彈性變形、塑性變形和液體滲流,實(shí)現(xiàn)了將混合耦合理論擴(kuò)展到彈塑性變形范圍。基于非平衡熱力學(xué)的混合耦合理論,確定塑性變形時(shí)土骨架結(jié)構(gòu)變化土顆粒重排產(chǎn)生的耗散熵,并利用Helmholtz自由能來連接宏觀尺度上的力學(xué)變形和納米尺度上的液體輸運(yùn)之間的相互作用,通過唯象定律進(jìn)一步發(fā)展了考慮固體塑性變形的拓展達(dá)西定律。然后通過與文獻(xiàn)中模型的比較,驗(yàn)證了該模型的有效性。最后,數(shù)值分析表明塑性變形對(duì)多孔介質(zhì)中的流固耦合行為具有比較顯著的影響。

附 錄:

對(duì)于彈塑性本構(gòu)模型,[Dep]的一般表達(dá)式

(A1)

按相容條件

(A2)

(A3)

(A4)

代入式(A1),有

([De]-[Dp]){dε}=[Dep]{dε}

(A5)

式中 [Dep]=[De]-[Dp]

(A6)

(A7)

當(dāng)彈塑性本構(gòu)模型選用修正Cam-Clay模型時(shí),模型中,土的狀態(tài)用平均有效應(yīng)力p′、偏應(yīng)力q和孔隙比e共同描述。彈塑性剛度系數(shù)[Dep](即D(Dijkl))的推導(dǎo)過程如下。

(1) 總應(yīng)變。多孔介質(zhì)的總應(yīng)變包括彈性應(yīng)變和塑性應(yīng)變兩部分,即

ε=εe+εp

(A8)

式中ε為多孔介質(zhì)的總應(yīng)變,εe為多孔介質(zhì)的彈性應(yīng)變,εp為多孔介質(zhì)的塑性應(yīng)變。

(2) 彈性應(yīng)變?cè)隽?/p>

(A9,A10)

(3) 塑性應(yīng)變?cè)隽?/p>

① 屈服方程

(A11)

② 流動(dòng)法則。修正劍橋模型假設(shè)相關(guān)聯(lián)塑性,因此,塑性勢(shì)方程g和屈服方程f一樣

(A12)

因此,描述塑性應(yīng)變的流動(dòng)法則是

(A13)

(A14,A15)

式中λ為正常固結(jié)曲線的斜率。

④ 塑性應(yīng)變和應(yīng)力的關(guān)系。根據(jù)上述方程,可以推導(dǎo)得到塑性應(yīng)變和應(yīng)力的關(guān)系為

(A16)

(4) 修正的剛度D(Dijkl)。將方程(A9,A10,A16)代入方程(A8),可計(jì)算出總應(yīng)變?yōu)?/p>

(A17)

(A18)