海底隧道管幕凍結法水熱耦合溫度場數值模擬分析

關鍵詞：滲流作用； 管幕凍結法； 有限元模擬； 水熱耦合； 凍結方案優化

中圖分類號：TU472 文獻標識碼：A DOI：10. 7525/j. issn. 1006-8023. 2024. 06. 022

0引言

隧道是廣泛應用于地層內的工程建筑物，是高效利用地下空間的一類形式。但在修建穿過飽和富水軟土地層的城市隧道時，由于城市建筑物密集，需在隧道施工過程中不影響周遭建筑的同時預防地下水滲漏的風險，如何滿足城市隧道施工的高要求是需要解決的難題［1-3］。管幕凍結法是在管幕法的基礎上結合人工地層凍結技術的一種新型地下暗挖工法。與常規管幕法不同的是，管幕凍結法利用制冷技術將土體凍結形成凍土帷幕，用以代替鋼管止水，不僅有效地提高結構的防水性，而且減小在復雜地層施工止水的難度，同時還保留管幕的大承載優點，能有效控制土體的變形，不破壞施工區域周遭的敏感建筑，該施工方法適用于修建具有淺埋大斷面特點且下穿富水飽和軟土地層的城市隧道［4-5］。

拱北隧道口岸段為港珠澳大橋珠海連接線工程最為特殊的地段，是世界最大暗挖斷面之一，該地段土層具有飽和水、高壓縮和高滲透性等特點，為保障施工人員安全、施工要求需達到超高的防滲水和擾動標準，同時不影響地上建筑和地下管線，為滿足這高難度的施工條件，拱北隧道在施工中采用管幕支撐與凍結封水結合的支護體系，是管幕凍結法在國內隧道暗挖中的首次成功應用［6-8］。過去十余年來，國內學者對管幕凍結法溫度場的研究已有相關成果，龍偉等［9］結合港珠澳大橋拱北隧道段實況，采用數值模擬研究異形凍結管開啟前后的管幕凍結法溫度場變化規律，其模擬結論滿足現場驗證。盧亦焱等［10］以港珠澳大橋珠海連接線拱北隧道管幕凍結工程為實際依托，使用有限元軟件分析其凍結管布置方案的溫度場分布和凍結帷幕厚度在鹽水凍結期間的變化。胡俊等［11］結合港珠澳大橋珠海連接線拱北隧道工程管幕凍結工法，提出將凍結管布置于鋼管之間，由鋼管與凍土咬合形成新的支護體系，同時采用有限元軟件模擬了該支護體系下的溫度場發展規律。以上的學者對管幕凍結法溫度場的研究成果大部分未考慮滲流條件，而在土體凍結期間，地下水的流動是影響凍結效果的主要因素之一，由于土中流體具有傳熱特性，土體在凍結期間會受滲流作用的影響產生熱量交換，凍土帷幕的形成將會受到地下水滲流的干擾［12-13］。對此，國內學者對凍土帷幕溫度場的水熱耦合機制也進行了相關研究，楊平等［14］將多孔介質熱傳導理論與達西定律理論結合，建立凍土帷幕在地下水條件下的數學模型，采用數值計算得出粉砂在凍結期間的溫度場與滲流場發展變化規律；白天麒等［15］在滲流場和溫度場耦合理論基礎上，建立數值模型研究圓形改橢圓形凍結管的凍土帷幕變化和溫度場分布規律；潘旭東等［16］結合太原地鐵某聯絡通道施工工況，建立數值模型研究多排凍結管的凍結效果，并改變模型水頭差分析該模型滲流場對溫度場的影響。

距今為止，國內工程采用管幕凍結法的實例仍較少，且地下水流動是凍土帷幕形成過程中的危害因素，管幕凍土共承載結構與常規的管幕結構在破壞形式上存在差異。三亞河口通道工程項目采用管幕凍結法施工，是繼大盾構地下施工之外最新的、技術含量高的地下工程項目，對于海南自由貿易島地下空間施工創新意義重大。為研究滲流作用對管幕凍土結構的不利影響，結合三亞河口通道海底隧道管幕凍結工程，運用COMSOL有限元軟件分析管幕凍結法溫度場在不同滲流流速下的變化規律，并設置觀測路徑進一步分析溫度場在上下游的分布差異；同時基于以上結果，對凍結方案進行優化分析，為之后類似管幕凍結工程的設計施工提供參考根據。

