寒區引水隧洞圍巖凍脹特性分析
——以新疆布倫口水電站引水隧洞工程為例

彭小麗，王海娟，楊豐春，姜海波，李唱唱

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

0 引言

地球上凍土區面積約占陸地面積的50%，我國多年凍土和季節凍土面積分別占國土面積的21.5%和53.5%[1]。隨著“一帶一路”倡議深入實施，越來越多的隧洞需要建立在我國西部寒區，因此凍巖問題已經成為近年來的研究熱點。

寒區隧洞普遍受低溫、滲流共同作用而產生凍脹力。低溫凍融環境下，圍巖因溫度梯度進行熱傳導而使得內部圍巖溫度場發生改變，圍巖裂隙中的液態水狀態也發生改變。當溫度達到冰點時，孔隙水凍結，低溫凍結作用會產生9%的體積膨脹。楊更社等[2]對軟巖材料進行了在一定溫度梯度下的水熱試驗，結果表明溫度梯度是巖體水分遷移的主要動力，巖石孔隙率與凍結時間成正相關。譚賢君等[3]建立含相變的水熱耦合模型，研究了隧道溫度場分布規律和凍融圈大小。夏才初等[4]基于穩態假定，積分推導了考慮圍巖、襯砌未凍水含量的最大凍結深度解析計算式，并通過溫度場反演獲取解析式中的參數從而求得最大凍結深度半解析解。

在凍土和凍巖的凍脹特性研究方面，國內外許多學者采用室內試驗、數值模擬和現場監測等方法相結合，取得了豐碩的成果。J.M.Konrad[5]提出了凍土分凝勢理論，并根據分凝勢理論研究了凍結方式對凍結特性的影響。賴遠明等[6]利用黏彈性原理推導出凍脹力求解相關公式，并采用數值逆變換方法求出寒區隧洞凍脹力，通過算例計算表明是否考慮凍脹因素對應力值大小存在著較大的差異。張德華等[7]、Gao等[8]、黃繼輝等[9]基于彈塑性理論，給出了不同應力準則下圍巖凍脹力的彈性力學解析解。康永水等[10]運用應變片法對低溫環境下飽和巖樣和干燥巖樣進行測試，得到巖樣的低溫應變特征，研究其凍脹融縮效應，并給出巖石凍脹變形規律。渠孟飛等[11]、嚴健等[12]為研究凍脹壓力和凍融圈厚度變化，分別進行了室內隧道模型試驗和現場原位測試試驗。耿珂[13]、吳劍等[14]依靠數值模擬方法對隧道圍巖凍脹力和凍脹特性進行研究。

寒區隧洞圍巖因水-冰相變的影響，在環境溫度、孔隙滲流以及應力三者之間相互影響、相互作用下表現出“凍脹”現象，因凍脹作用而產生的隧洞凍害嚴重時會影響隧洞的施工和運營。大多數學者在凍脹數值模擬中僅考慮了熱力耦合的過程，而未將孔隙水相變過程考慮在內。鑒于此，本文以新疆某水電站引水隧洞為研究對象，使用大型有限元模擬軟件，考慮水-冰相變現象，建立低溫凍結環境下的圍巖瞬態溫度傳導模型和溫度-滲流-應力三場耦合模型，計算凍結深度和凍脹力特性，并將凍脹發生前后的凍脹特性進行對比，以此分析凍脹作用對該引水隧洞的影響，以期為引水隧洞的安全施工和運營提供一定指導。

1 凍脹力彈塑性力學解析解

關于圍巖凍脹模型，存在多種研究理論[15]。本文采用圍巖整體凍脹模型，認為受溫度影響圍巖沿徑向范圍內形成凍融圈，當溫度降低時，圍巖孔隙中的水結冰使得凍融圈整體膨脹，因此產生了凍脹力。

寒區隧洞凍脹力彈塑性力學解析解求解模型如圖1所示。其中，ph為圍巖凍脹后凍結層內壁作用在襯砌上的圍巖壓力;pf為凍脹后凍結層外壁作用在未凍圍巖上的圍巖壓力;hf為凍結深度;r0、rⅠ、rⅡ分別為襯砌內徑、襯砌外徑(或凍結層內徑)、凍結層外徑;rⅢ為任意一處未凍結圍巖距圓心的距離。為簡化計算并使解析解的值更接近實際值，需對模型作出如下假設： 1)假定襯砌、已凍結巖、未凍結巖均為彈性體； 2)襯砌與已凍結圍巖為綁定接觸，凍融圈與原始圍巖為彈性接觸； 3)各材料符合各向同性假設； 4)不考慮風化破碎層以及孔隙水滲流的耦合作用。

