楊韜 郭穎 張程程 單煒

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

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

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

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

1 研究方法

1.1 耦合模型的建立

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

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

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

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

(1)

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

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

(2)

(3)

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

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

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

(4)

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

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

(5)

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

通過廣義胡克(Hooke)定律，建立應(yīng)力應(yīng)變兩者之間的關(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為導(dǎo)熱系數(shù)；Cpf為高溫凍土的恒壓比熱容；ρf為高溫凍土的密度；Cpw為水的恒壓比熱容；uj為高溫凍土形成的多孔介質(zhì)內(nèi)流體的流速。

1.1.2 耦合項(xiàng)

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

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

wu=α·|T|-β。

(9)

式中：α、β為試驗(yàn)系數(shù)；wu為未凍水含量。通過相關(guān)計(jì)算可得出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ù)計(jì)算方法

1.2.1 力學(xué)參數(shù)計(jì)算方法

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

E′=E(1-Dc)。

(11)

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

Dc=DM+DT-DM·DT。

(12)

式中：DM為應(yīng)力損傷因子；DT為溫度損傷因子。

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

(13)

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

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

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

(14)

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

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

(15)

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

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

DT=1-ET/ET0。

(16)

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

ET=β+γ|T|0.6。

(17)

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

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

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

k=CeD。

(18)

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

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

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

(19)

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

(20)

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

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

(21)

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

2 結(jié)果與分析

2.1 驗(yàn)證試驗(yàn)方法

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

圖1 粒徑級(jí)配累積曲線

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

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

2.2 模型參數(shù)計(jì)算

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

圖2 -2.0 ℃時(shí)三組圍壓條件對(duì)應(yīng)的應(yīng)力-應(yīng)變曲線

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

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

圖3 圍壓0.25 MPa時(shí)三組溫度條件對(duì)應(yīng)的應(yīng)力-應(yīng)變曲線

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

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

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

(22)

D=7.52e-0.25IL。

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

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

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

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

2.3 模型有效性檢驗(yàn)

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

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

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

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

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

圖7 試驗(yàn)固結(jié)度試驗(yàn)曲線與模擬曲線對(duì)比

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

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

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

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

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

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

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

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

3 結(jié)論

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

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

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