蔣中明，鄧自源，唐 棟 ,3，萬 發(fā)，歐陽鈺榕

(1.長沙理工大學(xué)水利工程學(xué)院，湖南 長沙 410114；2.長沙理工大學(xué)水沙科學(xué)與水災(zāi)害防治湖南省重點實驗室，湖南 長沙 410114；3.長沙理工大學(xué)洞庭湖水環(huán)境治理與生態(tài)修復(fù)湖南省重點實驗室，湖南 長沙 410114)

氣墊式調(diào)壓室作為一種性能優(yōu)越的水錘和涌浪控制設(shè)施，相較于常規(guī)調(diào)壓室，在交通、環(huán)境以及投資等方面更具有優(yōu)勢[1-2]。電站在變負荷工況下會發(fā)生水錘現(xiàn)象，調(diào)壓室水面上下波動，氣室氣體因此不斷膨脹、收縮，壓縮空氣溫度和壓力發(fā)生循環(huán)升降變化[2-3]。調(diào)壓室內(nèi)壓縮空氣壓力和溫度的變化不但影響到調(diào)壓室的水力特性[4]，同時也會改變圍巖的應(yīng)力變形狀態(tài)，進而影響圍巖的穩(wěn)定性。

目前國內(nèi)外關(guān)于氣墊式調(diào)壓室的研究重點主要集中在調(diào)壓室的水位變化規(guī)律和壓力變化過程等水力過渡過程特性方面[5-7]，對水位波動過程中調(diào)壓室內(nèi)壓縮空氣的溫度變化較少關(guān)注。實際上空氣的溫度也隨著壓力的改變而變化，并將影響到圍巖和密封層等結(jié)構(gòu)的受力和變形狀態(tài)[8-11]。為分析壓縮空氣熱力學(xué)過程對調(diào)壓室圍巖應(yīng)力和變形的影響，需要準確獲得調(diào)壓室水位變化條件下壓縮空氣的壓力及溫度的動態(tài)變化過程，為此，筆者以考慮傳熱條件下的熱力學(xué)控制方程為基礎(chǔ)，根據(jù)氣墊式調(diào)壓室水位波動變化特點，建立適用于氣墊式調(diào)壓室的壓縮空氣熱力學(xué)分析模型，并利用所建立的模型計算調(diào)壓室內(nèi)的空氣壓力和溫度變化過程，進而將空氣壓力和溫度作為調(diào)壓室圍巖應(yīng)力變形分析的邊界條件，對氣墊式調(diào)壓室圍巖在瞬態(tài)熱傳導(dǎo)和壓力變化工況條件下的受力特性進行研究。

1 氣墊式調(diào)壓室壓縮空氣熱力學(xué)模型

1.1 氣墊式調(diào)壓室水位波動方程

水電機組突然甩負荷工況下，其流量由Q0突變減至0，以調(diào)壓室初始水位為基準水位，建立氣墊式調(diào)壓室壓力變化分析的運動方程、連續(xù)性方程、氣體狀態(tài)方程[3]分別如下：

(1)

(2)

(3)

式中：L、f、v分別為引水隧洞的長度、斷面面積和斷面流速；F為調(diào)壓室水平斷面面積；p0、l0分別為機組在額定工況下，初始時刻t=0 s時氣墊式調(diào)壓室內(nèi)氣體初始壓力和氣室初始高度；p為任意時刻氣墊式調(diào)壓室內(nèi)氣體壓力；α、β分別為引水隧洞水頭損失系數(shù)和阻抗孔水頭損失系數(shù)；n為氣體多方指數(shù)；t為時間；z為水位；g為重力加速度。

氣墊式調(diào)壓室水位波動[12]建議方程如下：

z=Acosφ

(4)

(5)

φ=wt+φ0

(6)

(7)

(8)

式中：ρ為水的密度；v0為初始時刻引水隧洞中的初始流速；A、φ分別為水位偏差值和相位；A0、φ0分別為初始時刻水位偏差值和相位。