1凍結方案

1. 1工程概況

三亞河口通道海底隧道工程橫跨三亞河口，周邊環境敏感，地理位置特殊，河口段隧道采用上下疊層的單洞雙層隧道，采用管幕凍結法施工。建設隧道總長度3 118 m，地質環境復雜，工程穿過粉質黏土以及強透水性細砂等多種對施工不利的土層。

1. 2凍結管、鋼管布置

管幕段隧道環向布置28根鋼管，其中，8組咬合鋼管，4根獨立鋼管，鋼管切割處通過凍結支護，管幕內外圈設計有效凍土帷幕厚度大于等于4. 0m。凍結管總數為176根，分內外圈布置，外圈共布置96根凍結管，內圈布置80根凍結管。凍結區域按功能分為強凍結區與弱凍結區，強凍結區為鋼管間切割、連接焊板處內外側的區域，其凍土帷幕滿足承載和封水要求，弱凍結區為咬合管內外側的區域，其凍土帷幕僅滿足組合管在咬合焊接時的封水要求，強凍結區凍結管間距800mm，弱凍結區凍結管間距950 mm，強弱凍結區劃分與凍結管、鋼管布置如圖1所示。

3建立數值模型

3. 1計算基本假定

數值模型為水熱耦合機制，考慮多孔介質在滲流、凍結等條件下的作用機理復雜，參考相關依據［18］，基于以下的基本假定建立數值模型。

1）假定地層初始溫度為18 ℃，土體為多孔飽和狀態介質，土層均質分布且各項同性。

2）假定溫度荷載直接施加在凍結管外壁上，鹽水冷媒與土體換熱時無熱量損失。

3）假定凍結區域土體中滲流屬于層流，適用于達西定律，水流單向均勻。

4）假定凍結過程僅存在溫度場與滲流場耦合作用，忽略應力場影響。

5）假定土體溫度在低于-1 ℃時，土體開始產生凍結，-1 ℃等溫線包絡的面積為最大凍結區域，而-10 ℃等溫線包絡的面積為最小凍結區域。

3. 2模型參數選取

1）在基本假定的基礎上建立二維數值模型，實際工況管幕斷面尺寸約為21 m（橫）×21 m（縱），本研究需觀察熱量云圖以及考慮模型邊界誤差，將模型幾何尺寸設為30 m（橫）×30 m（縱），網格劃分與模型邊界條件如圖2所示。

2）本研究選取工程最不利土層并參考相關凍土試驗研究［19-20］，得到土體的導熱系數、滲透系數等，物理模型選用的參數值見表1。

在凍結期間，土體在低溫的作用下會變為凍土，土體的導熱系數、滲透系數、比熱、密度隨著溫度的改變而變化，未凍結溫度區間為-1～30 ℃，凍結溫度區間為-30～-1 ℃。將土體凍結狀態下視為完全不透水，凍土的滲透系數取至1. 91×10^-30m/s。

3）模型采用瞬態分析，凍結管管壁為熱荷載邊界，模型四周為絕熱邊界，凍結計劃時間步共40 d，每個時間步24 h，鹽水降溫計劃通過內插函數實現，在凍結5 d冷媒溫度降至-15 ℃，10～40 d冷媒溫度均為-28 ℃，見表2。

4模擬結果分析

將模型上下游邊界水頭差分別設置為ΔH=0、10、20、25、35 m，通過式（3）和式（4）計算得到相應的平均達西滲流流速，通過分析模型無滲流、低滲流流速、高滲流流速3種情形下的-1 ℃與-10 ℃等溫線在凍結期間的變化得到其凍土帷幕溫度場在滲流作用下的發展規律。

4. 1無滲流情形

當模型上下游邊界水頭差為0 m時，滲流速度V=0 m/d，此時模型無滲流流動。由圖3可知，凍結5 d時，-1 ℃與-10 ℃等溫線以凍結管為中心向四周擴散，呈同心圓狀；凍結10d時，位于外圈弱凍結區的-1 ℃與-10 ℃等溫線與其他凍結區域相比交圈速度較慢，與強凍結區還未相連，弱凍結區的凍結管間距較大，冷量較為分散，凍結效果較弱；凍結進行至21 d時，-1 ℃與-10℃等溫線均完成交圈，整體凍土帷幕已封閉；凍結40d時，隨著凍結的進行，凍土帷幕均勻地向四周進一步擴大。

