孟 堯，姜海波

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

0 引 言

高海拔寒區(qū)隧洞在自然通風條件下寒冷空氣與隧洞圍巖發(fā)生對流換熱，圍巖之間也存在著導熱，從而產(chǎn)生明顯的溫度差且出現(xiàn)溫度應力，使圍巖的應力分布及穩(wěn)定性發(fā)生變化。寒區(qū)隧洞的安全穩(wěn)定離不開隧洞溫度-應力耦合的研究分析，因此有必要基于對流-導熱耦合模型，對寒區(qū)水工隧洞圍巖在不同環(huán)境條件下圍巖不同位置的溫度-應力耦合進行分析。目前，關(guān)于寒區(qū)隧洞圍巖的溫度場及溫度-應力耦合的研究已取得很多成果。丁浩、孟堯等[1-2]研究發(fā)現(xiàn)，洞內(nèi)外大氣和圍巖發(fā)生熱傳遞，得到隧洞徑向和縱向溫度場的變化規(guī)律。嚴健等[3]研究得出，升溫對掌子面周圍的洞壁溫度影響較大。徐衍、劉江川等[4-5]研究得出，圍巖徑向距離增加，溫度影響逐漸消失，圍巖的溫度變化是影響圍巖受力狀態(tài)的重要因素。劉杰等[6]對比不同洞型開挖的應力場和塑性區(qū)分布，確定圓形斷面為最優(yōu)開挖斷面。李松磊等[7]研究發(fā)現(xiàn)，隧洞在無襯砌時通水，溫度荷載會產(chǎn)生較大的拉應力，影響工程的安全穩(wěn)定。劉春龍等[8]研究得出，變溫區(qū)溫度場與隧洞圍巖自重應力場的應力耦合關(guān)系曲線。王偉[9]研究發(fā)現(xiàn)，在對流荷載作用下，距離隧道圍巖越近圍巖溫度變化越大，對流荷載可使隧道開挖后圍巖應力減小，對隧道不同位置圍巖應力的影響不同。王睢等[10]推導出深埋含襯砌有壓圓形隧洞在施工和運行不同工況下的彈塑性解，驗證了模型的正確性。賀耕夫等[11]研究發(fā)現(xiàn)，塑性區(qū)半徑對圍巖單軸抗壓強度在一定范圍內(nèi)的變化非常敏感，甚至達到8～10倍洞徑。湯盛顯等[12]利用彈塑性交界面上應力連續(xù)原則，使用隧道壁面應變?yōu)檫吔鐥l件，求出塑性區(qū)半徑和對應的徑向應變。王明年等[13]得出隧道洞周徑向位移與強度比呈負相關(guān)，塑性區(qū)半徑、隧道洞周徑向位移均與開挖半徑呈線性正相關(guān)的結(jié)論。劉思杰等[14]對開挖期間圍巖破壞機理研究發(fā)現(xiàn)，在頂拱、拱肩多個部位出現(xiàn)應力集中區(qū)。任喜平[15]研究得出羊曲水電站泄洪洞圍巖位移變形及應力分布規(guī)律。馬康等[16]研究發(fā)現(xiàn)，隨著溫度荷載的增加，隧道產(chǎn)生徑向應力與環(huán)向應力，徑向應力在襯砌與圍巖交界處出現(xiàn)應力集中現(xiàn)象。陳張華等[17]研究軟巖力學參數(shù)的敏感性發(fā)現(xiàn)，降低軟巖變形模量后廠區(qū)的位移、應力狀態(tài)、塑性區(qū)均有變化。王平等[18]研究得出，深部高地應力軟巖巷道塑性區(qū)軸比與應力集中是惡性循環(huán)的，提出控制巷道塑性區(qū)最大半徑方向上的圍巖變形是關(guān)鍵。李志雨[19]得出地下隧洞只在溫度作用下以及在溫度和原始地應力共同作用下隧洞圍巖應力的數(shù)值解，并與解析解對比驗證，兩者結(jié)果吻合較好。以上研究分析了不同環(huán)境下寒區(qū)隧洞溫度場的變化規(guī)律，并且進一步研究了溫度荷載產(chǎn)生的溫度應力作用下隧洞圍巖發(fā)生的變化。

