寒區引水隧洞低溫相變溫度場-滲流場耦合模型及其數值模擬研究

姜海波,貊祖國,后雄斌,金 瑾

(石河子大學水利建筑工程學院,新疆石河子832000)

0 引 言

寒區隧洞圍巖凍脹問題涉及到巖體的應力場、溫度場、滲流場等多場耦合作用。國內外眾多學者都致力于低溫巖體“溫度場—滲流場”耦合問題的研究[1-5]?；谶B續介質或等效連續介質模型的裂隙巖體滲流場與溫度場耦合作用是目前該領域研究的熱點。巖體凍脹產生的根本原因是液態水在低溫下相變成固態冰,相變過程中水的體積膨脹引起的凍脹力造成了巖體的凍脹破壞。大量工程研究表明：巖體中的裂隙水低溫凍結膨脹會造成巖體的凍脹破壞,而凍結過程中水分不斷向凍結鋒面遷移、聚集并結晶成冰透鏡體也是巖體產生凍脹的重要原因[6- 8]。TAKEDA等[9-12]對水熱遷移機理及其相關問題進行了比較系統的分析和研究,還有學者提出了多個水、熱耦合計算模型[13-15],并對模型的有效性進行了驗證。然而,寒區巖體的凍脹是一個復雜的熱力耦合和水冰相變過程,巖體的凍脹特性及凍融劣化機理涉及低溫凍融環境下復雜的溫度場(T)、滲流場(H)和應力場(M)的耦合問題。裂隙巖體中的水分遷移機制、凍脹力的計算及其萌生消散機制、裂隙凍融開裂擴展機制和巖體多次凍融強度損傷及穩定性評價等都是研究凍巖過程中亟待解決的關鍵問題。在T—H—M三場耦合作用下巖體的凍脹特性和凍融劣化規律方面,前人取得了一系列研究成果[1,16-17],但總的來說,國內外關于低溫凍融循環作用下巖體的多場(溫度場、應力場、滲流場及化學場等)耦合研究較多,并取得了豐富的研究成果。而關于低溫相變條件下巖體溫度場—滲流場耦合問題的研究并不系統,并沒有考慮凍脹力的消散和低溫凍融循環對圍巖力學特性及其襯砌結構穩定性的影響。

鑒于此,本文基于連續介質力學、熱力學、滲流力學,考慮低溫相變對巖體溫度場和滲流場的影響,結合低溫凍結條件下圍巖的溫度場—滲流場耦合模型,采用有限元法對新疆某電站引水隧洞工程進行圍巖的凍結過程分析,研究低溫相變條件下隧洞圍巖溫度場和滲流特性的變化規律,對引水隧洞經受長時間持續凍結作用后的圍巖溫度場、滲流特征進行分析,定量分析長期低溫凍結作用對隧洞圍巖溫度場及其滲透特性的影響,為引水隧洞的安全運行提供理論依據。

1 低溫巖體溫度場和滲流場耦合數學模型

低溫凍結巖體是巖體骨架、孔隙水、冰以及空氣組成的多相介質。在凍結的過程中,巖體中的空氣對巖體相變的影響比較小,可以不考慮空氣影響,因此可以假定巖體處于飽和狀態。巖體的凍脹破壞是低溫情況下水冰相變引起的,巖體溫度的變化引起了水冰相變和水分遷移的發生。在低溫相變過程中,熱量的遷移受到了溫度梯度和巖體孔隙對流換熱的影響。低溫條件下巖體的水冰相變過程受到水力梯度和溫度梯度的作用,它們共同造成了水分不斷地向凍結鋒面遷移、聚集,最終加速了巖體的凍脹破壞。

根據以上分析,以及孔隙介質的對流換熱原理,低溫相變巖體介質熱傳導的微分控制方程為

(1)

低溫巖體的滲透方程為

(2)

式中,S為比貯水系數,取常數；p為滲透壓力；k為滲透率；η為水的粘滯系數,計算中取0.001 kg/(m·s)；Hg為重力水頭；SP0為分凝勢,在計算中溫度低于冰點時為正常數,在高于水的冰點溫度時取0；T為溫度；QH滲流場的源或者匯。求解上述微分控制方程,在移動邊界上必須滿足連續性條件和能量守恒條件[18]。

對于低溫巖體溫度場與滲流場耦合問題,求解時需要知道的邊界條件主要有兩種,一種為Dirichlet邊界條件(也稱固定邊界)[19]；另一種為Neumann邊界條件(也稱流量邊界)[20]。在有限元計算時,采用這兩類邊界條件進行計算。

(1) 溫度場邊界條件和初始條件。Dirichlet溫度邊界為T=Ta,該邊界條件可通過現場監測獲得計算邊界的溫度值,一般認為在足夠遠處該溫度值為一恒定的常數。Neumann滲流邊界為-K該邊界的取值,隨溫度梯度的變化而變化,在計算中,需要根據現場監測資料進行計算。

