馬林建，劉新宇，許宏發(fā)，楊石剛，王子甲

（解放軍理工大學(xué) 工程兵工程學(xué)院，南京 210007）

1 前 言

鹽巖體是國際公認(rèn)的能源（石油、天然氣）存儲的最理想介質(zhì)。自全世界第一座鹽巖地下儲庫在加拿大建成以來，國外共有鹽巖溶腔儲庫已逾 44座，并平穩(wěn)運(yùn)行40年之久。國外鹽巖介質(zhì)中能源存儲相關(guān)力學(xué)問題的研究經(jīng)歷了4個(gè)發(fā)展階段[1]：（1）簡單應(yīng)用階段（20世紀(jì)60年代初）；（2）模型研究與應(yīng)用階段（20世紀(jì)60年代中期至70年代初）；（3）數(shù)值模擬計(jì)算階段（20世紀(jì)70年代末至90年代中期）；（4）損傷模型設(shè)計(jì)階段（20世紀(jì)90年代中期至今）。隨著我國“西氣東輸”配套工程江蘇金壇儲庫建設(shè)工程起步及能源地下戰(zhàn)略儲備項(xiàng)目的興起，國內(nèi)學(xué)者掀起了對鹽巖地下溶腔能源儲庫的穩(wěn)定性分析和評價(jià)方面研究的熱潮[2－5]。

現(xiàn)有的研究工作大多是基于儲庫平穩(wěn)建設(shè)運(yùn)營狀態(tài)下鹽巖溶腔穩(wěn)定性分析，而極端災(zāi)害模式下的相關(guān)穩(wěn)定性研究則鮮有涉及。盡管鹽穴作為地下儲庫安全性比地面結(jié)構(gòu)優(yōu)良得多，但近30年來，國外災(zāi)儲庫泄漏、爆炸乃至鹽穴失穩(wěn)等災(zāi)難性事故時(shí)有發(fā)生。如2004年8月18日，美國德克薩斯州Moss Bluff儲氣庫因密封系統(tǒng)失效造成天然氣噴出，從而釀成持續(xù)4天的大火，燒烤半徑達(dá)400英尺，事故致使價(jià)值至少3600萬美元的天然氣損失。由此開展極端條件下（戰(zhàn)爭、地震、恐怖襲擊、壓力失控）儲庫安全穩(wěn)定性及災(zāi)害防護(hù)控制措施的研究成為國家能源戰(zhàn)略儲備安全的重大需求。本文通過數(shù)值模擬的手段，分析了建腔初期不同初始運(yùn)營壓力失控井噴條件下，儲庫溶腔腔周應(yīng)力狀態(tài)、損傷開展演化規(guī)律及變形收斂特征。

2 井噴失控工況解析

壓力失控的極端情形是由于儲氣庫密閉性失效、天然氣泄漏起火爆炸導(dǎo)致的井噴災(zāi)難。井噴狀態(tài)下井管中氣體為瞬態(tài)流動(dòng)，通過一維穩(wěn)態(tài)流動(dòng)假設(shè)可將問題簡化為有摩阻作用的絕熱管流問題。對于理想氣體絕熱管流運(yùn)動(dòng)適用摩阻作用歐拉方程為

式中：p為氣體絕對壓力；v為氣體截面流速；R為氣體常數(shù)；T為氣體熱力學(xué)溫度；λ為沿程阻力系數(shù)；A和D分別為井管截面面積和井管直徑。

將上式沿管底截面積分到管頂截面可解得絕熱管流質(zhì)量流量公式為

式中：L為井管管長；p1和p2分別對應(yīng)于井管底部和井管口部的氣體壓強(qiáng)。

根據(jù)范諾線理論可知，在等截面管道中，摩阻作用總是使氣流的速度趨于臨界音速。 因此，無論進(jìn)口截面（管底）氣流是亞音速還是超音速，它們的極限速度都是臨界速度，且臨界狀態(tài)只能出現(xiàn)在出口截面（管口）。將底部井口截面處氣體壓力取儲氣庫內(nèi)壓（p =p1），截面流速為

式中：a*為臨界音速。

將式（2）～（4）代入絕熱流動(dòng)能量方程，有

解得

式中：κ為氣體絕熱指數(shù)；N為κ的函數(shù)。