高海拔寒區(qū)隧洞貫通后在自然通風的低溫條件下，圍巖的力學性質(zhì)發(fā)生變化，從而影響水工隧洞圍巖的應力分布及穩(wěn)定性。因此，有必要研究不同自然通風溫度和不同風速下的水工隧洞洞口及洞中位置圍巖溫度-應力耦合及其演化規(guī)律。本文以新疆某寒區(qū)水工隧洞為依托，基于M-C本構(gòu)模型，采用有限元仿真計算，對不同自然通風溫度和不同風速下的水工隧洞洞口及洞中位置圍巖溫度-應力耦合及其演化規(guī)律進行系統(tǒng)分析，定量分析計算寒區(qū)隧洞溫度-應力耦合的分布規(guī)律及其變化規(guī)律，研究隧洞內(nèi)風溫、風速對圍巖溫度-應力耦合的影響，并分析其影響機制，得到圍巖瞬態(tài)的應力及塑性應變分布規(guī)律，可為寒區(qū)水工隧洞安全運行提供理論依據(jù)。

1 數(shù)值模型建立

1.1 工程實例

新疆某引水隧洞位于高寒半干旱氣候區(qū)，季節(jié)性溫差較大，最低氣溫-34.3 ℃，最高氣溫35.9 ℃，最冷月平均氣溫-16.66 ℃；多年平均風速2.1 m/s，風向西。根據(jù)現(xiàn)場監(jiān)測得知，引水隧洞周圍的圍巖實測溫度為9 ℃。洞址區(qū)為第四系地層，未見層理，洞室圍巖完整性較好，上覆蓋巖層厚260～300 m。通過地質(zhì)勘查，隧洞區(qū)只有少量上層滯水，圍巖及周圍巖體水分較少，可忽略水對應力的影響。

1.2 對流-導熱耦合模型構(gòu)建及參數(shù)

計算流體與固體壁面接觸換熱用牛頓冷卻公式表達，即

q=h(Ta-Tb)

(1)

式中，h為物體與周圍介質(zhì)的對流換熱系數(shù)；Ta為物體的壁面溫度；Tb為介質(zhì)的平均氣溫。熱傳導用傅立葉導熱基本定律表達

(2)

式中，q為單位傳熱面積熱流量；n為等溫面法線方向向量。對流-導熱公式表示為

(3)

隧洞直徑3 m，根據(jù)實際工程，采用有限元軟件建立模型，圍巖計算范圍取為21 m×21 m×1 000 m的立方體模型。因為溫度荷載對隧洞圍巖應力的影響是一個比較長的時間過程，所以分析某年寒冷季節(jié)11月到翌年3月共5個月的隧洞洞口部位與隧洞中間部位自然通風150 d的溫度-應力耦合分布規(guī)律，取每3 d的平均溫度為1個溫度荷載，共50個溫度荷載，作用在洞壁圍巖和進出口橫斷面上。模型上部施加地應力荷載且無約束，前后面y方向約束，左右面x方向約束，下部鉸鏈約束。模型網(wǎng)格采用C3D8T八結(jié)點熱耦合六面體單元，考慮到橫斷

面有圓形隧洞，給隧洞進行邊布種，使得模型網(wǎng)格劃分更合理。模型網(wǎng)格劃分見圖1所示。根據(jù)工程基本資料得到圍巖力學參數(shù)，見表1。

表1 圍巖力學參數(shù)

圖1 模型網(wǎng)格劃分

1.3 初始條件與邊界條件

本文中圍巖初始溫度9℃，隧洞為通風邊界，空氣和圍巖發(fā)生對流換熱，因此洞壁為對流換熱邊界條件，通風初始溫度4℃。由文獻[20]可知，對流換熱與流體的狀態(tài)、熱力學性質(zhì)及固體壁面的幾何尺寸、粗糙度有關(guān)，并且利用理想化的半圓拱形模型試驗擬合出對流換熱系數(shù)公式。本文研究隧洞不同位置圍巖應力受通風溫度、風速的影響，為滿足單一變量原則，模型中假定：①隧洞圍巖是連續(xù)均質(zhì)、各向同性的，自然風流為連續(xù)介質(zhì)，風速大小沿程不變。②洞壁圍巖的換熱條件是一致的，沿隧洞軸向和徑向不發(fā)生變化。