圖1 彈性力學求解模型

區域Ⅰ為圓形隧洞襯砌，將襯砌看作受外壓ph作用的圓筒結構，根據彈性力學[16]，該區域應力場拉密解可表示為

(1)

徑向位移

(2)

式中EⅠ、μⅠ分別為襯砌的彈性模量和泊松比。

當r=rⅠ時，襯砌外壁位移

(3)

區域Ⅱ為已凍結圍巖，已凍結區受ph、pf共同作用，該區域內應力拉密解為

(4)

其徑向位移

(5)

式中：rⅡ=rⅠ+hf；hf為凍結深度，m，計算時hf的取值通過對拱腰位置布置溫度測量元件進行圍巖內部徑向溫度場監測，確定為2.35 m；EⅡ、μⅡ分別為已凍結巖的彈性模量和泊松比。

當r=rⅠ時，凍結圈內壁發生的位移

(6)

當r=rⅡ時，凍結圈外壁發生的位移

(7)

區域Ⅲ為未凍結圍巖，未凍圍巖外半徑趨于無窮大時，令無窮遠處凍脹壓力p=0，該區域內的拉密解為

(8)

徑向位移

(9)

式中EⅢ、μⅢ分別為未凍結圍巖的彈性模量和泊松比。

當r=rⅡ時，未凍結圍巖內壁發生的位移

(10)

根據位移連續條件，襯砌與凍融圈內壁(r=rⅠ)、凍融圈外壁與未凍結圍巖內壁(r=rⅡ)需滿足的位移條件為[17]：

(11)

式中Δ1、Δ2分別為凍脹后已凍結圍巖內壁和外壁的變形量，表達式見式(2)。

(12)

式中ΔV為圍巖凍脹前后的體積差。

由幾何關系，徑向凍脹變形量可表示為

(13)

式中ηr為考慮凍結過程中巖石自身約束作用的飽和凍脹率[18]，表達式見式(14)。

ηr=2.17%n。

(14)

式中n為圍巖孔隙率。

聯立式(1)—(14)，求得凍結圍巖的凍脹力彈性力學解析解為：

(15)

(16)

本文考慮到圍巖和圍巖孔隙水在凍脹受凍脹過程中巖體骨架的約束作用，對徑向凍脹變形量式(13)中的凍脹率進行修正，并代入其他工程數據對修正后的凍脹力解析式進行對比以驗證公式的正確性，如表1所示。

表1 其他凍脹力解析計算值與本文解析計算值對比

通過將如表1所示的4個文獻中的參數代入到本文公式中進行計算，并與其他文獻中的計算值及實測值進行對比，得出利用本文所推導的公式計算出的凍脹力理論值接近其實測值。

利用式(15)—(16)，對新疆布倫口水電站引水隧洞進口50 m處斷面進行計算，得到凍結圈圍巖作用在襯砌上的凍脹力ph=10.856 MPa，作用在未凍結圍巖上的凍脹力為pf=4.002 MPa。將所求ph和pf值代回到式(5)、式(6)中，即可得到凍結圈圍巖受凍脹作用應力場的σr和σθ分別為：

(17)

2 基于水-冰相變的引水隧洞圍巖溫度場分析

2.1 工程地質概況

新疆布倫口水電站引水隧洞工程全長20.14 km，地質情況較為復雜。根據地質測繪，壩址區發育規模較大的斷裂主要為F1、F2、F5。引水隧洞位于高寒半干旱氣候區，季節性溫差較大，雨季和旱季明顯，年降水量不大，洞址區地下水多以基巖裂隙水為主。隧洞洞址處多年平均氣溫為0.7 ℃，年最低溫度為-34.3 ℃，年最高溫度為 35.9 ℃，最冷月平均氣溫為-13.5 ℃；多年平均降水量為127.5 mm，年平均蒸發量為2 297.9 mm；最大積雪厚度約160 mm，最大凍結深度為2.35 m。現場監測結果表明，沿隧洞軸線方向共計出現120多處結冰點，且其凍結程度不同。隧洞平導進口段部分斷面襯砌漏水產生凍脹，洞內凍結現象嚴重，并且凍結的時間較長，影響隧洞的施工和運營，不利于隧道的結構安全。寒區隧洞是否發生凍脹破壞的首要因素是隧洞所處的溫度環境，研究隧洞進口段的溫度場和圍巖凍脹力，是采取防凍害措施的重要依據。

