高溫凍土融化固結的水熱力耦合模型構建1)

楊韜 郭穎 張程程 單煒

(東北林業大學，哈爾濱，150040)

隨著全球氣候變暖，多年凍土地區的平均地溫升高，位于多年凍土南界及高海拔多年凍土外緣的不連續多年凍土分布區，普遍出現了退化現象[1-5]。溫度的改變，使冰在凍土中的存在形式發生了變化，高溫凍土中的未凍水含量的增加，使凍土由脆性逐漸向塑性轉變[6-7]。因此，由于年平均地溫的升高，凍土的穩定性變差，熱敏感性提高，其上的建筑物及植被受到下部土體條件變化的影響也更為明顯[8-9]。

研究人員最早結合太沙基(K. Terzaghi)一維固結理論，采用以0 ℃的溫度界限為求解域邊界的移動邊界方法，求解凍土的一維融沉固結問題，但是在求解二維及三維融化固結問題時適用性較差[10]，Sykes et al.[11]通過結合三維固結理論及引入相變潛熱作為溫度場熱源等方法，對其進行了改進。其后，Yao et al.[12]考慮到溫度場參數在混合相變區是一個隨溫度變化的變量，并結合大應變理論，建立了三維大應變固結模型。Na et al.[13]應用有限應變下的廣義硬化準則，結合線性運動平衡、物質平衡、能量平衡建立了水熱力三項耦合模型，用以分析相變和傳熱作用下彈塑性響應和臨界狀態的演化過程。目前，對于凍土融化固結的研究，大部分采用的是以溫度場求得的某一溫度值為界限，溫度場與流固耦合非同步的求解方式，只有少數研究考慮了凍土區變形的水熱力三場同求解域求解[13]。

本研究依據傳統的比奧(Biot)固結理論，以未凍水為耦合點，結合非穩態的對流傳熱方程，通過引入高溫凍土應力-溫度耦合損傷模型，建立水熱力三個物理場同求解域的耦合模型，依據該模型分析高溫凍土融化固結過程。

1 研究方法

1.1 耦合模型的建立

1.1.1 擴展的比奧大變形固結模型

比奧理論依據太沙基有效應力原理，建立了孔隙水壓消散與土骨架變形的流固耦合關系[14]。凍土融化固結過程中，當施加荷載較大時，其即時構型與參考構型不再滿足小應變狀態時的近似重合假設[12]。因此，依據總拉格朗日(Total Lagrangian)方法建立大變形應力場控制方程，使模型可以更為準確地分析凍土融化固結過程中的幾何非線性問題[15-16]。

在求解域內，張量形式的平衡方程可用式(1)表示。

?Tij/?Xj+F0i=0，i=1、2、3，j=1、2、3。

(1)

式中：Tij為第一類皮奧拉-基爾霍夫(Piola-Kirchhoff)應力；Xj為物質坐標；F0i為初始時刻參考構型之上的體積荷載。由于第一類皮奧拉-基爾霍夫應力張量為非對稱張量，在后續計算中存在一定的不便。因此引入基矢量完全位于參考構型之上的第二類皮奧拉-基爾霍夫應力張量(Slk)對平衡關系進行分析，并將其轉換為物質坐標的表達式。

(?/?Xk)[Slk(δil+?ui/?Xl)]+F0i=0。

(2)

(3)

引入不含空間坐標格林(Green)應變張量，建立宏觀變形與微觀應變之間的關系。

εij=(1/2)[(δki+?uk/?Xi)(δkj+?uk/?Xj)-δij]=

(1/2)[?ui/?Xj+?uj/?Xi+(?uk/?Xi)(?uk/?Xj)]。

(4)

εij表示的是宏觀變形在微觀應變上的總體映射，分析宏觀變形發生的成因。εij由兩部分組成，應力作用導致的土骨架變形(εeij)以及冰相變導致的體積收縮(εTij)，εTij可用式(5)表示。

εTij=(1/3)εVT=(1/3)(0.09/1.09)(dθice/dt)，i=j。

(5)

式中：εVT為冰相變引起的體積應變，作為一種等向變形，僅對應變球張量作出貢獻；θice為體積含冰率，為該耦合模型的耦合節點。

