軟黏土三維水-熱-力耦合凍融模型及應用

周 潔，李澤垚,，唐益群，田萬君

(1.同濟大學 土木工程學院，上海 200092；2.中國建筑第二工程局有限公司，北京 100074)

中國地鐵建設的過程中，人工凍結法是針對沿海地區飽和地層加固的一種有效工法。沿海地區人工凍土具有溫度變化速率快、凍結范圍較小、水文地質復雜的明顯特征[1-3]。人工凍土的凍結過程是一個復雜的不穩定導熱問題[4]，溫度場、滲流場和應力場相互耦合[5-7]。考慮滲透作用的人工凍土的研究有程樺等[8]基于Harlan水熱耦合模型，模擬了各向同性的飽和砂土在豎井凍結時凍結壁的擴展過程。攝宇[9]模擬了滲流作用下砂土單圈管凍結溫度場的凍結效果。崔灝[10]模擬了滲流作用下砂土和粉土凍結壁發育規律及對溫度場的影響。Vitel等[11]建立了一個完全熱力學一致的耦合熱-水-力數值模型，并模擬了滲流條件下飽和不可變形多孔介質的人工地層凍結。Ahmed等[12]建立了一種水-熱耦合有限元模型，可在滲流情況下通過尋找最佳位置來優化砂土地層凍結管的排布。目前，針對滲流影響的凍結效應研究集中在飽和無黏性土。相比之下，細粒的黏性土在溫度梯度作用下會產生可觀的水分遷移，從而形成分凝凍脹，且分凝凍脹比原位凍脹大很多[9]。

考慮水分遷移的分凝凍脹理論廣泛應用于天然正凍細粒土的模型計算[13-14]。與天然凍土不同，人工凍結區域范圍較小，溫度變化速率快，不能簡單抽象成一維的緩慢凍結，現有分凝凍脹模型的適用性較低。因此，在分凝凍脹理論的基礎上，通過對剛性冰模型的三維化，建立了考慮三維分凝凍脹的熱-水-力耦合作用模型，可以模擬計算沿海各向異性飽和地層使用人工凍結法應用時的溫度場、應力場、位移場。最后通過COMSOL有限元模擬上海楊樹浦路人工地層凍結法縮尺模型，并與縮尺試驗結果進行對比，驗證了模型的有效性。

1 考慮水分遷移的多場耦合模型建立

1.1 計算模型假定

軟黏土是一種細粒的黏性土，其顆粒表面具有一層吸附水，土顆粒通過吸附水傳遞土體應力[13]。吸附水層之外的水為重力水，產生孔隙水壓并驅動孔隙水流動。土體發生凍結時，首先重力水發生相變，隨著溫度的進一步降低，冰水相界面向吸附膜內侵入，分凝冰產生[14]。當冰相生長至完全隔絕分凝方向上土顆粒之間的聯系時，暖端一側的水分不能補給冷端一側，從而致使冷端一側的分凝冰停止生長[13]。由此，依據孔隙水的凍結狀態可將土體分為未凍結區、相變區和已凍結區。土體孔隙水的凍結狀態如圖1所示。

圖1 凍結土體凍結狀態Fig.1 Frozen state of frozen soil

人工凍結土多為飽和土，故假設:土體中的土骨架、孔隙水、孔隙冰皆為剛性且不可壓縮，土體壓縮僅與孔隙率的變化有關；水分遷移僅發生在未凍結區和相變區，已凍結區不發生水分遷移。

1.2 考慮三維分凝凍脹的水-熱-力耦合模型

飽和細粒黏性土的凍脹，可以看成由原位孔隙水相變膨脹所產生的原位凍脹與由水分遷移所產生的分凝凍脹的疊加。

1.2.1 原位凍脹模型

原位凍脹是由孔隙水相變時比容差所引起，即dV=Δνd[W/(1+W)]，其中，V為單位質量土體的體積，W為土體的未凍結水的含水率，Δν為H2O的比容差。

1.2.2 分凝凍脹模型

由于孔隙內冰、水處于相壓力平衡狀態，即冰、水兩相應符合Clapeyron方程[13]：

(1)

式中：pw為相界面水相壓強，L為相變潛熱。將方程積分并展開，取液態水在273.15 K、100 kPa時會發生的相變作為積分下限。

(2)

(3)