4. 2低滲流流速情形

將模型上下游邊界水頭差設為ΔH=10 m，模型平均滲流流速為V=2. 87 m/d，小于人工凍結設計要求的5 m/d范圍，視土體此時處于低滲流流速作用的情形，土體-1 ℃與-10 ℃等溫線在低滲流流速作用下產生相應的變化。

由圖4可知，圖4中箭頭指向表示滲流流向，凍結僅5 d，-1 ℃與-10 ℃等溫線已出現朝滲流的流向發展的趨勢，說明滲流作用在凍結早期就開始干擾凍土帷幕的形成；凍結10 d時，與無滲流情形相比，滲流作用下的上游外圈-1 ℃與-10 ℃等溫線的交圈速度慢于下游，下游處的弱凍結區等溫線已與強凍結區相連，說明在凍結的過程中，上游側凍結管釋放的冷量在滲流的作用下被傳遞至下游側，抑制上游側凍土帷幕的形成，促進下游側凍土帷幕的發展；凍結至22 d時，對比上下游-1 ℃與-10 ℃等溫線的交圈情況，下游的凍結區域比上游凍結范圍更大，凍土帷幕已明顯向下游一側偏移，上下游兩側的溫度場分布不均勻；凍結40d時，隨著凍結的進行，凍土帷幕區域的不均勻性進一步增大。

4. 3高滲流流速情形

為進一步分析滲流流速的變化對凍土帷幕發展的影響，將模型上下游邊界水頭差分別增大至ΔH=20、25m，滲流流速分別為5. 75、7. 17m/d，視土體處于高滲流流速作用的情形。

由圖5和圖6可知，在凍結5 d時，等溫線以凍結管為中心向滲流流向發展，與圖4情況類似，此時滲流速度的變化對凍土帷幕影響較小；凍結10d時，與低流速情形相比，高流速作用下整體的-1 ℃與-10 ℃等溫線交圈情況較差，滲流仍可穿過凍結區域，說明當滲流速度過大時，滲流作用加快了冷量在土體中的傳遞，雖使上下游的等溫線交圈的差異減小，但整體凍土帷幕的形成都受到了抑制，此時滲流速度的變化對凍土帷幕的形成造成顯著影響。

由圖3—圖6可知，ΔH=0、10、20、25 m4種水頭差下的-1 ℃與-10 ℃等溫線交圈所需的時間分別為21、22、26、31 d，隨著滲流流速的增大，凍土帷幕交圈所需的時間由21 d延長至31 d，在小滲流流速作用下的凍土帷幕交圈時間點僅延長1d，在最大滲流流速作用下的凍土帷幕交圈時間延長了10 d，距凍結完成僅剩9 d。低滲流流速作用對凍土帷幕完整交圈的時間影響較小，當滲流流速增大至5 m/d以上，凍土帷幕完整交圈時間點在高滲流流速作用下明顯延長。

對比3種情形的整體凍土帷幕在凍結40 d時的差異，在最大滲流流速作用下的凍土帷幕向下游一側偏移的程度最大，整體凍土帷幕的不均勻性最大。

將模型上下游邊界水頭差進一步增大至ΔH=35 m，此時滲流流速V=10.02 m/d，由圖7可知，凍結5、10 d的凍土帷幕發展情況與圖4類似；凍結20d，上游-1 ℃與-10 ℃等溫線交圈速度遠落后于下游，上游鋼管間的凍土帷幕仍未封閉，滲流仍可穿過凍結區域；凍結40 d時，上游弱凍結區的外圈凍土帷幕未交圈，滲流作用對上游的弱凍結區外圈凍土帷幕的形成造成了顯著影響。當滲流速度過大時將會導致凍土帷幕在凍結期間無法成功交圈，需在施工中加強上游弱凍結區的滲流流速的監測，必要時增加外圈弱凍結區的凍結管數量。

由于強凍結區土體凍結范圍大于弱凍結區，需分別討論強凍結區與弱凍結區凍土帷幕厚度，通過對比不同水頭差下溫度場云圖的變化，計算凍土帷幕外緣-10℃等溫線至凍結管外壁的距離，得出強弱凍結區在不同水頭差條件下的凍土帷幕厚度，如圖8所示，溫度場云圖如圖9所示，弱凍結區凍土帷幕厚度由4. 21 m減少至4. 03 m，強凍結區凍土帷幕厚度由4. 95m 減少至4.69m，隨著滲流速度的增大，強弱凍結區的凍土帷幕厚度均在減小，滲流將凍結管釋放的冷量帶至凍結區域外，雖擴大了凍結的范圍，但減弱了凍結效果，整體凍土帷幕厚度在減小。

