基于隨機介質理論的隧道凍結期地表凍脹預測方法

蔡海兵 ，彭立敏，鄭騰龍

(1. 安徽理工大學 土木建筑學院，安徽 淮南，232001；2. 中南大學 土木工程學院，湖南 長沙，410075)

20 世紀50 年代，Litwiniszyn[1]基于砂箱模型試驗，為研究地下采礦過程中巖層與地表移動等問題而提出了隨機介質理論，該理論將巖層視為一種“隨機介質”，將地下開采所引起的巖層與地表變形看成是一隨機過程，應用概率統(tǒng)計的方法可實現對巖層與地表變形的預估。隨機介質理論經過Liu 等[2-3]的發(fā)展和完善，其應用領域從最初的地下采礦地表變形分析，發(fā)展到露天采礦、近地表開挖及地層疏水所引起的地表變形分析等。陽軍生等[4]和楊小禮等[5]均針對城市隧道施工所引起的地表移動和變形問題，采用隨機介質理論進行分析討論，進一步拓展了隨機介質理論的應用領域。施成華等[6]考慮隧道開挖地層移動與變形的時間－空間發(fā)展過程，基于隨機介質理論，系統(tǒng)地對城市隧道礦山法及盾構法等施工引起的地層移動與變形相關問題進行了研究。目前該理論已成為我國城市地鐵隧道工程中地表變形預測的實用方法之一。如將地層凍結引起的地表凍脹變形看成是一隨機過程，該過程可以認為是巖土體開挖引起地表變形的逆過程，則隧道凍結施工引起的地表凍脹變形問題亦可以采用隨機介質理論進行預測分析。陽軍生等[4]、李方政[7]、陶德敬等[8]和周太全等[9]都采用該理論對地鐵隧道凍結法施工引起的地表凍脹變形進行了分析，并得出了一些有意義的結論，但共同的不足之處為在分析過程中均未考慮凍結過程，即凍結壁的形成過程。本文作者考慮凍結壁的形成過程，基于隨機介質理論，以圓形隧道全斷面水平凍結工程為例，以建立隧道凍結期地表凍脹的預測方法。

1 地表凍脹的隨機介質預測方法

1.1 凍結壁交圈前的預測方法

圖1 所示為凍結壁形成過程圖。開始凍結時，冷媒劑與凍結管周圍的地層產生劇烈的熱交換，在每個凍結管周圍很快形成近似圓柱形的凍結巖土柱，如圖1(a)所示；凍土柱進一步擴展，直至相鄰凍土柱相互交圈，在隧道周圍形成封閉的凍結壁，如圖1(b)所示；交圈后，原各凍結管的凍結鋒面連成向隧道內擴展的內凍結鋒面和向隧道外擴展的外凍結鋒面，且內、外凍結鋒面很快趨于平滑，如圖1(c)所示。

圖2 所示為凍結壁交圈前坐標系轉換圖。假定隧道采用單圈全斷面水平凍結，n 根凍結管均勻布置在同一圈徑上，各凍結管所處的土層性質相同，由凍結壁的形成過程可知，凍結壁交圈前，經t 時刻后，在隧道周邊形成n 個外半徑為R(t)的凍土柱，各凍土柱形成過程中發(fā)生的體積膨脹均為半徑從R(t)均勻膨脹至RΔ(t)，膨脹的各圓環(huán)柱截面積也均為Δ(t)，如圖2所示。由隨機介質理論可知[4]：對于平面問題，在εoη(xoz)直角坐標系下，單個凍土柱在t 時刻發(fā)生膨脹變形而引起的地表各點凍脹位移表達式為

圖1 凍結壁形成過程圖Fig.1 Formation processes of frozen wall

圖2 凍結壁交圈前坐標系轉換圖Fig.2 Coordinate system transformation diagram before circle-crossing period of frozen wall

式中：β 為巖土層的主要影響角。

記各凍結管的編號為1，2，…，n，第i 根凍結管中心在εoη(xoz)直角坐標系下的橫坐標為xi，縱坐標為zi，將εoη(xoz)直角坐標系變換到rioiθi極坐標系，如圖2 所示，則轉換公式為

式中：h 為隧道中心距地表的距離；Rd為凍結管布置圈半徑。

由式(1)可得：在平面問題條件下，第i 根凍結管水平凍結引起地表各點凍脹位移表達式為

如忽略相鄰凍結管凍結引起地表凍脹位移的相互影響，則在凍結壁交圈前，地表凍脹由多個凍土柱的疊加膨脹變形引起。由此，根據疊加原理，則有凍結壁交圈前的地表凍脹位移表達式為

