高溫凍土融化固結(jié)的水熱力耦合模型構(gòu)建1)

楊韜 郭穎 張程程 單煒

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

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

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

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

1 研究方法

1.1 耦合模型的建立

1.1.1 擴展的比奧大變形固結(jié)模型

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

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

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

(1)

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

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

(2)

(3)

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

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

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

(4)

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

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

(5)

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

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

(6)

式中：Dijlm為本構(gòu)張量。

依據(jù)質(zhì)量守恒，孔隙內(nèi)部的未凍水被排出，系統(tǒng)的質(zhì)量發(fā)生了變化，則水力場的控制方程為：

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

(7)

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

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

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

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

(8)

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

1.1.2 耦合項

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

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

wu=α·|T|-β。

(9)

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

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

(10)

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

1.2 模型參數(shù)計算方法

1.2.1 力學參數(shù)計算方法

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

E′=E(1-Dc)。

(11)

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

Dc=DM+DT-DM·DT。

(12)

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

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

(13)

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

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

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

(14)

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

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

(15)

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

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

DT=1-ET/ET0。

(16)

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

ET=β+γ|T|0.6。

(17)

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

1.2.2 水力及溫度場參數(shù)計算方法

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

k=CeD。

(18)

式中：C、D為材料參數(shù)[21]。

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

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

(19)

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

(20)

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

導熱系數(shù)(λpf)為：

(21)

式中：λf為穩(wěn)定凍結(jié)區(qū)的導熱系數(shù)；λu為完全融化區(qū)導熱系數(shù)。

2 結(jié)果與分析

2.1 驗證試驗方法

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

圖1 粒徑級配累積曲線

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

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

2.2 模型參數(shù)計算

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

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

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

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

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

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

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

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

(22)

D=7.52e-0.25IL。

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

圖4 圍壓0.25 MPa時彈性模量與溫度關(guān)系曲線

圖5 未凍水含量與溫度關(guān)系曲線

表1 水力場與溫度場相關(guān)參數(shù)

2.3 模型有效性檢驗

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

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

A段為變形段；B段為轉(zhuǎn)變段；C段為穩(wěn)定段。

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

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

圖7 試驗固結(jié)度試驗曲線與模擬曲線對比

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

2.4 融化固結(jié)過程分析

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

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

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

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

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

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

3 結(jié)論

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

考慮大應變狀態(tài)的水熱力耦合模型對融化變形的預測精度，優(yōu)于依據(jù)小應變假設的模型。

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