利用給定的邊界條件,便可求解低溫巖體溫度場和滲流場。由于巖體在低溫下其物理力學性質隨溫度變化而變化,且凍結緣帶的位置也是不斷變化的。因此,該問題為強非線性問題,不可能用解析方法求得其精確解,只有通過數值方法獲得近似解。本文采用有限元方法對該問題進行數值分析,以獲得較為精確的近似解。

2 寒區輸水隧洞溫度場與滲流場耦合有限元分析

2.1 工程概況

新疆某電站引水隧洞工程,地質情況較為復雜。一年凍結期達200 d,冰凍層平均厚度3～5 m,最厚達8 m。最冷月(1月)平均氣溫為-9.6 ℃,最熱月(7月)平均氣溫為9.0 ℃,極端最低氣溫-37.6 ℃?？紤]到該地區典型的大溫差環境和長時間的凍結環境,施工及其隧洞運行的影響,原有的巖體溫度場和滲流條件遭到破壞,可能出現的凍脹問題會影響隧洞的安全運行,因此需要對隧洞圍巖的溫度場和滲流場進行重點研究。結合當地歷史氣溫及其現場監測資料,選用-30 ℃～+20 ℃的溫度循環進行寒凍風化作用的疊加模擬,對引水隧洞經受長時間持續的凍結作用后圍巖的溫度場、滲流特征進行分析,從而定量分析長期低溫凍結作用對隧洞圍巖溫度場及其滲透特性的影響,為引水隧洞的安全運行提供理論依據。

2.2 有限元計算模型、材料參數及其邊界初始條件

該引水隧洞長2.5 km,為圓形隧洞,直徑5 m,平均埋深7.5 m。由于隧洞較長,可以近似視為平面問題進行分析。根據引水隧洞模型的對稱性,本文選用模型的1/2作為有限元計算區域,如圖1所示?？紤]到有限元計算的收斂性,當溫度在-30～+20 ℃之間變化時,圍巖的滲透系數隨溫度連續性變化而變化。考慮到低溫情況下圍巖孔隙水的相變,孔隙水的熱物理學參數通過現場試驗可獲得,具體如表1所示。

材料密度/kg·m-3比熱/J·(kg·K)-1導熱系數/W·(m·K)-1相變潛熱/J·kg-1孔隙率/%圍巖36509001.8—0.37水100041800.56333000—冰917.521002.23——

熱力學邊界條件：AB、CD為對稱邊界,在該邊界上滿足-K；AF為固定邊界,TAF=-15 ℃；DE為滲流邊界,在該邊界上滿足-K在計算中,認為巖體的初始溫度場均勻分布,即滿足T|t=0=3.5 ℃。

隧洞圍巖的計算參數如表1所示。計算模型取時間步長為1 d,計算總時長取200 d。

2.3 計算結果分析

通過計算得到了不同時期隧洞圍巖溫度場的分布,如圖2所示,圖中“淺色”區域表示凍結區,“深色”區域代表為未凍結區。從圖2可以看出,在給定的邊界及其初始條件下,圍巖在最初的10 d內,凍結速率非?？?隧洞斷面凍深達到了21.25 m；在10～100 d,凍結深度呈現穩定增長的趨勢,表現為凍結區范圍擴大,未凍結區范圍縮小；當計算時間為200 d時,溫度場的分布趨于穩定。隨著凍結時間持續,從計算結果中可以看出,凍結寬度隨時間的推移不斷增加。

分析圖2可知,在最初的時刻,隧洞襯砌結構內的溫度為-16 ℃,隨著時間的持續,隧洞襯砌結構內的溫度變化不大,在-16～-12 ℃之間變化。在計算區域內,受到巖體導熱系數、相變潛熱等因素的影響,隨著水平深度的增大,凍結溫度逐漸回升。凍結區與未凍結區的界面為凍結鋒面,隨著凍結時間的持續,凍結鋒面向圍巖內部推移。

含水巖體在凍結過程中,由于分凝勢的作用,水分由未凍結區向凍結區遷移,如圖3所示。根據有限元計算的結果,由于凍結溫度的作用,孔隙中的原位水由未凍結區向凍結區遷移,滲流速度在凍結緣處達到最大,為1.359×10-4m/s,在未凍結區滲流速度較小,為4.1×10-6m/s。

3 結 論

(1)根據連續介質力學、熱力學、滲流力學等理論,結合低溫凍結條件下巖體的溫度—滲流耦合模型,并考慮水冰相變和長期凍結作用對巖體溫度場和滲流場的影響,分析了低溫相變條件下隧洞圍巖溫度場和滲流特性的變化規律,定量分析了長期低溫凍結作用對隧洞圍巖溫度場及其滲透特性的影響,模型可靠,計算結果具有一定的合理性。