1.2 凍結壁交圈后的預測方法

圖3 所示為凍結壁交圈后坐標系轉換圖。隨著凍結時間的延長，隧道周邊各凍土柱進一步擴展，相鄰凍土柱相互交圈后即在隧道周邊形成圓環(huán)柱狀凍結壁，如圖3 所示，凍結壁外鋒面半徑為R(t)，凍結壁擴展過程中的體積膨脹假定為從半徑R(t)均勻膨脹至RΔ(t)，凍脹區(qū)域面積為Δ(t)。凍結壁交圈后，地表凍脹則由整個凍結壁膨脹變形引起。

圖3 凍結壁交圈后坐標系轉換圖Fig.3 Coordinate system transformation after circle-crossing period of frozen wall

將εoη(xoz)直角坐標系轉換成ro′θ 極坐標系，則該兩坐標系之間的轉換公式為

由式(1)可得在平面問題條件下，凍結壁交圈后的地表凍脹位移表達式為

當隧道采用全斷面水平凍結法施工時，考慮凍結壁的形成過程，凍結壁交圈前采用式(5)預測地表各點的凍脹位移，凍結壁交圈后則采用式(8)預測地表各點的凍脹位移。

2 凍結外鋒面半徑的確定方法

在凍結壁交圈前，凍結外鋒面半徑R(t)為單個凍土柱的半徑，在凍結壁交圈后，則為凍結壁外鋒面半徑。凍結壁交圈前、后凍結外鋒面半徑的確定問題即為單管凍結和多管凍結條件下溫度場分布規(guī)律的求解問題。國內外眾多科研學者對此問題進行了相關研究，總的來說，現有凍結溫度場的計算理論主要包括平板凍結理論和單管凍結理論[10-11]，其中平板凍結理論為一維半無限域熱傳導問題，而單管凍結理論為二維熱傳導問題。凍結壁交圈后的瞬態(tài)溫度場可由平板凍結理論近似求解，凍結壁交圈前的瞬態(tài)溫度場則可由單管凍結理論近似求解。

2.1 平板凍結理論

圖4 平板凍結理論示意圖Fig.4 Schematic diagram of flat-panel freezing theory

圖4 所示為平板凍結理論示意圖。將薄板垂直置入土中，假定薄板僅在其側面與土體發(fā)生熱交換，其頂、底面均與土體不發(fā)生熱交換，薄板冷卻面溫度為θc且保持恒定，凍結區(qū)土體溫度為θf，未凍區(qū)土體溫度為θu，均為時間t 和坐標x 的函數。凍結鋒面為一移動的相變邊界，其與冷卻面的距離為X(t)，且在凍結鋒面處，土體凍結溫度θd=0 ℃。

對于該一維半無限域熱傳導問題，凍結區(qū)和未凍區(qū)微分方程可寫為

式中：αf和αu分別為凍結區(qū)、未凍區(qū)土體的熱擴散系數，m2/d；kf和ku分別為凍土和未凍土的導熱系數，kJ/(m·d·℃)；cf和cu分別為凍土和未凍土的比熱容，kJ/(kg·℃)；ρf和ρu分別為凍土和未凍土的飽和密度，kg/m3。

微分方程的初始條件為

式中：θ0為土體初始溫度。

微分方程的邊界條件為：

除上述邊界條件外，在凍結鋒面，即相變邊界面X(t)處的熱平衡方程為

式中：L 為單位容積土體的相變潛熱，kJ/m3，由下式給出：

式中：Lw為單位質量水的相變潛熱；ρd為土體的干密度，kg/m3；w0為土體的初始含水量；wu為凍土中未凍水含量。

結合上述初始條件和邊界條件，對微分方程進行求解，則可得凍結鋒面X(t)與凍結時間的平方根成正比，即：

式中：A 為待定常數，需滿足下式

式中：Φ(y)為高斯誤差函數，且有：

實際的凍結工程一般都依靠多根凍結管來形成凍結壁，凍結壁交圈后的溫度場可由平板凍結理論近似求解。當未凍土和凍土的導熱系數、比熱容、密度等熱物理參數由室內試驗獲取后，凍結鋒面位置可采用式(18)和(19)計算。若已知凍結管布置圈半徑Rd，則凍結外鋒面半徑R(t)、內鋒面半徑R1(t)可按下式計算：

2.2 單管凍結理論