5觀測點分析

5. 1溫度觀測點

土體的凍結溫度是判斷凍土帷幕是否形成的重要依據，由于地下水的流動會干擾凍結的正常進行，使凍土帷幕發展不均勻，為進一步分析溫度場分布在上下游的差異，在模型上選取2條路徑，滲流流速選取V=0 m/d和V=7. 17 m/d，路徑1與路徑2長度均為4 m，每個觀測點間隔1 m，且觀測點在上下游呈對稱分布，如圖10所示。

由圖11 可知，由于5 號點與4 號點與凍結管間距較小，降溫速率最大，凍結僅7 d，溫度就已降至-4. 90、-4. 96 ℃ ，凍結至40 d，溫度約降至-27. 32、-27. 33 ℃；2號和7號2點的降溫速率僅次于4號點和5號點，4個點的降溫曲線基本重合；1號與8號點位于凍結區域兩側，距凍結管最遠，在凍結22 d時，溫度約為-1. 01、-1. 48 ℃，凍結40 d，溫度約為-7. 84 、-8. 11 ℃，最終溫度較高，凍結效果較差；3 號點與6 號點位于鋼管之間，均在凍結17 d時溫度降至-1 ℃以下，凍結40 d，溫度分別降至-24. 92、-24. 93 ℃。

無滲流作用下的凍土帷幕上下游均勻發展，路徑1與路徑2的觀測點降溫曲線走向基本一致，距離凍結管越近的點開始凍結的時間越早，降溫速率越大，最終溫度越低。

路徑1（1—4號點）位于上游側，路徑2（5—8號點）位于下游側，由圖12降溫曲線可知，1號點在凍結進行至16 d時，溫度仍為原始地溫18 ℃，凍結40 d時，該點溫度仍高于0 ℃，約為1. 52 ℃，滲流作用對該點的凍結造成了顯著的影響，體現了滲流對上游側凍土帷幕發展的抑制作用，8號點在凍結17 d時溫度已降至約-1. 38 ℃ ，凍結40 d 時，溫度約為-14. 50 ℃，1號點未形成凍土帷幕，8號點形成穩定的凍土帷幕，2點的凍結效果有明顯差異，最終溫度相差較大；2號點于凍結18 d降至-1 ℃以下，7號點凍結16 d溫度約為-1. 18 ℃，表示凍結開始進行，凍結40 d 時，2 號點與7 號點溫度約為-26. 85、-27. 32 ℃，7號點開始凍結的時間點提前了2 d，且最終溫度略低于2號點；3號點和6號點、4號點和5號點的降溫曲線對比情況與2號點和7號點類似，6號點開始凍結的時間點相較于3號點提前了3 d，5號點相較于4號點提前1d，凍結40d時，6號點最終溫度低于3號點，約為-24. 36 ℃，5號點溫度低于4號點，約為-27. 33 ℃。上游側凍結管釋放的冷量隨滲流的流動被帶至下游側，增強了下游側土體的降溫效果，使得下游側土體開始凍結的時間點早于上游側，且最終溫度低于上游側。

滲流作用下的上下游路徑點的降溫曲線走向有較大區別，體現了上下游凍土帷幕的不均勻發展，位于凍結區域兩側的溫度測點在滲流作用下變化最大。

5. 2滲流流速觀測點

由以上分析可知，上游側外圈弱凍結區凍結管的交圈受滲流影響較大，在外圈弱凍結區中選取一凍結管，在其四周布置4個流速測點，如圖13所示，模型平均滲流流速設置為2. 87m/d，通過分析凍結期間滲流流速的變化得到相關溫度場發展規律。

由圖14知，測點1位于上游側，流速降低最慢，凍結27 d，滲流速度才降低至0，此時測點1剛與封閉的凍土帷幕相連，滲流已無法通過該點；測點2位于下游側，流速降低最快，在凍結的12 d內，滲流速度較為穩定地降低至0，與測點1相比，形成封閉凍土帷幕的時間點提前了14 d，兩點流速變化的對比體現了凍土帷幕隨滲流流向發展的特點。測點3與測點4的滲流速度分別在凍結18 d和凍結13 d減少至0，測點4形成封閉凍土帷幕的時間點早于測點3。除測點2外，其他3個測點在凍結一段時間后滲流流速驟降至0，且4個測點的流速在不同的時間點降至0，體現了滲流作用下凍土帷幕的不規則發展。