依據(jù)文獻[20]擬合出的計算公式求對流換熱系數(shù)h。風速為1 m/s時，h=39.96 W/(m2·℃)；風速為2 m/s時，h=59.61 W/(m2·℃)；風速為3 m/s 時，h=88.93 W/(m2·℃)；風速為4 m/s時，h=132.67 W/(m2·℃)。隧洞埋深均取300 m，因為隧洞不同位置埋深相差不大，取最大極限有利于安全穩(wěn)定，還能避免埋深變化的影響。

1.4 溫度-應力耦合機理及求解方法

由文獻[21]可知，溫度-應力耦合有2個方面：一方面，是溫度場對應力場的影響，圍巖溫度變化由溫差產(chǎn)生的溫度應力對圍巖應力場的影響以及溫度變化對圍巖熱力學參數(shù)的影響；另一方面，是應力場對溫度場的影響，圍巖應力的變化使圍巖的內(nèi)部結(jié)構(gòu)發(fā)生變化，影響圍巖的導熱系數(shù)，進而影響圍巖的溫度場，與此同時，圍巖內(nèi)部變形會產(chǎn)生一定熱量也會影響圍巖溫度場。

力學耦合機理：溫度場和應力場之間通過某種力學作用相互耦合；參數(shù)耦合機理：溫度場對巖體的物理性質(zhì)參數(shù)的影響及巖體在不同風溫、不同風速下熱物理性質(zhì)的改變影響溫度場而表現(xiàn)出耦合作用。由文獻[22]可知，溫度場的變化對應力場的影響較大，而應力場在彈性變形范圍內(nèi)對溫度場的影響非常小，從而可以忽略應力場對溫度場的影響。本文為簡化耦合計算，只考慮溫度場對應力場的影響，忽略應力場改變對溫度場的影響。

2 基于對流-導熱耦合模型的圍巖溫度-應力分析

2.1 隧洞洞口和洞中位置溫度模擬結(jié)果

圖2、3分別為基于對流-導熱耦合模型通風風速v為2 m/s時，距離隧洞洞口1 m和隧洞洞中位置寒冷季節(jié)不同天數(shù)的溫度-應力耦合溫度云圖。表2為風速2 m/s時，自然通風情況下距離洞口1 m和隧洞洞中位置洞壁圍巖的溫度值。

表2 不同通風天數(shù)不同位置洞壁圍巖溫度 ℃

從圖2、3可知，自然通風時間增加，對流-導熱對隧洞周圍圍巖溫度的影響范圍逐漸增大，最終穩(wěn)定在以隧洞中心為圓心半徑10 m左右的范圍，距離洞壁圍巖2 m范圍內(nèi)的圍巖溫度變化最大，距離洞壁圍巖越遠，圍巖的溫度變化越小。從表2可知，隧洞60 d與30 d洞口的溫差為1.15 ℃，洞中為1.31 ℃；隧洞90 d與30 d洞口的溫差為6.26 ℃，洞中為6.48 ℃；隧洞120 d與30 d洞口的溫差為1.88 ℃，洞中為2.48 ℃；隧洞150 d與30 d洞口的溫差為9.01 ℃，洞中為9.425 ℃。隧洞在自然通風時，隧洞洞中位置產(chǎn)生的溫差大于洞口。

圖2 不同通風天數(shù)洞口位置圍巖溫度-應力耦合下溫度分布(單位：℃)

2.2 隧洞洞口位置應力及塑性模擬結(jié)果

圖4、5為通風風速2 m/s時，距離隧洞洞口1 m 的位置寒冷季節(jié)不同天數(shù)的溫度-應力耦合最大主應力和最小主應力分布。表3為風速2 m/s時，距離洞口1 m位置圍巖不同部位的主應力(拉應力為正，壓應力為負)。

從圖4、5和表3可知，頂拱、腰拱、底拱部位隨著時間的增加主應力逐漸減小，主應力經(jīng)過削減期到達穩(wěn)定期。主要原因是圍巖發(fā)生塑性變形使圍巖內(nèi)的能量釋放，次要原因是低溫下由溫度差產(chǎn)生的溫度應力(圍巖單元受低溫影響體積收縮產(chǎn)生拉應力)抵消了一部分圍巖壓應力。最大主應力與最小主應力差值在腰拱部位最大，并且隨著時間的增加差值逐漸減小，最大差值為15.05 MPa。

圖3 不同通風天數(shù)洞中位置圍巖溫度-應力耦合下溫度分布(單位：℃)

圖4 不同通風天數(shù)洞口位置溫度-應力耦合下最大主應力分布(單位：Pa)