通過廣義胡克(Hooke)定律，建立應力應變兩者之間的關系方程，即本構方程。

(6)

式中：Dijlm為本構張量。

依據質量守恒，孔隙內部的未凍水被排出，系統的質量發生了變化，則水力場的控制方程為：

-(k/yw)Δp=dεe/dt。

(7)

式中：k為滲透系數；Δ為二階偏微分算子。

依據能量守恒，高溫凍土單元溫度場的控制方程如下：

Liρi(dθice/dt)+λpfΔT=Cpfρf(dT/dt)+

Cpwρw∑[?(Tuj)/?xj]。

(8)

式中：Li為冰的相變潛熱；ρi為冰密度；λpf為導熱系數；Cpf為高溫凍土的恒壓比熱容；ρf為高溫凍土的密度；Cpw為水的恒壓比熱容；uj為高溫凍土形成的多孔介質內流體的流速。

1.1.2 耦合項

體積含冰率(θice)不僅對各物理場的因變量產生影響，而且成為了三個物理場共有的物理量，可由未凍水含量推導得出。

由于土顆粒表面能的存在，未凍水的含量與溫度是一個非線性的關系，未凍水與溫度的經驗公式為[17]：

wu=α·|T|-β。

(9)

式中：α、β為試驗系數；wu為未凍水含量。通過相關計算可得出wu與θice的關系：

θice=ρsρw(w0-wu)/[ρiρw+(ρi-ρw)wu+ρsρww0]。

(10)

式中：w0為初始含水率；聯立式(9)、(10)即可得出θice與溫度T的關系式，消除了控制方程的多余未知量。

1.2 模型參數計算方法

1.2.1 力學參數計算方法

依據損傷力學應變等價理論[17]，建立了高溫凍土應力-溫度耦合損傷本構模型，普通的楊氏模量被損傷模量所替代。

E′=E(1-Dc)。

(11)

式中：E為材料初始狀態下的楊氏模量；Dc為復合損傷因子，可通過式(12)計算[18]：

Dc=DM+DT-DM·DT。

(12)

式中：DM為應力損傷因子；DT為溫度損傷因子。

凍土受力劣化，應力損傷因子(DM)可通過某一時刻材料損傷單元數與初始總單元數的比值確定，因此，在某一應力水平下，損傷因子(DM)可表示為：

(13)

式中：Nt為損傷單元數；N為總單元數；η、F0是韋伯(Weibull)分布的參數。

依據莫爾-庫倫(Mohr-Coulomb)強度準則，定義三維軸對稱應力狀態下的應力損傷表征量(F)[19]：

F=Eε1[σ1(1-sinφ)-σ3(1+sinφ)]/(σ1-2vσ3)。

(14)

式(13)中的2個試驗參數(η、F0)，可通過直線擬合試驗參數求解。由三軸壓縮試驗應力應變關系可得：

ln{-ln[(σ1-2vσ3)/Eε1]}=ηlnF-ηlnF0。

(15)

由上式通過擬合應力應變曲線，可直接得出參數η、F0的值。

隨著溫度升高，凍土內部冰相強度降低，膠結能力減弱。該種溫度劣化效果，反映在凍土宏觀層面即為凍土剛度降低，由此引入溫度損傷因子(DT)進行分析。

DT=1-ET/ET0。

(16)

式中：ET為某一溫度條件下對應的初始彈性模量；ET0為設無損狀態下的初始彈性模量。彈性模量與溫度的關系，可由式(17)表示[20]：

ET=β+γ|T|0.6。

(17)

聯立式(11)、(12)、(13)、(16)，可得應力-溫度耦合損傷模型。

1.2.2 水力及溫度場參數計算方法

在水力場控制方程中，采用了表征土骨架透水能力的系數滲透率(k)。粉質黏土滲透率與孔隙比(e)的關系為：

k=CeD。

(18)

式中：C、D為材料參數[21]。

由于高溫凍土中存在的冰相會對水分遷移形成阻滯作用，因此，凍結狀態下高溫凍土滲透率(kf)表達式為[22]：