在凍結區，為保證任意一微分凝冰層顆粒處于靜力平衡狀態，且能完成應力傳遞(如圖2所示)，需要在分凝冰暖端的相界面水相壓力與外荷載之間引入一個平衡項uw，即

圖2 土體分凝狀態Fig.2 Soil segregation state

uw+pw=σl

(4)

其中σl為凍結區域土體的主應力。

Miller等[13]認為平衡項uw即為相變區致使未凍結水遷移的驅動壓力，從而有

duw+dpw=0

(5)

即相變區驅動孔隙水流動的水力梯度應為

(6)

冷端需要為土體中提供足夠的冷量才足以產生分凝冰。當該位置土體冷端吸熱能力使得凍結鋒面處水分遷移過來的孔隙水完全發生相變時，分凝冰發育迅速。土體體積應變速率應為水分遷移所引起的體積變化，即

(7)

針對某一個微單元體，倘若冷端沒有的冷量以支持孔隙水發生相變，該微單元體不產生分凝冰。倘若微單元體冷端傳遞過來的冷量不足以支持分凝所需，該微單元體就不能為凍結方向上后一個微單元體傳遞冷量，下一個微單元體就不產生分凝。這也表明凍結方向上的前一個微單元體有足夠的冷量支持孔隙水發生相變。在這個凍結方向上，這種微單元體有且僅有一個，不能形成有效凍結帶，故而可以忽略此類型的單元體。

溫度梯度所引起的水分遷移是土體產生凍脹的原因，水分遷移的方向是溫度梯度方向，垂直于溫度梯度方向可以近似地視為等溫線，不具備能使水分發生遷移的動力因素。但是凍結鋒面擴展前后都是一個連續的曲面。在溫度梯度方向上遷移動力因素的影響下，為保證凍結鋒面前后的連續性，垂直于溫度梯度方向上會獲得一個補償應變量，如圖3所示。

圖3 土體分凝凍脹量示意Fig.3 Segregation frost heave of soil

單元土體由于水分遷移所產生的總分凝凍脹量與凍結垂向和徑向分凝凍脹關系應滿足：

dε=dεH+dεHdεV1+dεHdεV2

(8)

式中dεV1、dεV2分別為凍結鋒面上徑向的應變量和緯向上的應變量。

對于單元土體，dεH相對較小，忽略去高階無窮量dεHdεV1+dεHdεV2的影響，得dε=dεH。則該位置所產生凍結方向上的分凝凍脹速率為

(9)

垂直于凍結方向上補償應變量應與凍結方向上的應變量和凍結鋒面曲率相關。該位置所產生垂直于凍結方向上的補償應變速率為

(10)

式中Fs1、Fs2為凍結鋒面上的徑向的曲率半徑與緯向的曲率半徑，可根據微分幾何學公式進行計算(各向凍脹計算見圖4)。

圖4 各向凍脹計算示意Fig.4 Schematic of calculation of frost heave in all directions

其中

1.2.3 軟黏土的融沉模型

土體融化時相變界面水相壓力差依然存在，但遷移的水分無法有效被消耗，即水分遷移路徑上沒有泄水出口，從而在根本上抑制了水相壓力差引起的水分遷移。因此，不考慮融化時水熱耦合關系引起的水分遷移。

土體融沉的過程由原位融化和融化后土體的固結沉降構成。原位融化與原位凍結相似，都是由土體孔隙冰相變所造成的體積變化引起，不同的是原位融化的水(冰)量不僅有凍結之前土體孔隙中的水分，還包括凍結過程中所遷移過來的水分量[14]。

排水固結的過程使用太沙基-倫杜利克三維固結微分方程：

(11)

式中：u為超孔隙水壓力，n為融化鋒面的法向，σ為融化時土體融化方向上的主應力，Cv為固結系數。

飽和土體融化相變時孔隙水承擔所有的應力傳遞，融化后的孔隙水壓力應與融化時土體凍脹力相同，計算時與土體一次性賦予孔隙水壓力初值。隨著融化的進一步進行，凍脹力消散，超孔隙水壓力產生。當融化鋒面在Δt時間內穿過一個長度增量時，在凍土層中排出的孔隙水體積為

(12)

式中：A為土體單位的橫截面面積，?u/?n為融化界面上的孔隙壓力梯度，γw為水的重度。

水流ΔV等于厚度為Δn土層體積的變化量，則體積應變為

(13)

1.3 模型參數的計算取值

1)土體含冰率：I=n0-W。

2)土體骨架比率：S=1-I-W。