(2)在長期的凍結過程中,圍巖內部形成了凍結區與未凍結區；隨著凍結時間的持續,凍結鋒面向圍巖內部推移。

(3)隨著凍結時間持續,凍結寬度和深度隨時間的推移不斷增加,這將進一步加劇圍巖的凍結范圍和凍融損傷程度。需要說明的是,圍巖凍脹問題涉及到巖體溫度場、滲流場、應力場以及凍融損傷多場耦合作用,考慮低溫相變的影響,圍巖溫度場和滲流特性對圍巖凍脹損傷的影響還需更進一步的研究。

圖2 隧洞溫度場分布(單位：℃)

圖3 200 d時滲流速度分布及其矢量方向

[1] NEAUPANE K M, YAMABE T, YOSHINAKA R. Simulation of a fully coupled thermo-hydro-mechanical system in freezing and thawing rock[J]. International Journal of Rock Mechanics and Mining Sciences, 1999, 36(5): 563- 580．

[2] 徐光苗, 劉泉聲, 張秀麗. 凍結溫度下巖體THM完全耦合的理論初步分析[J]. 巖石力學與工程學報, 2004, 23(21): 3709- 3713．

[3] KANG Y, LIU Q, HUANG S. A fully coupled thermo-hydro-mechanical model for rock mass under freezing/thawing condition[J]. Cold Regions Science and Technology, 2013, 95(11): 19- 26．

[4] 張國新, 沙莎. 混凝土壩全過程多場耦合仿真分析[J]. 水利水電技術, 2015, 46(6): 87- 93．

[5] AKAGAWA S, FUKUDA M. Frost heave mechanism in welded tuff[J]. Permafrost and Periglacial Processes, 1991, 2(4): 301-309．

[6] SETZER M J. Development of the micro-ice-lens model[C]∥International RILEM Workshop on Frost Resistance of Concrete. Bagneux: RILEM Publications SARL, 2002: 133- 145．

[7] MURTON J B, PETERSON R, OZOUF J C. Bedrock fracture by ice segregation in cold regions[J]. Science, 2006(314): 1127- 1129．

[8] OZAWA H. Thermodynamics of frost heaving: a thermodynamic proposition for dynamic phenomena[J]. Physical Review E, 1997, 56(3): 2811- 2816．

[9] TAKEDA K, OKAMURA A. Microstructure of freezing front in freezing soils[C]∥Proceedings of the 1997 International Symposium on Ground Freezing and Frost Action in Soils. Rotterdam: A A Balkema, 1997: 171- 178．

[10] KONRAD J M,DUQUENNOI C A. A model for water transport and ice lensing in freezing soil[J]. Water Resources Research, 1993, 29(9): 3109- 3124.

[11] LAI Y M, WU Z, ZHU Y, et al. Nonlinear analysis for the coupled problem of temperature and seepage fields in cold region tunnels[J]. Cold Regions Science and Technology, 1999, 29(1): 89- 96．

[12] HARLAN R L. Analysis of coupled heat-fluid transport in partially frozen soil[J]. Water Resources Research, 1973, 9(5): 1314- 1323．

[13] SHEPPARD M I, KAY B D, LOCH J P G. Development and testing of a computer model for heat and mass flow in freezing soils[C]∥The 3rd International Conference on Permafrost. Ottawa, Ontario: National Research Council of Canada, 1978: 76- 81．

[14] TAYLOR G S, LUTHIN J N. A model for coupled heat and moisture transfer during soil freezing[J]. Canadian Geotechnical Journal, 1978, 15(4): 548- 555.

[15] FUKUDA M. Experimental studies of coupled heat and moisture transfer in soils during freezing[R]. Institute of Low Temperature Science, Hokkaido University, Sapporo, Japan: 1982: 35- 91．

[16] NEAUPANE K M, YAMABE T. A fully coupled thermo-hydro-mechanical nonlinear model for a frozen medium[J]. Computers and Geotechnics, 2001, 28(8): 613- 637．

[17] LI N, CHEN B, CHEN F X. The coupled heat-moisture-mechanic model of the frozen soil[J]. Cold Regions Science and Technology, 2000, 31(3): 199- 205．

[18] 賴遠明. 寒區隧道溫度場、滲流場和應力場耦合問題的非線性分析[D]. 蘭州： 中國科學院蘭州冰川凍土研究所, 1999.

[19] KONRAD J M, MORGENSTERN N R. The segregation potential of a freezing soil[J]. Canadian Geotechnical Journal, 1981, 18(4): 482- 491．

[20] DRNEVICH V P, KONRAD J M. Procedure for determining the segregation potential of freezing soils[J]. Geotechnical Testing Journal, 1987, 10(2): 51- 58．