kf=k/I=k/1010θi。

(19)

由于本研究的控制方程是依據高溫凍土的多孔介質假設建立的，因此，使用復合材料的理論計算模型更為貼切，高溫凍土的恒壓比熱(Cpf)為[12]：

(20)

式中：Cf為穩定凍結區的比熱容；Cu為完全融化區的比熱容；Tp、Tb分別為相變區的上下限溫度(分別為1、-17 ℃)。

導熱系數(λpf)為：

(21)

式中：λf為穩定凍結區的導熱系數；λu為完全融化區導熱系數。

2 結果與分析

2.1 驗證試驗方法

驗證試驗選取土樣為哈爾濱市近郊粉質黏土(見圖1)，相關物理力學參數：最大干密度1.86 g·cm-3、最佳含水率15.50%、液限32.70%、塑限19.50%。

圖1 粒徑級配累積曲線

采用該低液限粉質黏土，壓實制備直徑為10 cm、高為20 cm、含水率為21%的圓柱體試件，干密度為1.75 g/cm3。為減少試件凍結時產生的水分遷移量，將其置于低溫試驗箱內-25 ℃快速凍結，并恒定48 h。試驗前24 h，置于恒溫箱內-2 ℃，使試件整體溫度分布均勻；隨后置于低溫三軸試驗儀低溫艙內，以-2 ℃恒溫2 h，使系統溫度恒定。

施加圍壓至0.25 MPa，設定目標軸向應力至1.57 kN；施加圍壓及軸壓過程中，記錄總軸向變形量。待達到目標軸壓，且圍壓穩定后，以0.1 ℃/min速度均勻升高圍壓室溫度；升溫過程中，打開試件底部排水通道，即可實現對三向應力狀態下高溫凍結粉質黏土的融沉固結模擬試驗。

2.2 模型參數計算

以-2 ℃為無損溫度狀態，通過0.05、0.15、0.25 MPa三組圍壓條件的低溫凍土三軸壓縮試驗，確定應力損傷因子，試驗曲線見圖2。

圖2 -2.0 ℃時三組圍壓條件對應的應力-應變曲線

依據力學參數計算方法對圖2數據進行處理，可得應力損傷因子相關數據，圍壓0.25 MPa、E為141.64 MPa、F0為2.151 6、η為0.443 6、φ為14.66°。通過這些參數，即可得到圍壓為0.25 MPa時，應力-溫度耦合損傷模型應力損傷因子。

在0.25 MPa圍壓時，-0.5、-1.0、-2.0 ℃三組溫度條件對應的應力-應變曲線對模型溫度損傷因子進行求解，試驗曲線見圖3。依據圖3數據，擬合求解式(17)，最終結果見圖4。融化狀態下，則采用壓縮試驗所得壓縮模量換算為楊氏模量予以確定，壓縮模量(Es)為2.39 MPa。

圖3 圍壓0.25 MPa時三組溫度條件對應的應力-應變曲線

本研究粉質黏土融化階段，存在一定的滯后性，其相變區上限為1 ℃[23]。因此，將式(9)進行變形，最終所得融化過程中未凍水含量與溫度變化關系曲線見圖5。

依據水力及溫度場參數計算方法，確定水力及溫度場參數。式(18)中C、D兩參數，與材料的塑性指數、液性指數有關，表達式如下：

C=e-5.51-4ln(Ip)。

(22)

D=7.52e-0.25IL。

(23)式中：IP為塑性指數；IL為液性指數[21]。經計算，獲得水力場、溫度場相關參數(見表1)。

圖4 圍壓0.25 MPa時彈性模量與溫度關系曲線

圖5 未凍水含量與溫度關系曲線

表1 水力場與溫度場相關參數

2.3 模型有效性檢驗

將計算獲得的各項物理場參數帶入水熱力耦合(THM)模型控制方程進行計算，對比試驗結果，驗證模型有效性(見圖6)。

由圖6可見：A段內，高溫凍土變形經歷了緩慢變形及變形加速兩個階段，考慮幾何非線性的THM模型計算結果，與試驗曲線吻合度較高。B段為轉變段，該處變形速率逐漸減小，模型計算結果與試驗曲線出現少量偏離。C段為穩定段，此時模型計算結果與試驗曲線在B段產生的偏離并未進一步擴大，兩條變形曲線近乎平行變化。