1.2 氣墊式調(diào)壓室空氣壓力及溫度計算方程

利用式(1)至式(4)得到調(diào)壓室內(nèi)空氣壓力時，沒有考慮空氣壓縮引起的溫度變化對壓力的影響。研究表明調(diào)壓室內(nèi)的空氣溫度和壓強的變化過程具有高度復(fù)雜的相關(guān)性[13-16]。為獲得調(diào)壓室內(nèi)水位變化條件下空氣壓力和溫度計算的新方法，參考壓氣儲能地下儲氣庫內(nèi)壓縮空氣熱力學(xué)過程的分析方法[14]，建立考慮空氣與圍巖和水體間的熱量相互傳遞影響的氣墊式調(diào)壓室的壓縮空氣溫度和壓力計算模型。根據(jù)熱力學(xué)第一定律，調(diào)壓室內(nèi)部空氣的能量守恒方程為

(9)

單位質(zhì)量空氣內(nèi)能可用下式表示：

(10)

將式(10)兩邊乘以m，然后對時間t求導(dǎo)，可得

(11)

利用體積做功關(guān)系式dW=-pdV和式(11)，可將調(diào)壓室內(nèi)空氣能量守恒方程(9)改寫為

(12)

空氣與接觸界面(壁面和水面)之間的對流換熱速率計算表達式為

(13)

(14)

式中：hcw、hcr分別為空氣與圍巖和水體間的熱交換系數(shù)；Acw、Acr分別為氣體與圍巖、水間的接觸面積；Tsw、Tsr分別為圍巖和水的溫度。

調(diào)壓室中的氣體狀態(tài)方程為

pV=ZmRT

(15)

式中：Z為空氣壓縮因子，假定調(diào)壓室氣體為理想氣體，取壓縮因子Z=1；R為氣體常數(shù)。將式(15)左右兩邊分別對t求導(dǎo)，有

(16)

由(14)和式(16)得壓力對時間的導(dǎo)數(shù)表達式為

(17)

式(14)和式(17)即為調(diào)壓室內(nèi)壓縮空氣的壓力和溫度動態(tài)變化方程。由此可知，調(diào)壓室內(nèi)空氣壓力變化與調(diào)壓室內(nèi)空氣占據(jù)的體積空間V(可由水位z和調(diào)壓室截面積計算求得)和溫度密切相關(guān)。計算體積V所需的z值可以由式(4)得到。利用式(4)計算水位z時所需的壓力由初始條件和式(17)提供，而不是利用式(3)計算氣室內(nèi)的氣體壓力。

在對調(diào)壓室圍巖應(yīng)力及變形特性進行分析時，需將式(14)和式(17)計算得到的壓力和溫度值分別作為壓力邊界和對流換熱邊界條件直接施加到調(diào)壓室的洞壁上作為計算邊界，從而可實現(xiàn)對圍巖的單向熱力耦合變形與受力特性分析。

2 模型合理性分析

為驗證式(14)和式(17)的合理性，假定氣墊式調(diào)壓室不發(fā)生氣體滲漏，即m為常數(shù)的條件下，以小天都氣墊式調(diào)壓室為研究對象，以蔣中明等[17]建議的差分法對式(14)(17)進行求解及合理性分析。表1為小天都水電站氣墊式調(diào)壓室的運行參數(shù)[18]，熱學(xué)計算參數(shù)取值如表2所示。

表1 氣墊式調(diào)壓室計算參數(shù)

表2 熱學(xué)計算參數(shù)

