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

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

1 前 言

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

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

2 井噴失控工況解析

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

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

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

式中：L為井管管長；p1和p2分別對應于井管底部和井管口部的氣體壓強。

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

式中：a*為臨界音速。

將式（2）～（4）代入絕熱流動能量方程，有

解得

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

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

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

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

圖1 計算工況Fig.1 Operating mode for simulation

3 模型建立

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

依據(jù)江蘇金壇鹽礦典型地層鹽巖及夾層分布的鉆井資料[6]，計算選取豎直方向高程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）。儲庫溶腔采用上下兩個半橢球，中間接一理想圓柱體，高為98 m，半徑為26 m，容積約為17×104m3。鑒于結構形式、荷載分布及邊界條件的對稱性，取1/4對稱模型，采用六面體單元，網(wǎng)格在鹽巖溶腔和夾層地段密集布置，精細到1 m 1個，共劃分為317 920個計算單元。建模過程中，通過FLAC3D內(nèi)置Fish語言調(diào)整軟件內(nèi)置基本形狀網(wǎng)格節(jié)點坐標，實現(xiàn)模型網(wǎng)格的精確匹配。模型尺寸及網(wǎng)格劃分見圖2。

圖2 計算模型尺寸及網(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)典鹽巖本構模型對同一鹽巖試驗數(shù)據(jù)擬合對比，認為在不考慮溫度的影響下，所有模型均能較好地匹配試驗數(shù)據(jù)，相比較Munson-Dawson模型具有17個材料參數(shù)，Norton Power僅有2個材料參數(shù)，易從試驗數(shù)據(jù)求得。因此，本文針對實際工程，選用簡單實用的 Norton Power 模型描述鹽巖的蠕變行為。

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

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

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

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

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

式中：σt為鹽巖的抗拉強度。

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

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

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

圖3 Mohr-Coulomb屈服準則示意圖Fig.3 Sketch of Mohr-Coulomb failure criterion

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

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

4 計算結果及分析

4.1 應力場分析

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

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

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

4.3 腔體蠕變變形

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

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

從關鍵部位位移監(jiān)測時程曲線(圖 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)測點位移不斷擴大，腔體加速收縮。據(jù)圖10可知，最高運營壓力下（17 MPa）失控卸壓階段儲庫體積減小量為 0.21%，平均體積減小速率為0.015 %/d；12 MPa下自由卸壓使得儲庫體積縮小0.71%，平均速率為0.0358 %/d；7 MPa下卸壓造成儲庫體積損失5.02%，則每天達0.142%是12 MPa下自由卸壓的近4倍。

圖9 腔周關鍵測點位移-時程曲線Fig.9 Curves of key measuring points displacements vs. time

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

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

4.4 圍巖應變速率分析

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

圖11 腔頂單元應變速率時程曲線Fig.11 Top zone strain rate curves of storage cavern

圖12 腔腰單元應變速率時程曲線Fig.12 Middle zone strain rate curves of storage cavern

圖13 腔底單元應變速率時程曲線Fig.13 Bottom zone strain rate curves of storage cavern

5 結 論

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

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

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

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

[1]吳文, 楊春和, 侯正猛. 鹽巖中能源（石油和天然氣）地下存儲力學問題研究現(xiàn)狀及其發(fā)展[J]. 巖石力學與工程學報, 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]. 巖石力學與工程學報, 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].巖土力學, 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]. 巖土力學, 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)忠, 楊春和, 等. 金壇鹽巖儲氣庫運營模型試驗研究[J]. 巖土力學, 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]李銀平, 尹雪英, 高敬. 中鹽金壇壓氣蓄能電站工程初步可行性研究項目（地質(zhì)）報告[R]. 武漢: 中國科學院武漢巖土力學研究所, 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)定性評價技術服務報告[R]. 德國: 克勞斯塔大學, 2004: 18－25.

[10]陳鋒, 楊春和, 白世偉. 鹽巖儲氣庫最佳采氣速率數(shù)值模擬研究[J]. 巖土力學, 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]廖紅建, 蒲武川, 殷建華. 軟巖的應變速率效應研究[J].巖石力學與工程學報, 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)國, 徐素國, 莫江, 等. 鹽巖力學特性應變率效應的試驗研究[J]. 巖石力學與工程學報, 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.