3)未凍土孔隙率與土體應變的關系：n=n0-(εx+εy+εz)(1+n0)。其中,εx，εy，εz為土體各方向的應變；凍結區與相變區的孔隙率：n=n0-W。

6)土體密度：ρ=ρwW+ρiI+ρsS，其中,ρw、ρi、ρs分別為水、冰、土顆粒的密度；土體熱傳導系數：λ=λwW+λiI+λsS，其中,λw、λi、λs分別為水、冰、土顆粒的熱傳導系數；土體比熱：C=CwW+CiI+CsS，其中,Cw、Ci、Cs分別為水、冰、土顆粒的比熱。

7)當土體完全凍結后，由于分凝作用土體出現了垂直于凍結方向與平行于凍結方向的分凝冰帶，層狀的分凝冰帶將隨機出現在凍結體中，宏觀可以體現一定的均勻性。將3個凍結方向上分凝冰帶、原狀土比擬成相互穿插的彈簧模型，從而對分凝凍結土的壓縮模量進行預測(計算模型示意如圖5)：

圖5 計算模型示意Fig.5 Schematic of the calculation model

(14)

(15)

(16)

式中：Ec為土體在封閉系統內凍結完成后的壓縮模量，可以由實驗室測定得到；EH、EV1、EV2分別為凍結方向上的變形模量、垂直于凍結方向徑向上的變形模量、垂直于凍結方向緯向上的變形模量。凍結前土體、融化后土體的彈性模量可以由試驗測定得到。

1.4 本模型的優點

使用未凍土孔隙率、相變界面的壓強、相變溫度、未凍結水含水率等作為模型的基本參量，對分凝凍脹理論分凝冰模型進行三維化。相較廣泛應用的天然凍土分凝凍脹理論剛性冰模型，本模型更適用于凍結范圍小、溫度變化速率較快、土體各向異性明顯、水文條件復雜的人工凍結的應用。具體為

1)能為凍結鋒面曲率半徑小、凍結溫度變化迅速的工況更準確地預測(溫度變化、應力應變變化等)；

2)可以考慮凍結時具有強烈地下水滲流的工況，能夠對具有滲透各向異性地層不規則形狀的凍結進行模擬預測；

3)能夠預測各相異性凍結土的彈性模量，進而可以模擬凍結帷幕的變形特征、受力特征。

2 數值模型計算與結果分析

基于上述建立的三維分凝凍脹模型，通過COMSOL有限元模擬上海12號線楊樹浦路站人工地層凍結法縮尺模型，并與縮尺試驗結果進行對比，對模型的有效性進行驗證。

2.1 模擬工況

2.1.1 工程背景

上海市地鐵12號線楊樹浦路地鐵站是越江站(黃浦江)，隧道埋深15 m，地鐵聯絡通道建設在淤泥質軟黏土中。由于臨近黃浦江泄水通道，聯絡通道下伏著流速相對較高的微承壓水粉細砂層，平均流速為0.41 m/d。地鐵聯絡通道凍結管長為12 m，直徑為89 mm，凍結帷幕一排設置7根凍結管，凍結管中冷鹽水流量為3 m3/h，溫度為-30 ℃。對地鐵聯絡通道凍結工程有影響的土層為第④層淤泥質黏土、第⑤層粉細砂，其巖土、水文特性如表1、2所示。

表1 淤泥質黏土層特性Tab.1 Characteristics of silty clay layer

表2 粉細砂層特性Tab.2 Characteristics of silt fine sand layer

2.1.2 模擬試驗