聯(lián)合式（6）與理想氣體狀態(tài)微分方程，對于特定的初始條件p(0)= p0，解得

式中：p0為腔內(nèi)氣體初始壓強(qiáng)，解析得到鹽巖溶腔內(nèi)壓力下降隨時(shí)間呈指數(shù)型降低。

結(jié)合我國金壇鹽巖儲庫建設(shè)運(yùn)營的實(shí)際工況，數(shù)值分析綜合考慮水溶建腔(Ⅰ)、注氣排鹵(Ⅱ)和壓力失控(Ⅲ)3個(gè)階段。其中，建腔階段的分析假定溶腔內(nèi)壓從初始地應(yīng)力值線性減小至相應(yīng)深度的飽和鹵水壓力，建腔周期為2年。建腔后的注氣準(zhǔn)備期為9個(gè)月，注氣排鹵階段溶腔從飽和鹵水壓力線性上升至最大運(yùn)營內(nèi)壓，周期為3個(gè)月。壓力失控階段根據(jù)式（7）分別考慮從17、12、7 MPa卸載至標(biāo)準(zhǔn)大氣壓的情況，具體工況如圖1所示。

圖1 計(jì)算工況Fig.1 Operating mode for simulation

3 模型建立

3.1 地質(zhì)模型的建立

依據(jù)江蘇金壇鹽礦典型地層鹽巖及夾層分布的鉆井資料[6]，計(jì)算選取豎直方向高程1 333.2 m，從上至下依次為地表沉積層（433.2 m）、上覆泥巖層（400 m）、層狀鹽巖地層（158.4 m）和下臥泥巖層（341.6 m），鹽巖層中泥巖夾層厚度分別為1.5、1.3、2.0、2.1 m。水平方向取約6倍溶腔直徑范圍（-300～300 m）。儲庫溶腔采用上下兩個(gè)半橢球，中間接一理想圓柱體，高為98 m，半徑為26 m，容積約為17×104m3。鑒于結(jié)構(gòu)形式、荷載分布及邊界條件的對稱性，取1/4對稱模型，采用六面體單元，網(wǎng)格在鹽巖溶腔和夾層地段密集布置，精細(xì)到1 m 1個(gè)，共劃分為317 920個(gè)計(jì)算單元。建模過程中，通過FLAC3D內(nèi)置Fish語言調(diào)整軟件內(nèi)置基本形狀網(wǎng)格節(jié)點(diǎn)坐標(biāo)，實(shí)現(xiàn)模型網(wǎng)格的精確匹配。模型尺寸及網(wǎng)格劃分見圖2。

圖2 計(jì)算模型尺寸及網(wǎng)格劃分Fig.2 Model size and mesh division

3.2 蠕變模型及參數(shù)

DURUP[7]等通過對 Lemaitre、Norton-Hoff Power、Power、Transient Power、Munson-Dawson等經(jīng)典鹽巖本構(gòu)模型對同一鹽巖試驗(yàn)數(shù)據(jù)擬合對比，認(rèn)為在不考慮溫度的影響下，所有模型均能較好地匹配試驗(yàn)數(shù)據(jù)，相比較Munson-Dawson模型具有17個(gè)材料參數(shù)，Norton Power僅有2個(gè)材料參數(shù)，易從試驗(yàn)數(shù)據(jù)求得。因此，本文針對實(shí)際工程，選用簡單實(shí)用的 Norton Power 模型描述鹽巖的蠕變行為。

Norton材料模型(Norton)一般用于模擬鹽的蠕變特性，標(biāo)準(zhǔn)形式為[8]

根據(jù)Mises應(yīng)力定義及蠕變流動(dòng)法則，鹽巖穩(wěn)態(tài)蠕變率的張量形式為

巖層的塑性行為可用Mohr-Coulomb屈服準(zhǔn)則來判斷（見圖3）。以壓應(yīng)力為負(fù)，則剪切屈服破壞函數(shù)為

式中：N?=(1+sin?)/(1-sin?)；σ1、σ3分別為最大和最小主應(yīng)力；c、?分別為巖石的黏聚力和內(nèi)摩擦角。

考慮到層狀鹽巖張拉破壞特征，拉伸破壞屈服函數(shù)為

式中：σt為鹽巖的抗拉強(qiáng)度。