圖1(a)為利用式(17)與式(3)得到的調(diào)壓室內(nèi)氣體壓力變化過程對比，在同一水位波動條件下，采用兩種方法得到的氣體壓力變化規(guī)律基本相同，壓力最大值分別為4.57 MPa、4.55 MPa，最小值分別為3.14 MPa、3.11 MPa。雖然兩者在數(shù)值上存在一定的差異，但偏差較小，均方差僅為0.059 1，且均與杜鵬俠等[18-19]監(jiān)測得到的壓力變化范圍3.28～4.44 MPa較為接近，表明筆者提出的調(diào)壓室內(nèi)壓縮空氣壓力的計算方法是合理的。造成壓力值偏差的原因是在利用式(3)進行求解計算時僅依據(jù)氣體多方指數(shù)的取值來反映壓力的變化規(guī)律，沒有完全考慮氣體溫度變化對壓力的變化的影響；而實際運行過程中壓力與溫度之間是相互影響的。正是由于筆者提出的壓縮空氣壓力計算模型考慮了溫度變化對壓力的影響，從而使得兩種方法得到的計算結(jié)果之間存在一定偏差。

圖1 氣體狀態(tài)演化過程計算結(jié)果分析

圖1(b)為根據(jù)式(14)計算得到的空氣溫度變化過程圖，調(diào)壓室內(nèi)氣體的溫度變化趨勢與壓力變化趨勢基本相同，其中空氣溫度的最高、最低值分別為33.3℃、1.2℃。由圖1(a)可知，調(diào)壓室內(nèi)空氣的峰值壓力在初始值3.68 MPa基礎(chǔ)上增加了約0.89 MPa，空氣的峰值溫度在初始值15℃基礎(chǔ)上增加了約18.3℃，與Lux[19]對地下洞室壓縮空氣的溫度變化研究結(jié)果(即氣體壓力每增加1MPa氣體溫度即升高約13℃)基本一致，表明筆者提出的調(diào)壓室內(nèi)壓縮空氣壓力和溫度變化分析的熱力學(xué)方程是合理的。圖1(b)還表明，在運行過程中氣墊式調(diào)壓室的氣體溫度變化幅度達到了30℃以上，同時最低溫度接近于0℃，存在溫度較低的現(xiàn)象，因此在評價調(diào)壓室結(jié)構(gòu)性能與安全時，不可忽略溫度變化對結(jié)構(gòu)應(yīng)力的不利影響。

3 熱力耦合作用下調(diào)壓室圍巖受力特性分析

為分析氣墊式調(diào)壓室圍巖在調(diào)壓過程中的受力及變形演化特性，基于3DEC軟件平臺，編寫調(diào)壓室內(nèi)氣體壓力和溫度計算的FISH程序，采用熱力耦合數(shù)值分析方法對小天都氣墊式調(diào)壓室圍巖受力和變形特性進行數(shù)值計算。

3.1 數(shù)值模型及計算條件

小天都氣墊式調(diào)壓室為圓拱直墻形，尺寸為80 m×16 m×16.5 m，圍巖以II～III類為主。在模型建立時，考慮部分主要裂隙的影響，其空間分布如圖2(a)所示，裂隙總數(shù)為35條，呈橢圓盤狀；數(shù)值計算模型尺寸為300 m× 200 m×200 m，網(wǎng)格計算單元數(shù)量為670 970個，如圖2(b)所示。為便于分析，在調(diào)壓室周圍設(shè)置若干監(jiān)測點，位置如圖3所示。

圖2 三維數(shù)值模型

圖3 洞室圍巖測點分布

計算參數(shù)采用小天都電站相關(guān)地質(zhì)勘測資料推薦參數(shù)[19-21]；熱力學(xué)計算參數(shù)參考Zhou等[22-23]論文數(shù)據(jù)，具體取值見表3、表4。圍巖的力學(xué)本構(gòu)模型采用摩爾庫倫彈塑性本構(gòu)模型，熱力學(xué)本構(gòu)模型采用各向同性熱傳導(dǎo)模型。模型頂部為壓力邊界，其余邊界均設(shè)置位移約束；調(diào)壓室水位以上部分內(nèi)表面壓力邊界和對流換熱邊界。按自重應(yīng)力計算初始應(yīng)力場，按15℃計算初始溫度場。

表3 物理力學(xué)參數(shù)

表4 熱學(xué)計算參數(shù)