縮尺模型試驗以現場凍結工況為基礎，采用考慮了材料、尺寸、滲流水頭、土體傳熱、泵送鹽水流量相似性的模型試驗系統[17]進行試驗。模型試驗系統相似比為1∶30，凍結管埋深50 cm，凍結管直徑8 mm，數量為7根，凍結鹽水泵流量經熱流密度相似計算需為30 L/min，凍結冷鹽水溫度為-30 ℃。試驗裝置如圖6所示。試驗箱中填置重塑淤泥質黏土與粉細砂以模擬現場施工地層，淤泥質黏土的厚度為46 cm，下層砂土的厚度為40 cm。箱體兩側設置水槽，試驗通過抽水泵為土體兩側施加水頭差以模擬臨江地下水的流動。設置滲流速度為0、0.25、0.50、0.75、1.2、3.0 m/d的6種試驗工況。通過低溫恒溫箱為凍結管進行制冷以模擬現場凍結。試驗首先為重塑土進行預固結，對模型土體進行制冷后實時監測土體內部各傳感器的監測指標變化規律。溫度傳感器在相同深度沿水平滲流方向以8 cm的間距(凍結管附近間距10 cm)串聯布置成一排測線，共布設了4排，分別位于凍結管、凍結帷幕上下邊緣及砂層處(縱向間距6 cm)。土壓力傳感器在凍結管、凍結帷幕上緣及砂層各位置的中心位置各布置1個，凍結帷幕下沿處由于滲流影響較大，以10 cm的間距沿滲流方向布置了3個。位移傳感器在土層上表面以8 cm的間距沿滲流方向均勻布置了5個。成排分布的土壓力傳感器、溫度傳感器使用細鐵絲進行固定，以使試驗過程中傳感器的位置不發生較大的偏移。傳感器排布位置如圖6(b)所示。溫度傳感器精度為0.01 ℃、土壓傳感器精度為0.2 kPa、位移傳感器精度為0.01 mm。傳感器精度較高且工作穩定性較好、尺寸效應不明顯、監測頻率滿足試驗設計要求。凍結帷幕擴展到規定厚度或當滲流速度特別大時，即使凍結足夠長時間，凍結傳感器監測的溫度仍然達不到-10 ℃(1 h內溫度下降不超過0.1 ℃)，關閉制冷系統，凍結結束。當溫度全部回升到17 ℃左右且地表沉降不再增加時，試驗結束。

圖6 縮尺試驗裝置Fig.6 Reduced-scale test device

數值模擬的模型工況與縮尺模型試驗工況完全相同，其等尺寸建立的模型如圖7所示。模擬凍結管溫度為-30 ℃，數值監測點位置、地下水流速、凍結厚度控制、凍結溫度控制等皆與縮尺模型試驗相同。試驗土體物理力學參數[18]如表3所示。

圖7 幾何模型示意Fig.7 Geometric model diagram

表3 模擬參數Tab.3 Simulation parameters

2.2 物理場控制方程

使用COMSOL有限元的PDE模塊對物理場方程進行編譯。滲流在未凍結區連續性方程為

(17)

滲流在相變區的連續性方程為

(18)

傳熱方程：在一個土體單元內，質量累積量的改變量、熱傳導量、熱對流量與潛熱釋放量的總和為0，即

(19)

土體應力場采用彈塑性本構模型。

2.3 結果分析

模擬結果表明，滲流對凍結效果存在臨界點，即當滲流小于1.2 m/d時，凍結帷幕才能凍結完成。

2.3.1 組合地層滲流對凍結溫度場的影響

凍結穩定時數值模擬區域等溫線分布如圖8所示?？梢钥闯觯瑳]有滲流作用時模型土體的等溫線相對水平；有滲流存在時，等溫線呈一定程度的不均勻化偏移，上游等溫線密集，下游等溫線稀疏。滲流對凍結帷幕的厚度整體都有削弱作用，對上游削弱最為明顯。

圖8 不同滲流速度下等溫線Fig.8 Isotherm diagram under different seepage velocities

2.3.2 組合地層滲流對軟黏土凍脹力的影響

土體凍結過程發生凍脹后，由于膨脹受限從而會產生向四面擴張的內應力。凍結工程中將土體凍脹受限所產生的內應力稱為凍脹力。

不同滲流作用下土體凍結穩定時，凍結管中心位置土體所產生凍脹力的切片云圖見圖9?？梢钥闯?，凍結管附近凍脹力集中效應明顯，凍結帷幕下沿下游凍脹力明顯大于上游，帷幕上沿變化不明顯。凍結帷幕上沿的凍脹力整體均大于凍結帷幕下沿。

圖9 不同滲流速度下凍脹力云圖Fig.9 Cloud diagram of frost heave force under different seepage velocities

針對能夠凍結完成的工況，凍結帷幕下上沿監測點A、B在不同滲流速度下凍脹力的增長如圖10、11所示。結果表明，監測點A、B凍脹力發展模式存在較大不同。監測點A凍結階段的凍脹力開始時增長迅速，后期增長緩慢。對于監測點B，凍脹力在整個凍結階段都在緩慢增長，凍脹力平穩不變的階段不明顯。在融沉階段，滲流速度大的工況凍脹力消散得更快。