凍結壁交圈前的溫度場則由單管凍結理論進行求解，如圖5 所示，根據凍結鋒面將溫度場分成凍結區(qū)和未凍結區(qū)，凍結區(qū)土體溫度為θf，未凍區(qū)土體溫度為θu，均為時間t 和徑向坐標r 的函數，凍結鋒面半徑為R(t)，且凍結鋒面處土體凍結溫度θd=0 ℃。凍結管外半徑為Rp，凍結管熱流密度為Q(kJ/(m·d))。

圖5 單管凍結理論示意圖Fig.5 Schematic diagram of single pipe freezing theory

對于該二維熱傳導問題，凍結區(qū)和未凍區(qū)微分方程可寫為

微分方程的初始條件為

微分方程的邊界條件為

在相變邊界面R(t)處的熱平衡方程為

這種模式不僅能保障購房者的基本住房需求，更能有效化解資金流動性風險，具體表現在兩個方面。一是維護公積金制度設立初衷，有效保障中低收入群體的購房貸款需求。這種模式對中低收入群體是一種保護，尤其對效益不好的單位的職工而言，如果不實行該模式，他們會因單位繳存基數下降而無法享受公積金購房政策。二是可以盤活資金使用效率，減少資金使用風險。在這種模式下，通過分析職工的購房面積和購房實際需求，來調整公積金額度，既能增加市場上個人貸款資金的供應總量，也提高了公積金貸款對區(qū)域居民的覆蓋率。

結合上述初始條件和邊界條件，對微分方程進行求解，則可得凍結鋒面半徑R(t)與凍結時間的平方根也成正比，即：

式中：B 為待定常數，需滿足下式

式中：E(y)為指數積分函數，且有

凍結壁交圈前的溫度場可由單管凍結理論近似求解，由單管凍結理論可知，當未凍土和凍土的熱物理參數確定后，單管凍結條件下凍結外鋒面半徑可采用式(30)和(31)進行計算，凍結內鋒面半徑即為凍結管外半徑Rp。

相鄰凍土柱的凍結鋒面開始重疊，即意味著凍結壁剛好交圈，則凍結壁交圈時間可由下式計算：

式中：tj為凍結壁交圈時間；l 為相鄰凍結管間距。

3 凍脹區(qū)域外半徑的確定方法

工程中反映土體凍脹強弱的指標一般采用凍脹率，它是指土體試樣在無側向變形無縱向荷載條件下，經單向凍結，其縱向的高度增量與試樣原高度的比

式中：εf為土體凍脹率；Δhf為試樣縱向凍脹量；h 為試樣原高度。

據其原理，假定凍結壁向外均勻膨脹，在凍結壁交圈前，凍脹區(qū)域外半徑由下式計算：

在凍結壁交圈后，凍脹區(qū)域外半徑由下式計算：

4 工程案例計算及分析

廣州地鐵3 號線天河客運站折返線隧道斷面為馬蹄形，長138.8 m，寬11.4 m，高9.146 m，埋深8～10 m，屬于淺埋大斷面隧道。由于該隧道斜下穿廣汕公路和沙河立交橋，不能采用明挖法施工，最終采用選擇全斷面水平凍結暗挖法施工。

隧道水平凍結長度為138.8 m，分南、北兩段施工，南、北段凍結長度均為71.9 m，末端搭接凍結長度為5 m。南、北端均在距隧道開挖邊界1 m 處布置一圈凍結孔，凍結孔孔數為46 個，凍結孔間距為700～950 mm，南、北端凍結孔錯位布置。凍結壁有效厚度設計值為2.5 m，凍結壁設計平均溫度為-8 ℃，積極凍結期約為150 d[13]。隧道北端凍結孔布置設計如圖6所示。

圖6 隧道北端凍結孔布置圖Fig.6 Freezing hole design at northern end of tunnel

隧道斷面及凍結管布置圈均為五心圓拱線形，為計算簡便起見，采用等代圓方法將其簡化成圓形，經簡化后的隧道外半徑為5.137 m，凍結管圓形布置圈半徑為6.137 m，相鄰管間距為0.838 m。

根據折返線隧道巖土層分布情況，結合隧道及凍結壁的所處位置，在現場針對具有代表性的土層鉆取土樣，然后進行人工凍土物理力學性能試驗，鑒于對地表凍脹位移起決定作用的凍結管主要位于可塑砂質黏土層中[13-14]，因此，本次預測分析中所需輸入的熱物理力學參數按該土層進行取值，如表1 所示。

表1 土體熱物理力學參數Table 1 Thermophysical parameters of soil

由于該隧道相鄰凍結管平均間距僅為0.838 m，凍結壁交圈時間小于10 d，因此僅采用凍結壁交圈后的凍脹預測方法進行分析。