工程實踐表明：當(dāng)水電站突然完全甩負荷時，水流涌向調(diào)壓室，一般按水庫最高設(shè)計水位計算調(diào)壓室水位最高涌波水位；當(dāng)水電站突然增加全負荷時，計算調(diào)壓室的最低涌波水位。因此，計算工況采用完全甩負荷工況，此時調(diào)壓室內(nèi)的水位最高，空氣壓力和溫度最高，對結(jié)構(gòu)影響最為不利。

3.2 調(diào)壓室圍巖溫度時空分布特性

圖4為洞壁處測點的溫度變化過程，其中D組測點位于水下，不受調(diào)壓室內(nèi)溫度的影響。由圖4可知，隨著調(diào)壓室內(nèi)水位波動變化，圍巖溫度也隨著氣體溫度的變化而變化；盡管洞壁測點所處的位置不同，但其溫度變化規(guī)律卻基本相同；不同部位測點的溫度數(shù)值存在一定的差異，其中左側(cè)洞壁A1點的溫度相對于其他測點的溫度更高，極大值達到16.4℃。其原因是距A1點不遠處存在節(jié)理裂隙交匯點，而溫度在裂隙中傳播速度比在壁面中的傳播速度更快[24]，故靠近裂隙的A1點呈現(xiàn)出較高的溫度。此外，B1和C1測點的溫度在工況發(fā)生變化后的第5min達到極大值15.95℃左右，且將在該溫度狀態(tài)維持一定的時間。

圖4 30 min內(nèi)調(diào)壓室各組測點洞壁處溫度變化過程

圖5為調(diào)壓5 min末期時調(diào)壓室橫剖面圍巖的溫度分布。隨著氣體溫度的升高，調(diào)壓室內(nèi)壁直接與氣體進行熱交換，導(dǎo)致圍巖內(nèi)壁的溫度升高至16.0℃左右。由于節(jié)理裂隙的影響，圍巖溫度分布并不均勻，局部區(qū)域特別是節(jié)理裂隙交匯點處的最高溫度升高到16.64℃，見圖中N′、N處。同時，在距洞壁一定范圍0.5 m內(nèi)圍巖的溫度變化梯度較小，與內(nèi)壁溫度基本保持一致，但超過0.5 m后，熱能向圍巖深處傳遞的阻力越來越大，圍巖的溫度變化梯度先增大后減小，溫度逐漸降低至與遠處圍巖溫度一致。由于水的存在，調(diào)壓室底部圍巖沒有與氣體接觸進行熱交換，故其溫度受到氣體溫度變化的影響較弱。

圖5 圍巖的溫度分布

3.3 圍巖應(yīng)力變形特性

圖6給出了調(diào)壓室工況改變后第5 min橫剖面上的第一主應(yīng)力和第三主應(yīng)力分布。由圖6(a)可知，調(diào)壓室內(nèi)壁附近區(qū)域的第一主應(yīng)力基本維持在-6～-4 MPa，且在調(diào)壓室頂部左側(cè)節(jié)理附近及底部拐角區(qū)域出現(xiàn)了點狀區(qū)域的應(yīng)力集中的情況，局部甚至達到-10 MPa，但其分布范圍較小。洞室頂部右側(cè)節(jié)理處出現(xiàn)了高達1.26 MPa的拉應(yīng)力。調(diào)壓室邊墻處的第三主應(yīng)力均表現(xiàn)為壓應(yīng)力，其壓應(yīng)力最大值為-2.38 MPa。由于上部圍巖及底部圍巖受開挖卸荷作用的影響，調(diào)壓室頂拱和底板附近的圍巖均出現(xiàn)了一定程度的拉應(yīng)力，且在洞室頂部節(jié)理處存在拉應(yīng)力集中的情況，其最大值達到1.26 MPa。結(jié)果表明洞室圍巖的節(jié)理裂隙對圍巖的應(yīng)力狀態(tài)有較大的影響，在進行地下洞室設(shè)計時應(yīng)考慮節(jié)理裂隙對洞室穩(wěn)定性的影響。

圖6 圍巖應(yīng)力分布