2.2 相變模擬方法

寒區隧洞由于環境溫度與圍巖自身溫度之間存在一定的溫差，從而隨時間產生瞬態溫度傳導，加之隧洞圍巖中裂隙等結構面中水的存在會伴隨溫度的降低而發生相變。水冰相變影響熱傳導參數，設定水的相變域為(-0.5 ℃，0.5 ℃)，使用顯熱容法，將水-冰相變潛熱的影響等效為圍巖比熱容的變化，可構造等價比熱容和熱傳導系數[19]如下：

(18)

(19)

式(18)—(19)中：Cf、Cu分別為圍巖中水融化和凍結時圍巖比熱容；λf、λu分別為融化和凍結時熱傳導系數；L為相變潛熱；t為圍巖溫度。

2.3 有限元計算模型

引水隧洞洞徑D=5 m，越靠近進水口其圍巖破碎程度越高。對該引水隧洞現場溫度監測分析可知，越往隧洞內部越不易發生凍脹[20]，故選取隧洞進口段洞深50 m處的截面作為計算斷面，對圍巖溫度場200 d的變化情況進行分析。

本文采用有限元模擬軟件先對溫度場傳熱進行模擬分析。在建立數值模型時根據圓形截面隧洞的對稱性，選用隧洞斷面的一半作為計算區域，圍巖計算范圍為22.5 m×40 m，斷面總共劃分748個單元，共計803個節點，有限元計算模型網格劃分及尺寸如圖2所示，圖中a、b、c3點分別為拱頂、拱腰、拱底數據提取處。溫度場傳熱過程為瞬態傳熱，計算總時長為200 d，每4 h設置1個增量步，共計800增量步。

(a) 模型網格 (b) 模型尺寸(單位： m)

根據工程地質資料及現場數據，考慮低溫相變，圍巖初始溫度場為3.5 ℃，AE邊界及隧洞洞腔BC邊界的溫度取10月首次出現負溫至次年4月期間200 d內各月的平均溫度，隧洞進口段月平均溫度見表2，熱流邊界DF的熱流密度為0.06 W/m2，其余熱物理學參數見表3。

表2 隧洞進口段月平均氣溫

表3 圍巖熱物理參數

2.4 溫度場結果分析

瞬態溫度傳導200 d后隧洞圍巖溫度發展規律如圖3和圖4所示。由圖3可知，隧洞圍巖溫度場受洞腔內低溫環境影響顯著，圍巖孔隙水受凍結發生相變，凍結鋒面隨著凍結時間的推移向圍巖內部逐漸推移，表現為凍結區逐漸擴大，未凍區逐漸縮小。通過路徑b—d(見圖2)，從溫度場分布云圖中提取凍結深度結果，觀察凍結深度隨時間的變化規律(見圖4)，自傳熱開始到傳熱時間為130 d時，徑向凍結深度隨凍結時間的增長而增加，最大凍結深度達2.04 m，與實際測得的凍結深度2.35 m存在一定誤差，這是由于采用有限元數值模擬水冰相變瞬態傳熱過程時將圍巖默認為各向同性介質，未考慮圍巖內部實際存在的裂隙、節理結構，熱傳導系數會因結構面的存在而受到一定影響。隨著傳熱過程進行至130 d左右，由于環境溫度的升高，隧洞邊界同凍結鋒面溫差Δt減小，加之已凍區和未凍區圍巖導熱系數不同，凍結面開始往外部移動，圍巖已凍區逐漸減小，凍結深度開始減小。至180 d時，邊界出現正溫，圍巖導熱系數增大，溫度變化加快，圍巖溫度全部變為正溫區。

(a) 50 d (b) 100 d

(c) 150 d (d) 200 d

圖4 凍結深度隨時間的變化曲線

3 基于水-冰相變的圍巖凍脹特性分析

巖土體產生凍脹需要滿足以下3個條件： 1)凍脹敏感巖石； 2)凍結溫度和時間； 3)初始水分或補給水分。因此，在計算凍脹力時溫度傳導和孔隙滲流也需要考慮在內。考慮到工程實際，圍巖孔隙水在低溫凍結條件下發生相變產生凍脹，凍結鋒面在圍巖內部移動，并且隧洞凍脹實際邊界條件復雜，是一種強非線性問題。式(17)求得的凍脹作用下的應力場σr和σθ，可通過坐標轉換求得σx和σy，進而解得主應力。但基于一定假設基礎的理論解計算存在一定誤差，且無法直接反映出圍巖的應力分布狀態。為反映圍巖的應力狀態并求得更接近于真實值的近似解，采用有限元軟件對引水隧洞進行凍脹模擬，分析圍巖凍脹特性。