剪切破壞勢函數(shù)采用不相關(guān)聯(lián)流動(dòng)法則，其表達(dá)式為

式中：Nψ=(1+sinψ)/(1-sinψ)，ψ為膨脹角。

拉伸破壞勢函數(shù)采用相關(guān)聯(lián)流動(dòng)法則，其表達(dá)式為

圖3 Mohr-Coulomb屈服準(zhǔn)則示意圖Fig.3 Sketch of Mohr-Coulomb failure criterion

根據(jù)金壇鹽巖試樣（純鹽巖、含夾層鹽巖和純夾層）巴西劈裂、單軸三軸壓縮和單軸三軸蠕變試驗(yàn)結(jié)果，參考國外鹽巖儲庫計(jì)算所用參數(shù)，并考慮鹽巖層、含鹽夾層與泥巖層蠕變特性的差異，確定本文穩(wěn)定性分析的參數(shù)如表1所示。

表1 儲庫穩(wěn)定性計(jì)算參數(shù)Table 1 Parameters for storage cavern stability

4 計(jì)算結(jié)果及分析

4.1 應(yīng)力場分析

圖4給出了不同階段（工況3）最大、最小主應(yīng)力變化圖。溶腔過程考慮腔體的流變，腔周未發(fā)現(xiàn)明顯應(yīng)力集中。同一區(qū)域最大、最小主應(yīng)力差值不大，2倍洞徑之外與初始地應(yīng)力場趨于一致。但夾層位置產(chǎn)生一定集中度較高的拉應(yīng)力區(qū)（以拉為正），腔周主應(yīng)力等值線圖被夾層隔斷，這在最小主應(yīng)力云圖上更為顯著。相同豎直位置夾層中最小主應(yīng)力絕對值要低于鹽巖層中的，而夾層中最大主應(yīng)力絕對值則要高于鹽巖層，即夾層位置最大剪應(yīng)力要普遍高于鹽巖層。注氣排鹵后應(yīng)力場分布較建腔后變化不大。卸壓階段，腔周出現(xiàn)明顯的應(yīng)力集中，且隨時(shí)間的發(fā)展不斷加強(qiáng)；夾層位置拉應(yīng)力區(qū)逐漸擴(kuò)展，且應(yīng)力集中程度加大。對比3種卸壓工況，這種腔體較大范圍的應(yīng)力不均衡分布在低壓卸荷時(shí)更為突出，因而對儲庫穩(wěn)定性較為不利。

4.2 塑性區(qū)分布擴(kuò)展

圖5～7為3工況（17、12、7 MPa）卸壓初始時(shí)刻(圖 5(a)～圖 7(a))、卸壓中間時(shí)刻（圖 5(b)～圖7(b)）和卸壓后（圖5(c)～圖7(c)）塑性區(qū)分布開展三維視圖（視角與圖8一致）。從圖5(a)～圖7(a)可清楚看到塑性區(qū)最先出現(xiàn)在夾層位置，屈服形式以拉伸屈服為主，并逐漸擴(kuò)展成3條明顯的屈服帶。失控卸壓階段儲庫內(nèi)腔塑性區(qū)向上下延伸直至連成一片，隨后出現(xiàn)在腔頂和腔底部位，此時(shí)屈服以剪切屈服為主。卸壓后塑性區(qū)布滿全腔并向徑向擴(kuò)展一定厚度。腔的穩(wěn)定性是不利的。

4.3 腔體蠕變變形

為分析失控卸壓狀態(tài)下儲庫變形特性及破壞規(guī)律，沿溶腔內(nèi)表面腔頂、腔身和腔底部位布置了 7個(gè)監(jiān)測點(diǎn)，如圖8。

圖8 鹽巖溶腔監(jiān)測點(diǎn)布置示意圖Fig.8 Locations of measuring points of rock salt cavern