圖10 監測點A凍脹力監測圖Fig.10 Monitoring results of frost heave force at point A

圖11 監測點B凍脹力監測圖Fig.11 Monitoring results of frost heave force at point B

產生該現象的原因是：靠近滲流邊界的地方，滲流抑制凍結方向上的溫度變化，從而抑制了孔隙冰的生成速率，致使凍結帷幕下沿凍脹力增長速率緩慢。同時，滲流減少了靠近滲流邊界地區凍結方向上分凝冰生成的質量，從而減小了凍脹力的最終值。

在實際工程施工中，凍脹力是監測預警的關鍵指標。對于凍結帷幕上沿部分，凍結前期就應重點關注凍脹力的突變增長，采取預警及措施；而凍結帷幕下沿更應該關注的是凍結后期，因為前期凍脹力增長較緩慢且數值較小。在融沉初始階段，凍結帷幕上沿凍脹力的快速消散應是工程注意的重點。同時，滲流速度增大時，土體產生的凍脹力較小，特別是靠近滲流邊界的凍結帷幕下沿。凍脹力的降低可以減小土體對隧道產生的壓力，增強工程的安全性。

2.3.3 組合地層滲流對地表沉降的影響

針對能夠凍結完成的工況，不同滲流速度地表最終沉降量如圖12所示。可以看出，沿著滲流方向，地表最終沉降越來越大，且滲流速度越大，最終沉降越低。因此，當存在下覆較大滲流砂層情況的軟黏土凍結，不均勻沉降也是需要關注的重要方面。

圖12 不同滲流速度地表位移Fig.12 Surface displacement at different seepage velocities

產生該現象的原因是：下游區域溫度較低，該區域生成了更多的孔隙冰從而最終沉降量更大。滲流可以抑制分凝冰的生成，從而使滲流大的工況分凝冰含量低，也就產生了滲流速度越大，凍脹、沉降量越低的現象。

3 模型驗證

選取監測點A在0.25 m/d滲流速度下的溫度、凍脹力及位移分別與模型試驗結果進行對比，結果如圖13所示。可以看出，溫度場與應力場的數值模擬結果與實測結果接近度較大。而位移的計算結果發展規律基本一致，凍脹位移存在一定誤差，融沉位移相對接近。主要原因為數值模型模擬的是完全考慮了三維水分遷移的理想開放系統，而模型試驗需要對滲流進行控制，并沒有做到在各個方向上給軟黏土同時補水,出現了試驗測定的凍脹量要小于數值模擬估計值的現象。

圖13 數值模擬與模型試驗結果對比Fig.13 Comparison of numerical simulation and model test

4 組合地層滲流對凍結法安全影響分級

以上海市地鐵12號線楊樹浦路地鐵站為工程背景，針對不同滲流速度的凍結工況進行數值模擬，對照現行隧道安全體系，分別對隧道位移、隧道附加凍脹應力、地表隆起、地表沉降進行預測并分級。各指標的管理等級可分為：等級Ⅰ為所有滲流工況下，物理指標處于極限值的2/3～極限值，是安全隱患較高的等級階段；等級Ⅱ為所有滲流工況下，物理指標處于極限值的1/3～2/3，是安全隱患發展的等級階段;等級Ⅲ為所有滲流工況下，物理指標小于極限值的1/3，是不利因素產生的初始階段，相對安全。分級結果如表4、5。

表4 凍結階段安全控制Tab.4 Safety control during freezing phase

表5 融沉階段安全控制Tab.5 Safety control during thawing phase

分級預警可以提前對工程地質不利因素的嚴重程度進行預先評估，從而指導工程及時應對工程地質不利因素，增強凍結施工的安全性，降低地質災害的影響程度。

5 結 論

1)在組合地層條件下，滲流會削減凍結穩定時凍結帷幕的厚度，且上游削弱最明顯；當速度大于1.2 m/d時，凍結工程將不能達到設計厚度。

2)在組合地層條件下，凍結工程靠近滲流邊界的一側凍脹力增長緩慢、最大值較低。

3)滲流會致使地表產生不均勻沉降，下游區域的沉降量更大。對于實際工程而言，地表的不均勻沉降也是組合地層滲流凍結法需要關注的危害。

4)通過建立的熱-水-力三場耦合模型對組合地層滲流作用下凍結法施工及周圍環境的影響進行參數分析，提出了結合組合地層滲流因素的隧道安全凍結法分級標準，可以為工程建設提供參考依據。