3.1 計算模型及參數

建立考慮低溫水冰相變的水熱力三場耦合的凍脹特性計算模型(見圖2)，對模型左右邊界進行水平約束，下部邊界進行豎向約束，上邊界無約束，并在此基礎上添加力學和滲流參數，進行多場耦合計算。但由于本文使用的有限元分析軟件不能建立直接的溫度-滲流-應力耦合關系，故采用溫度場對滲流-應力耦合的影響來建立間接瞬態耦合模型進行凍脹計算。考慮低溫水-冰相變，水的相變域為(-0.5 ℃，0.5 ℃)，相變潛熱Lf=333.56 kJ/kg，熱膨脹系數為8×10-6。假定巖體為彈塑性材料，服從Mohr-Coulumb屈服準則，根據相關規范和現場資料，計算所需的其他相關參數如表4所示。

表4 隧洞圍巖力學計算參數

3.2 計算結果分析

圍巖受凍脹200 d過程中的應力及位移變形如圖5和圖6所示，從數值分析軟件中提取受凍脹前后隧洞拱頂、拱腰、拱底的應力及位移，對比變化如表3所示。觀察圖5可知，隧洞受凍脹后，拱腰周圍圍巖出現應力最大值，200 d內最大拱腰應力達0.635 MPa，拱頂和拱底位置應力較小。對比凍脹前后各部位應力值，拱頂應力由凍脹前的0.085 MPa增大到0.392 MPa，拱腰應力從0.604 MPa減小到0.527 MPa，拱底應力由0.115 MPa增加到0.399 MPa。凍脹后圍巖拱頂、拱底應力值呈增大趨勢，拱腰呈減小趨勢，各部位變化情況見表5。

(a) 50 d (b) 100 d (c) 150 d (d) 200 d

(a) 凍脹前水平位移 (b) 200 d水平位移 (c) 凍脹前豎直位移 (d) 200 d豎直位移

表5 凍脹前后圍巖應力及變形對比

觀察圖6(a)—(b)可知，受凍脹影響，隧洞拱頂和拱底處水平位移因受到水平方向的位移約束無明顯變化；從有限元分析軟件后處理中提取凍脹前后隧洞拱腰周圍的水平位移，呈水平收縮趨勢，變為-0.36 mm。由圖6(c)—(d)可知，圍巖拱頂豎直方向上位移發生沉降，位移由凍脹前的-9.59 mm變為凍脹后的-8.41 mm，豎向位移變化量為1.18 mm。

綜上，受凍脹影響，圍巖拱頂、拱底應力呈增大的趨勢，拱腰應力逐漸減小，洞壁變形呈收縮趨勢。隨凍結時間的增加，凍脹引發的凍深變化、應力變化和位移變形會進一步使得圍巖結構發生凍害，進而威脅圍巖穩定性。

4 結論與討論

利用彈性力學理論對圍巖凍脹力解析解進行計算，基于有限元分析軟件，考慮水-冰相變下的溫度場、滲流場、應力場相互作用，建立有限元模型對凍脹特性進行數值模擬分析計算，得出如下結論：

1)引水隧洞圍巖在長時間的低溫作用下，孔隙水凍結發生相變，隨時間的推移，圍巖凍結區增大，凍結深度隨之增加，低溫凍結200 d內最大凍結深度為2.04 m。由于凍結深度變化速度受圍巖溫度梯度及環境溫度影響，130 d后凍結深度逐漸減小直至全部變為正溫區。

2)在相變引起的凍脹作用下，拱頂應力由凍脹前的0.085 MPa增大到0.392 MPa，拱腰應力從0.604 MPa減小到0.527 MPa，拱底應力由0.115 MPa增加到0.399 MPa，對比凍結前后應力狀態，拱頂、拱腰、拱底的變化幅度分別為361.2%、-12.7%、247.0%。拱腰處水平位移變為-0.36 mm，表現為水平收縮，圍巖拱頂豎直方向發生沉降，產生1.18 mm變化量。

寒區水工隧洞在低溫水-冰相變影響下凍深逐漸增大，圍巖的應力和位移受低溫凍結作用影響發生明顯改變，為保證圍巖穩定性和隧洞的安全運營，應及時對隧洞施加合理的保溫和支護措施。因此，在本文對圍巖凍脹特性研究基礎上，后續可通過數值模擬、現場監測等方法對保溫層和襯砌進行進一步研究。