從關(guān)鍵部位位移監(jiān)測時(shí)程曲線(圖 9)及溶腔體積損失曲線(圖10)可判定，儲庫處于整體流變收斂狀態(tài)。建腔和排鹵階段溶腔收縮比較平緩，腔頂、腔腰和腔底最大位移為-0.497、-0.653、0.904 m，平均位移速率分別為-0.45、-0.6、0.82 mm/d；體積縮小6.06%，收斂速率V˙/V約為-0.005 %/d。卸壓階段溶腔監(jiān)測點(diǎn)位移不斷擴(kuò)大，腔體加速收縮。據(jù)圖10可知，最高運(yùn)營壓力下（17 MPa）失控卸壓階段儲庫體積減小量為 0.21%，平均體積減小速率為0.015 %/d；12 MPa下自由卸壓使得儲庫體積縮小0.71%，平均速率為0.0358 %/d；7 MPa下卸壓造成儲庫體積損失5.02%，則每天達(dá)0.142%是12 MPa下自由卸壓的近4倍。

圖9 腔周關(guān)鍵測點(diǎn)位移-時(shí)程曲線Fig.9 Curves of key measuring points displacements vs. time

圖10 腔體體積損失率曲線Fig.10 Cavern volume convergence curves

德國克勞斯塔大學(xué)研究建議[9]，一年循環(huán)周期下平均每天的體積減少應(yīng)小于0.008 2%來保證正常運(yùn)營階段儲氣庫的長期使用壽命。陳鋒[10]通過數(shù)值計(jì)算得到恒定最低內(nèi)壓6.5 MPa時(shí)的體積減小速率約為0.004 5%/d，并假定降壓時(shí)間占一個(gè)運(yùn)行周期的1/3時(shí)間算得的降壓階段體積變化速率上限值為0.015 %/d。據(jù)此，只有最高運(yùn)營壓力下的自由降壓才能保能儲庫的穩(wěn)定。需要指出的是文獻(xiàn)[10]認(rèn)為體積收斂速率與降壓速率呈線性關(guān)系，即體積損失速率隨壓降速率增加而線性增加，而本文得到的結(jié)論與此相反。分析原因，區(qū)別在于文獻(xiàn)[10]所模擬的工況為正常運(yùn)營壓力范圍（6.5～14.5 MPa）內(nèi)，人為可控的不同降壓速率對儲庫穩(wěn)定性影響，相比本文研究的則是井噴自由卸壓至標(biāo)準(zhǔn)大氣壓條件下儲庫的穩(wěn)定性。因而，雖然高壓井噴卸荷速率要高于低壓卸荷情形，但溶腔內(nèi)壓按指數(shù)規(guī)律遞減，卸壓初始壓力越低儲庫在低壓維持的時(shí)間就越長（見圖1），從而加速了儲庫的收斂變形?；谝陨戏治觯Э刈杂尚秹弘A段儲庫穩(wěn)定的關(guān)鍵因素為失控初始內(nèi)壓。

4.4 圍巖應(yīng)變速率分析

圖11～13揭示了不同工況下，腔周單元應(yīng)變速率的變化規(guī)律。水溶建腔階段，腔頂、腔腰和腔底拉伸應(yīng)變速率（以拉為正）分別在 4×10-3、1.0×10-3、3×10-3s-1以內(nèi)；注氣排鹵后期，腔內(nèi)巖體由拉伸轉(zhuǎn)為壓縮，最大應(yīng)變速率各是-2×10-3、-1.48×10-3、-1.53×10-3s-1；失控卸壓階段，腔體內(nèi)壁拉伸應(yīng)變速率躍升至10-2s-1數(shù)量級。一般地，第Ⅲ階段降壓速率越快（高壓卸荷）所獲得的應(yīng)變速率越高；同一降壓速率之下腔頂和腔腰的應(yīng)變速率要高于腔底，如12 MPa初始卸壓下，腔頂、腔腰最大應(yīng)變速率為 6.4×10-2、4.9×10-2s-1，腔底為 4.7×10-2s-1。通過最高與最低運(yùn)營壓力失控工況計(jì)算最終界定的腔周圍巖應(yīng)變速率的范圍為：腔頂 1.0×10-2～7.2×10-2s-1，腔腰 1.0×10-2～7.5×10-2s-1，腔底1.0×10-2～5.3×10-2s-1。鹽巖是具有較強(qiáng)流變特性的致密軟巖。廖紅建等[11]對硅藻質(zhì)軟巖進(jìn)行了7.3×10-7～2.9×10-4s-1加載速率范圍內(nèi)的應(yīng)變和應(yīng)力控制式固結(jié)不排水三軸試驗(yàn)，結(jié)果表明，該軟巖具有明顯應(yīng)變速率效應(yīng)，其峰值強(qiáng)度和殘留強(qiáng)度皆隨著應(yīng)變速率的增大而提高，軸向應(yīng)變隨時(shí)間的對數(shù)坐標(biāo)近視呈線性增長趨勢。此外，梁衛(wèi)國等[12]針對層狀鹽巖礦床中 NaCl巖鹽與無水芒硝鹽巖，進(jìn)行了10-5～10-3s-1范圍內(nèi)單軸壓縮強(qiáng)度與變形特性的應(yīng)變速率效應(yīng)研究，研究表明，鹽巖泊松比、峰值應(yīng)變隨加載應(yīng)變速率增加而減小，其變形模量與應(yīng)變速率呈對數(shù)關(guān)系，但強(qiáng)度和彈性模量基本不變。根據(jù)這一結(jié)論結(jié)合正常運(yùn)營條件下鹽巖儲庫應(yīng)變率范圍（10-5～10-3s-1），儲庫腔體穩(wěn)定性是有所保證的。但失控階段的腔體應(yīng)變率相比穩(wěn)定注采運(yùn)行應(yīng)變率增大了幾十倍，10-3～10-1s-1應(yīng)變率范圍內(nèi)鹽巖的強(qiáng)度及變形特性對于失控井噴狀態(tài)下儲庫的安全穩(wěn)定至關(guān)重要，亟待開展進(jìn)一步試驗(yàn)研究。