表3 不同天數(shù)圍巖不同部位的主應力 MPa

提取圖5中的數(shù)據(jù)，作出腰拱部位不同時間不同深度的最小主應力分布，見圖6。從圖6可知，腰拱部位隨著時間的增加最小主應力逐漸減小，隨著圍巖深度的增加最小主應力逐漸減小至穩(wěn)定，表明隧洞圍巖經(jīng)過應力削弱達到應力穩(wěn)定。隧洞壁面圍巖的最小主應力最大，這是由于洞壁處無襯砌結(jié)構(gòu)并且出現(xiàn)應力集中現(xiàn)象造成的。

圖5 不同通風天數(shù)洞口位置溫度-應力耦合下最小主應力分布(單位：Pa)

圖6 腰拱部位在不同天數(shù)不同深度的最小主應力變化

由文獻[23]可知，在非均勻應力場測壓力系數(shù)小于1時，水平方向的塑性區(qū)半徑大于垂直方向的塑性區(qū)，塑性區(qū)形狀近似為橫蝶形。本文側(cè)壓力系數(shù)通過泊松比計算為0.43。圖7為通風風速2 m/s時，距離隧洞洞口1 m的位置寒冷季節(jié)不同天數(shù)溫度-應力耦合塑性應變值云圖。從圖7可知，塑性區(qū)形狀大致為橫蝶形，塑性應變最先出現(xiàn)在腰拱位置，并且隨著通風時間的增加塑性應變值逐漸增大，塑性區(qū)逐漸由腰拱向上下兩端延伸，塑性區(qū)范圍逐漸增大。通風30 d時，塑性應變值最大出現(xiàn)在腰拱位置，為7.906×10-4；通風150 d，圍巖塑性應變值最大出現(xiàn)在腰拱，為0.019 14。通風30 d時，塑性區(qū)半徑約為1倍洞徑，厚度約為0.69 m。隨著通風時間的增加，塑性區(qū)向上下兩端延伸，150 d時塑性區(qū)半徑有微小變化，厚度變?yōu)?.08 m。

圖7 不同通風天數(shù)洞口位置溫度-應力耦合圍巖塑性應變云圖

2.3 隧洞洞中位置應力及塑性模擬結(jié)果

圖8、9為通風風速2 m/s時隧洞洞中部位寒冷季節(jié)不同天數(shù)的溫度-應力耦合最大主應力和最小主應力分布。通風150 d的最大主應力和最小主應力見表4。從圖8、9可知，隧洞洞中部位的最大主應力與最小主應力變化規(guī)律和隧洞洞口位置基本一致。表4中應力值有明顯的變化，洞口位置與洞中位置壓應力值最大相差0.01 MPa。圖2、3溫度場結(jié)果可以很好解釋這種變化規(guī)律，隧洞洞中位置溫度高于洞口，自然通風時洞口位置的溫差小于洞中位置的溫差，所以洞口位置圍巖體積收縮產(chǎn)生的拉應力小于洞中位置，低溫產(chǎn)生的溫度應力抵消了一部分圍巖的壓應力，使洞中位置主應力的壓應力小于洞口位置。

圖8 不同通風天數(shù)洞中位置溫度-應力耦合下最大主應力分布(單位：Pa)

圖9 不同通風天數(shù)洞中位置溫度-應力耦合下最小主應力分布(單位：Pa)

表4 通風150 d不同位置的主應力 MPa

圖10為通風風速2 m/s時隧洞洞中位置寒冷季節(jié)不同天數(shù)的溫度-應力耦合塑性應變云圖。不同通風時間隧洞洞口與隧洞洞中的塑性應變值見表5。從圖10可知，通風150 d，隧洞洞中與隧洞洞口的塑性應變云圖變化規(guī)律基本一致，塑性應變都是30 d時最先出現(xiàn)在腰拱位置，并且隨著通風時間的增加塑性應變值逐漸增大，塑性區(qū)逐漸由腰拱向上下兩端延伸，厚度逐漸增加。從表5可知，隧洞洞口位置的塑性應變大于隧洞洞中位置，一是，洞中位置的溫差大，產(chǎn)生的溫度應力大于洞口位置；二是，洞口和洞中位置的埋深一樣，洞中位置的壓應力小于洞口，故洞中位置的塑性應變也小。