圖7給出了氣墊式調(diào)壓室工況改變后5min時各組測點的位移與距洞壁距離的關(guān)系。調(diào)壓室圍巖的位移變化的總體規(guī)律為：距離洞壁越遠，位移越??；其中底部圍巖(D組)位移相較于其他部位的位移值更大，而頂部圍巖(B組)的位移值則較??；兩側(cè)洞身圍巖的位移值大小基本一致。

圖7 測點位移與洞壁距離的關(guān)系

3.4 測點溫度、應(yīng)力及變形的變化過程

圖8為B組測點中的部分測點溫度變化過程。隨著調(diào)壓室內(nèi)壓縮空氣溫度和壓力的循環(huán)波動變化，圍巖的溫度也呈現(xiàn)出相應(yīng)的波動變化特征。調(diào)壓室洞壁附近區(qū)域(0.5 m范圍內(nèi))的溫度受到調(diào)壓室空氣溫度的影響較為顯著，如B1、B3點；距離洞壁較遠區(qū)域的溫度受到的影響較小，溫度上升較為緩慢；當(dāng)距離達到5 m后，由于距離遠、熱阻大等原因?qū)е聡鷰r溫度基本不受空氣溫度的影響，如B9點?？傮w上看，調(diào)壓室內(nèi)壓縮空氣溫度的顯著變化對洞周圍巖的溫度絕對升溫幅度影響較小。

圖8 圍巖測點的溫度變化過程

圖9～11分別給出了調(diào)壓室B組部分測點在考慮和不考慮溫度荷載條件下的位移、應(yīng)力變化過程的對比，圖中符號M表示僅考慮力學(xué)荷載的計算成果，TM表示考慮溫度荷載和力學(xué)荷載共同作用的計算成果。由圖9可知，兩種情況下測點處圍巖的位移值基本一致，考慮溫度影響條件下的圍巖位移僅比不考慮溫度時的略大0.01 mm，表明壓縮空氣溫度對圍巖的位移值影響不大。圖10和圖11表明靠近洞壁附近的圍巖應(yīng)力受溫度荷載影響比較顯著。在考慮溫度荷載條件下，圍巖內(nèi)壁附近第一主應(yīng)力較不考慮溫度荷載時增加了-0.14 MPa，為壓應(yīng)力。相較于不考慮溫度荷載的情況，考慮溫度荷載條件下的洞壁附近圍巖第三主應(yīng)力在調(diào)壓過程中先降后增的變化趨勢，其第三主應(yīng)力最大差值較不考慮溫度荷載時的應(yīng)力值減少約0.14 MPa，為拉應(yīng)力。第三主應(yīng)力在4.5 min時達到極小值后又呈現(xiàn)出波動起伏且上升的變化趨勢，在30 min后達到相對穩(wěn)定，并接近于不考慮溫度荷載下的應(yīng)力值。

圖9 典型測點位移變化過程對比

圖10 典型測點第一主應(yīng)力變化過程對比

圖11 典型測點第三主應(yīng)力變化過程對比

4 結(jié) 論

a.氣墊式調(diào)壓室運行工況發(fā)生改變時，氣體溫度可能出現(xiàn)大幅度的變化現(xiàn)象；溫度的大幅度變化將對調(diào)壓室結(jié)構(gòu)的應(yīng)力應(yīng)變狀態(tài)產(chǎn)生不利影響，故分析調(diào)壓室密封性能和圍巖穩(wěn)定性時應(yīng)考慮包含氣體溫度和壓力的熱力耦合效應(yīng)影響。

b.氣體溫度的變化是引起圍巖溫度變化的直接原因，圍巖各區(qū)域溫度分布呈現(xiàn)出顯著的不均勻性，節(jié)理裂隙處的溫度相對其他區(qū)域更高；圍巖溫度影響范圍一般在距內(nèi)壁0.5 m范圍內(nèi)。

c.相較于不考慮溫度荷載的影響工況，溫度荷載對圍巖的位移值影響不大，但對圍巖應(yīng)力(特別是第三主應(yīng)力)有較大的影響。