圖11 腔頂單元應(yīng)變速率時(shí)程曲線Fig.11 Top zone strain rate curves of storage cavern

圖12 腔腰單元應(yīng)變速率時(shí)程曲線Fig.12 Middle zone strain rate curves of storage cavern

圖13 腔底單元應(yīng)變速率時(shí)程曲線Fig.13 Bottom zone strain rate curves of storage cavern

5 結(jié) 論

（1）建腔和排鹵階段，由于應(yīng)變不協(xié)調(diào)硬夾層先行破損；失控卸壓階段，屈服區(qū)由夾層屈服帶擴(kuò)展直至布滿整個(gè)腔體，卸壓初始內(nèi)壓越小塑性區(qū)體積越大。對應(yīng)的應(yīng)力場變化表現(xiàn)為夾層位置出現(xiàn)拉應(yīng)力區(qū)，卸壓階段腔周出現(xiàn)明顯的應(yīng)力集中，且隨時(shí)間的發(fā)展不斷加強(qiáng)；對比3種卸壓工況，腔體較大范圍的應(yīng)力集中在低壓卸荷時(shí)更為突出。

（2）硬夾層的存在一方面對儲庫的收斂變形起到了限制和加筋作用，有利于儲庫單腔的穩(wěn)定；另一方面夾層部位的狹長破損區(qū)易于形成儲庫間的貫通造成儲庫密閉性實(shí)效，則對于儲庫腔群穩(wěn)定性不利。

（3）17 MPa初始內(nèi)壓下失控卸壓，儲庫平均體積損失率為0.015 %/d，12 MPa自由卸壓儲庫體積損失率為0.035 8 %/d，7 MPa自由卸壓儲庫體積損失率則高達(dá)0.142 %/d；根據(jù)克勞斯塔大學(xué)的研究建議，失控卸壓僅發(fā)生在最高運(yùn)營壓力下才能保證儲庫的穩(wěn)定性。

（4）相比建腔和排鹵階段，失控階段腔體內(nèi)壁圍巖應(yīng)變速率從10-3s-1級以內(nèi)增大至10-2s-1級，屬準(zhǔn)靜態(tài)范疇。10-3～10-1s-1應(yīng)變率范圍內(nèi)鹽巖的力學(xué)特性對于研究井噴失控狀態(tài)下儲庫的穩(wěn)定性尤為重要。

[1]吳文, 楊春和, 侯正猛. 鹽巖中能源（石油和天然氣）地下存儲力學(xué)問題研究現(xiàn)狀及其發(fā)展[J]. 巖石力學(xué)與工程學(xué)報(bào), 2005, 29(增刊2): 5561－5568.WU Wen, YANG Chun-he, HOU Zheng-meng.Investigations on studied situations associated with mechanical aspects and development for underground storage of petroleum and natural gas in rock salt[J].Chinese Journal of Rock Mechanics and Engineering,2005, 29(Supp. 2): 5561－5568.