A段為變形段；B段為轉變段；C段為穩定段。

同時對比考慮幾何非線性的大變形THM模型與不考慮幾何非線性的THM模型，由圖6可見：A段兩者計算結果近乎沒有差異；在B段，變形大于12 mm時，兩者的差異出現。大變形模型的變形速率衰減速度明顯小于小應變模型；但在C段，B段差異同樣沒有繼續擴大，兩者變形速率近乎相等。

融化固結為一個與時間相關的物理過程，為了分析變形速率的變化，定義試驗固結度(U(t))，U(t)=uz(t)/uz(t1)，式中uz為一點豎直方向的位移、t1為最終試驗結束時間，依此計算公式對沉降數據進行處理(見圖7)。

圖7 試驗固結度試驗曲線與模擬曲線對比

由圖7可見：無論是考慮幾何非線性的大變形THM模型，還是小應變假設下的THM模型，其固結進程實質上與試驗曲線差別并不大。模型計算結果與試驗結果固結曲線主要差別在于變形急速增加段后半部分，原因在于模型計算時，融土彈性模量采用了壓縮模量換算而來的楊氏模量，其與凍結狀態下的彈性模量存在不連續變化，造成了固結變形曲線的非平滑過渡。

2.4 融化固結過程分析

沿軸線方向取一垂直于上頂面與下底面的截面，選取0.9 h對比該截面上的溫度及孔隙水壓力(見圖8)。該時刻試件溫度處于高溫凍土相變區域，可較好反映融化狀態凍土內部物理場之間的相關關系。由圖8可見：由于邊緣先于中間部分融化，融化狀態下的土體排水能力優于凍結狀態，邊緣土體形成排水通道，孔隙水壓力先于中間部分消散；同時，孔隙水流動過程中，高溫流體經過低溫土骨架會交換部分熱量，加速中間部分土體融化。

進一步對比應力場與溫度場間的相關關系，同樣選取0.9 h時的溫度及范式等效(Von Mises)應力(見圖9)。由圖9可見：中部凍結部分土體所受應力高于邊緣融化部分；由于凍土剛度高于融土，而上部施加軸向荷載的加壓模塊，迫使土體軸向發生整體變形，屬于均勻變形，因此剛度較大的凍結部分承擔了更多的應力，即為混合狀態的土柱提供了更多剛度。

試驗時間0.9 h；圖面為孔隙水壓力等值線圖(云圖)；圖面數據為溫度，單位為℃。

試驗時間0.9 h；圖面為范式等效應力等值線圖(云圖)；圖面數據為溫度，單位為℃。

最后對應力場及水力場間的相互關系進行分析，對比0.9 h時孔隙水壓力和范式等效應力間的相關關系(見圖10)。由圖10可見：中部凍結部分應力最大，但是該處孔隙水壓力并非最大；由于排水面位于凍土土柱下部，而位于頂端的區域孔隙水排出通道最長，因此孔隙水壓力最大。依據有效應力原理，該處孔隙水壓力承擔了部分外荷載所造成的應力，則土骨架所受應力分布在該區域與孔隙水壓力分布具有顯著的相關性，隨著孔隙水壓力由邊緣向中部逐漸增大，應力水平逐漸降低。

試驗時間0.9 h；圖面為范式等效應力等值線圖(云圖)；圖面數據為孔隙水壓力，單位為kPa。

3 結論

該模型通過引入高溫凍土應力-溫度耦合損傷本構模型，其凍結狀態初始融化階段預測結果與試驗曲線較為吻合。在變形加速至緩慢變形過渡階段(1.5～2.5 h)，模擬結果與試驗曲線存在一定差異，但整體可有效預測高溫凍結粉質黏土的融化固結過程。

考慮大應變狀態的水熱力耦合模型對融化變形的預測精度，優于依據小應變假設的模型。

高溫凍土外部升溫狀態下，邊緣融化部分形成排水通道，加速上部土體水分排出。