在分析中，考慮到隧道斷面大，地表凍脹主要影響范圍取為45 m，隧道軸線埋深為14.728 m，則土層主要影響角度的正切值tan β 為0.327。

該淺埋大斷面地鐵隧道實際凍結施工中，共布置了28 個地表變形測點，其中15 和10-1 號測點均位于隧道中心軸線處[15]，16，17 號測點與15 號測點位于同一斷面上，與隧道中心水平距離分別約為8 和15 m。

采用自行編制的地表凍脹Maple 計算程序，地表計算區(qū)域選為0≤x≤40 m，得到隧道積極凍結150 d后，地表凍脹位移分布規(guī)律的理論分析與現場實測結果對比如圖7 所示。

由圖7 可知：隧道積極凍結150 d 后，地表凍脹位移分布規(guī)律的理論分析結果較現場實測結果(15 號、16 號和17 號測點)稍大一些，但可大致反應實際情況，如對比隧道中心軸線處地表凍脹位移(距隧道中心水平距離為0 m)，現場實測結果(15 號測點)為317 mm，理論分析結果為350 mm，二者相差33 mm，與凍脹位移相比，該差尚可忽略。地表凍脹位移隨機介質理論分析中，對隧道斷面及巖土分層均在一定程度上進行了簡化，特別是熱物理參數主要參考可塑砂質黏土層進行取值，相對于其他土層，該土層的凍脹率最大，這也是地表凍脹位移理論分析結果較現場實測略大的主要原因。而同樣位于隧道中心軸線處地表的10-1 號測點，其實測的凍脹位移達到近400 mm，這是因為該測點剛好布置在折返線隧道中部處地表，且在隧道中部，南、北端凍結管有約5 m 長的搭接范圍，從而造成 10-1 號測點的凍脹位移在所有測點中最大[12,14]。

圖7 凍結150 d 后地表凍脹位移分布Fig.7 Surface frost heave displacements distribution after freezing 150 d

積極凍結期內，隧道中心軸線處地表歷時凍脹位移的理論分析與現場實測結果對比如圖8 所示。

圖8 表明，在地表凍脹位移隨凍結時間的變化方面，理論分析與現場實測結果(15 號測點) 吻合，總體表現為，在凍結初期，凍脹位移增長速度較快，而到凍結后期，凍脹位移隨時間呈相對緩慢增長趨勢。10-1號測點因其所處的特殊位置，其凍脹速度最快。

地表凍脹位移的對比分析結果驗證了本文所提出的隨機介質預測方法的可靠性，且該方法用于工程實際均具有較高的精確性。

圖8 隧道中心軸線處地表歷時凍脹位移Fig.8 Surface frost heave displacements with time at central axis of tunnel

5 結論

1) 在隧道水平凍結期，凍結壁交圈前，地表凍脹位移由多個獨立凍土柱的疊加凍脹效應引起，凍結壁交圈后，地表凍脹位移則由整個凍結壁的凍脹效應引起。據此，考慮凍結壁的形成過程，基于隨機介質理論，推導了地表凍脹位移的計算表達式，建立了地表凍脹的預測方法。

2) 對凍結外鋒面半徑和凍脹區(qū)域外半徑這2 個關鍵參數的取值方法進行了探討。提出凍結壁交圈前的瞬態(tài)溫度場則可由單管凍結理論近似求解，凍結壁交圈后的瞬態(tài)溫度場可由平板凍結理論近似求解，得出多管凍結交圈條件下凍結壁的擴展厚度(單管凍結條件下的凍土柱半徑)與凍結時間的平方根成正比，并基于平板凍結理論和單管凍結理論，確定了凍結外鋒面半徑的計算方法，并在此基礎上，導出了根據土體凍脹率、凍結內、外鋒面半徑計算凍脹區(qū)域外半徑的解析方法。

3) 將所提出的地表凍脹隨機介質預測方法應用于實際工程案例中，通過與現場實測結果的對比分析，驗證了該預測方法的可靠性。

[1] Litwiniszyn J. Fundamental principles of the mechanics of stochastic medium[C]//Proceedings of the 3rd Conference on Theoretical Applied Mechanics. Bangalore, India, 1957: 18-26.

[2] LIU Baochen. Ground surface movement due to underground excavation in China[C]//Comprehensive Rock Engineering. New York: Pergaman Press, 1993: 781-817.