[2]陳衛(wèi)忠, 伍國軍, 戴永浩, 等. 廢棄鹽穴地下儲氣庫穩(wěn)定性研究[J]. 巖石力學(xué)與工程學(xué)報(bào), 2006, 25(4): 848－854.CHEN Wei-zhong, WU Guo-jun, DAI Yong-hao, et al.Stability analysis of abandoned salt caverns used for underground gas storage[J]. Chinese Journal of Rock Mechanics and Engineering, 2006, 25(4): 848－854.

[3]陳鋒, 楊春和, 白世偉. 鹽巖儲氣庫蠕變損傷分析[J].巖土力學(xué), 2006, 27(6): 945－949.CHEN Feng, YANG Chun-he, BAI Shi-wei. Investigation on creep damage of natural gas storage in salt rock layer[J]. Rock and Soil Mechanics, 2006, 27(6): 945－949.

[4]尹雪英, 楊春和, 李銀平. 泥巖夾層對層狀鹽巖體中儲庫穩(wěn)定性影響[J]. 巖土力學(xué), 2006, 27(增刊): 344－348.YIN Xue-ying, YANG Chun-he, LI Ying-ping. Influence of Mudstone Interlayers on stability of oil storage in laminates salt rock[J]. Rock and Soil Mechanics, 2006,27 (Supp.): 344－348.

[5]戴永浩, 陳衛(wèi)忠, 楊春和, 等. 金壇鹽巖儲氣庫運(yùn)營模型試驗(yàn)研究[J]. 巖土力學(xué), 2009, 30(12): 3574－3580.DAI Yong-hao, CHEN Wei-guo, YANG Chun-he, et al. A study of model test of Jintan rock salt gas storage’s operation[J]. Rock and Soil Mechanics, 2009, 30(12):3574－3580.

[6]李銀平, 尹雪英, 高敬. 中鹽金壇壓氣蓄能電站工程初步可行性研究項(xiàng)目（地質(zhì)）報(bào)告[R]. 武漢: 中國科學(xué)院武漢巖土力學(xué)研究所, 2007: 3－30.

[7]DURUP G G, XU J. Comparative study of certain constitutive laws used to describe the rheological deformation of salts[C]//The Mechanical Behavior of Salt Proceedings of the Third Conference. France: Trans Tech Publications, 1993: 75－83.

[8]Itasca Consulting Group Inc. FLAC3D(Version 2.1)users manual[M]. USA: Itasca Consulting Group Inc., 2002.

[9]侯正猛. 金壇地下儲氣庫15口采鹵溶腔穩(wěn)定性評價(jià)技術(shù)服務(wù)報(bào)告[R]. 德國: 克勞斯塔大學(xué), 2004: 18－25.

[10]陳鋒, 楊春和, 白世偉. 鹽巖儲氣庫最佳采氣速率數(shù)值模擬研究[J]. 巖土力學(xué), 2007, 28(1): 57－62.CHEN Feng, YANG Chun-he, BAI Shi-wei. Investigation on optimized gas recovery velocity of natural gas storage in salt layer by numerical simulation[J]. Rock and Soil Mechanics, 2007, 28(1): 57－62.

[11]廖紅建, 蒲武川, 殷建華. 軟巖的應(yīng)變速率效應(yīng)研究[J].巖石力學(xué)與工程學(xué)報(bào), 2005, 24(18): 3218－3223.LIAO Hong-jian, PU Wu-chuan, YIN Jian-hua. Study on strain rate effect of soft rock[J]. Chinese Journal of Rock Mechanics and Engineering, 2005, 24(18): 3218－3223.

[12]梁衛(wèi)國, 徐素國, 莫江, 等. 鹽巖力學(xué)特性應(yīng)變率效應(yīng)的試驗(yàn)研究[J]. 巖石力學(xué)與工程學(xué)報(bào), 2010, 29(1): 43－50.LIANG Wei-guo, XU Su-guo, MO Jiang, et al. Test study of strain rate effects on mechanical performances of salt rock[J]. Chinese Journal of Rock Mechanics and Engineering, 2010, 29(1): 43－50.