圖10 不同通風天數(shù)洞中位置溫度-應力耦合圍巖塑性應變云圖

表5 隧洞洞口與洞中位置塑性應變 10-3

2.4 隧洞洞中位置不同風速下應力及塑性模擬結(jié)果

圖11為隧洞洞中部位不同通風風速下經(jīng)過寒冷季節(jié)150 d的溫度-應力耦合最小主應力分布。表6為不同情況下溫度-應力耦合的最小主應力值。從圖11和表6可知，對照組(v=1 m/s且無溫度荷載)與v=1 m/s的最大壓應力差為0.16 MPa，最小壓應力差為0.14 MPa，由此可以得出，通風150 d溫度荷載在腰拱位置產(chǎn)生的最大溫度應力(拉應力)為0.16 MPa，并且抵消了部分圍巖壓應力。隨著通風風速的增加，最小壓應力在緩慢減小，最大壓應力也有緩慢減小，但表6中沒體現(xiàn)(云圖保留3位小數(shù))。腰拱位置的壓應力比溫度應力大得多，故腰拱位置的應力受風速影響比頂拱和底拱受風速的影響小，文獻[2]溫度場結(jié)果可以很好地解釋這種變化規(guī)律。當風速增加時，對流換熱系數(shù)增加，隧洞洞壁溫度逐漸降低，空氣對流產(chǎn)生的溫差逐漸增大，產(chǎn)生的溫度應力(拉應力)逐漸增加從而抵消了一部分圍巖壓應力，因而風速增加圍巖壓應力在減小。

表6 隧洞洞中不同情況下最小主應力 MPa

圖11 隧洞洞中不同風速作用下溫度-應力耦合最小主應力分布(單位：Pa)

圖12為隧洞洞中部位不同通風風速下經(jīng)過寒冷季節(jié)150 d的溫度-應力耦合塑性應變云圖。表7為不同風速作用下溫度-應力耦合的最大塑性應變值。從圖12和表7可知，對照組(無溫度荷載)的塑性應變值最大，為0.019 15，與v=1 m/s的塑性應變值相差0.000 41；隨著通風風速的增加，對流換熱系數(shù)增加，溫度應力逐漸增大，抵消了圍巖壓應力，使圍巖的塑性應變逐漸減小。由此可以得出，溫度荷載對隧洞圍巖的塑性應變影響較大，通風風速增大也可以影響隧洞圍巖的塑性應變值，只是影響較小。

圖12 隧洞洞中不同風速作用下溫度-應力耦合圍巖塑性應變云圖

表7 隧洞洞中不同情況下最大塑性應變 10-3

3 結(jié) 語

本文基于對流-導熱耦合模型，對寒區(qū)水工隧洞圍巖溫度-應力耦合場進行了研究，得出以下結(jié)論：

(1)自然通風150 d，對流-導熱對隧洞周圍半徑10 m的圍巖溫度場都有明顯的影響，距離洞壁2 m范圍內(nèi)的圍巖溫度場變化最大。隧洞洞中位置溫差大，對流產(chǎn)生的溫度應力大于洞口位置。

(2)隧洞圍巖頂拱、腰拱、底拱部位隨著時間的增加主應力逐漸減小。一方面，頂拱、腰拱、底拱位置產(chǎn)生的溫度應力(拉應力)抵消了部分圍巖的壓應力；另一方面，圍巖主應力經(jīng)過消減期到達穩(wěn)定期。在隧洞洞壁處出現(xiàn)應力集中現(xiàn)象，主應力值最大，距離圍巖深度越深，圍巖的主應力值逐漸減小到初始應力值。

(3)隧洞洞口與洞中位置的塑性應變、塑性區(qū)厚度都隨著時間增加而增加，經(jīng)過150d最終趨于穩(wěn)定不變，最大塑性應變出現(xiàn)在隧洞洞口腰拱位置，為0.0191 4，塑性區(qū)半徑約為1倍洞徑，塑性區(qū)厚度約為2.08 m。

(4)經(jīng)過不同通風風速150 d的作用，隧洞洞中位置的最小主應力隨風速的增加逐漸減小，但變化較?。凰矶炊粗形恢玫乃苄詰円仓饾u減??；風速改變時溫度差變化較小，對隧洞圍巖的主應力及塑性應變值都有影響但較小。溫差產(chǎn)生的溫度應力最大值在腰拱位置，為0.16 MPa。