[3] 劉寶琛, 張家生. 近地表開挖引起的地表沉降的隨機介質方法[J]. 巖石力學與工程學報, 1995, 14(4): 289-296.LIU Baochen, ZHANG Jiasheng. Stochastic medium calculating method of surface settlement caused by near surface excavation[J]. Chinese Journal of Rock Mechanics and Engineering, 1995, 14(4): 289-296.

[4] 陽軍生, 劉寶琛. 城市隧道施工引起的地表移動及變形[M].北京: 中國鐵道出版社, 2002: 10-19.YANG Junsheng, LIU Baochen. Ground movement and deformation induced by urban tunnel construction[M]. Beijing:China Railway Publishing House, 2002: 10-19.

[5] YANG Xiaoli, WANG Jinming. Ground movement prediction for tunnels using simplified procedure[J]. Tunnelling and Underground Space Technology, 2011, 26(3): 462-471.

[6] 施成華, 彭立敏, 雷明鋒, 等. 淺埋隧道施工地層變形時空統(tǒng)一預測理論與應用[M]. 北京: 科學出版社, 2010: 18-30.SHI Chenghua, PENG Limin, LEI Mingfeng, et al. Time space united prediction theory and application of stratum deformation caused by shallow tunnel excavation[M]. Beijing: Science Press,2010: 18-30.

[7] 李方政. 土體凍脹與地基梁相互作用的疊加法研究[J]. 巖土力學, 2009, 30(1): 79-85.LI Fangzheng. Research on superposition method of interaction between soil frost heave and beams on foundation[J]. Rock and Soil Mechanics, 2009, 30(1): 79-85.

[8] 陶德敬, 王明年, 劉大剛. 凍結法隧道施工引起的地表移動及變形預測[J]. 現代隧道技術, 2006, 43(6): 45-50.TAO Dejing, WANG Mingnian, LIU Dagang. Prediction of surface movement and deformation caused by freezing method construction of tunnels[J]. Modern Tunnelling Technology, 2006,43(6): 45-50.

[9] 周太全, 華淵, 呂寶華. 某地鐵區(qū)間隧道凍結施工凍脹效應隨機預測分析[J]. 冰川凍土, 2009, 31(1): 119-123.ZHOU Taiquan, HUA Yuan, Lü Baohua. Stochastic prediction for frost heave effect of a subway tunnel[J]. Journal of Glaciology and Geocryology, 2009, 31(1): 119-123.

[10] 蔣斌松, 王金鴿, 周國慶. 單管凍結溫度場解析計算[J]. 中國礦業(yè)大學學報, 2009, 38(4): 463-466.JIANG Binsong, WANG Jinge, ZHOU Guoqing. Analytical calculation of temperature field around a single freezing pipe[J].Journal of China University of Mining & Technology, 2009,38(4): 463-466.

[11] Lunardini V J, Varotta R. Approximate solution to Neumann problem for soil systems[J]. Journal of Energy Resources Technology, 1981, 103(1): 76-82.

[12] 劉鴻緒. 對土凍結過程中若干凍脹力學問題的商榷[J]. 冰川凍土, 1990, 12(3): 269-280.LIU Hongxu. Discussion of frost heave mechanics problems in soil freezing process[J]. Journal of Glaciology and Geocryology,1990, 12(3): 269-280.

[13] 蔡海兵, 彭立敏, 鄭騰龍. 隧道水平凍結施工引起地表凍脹的歷時預測模型[J]. 巖土力學, 2012, 33(6): 1761-1768.CAI Haibing, PENG Limin, ZHENG Tenglong. A duration prediction model of surface frost heave induced by tunnelling with horizontal freezing method[J]. Rock and Soil Mechanics,2012, 33(6): 1761-1768.

[14] 姚直書. 蔡海兵, 程樺, 等. 采用長距離水平凍結暗挖法的淺埋大斷面地鐵隧道施工技術[J]. 中國鐵道科學, 2011, 32(1):75-80.YAO Zhishu, CAI Haibing, CHENG Hua, et al. Multiparameter monitoring of frozen soil structure with super- long freezing hole drilling in Guangzhou metro[J]. China Railway Science, 2011, 32(1): 75-80.

[15] 姜耀東, 趙毅鑫, 周罡, 等. 廣州地鐵超長水平凍結多參量監(jiān)測分析[J]. 巖土力學, 2010, 31(1): 158-164, 173.JIANG Yaodong, ZHAO Yixin, ZHOU Gang, et al.Multi-parameter monitoring of frozen soil structure with super-long freezing hole drilling in Guangzhou metro[J]. Rock and Soil Mechanics, 2010, 31(1): 158-164, 173.