6凍結方案優化分析

由以上分析可知，當模型平均滲流速度約為10 m/d時，凍土帷幕才出現局部不交圈的情況，原凍結方案設計偏于保守，為提升工程經濟效益，在保證凍土帷幕厚度滿足設計要求的前提下，設置模型平均滲流速度為2. 87 m/d，對原凍結方案的凍結管布置進行調整：外圈96根凍結管布置不變，且不改變內圈弱凍結區凍結管布置，僅將內圈強凍結區靠近鋼管切割處的兩排凍結管改為一排，其數量由4根減少至2根，凍結管間距仍為800 mm，內圈凍結管總數由80根減少至56根。

由圖15可知，當凍結7 d時，凍結區域上下兩側的-1 ℃等溫線已相連，而左右兩側-1 ℃等溫線仍呈同心圓狀，且下游一側交圈速度快于上游一側，此時凍土帷幕還未形成，滲流從鋼管間直接穿過；當凍結12 d時，-1 ℃等溫線已在凍結區域內呈片狀相連，說明凍結區域大部分的土體已低于-1 ℃，處在穩定凍結階段，且部分區域內-10 ℃等溫線已交圈，局部穩定凍土帷幕已形成，開始發揮止水作用，滲流僅能穿過-10 ℃等溫線還未交圈的區域；當凍結17 d時，完整凍土帷幕輪廓已出現，下游一側的內外圈-10 ℃等溫線已交圈成功，上游一側的外圈弱凍結區的-10 ℃等溫線還未交圈，上游側的外圈弱凍結區受滲流影響較大；當凍結40 d時，完整凍土帷幕已形成，凍土帷幕朝下游一側發展。

兩凍結方案在凍土帷幕厚度與交圈時間上存在差異，由圖4和圖15可知，原方案在凍結10d時，鋼管間的凍土帷幕已經封閉，滲流已無法從鋼管間穿過；而優化方案在凍結12d時，鋼管間的凍土帷幕仍未封閉，在凍結27d時，-1℃和-10℃等溫線交圈成功，凍土帷幕完整封閉的時間點相較于原方案延長了5 d。計算凍土帷幕外緣-10℃等溫線至凍結管外壁的距離，由圖9和圖16可知，原方案弱凍結區與強凍結區凍土帷幕厚度分別約為4.14 、4.81m，平均厚度約為4.475m，優化后弱凍結區與強凍結區凍土帷幕厚度分別約為4.06、4. 44m，平均厚度約為4.28m，優化方案弱凍結區厚度減少了0.08m，強凍結區減少了0.37m，平均厚度減少了0.195m，優化方案適量減少了內圈凍結管的數量，雖使凍土帷幕完整交圈的時間延長，但最終凍土帷幕厚度減少較小，仍滿足凍結設計要求。

7結束語

本研究運用COMSOL有限元軟件分析管幕凍土溫度場在滲流作用下的發展規律，并基于分析的結果對原凍結方案進行優化，主要得出以下結論，

1）無滲流作用時，凍土帷幕在上下游兩側呈對稱分布，當滲流流速約為2. 87 m/d時，上下游的凍土帷幕交圈情況有明顯差異，但對整體凍土帷幕交圈時間影響較小，當滲流流速大于5 m/d時，上下游的凍土帷幕交圈差異減小，高滲流速度作用抑制了整體凍土帷幕的發展。

2）滲流流速越大，凍土帷幕向下游一側偏移的程度越大，完整交圈所需的時間越長，且厚度越小。

3）相較于上游側土體，下游側土體開始凍結的時間點更早，滲流流速降低速率更大，且在凍結期結束時溫度更低。

4） 原凍結方案布置下的凍土帷幕在滲流速度約為10m/d時才出現局部不交圈情況，凍結止水效果良好，在實際的施工中安全性強、可行性高，但凍結管數量使用過多，會提高施工的經濟成本，因此，適量減少了內圈凍結管的數量，雖會略微延長凍土帷幕交圈的時間，但整體凍土帷幕厚度相較于原方案相差較少。