層狀鹽穴難溶夾層爆破效應數值模擬研究

潘 越， 李二兵， 尤業超， 濮仕坤， 段建立，譚躍虎

(1.中國人民解放軍陸軍工程大學 國防工程學院，南京 210007；2.蘇交科集團股份有限公司，南京 210017)

由于鹽巖具有良好的低滲透性、蠕變性、損傷自愈性、易溶性等特性，深部鹽礦水溶開采所形成的鹽穴是國際上公認的石油、天然氣等能源地下儲備首選方式，已在美國、德國、法國、加拿大等國家廣泛應用[1]，我國目前也在江蘇金壇、江蘇淮安、湖北應城、湖北潛江、河南平頂山等地規劃了大量鹽穴油氣地下儲庫。我國鹽巖礦床大多屬于湖相沉積所形成的層狀結構，硬石膏、泥巖和鈣芒硝等難溶夾層較多[2]。在深部鹽礦水溶建腔過程中，難溶夾層的存在不僅影響溶腔速度，導致腔體形狀不規則，甚至突然垮塌給水溶建腔造成嚴重的安全隱患。研究表明，鹽井聚能爆破技術不失為促使難溶夾層可控垮塌和消除建腔安全隱患的有效辦法[3]。但由于鹽穴地下儲庫需要滿足穩定性和密閉性要求，鹽井爆破一方面要促使難溶夾層可控垮塌，另一方面不能因爆破效應過大而影響整個儲庫密閉性，因此在開展現場先導性試驗之前，進行爆破效應數值模擬分析顯得十分必要和重要，更具有直接的工程指導意義。

目前，國內外針對巖體爆破數值模擬開展了廣泛研究，肖正學等[4]根據室內試驗和國內外礦山爆破的案例研究了不同初始應力對爆破作用的影響；謝源[5]研究了高地應力場對爆破裂紋的影響；袁青[6]基于不確定理論研究了花崗巖地下儲氣庫洞室圍巖爆破特征參量，獲取了最佳特征參量預測模型；Chen等[7]研究了深埋隧道爆破開挖圍巖損傷區與地應力關系；夏祥等[8]通過數值模擬研究了巖體爆破裂紋的產生和擴展機制，得到了粉碎區和裂隙區范圍；孫金山等[9]在Flac3D中模擬了爆破應力波對鄰近隧洞的動力擾動；郭文章等[10]基于DYNA-2D模擬節理巖體爆破過程，量化了節理巖體在爆炸荷載作用下的應力場；高科[11]進行了巖石SHPB試驗模擬研究，證明了半正弦波加載方式是用于巖石類脆性材料大直徑SHPB測試的理想形式；Zhu等[12]利用數值分析方法研究了動靜組合加載SHPB試驗過程中巖體的破壞和損傷特性，得到復合靜力作用下巖石動強度增加機理；Zhu等[13]基于模擬分析手段探討了軸壓、圍壓等因素對巖石動態特性的影響；詹金武等[14]通過LS-DYNA模擬了破碎泥巖注漿結石體SHPB試驗，驗證了與試驗結果吻合較好的HJC模型及其參數；李瑩[15]借助LS-DYNA分析不同初始地應力場條件對爆破的影響，結果表明初始地應力對爆破作用影響在近爆區域較小而在遠爆區域較大；張鳳鵬等[16]借助有限元軟件研究節理條件下巖體爆破機理，發現節理面產生的反射拉伸應力波加強了破巖效果。陳江海等[17]對洋山港水下爆破震動進行測試分析；陳春歌等[18]對水下爆破進行數值模擬，并提出了人工干預措施；王宏[19]模擬分析了水下鉆孔爆破得到了水擊波的傳播規律。

上述分析可知，爆破效應數值模擬研究主要集中在隧道與地下工程、礦山邊坡等領域。關于動載作用下鹽穴穩定性研究分析很少[20]，且由于層狀鹽穴地質構造特殊性，鹽穴內部注滿鹵水，使得上述研究不能直接運用。為此，本文首先進行三軸壓縮狀態下泥巖的SHPB室內試驗，獲得難溶夾層塑性隨動硬化材料模型參數，利用LS-DYNA建立鹽穴實例模型模擬分析泥巖夾層和鹽穴腔體爆破振動等效應，以期為層狀鹽穴難溶夾層爆破方案設計和穩定性控制提供依據。

1 數值模擬模型建立

1.1 鹽穴實例模擬背景

數值模擬實例為湖北應城某埋深1 000 m的層狀鹽穴，其存在厚度8 m的難溶夾層(主要成分是泥巖)，夾層上部鹽巖層厚度約為90 m，夾層下部鹽巖層厚度約為20 m。建腔過程中不斷提升導管，使得溶腔由下而上逐層溶蝕，并在完成夾層下部鹽腔溶蝕后即裝藥進行爆破(見圖1(a)) 促使難溶夾層可控垮塌。保證形成的夾層通道暢通，加快建腔速度，并最終建成如圖1(b)所示的鹽穴儲庫。為簡化溶腔建模，近似看作以距泥巖夾層底部上方21 m的點為圓心作半徑29 m的球體并由泥巖夾層底面剪裁得到的下半部分。

圖1 層狀鹽穴爆破建腔示意圖Fig.1 Schematic diagram of blastingconstruction of bedded salt cavern

1.2 動態力學試驗裝置

考慮到層狀鹽穴難溶夾層埋深約1 000 m，所處地應力約25 MPa。方秦等[21]通過三向壓力下SHPB試驗得出了鹽巖的DIF隨應變率變化曲線，得到了25 MPa下鹽巖動態力學特性參數。因此，本文主要針對夾層泥巖進行圍壓25 MPa下的SHPB試驗。

本次試驗采用Φ75 mm帶主動圍壓的SHPB試驗裝置(見圖2)。裝置主要由動力系統(提供初速度的氮氣瓶和彈膛)、彈性桿件(子彈、入射桿和透射桿)、傳感器(應變片、動態測速儀等)，入射桿和透射桿組成的壓桿系統，長度分別為300 mm，3 500 mm和2 000 mm，材質均為高強度不銹鋼，密度為7 850 kg/m3，彈性模量210 GPa。

1-入射桿端部液壓油缸；2-試件環向圍壓裝置；3-透射桿端部液壓油缸；4-單作用液壓油缸圖2 帶主動圍壓SHPB試驗裝置Fig.2 SHPB experimental device with confininge pressure system

1.3 試樣制作及試驗過程

本試驗試件均采自湖北應城鹽礦泥巖夾層，其主要成分為灰質泥巖平均含量為40%～50%，次要成分為CaSO4、Na2SO4和NaCl。

為了降低摩擦效應(試件與壓桿之間的界面摩擦力)、慣性效應，使試驗過程中試件能較好的滿足均勻性假定，試件的長徑比宜選為0.5。本試驗選用Φ75 mm的壓桿系統。為了使試件截面完全受壓，在加工試件時，直徑略小于壓桿的直徑(直徑差不超過3 mm，否則導致數據失真)，因此，試驗前將試件加工成標準尺寸Φ74 mm×37 mm，打磨過程實時進行平整度檢查(見圖3)，完成后，刷一層油漆，以防空氣中的水分影響試件力學性能。選取5枚理想試樣如圖4所示，試樣各物理特性及巖性特征見表1。

試驗步驟：架設SHPB壓桿系統，在壓桿兩側架設提供反力的拉桿，沖擊桿對中→架設圍壓系統→粘貼應變片、連接線路→在接觸斷面涂抹凡士林進行空打試驗→加軸壓空打試驗→試樣準備→三向壓力下沖擊試驗，沖擊氣壓分別設定為0.7 MPa，0.9 MPa，1.1 MPa→波形信號分析→數據分析處理。

(a)泥巖巖芯

(b)取芯及巖樣打磨圖3 泥巖試樣制作Fig.3 Mudstone specimen making

圖4 試樣成品Fig.4 Specimens of the pre-test

編號厚度L/mm直徑d/mm質量/g密度/(g·cm-3)試樣巖性18-435.4271.92327.722.276灰質泥巖(純泥巖)4-12-336.3073.50373.582.424灰質泥巖(含鹽)5-1-134.2473.57357.52.455灰質泥巖(純泥巖)137-435.4573.47345.532.298灰質泥巖(純泥巖)30-6-535.5873.27351.572.343灰質泥巖(含鹽)

1.4 試驗結果分析

將試驗獲得的電信號波形圖數據導入SHPB數據處理軟件，在平衡圖中比較波與透射波近似重疊，說明試驗過程滿足SHPB試驗應力均勻性假設，結果具有準確性。利用三波法處理數據并導出動態應力應變曲線，如圖5所示，試驗處理結果見表2。從圖5和表2可知，試樣具有明顯的應變率效應，隨著應變率的增加，峰值應力也普遍增加。

圖5 泥巖夾層動態應力應變曲線Fig.5 Dynamic stress-strain curve of mudstone interlayer

編號圍壓/MPa沖擊氣壓/MPa子彈速度/(m·s-2)應變率/s-1峰值應力/MPa18-4251.124.02435.3193.734-12-3251.117.97345.32224.275-1-1250.924.7331.87183.69137-4250.721.24219.84168.2830-6-5250.718.05272.28146.06

1.5 模型參數驗證

王健等[22]通過對不同圍壓下泥巖靜態強度的擬合，獲得25 MPa圍壓下泥巖夾層峰值應力為116.12 MPa。結合表2各應變率下峰值強度計算得到25 MPa圍壓下動力放大系數DIF(見表3)。

表3 25 MPa圍壓下泥巖夾層DIF

在LS-DYNA中采用非線性塑形材料模型(PLSTIC-KINEMATIC)可以反應應變率效應，并且模型簡單，材料模型參數較易獲取。巖石類材料的塑形屈服點與破壞點近似，故不考慮硬化，模型簡化為Cowper-Symonds關系式

(1)

式中：σy為動態條件下的屈服應力(取峰值應力)；σ0為靜態下的初始屈服應力；C和p為應變率效應提高系數。用Cowper-Symonds關系式對DIF進行擬合得到C和p參數值分別為254.4和1.7。根據應城泥巖夾層靜態三軸試驗結果[23]擬合得到25 MPa圍壓下靜態力學參數，綜合得到25 MPa圍壓下泥巖材料參數如表4所示。

表4 25 MPa圍壓下泥巖材料參數

在LS-DYNA中建立與SHPB實體試驗相同尺寸的壓桿模型(入射桿、透射桿長度分別為3 500 mm,2 000 mm；直徑均為75 mm)，試件采用量取的尺寸。所有模型均采用Solid164單元。采用映射網格劃分方法將入射桿沿長度方向劃分350段，透射桿沿長度方向劃分200段，沿桿半徑方向劃分32段。試件軸向分成48段，徑向分成48段。整個系統為對稱結構，為減少計算時間，建立1/4的模型進行計算(見圖6)。為在接觸面單元發生材料失效時剩余單元保持接觸狀態，壓桿與試件之間采用面面侵蝕接觸。另外，考慮到透射桿比較短，而波相對較長，在透射桿的端部加入非反射邊界條件。對稱面采用對稱邊界條件，入射桿自由端端部為應力邊界(加載面)。

圖6 SHPB模型及網格劃分示意圖Fig.6 SHPB model and grid partition schematic diagram

將試驗測得的入射波作為沖擊荷載施加在入射桿的端部可以避免入射波波形對計算結果產生影響。完成荷載施加后，設置接觸、分析時間和能量控制等，通過調節材料模型參數獲得模擬工況下的波形圖，三波法處理得到模擬條件下的應力應變曲線。將數值模擬結果與SHPB實體試驗結果比對(見圖7)，大部分曲線與試驗所獲曲線較為接近，說明塑形隨動強化模型可以較好地擬合泥巖動力特性，材料參數合理并可為后續爆破模擬提供參考。

圖7 圍壓25 MPa下模擬應力應變曲線與試驗應力應變曲線Fig.7 Simulated stress-strain curves and experimental stress-strain curves under confining pressure of 25 MPa

1.6 材料參數確定

巖體材料模型均采用PLASTIC-KINEMATIC模型。上文對泥巖動力特性進行了研究和驗證，用25 MPa圍壓下泥巖的塑性隨動硬化模型參數填充泥巖夾層和下臥層材料；通過方秦等的研究，將25 MPa圍壓下鹽巖的塑形隨動硬化模型參數填充鹽巖材料，如表5所示。

表5 巖體材料參數

炸藥采用MATHIGH_EXPLOSIVE_BURN和EOS_JWL來描述高性能炸藥的材料模型和狀態方程。根據文獻[24]，獲得某高性能TNT炸藥材料參數和狀態方程參數，如表6所示。

鹵水通常以低濃度注入，為了簡化計算，鹵水狀態方程采用純水的參數，而密度略大于水的密度，取值1.025 g/cm3。

表6 炸藥材料參數

1.7 計算模型建立

在LS-DYNA中選用半徑的3.5倍作為模擬的計算模型范圍(模型半徑90 m，高度100 m，見圖8)，由于模型較大，內部結構并非線性尺寸，為了保證網格長短邊尺寸之比滿足計算要求同時達到計算精度要求，劃分80萬個單元，其中最小單元尺寸為0.5 m(指Lagrangian算法網格，ALE網格最小單元尺寸為0.2 m)。模型邊界條件設為25 MPa應力邊界，對稱面設為對稱邊界條件。采用多物質流固耦合法(流體和固體不共節點)進行爆炸數值模擬，避免了Lagrange方法中因水和炸藥產生大變形而導致的計算終止。

圖8 數值模擬模型Fig.8 Numerical simulation model

2 數值模擬結果分析

2.1 夾層爆破試算分析

為了確定爆破的效果以及為后期裝藥量研究設置合理的計算梯度，先進行爆破試算。采用柱狀裝藥，假設炸藥半徑為7.5 cm，高度6 m，重量為43.2 kg(為四分之一藥柱的重量)，考慮到藥柱下部為鹽腔，為了降低端部爆破效應，同時充分利用炸藥，炸藥端部離夾層底部約1 m的位置。由于模型中液體壓縮產生的沖擊波會對爆破作用分析產生干擾，待計算到0.115 s后所有波形均耗散，整個模型應力趨于平穩。故在LS-DYNA中利用重啟動將炸藥起爆時間更改為0.15 s，研究0.15～0.2 s過程中的爆破波形圖。依據文獻[25]將振動峰值速度和等效應力作為兩項重要的判據。

2.1.1 爆破作用下腔體圍巖分析

考慮到鹽腔圍巖損傷過大會影響儲庫的密閉性，在鹽腔圍巖內壁上均勻選取平均弧長5.23 m的5個點A～E(見圖9)。圖10為各個點的有效應力和合速度的時程曲線。從圖中可以看出，鹽腔圍巖各點有效應力變化較小、受沖擊波影響小，故鹽腔圍巖的安全主要受振動速度控制。對比兩圖，可以明顯看出E點距離藥柱間距最大，但是最先受到爆破沖擊的影響產生應力的變化和振動，而A點距離藥柱端部僅9 m，卻是最后受到爆破沖擊的影響。分析可知，沖擊波傳遞到鹽腔圍巖主要通過兩種路徑，一是通過鉆井中鹵水傳遞，而作用于鹽腔圍巖，另一種則是通過夾層巖體傳遞。應力波在巖體沖傳播的速度明顯大于在水中傳播的速度。故出現了A點在爆破作用下的響應滯后于E點。從振動的絕對值來看，在該當量下，A點的振動峰值最大，為15.31 cm/s，其次是E點，為13.61 cm/s。其余各點受爆破作用影響較小，振動速度均不超過 5 cm/s。

2.1.2 爆破作用下夾層動力分析

為了對爆破作用下夾層動力響應有全面認知，在夾層下部和距離夾層底部高度4.5 m處分別布置觀察單元F～I和J～M，上下對應的各點距藥柱間距分別

圖9 鹽腔圍巖觀察單元位置示意圖Fig.9 Schematic diagram of the location of the observation unit of the surrounding rock of a salt cavity

(a) 有效應力時程曲線

(b) 合速度時程曲線圖10 鹽腔圍巖觀察單元時程曲線Fig.10 The time history curve of the observation unit of the surrounding rock of the salt cavity

約為4.6 m，9.4 m，14.1 m和18.9 m。圖11為夾層觀察單元位置示意圖。圖12為各觀察單元的有效應力和合速度的時程曲線。從圖中看出隨著與藥柱間距的增大，應力波峰值和振動合速度峰值迅速降低。

圖11 夾層中部觀察單元位置示意圖Fig.11 Schematic diagram of the position of the observation unit in the middle of the interlayer

圖12 夾層觀察單元時程曲線Fig.12 Time history curve of interlayer observation unit

2.2 炸藥當量對爆破作用的影響

為了探究炸藥當量對夾層爆破產生的影響，依次選用4.8 kg，19.2 kg，49.2 kg和79.2 kg 4個當量(為四分之一藥柱的重量)的炸藥模擬對破壞效果的影響。根據3.1節所選的觀察單元，提取爆破作用下各點的峰值有效應力和峰值振動合速度，并分析其與炸藥當量的關系。圖13為各觀察單元峰值應力與炸藥當量的關系。

圖13 觀察單元峰值應力與炸藥當量的關系Fig.13 The relationship between the peak stress of the observation unit and the equivalent of the explosive

從圖中可以看出，隨著炸藥當量的增加，各觀察單元的有效應力峰值一般都增加，但由于各測點的位置不同，其變化的規律也不盡相同。如圖13中的A～E點，隨著炸藥當量變化，其有效應力變化很小，故而說明，爆破作用下炸藥對鹽腔的直接力學破壞作用小，可以忽略。另外，從圖中可以看出，在離炸藥較近的觀察單元(如F點、J點)，其峰值有效應力受到爆破的作用影響較大，受到炸藥當量的影響比較明顯，而距離藥柱越遠，其峰值應力受到爆破炸藥作用越小。

圖14是觀察單元振動合速度與炸藥當量的關系，從圖中可以看出，隨著炸藥當量的增加，同一個點的振動和速度也增加，距離藥柱較近的位置速度增加的幅度大，較明顯，而距離藥柱遠的位置，振動速度增加的幅度較小。圖14(a)可以看出，鹽腔圍巖B，C和D振動速度受炸藥當量影響較小，而A點和E點受炸藥當量影響較大。圖14(b)和圖14(c)中反應距離效應比較明顯，距離近的觀察點隨藥柱當量的曲線在距離遠的曲線之上，而在圖14(c)中，A點和E點沒有這個規律，其原因在于，巖石傳播的波速及阻抗與水中傳播的波速及阻抗不同。

圖14 觀察單元振動合速度與炸藥當量的關系Fig.14 Observe the relationship between unit vibration speed and explosive equivalence

2.3 炸藥位置對爆破作用的影響

基于上部分計算結果，振動速度是爆破損傷范圍的主要控制指標。選用半徑7.5 cm，高度為6 m的炸藥，藥柱底端到夾層底部的距離直接關系到鹽腔圍巖的振動速度。為此，根據藥柱下端面與夾層板底部的距離(分別為0 m，0.5 m，1.0 m和1.5 m)設置為4個工況。分別討論藥柱端部對鹽腔圍巖及夾層振動的影響。

圖15為A～I觀察單元圍巖振動情況，從圖中可以看出，藥柱下端部與夾層的間距對圍巖振動影響較大。隨著藥柱下端部與夾層的間距逐漸增大，A～D點的峰值振動合速度逐漸降低，主要是由于，藥柱距離下端部的距離越小，藥柱端部距離鹽腔越近，沖擊擾動就越大。另外，E點位于鹽腔邊緣部位，其峰值應力主要受到巖石中應力波傳播的影響，而鹵水傳遞的沖擊波對其影響較小，故而隨著間距的變化，其值變化相對較小。圖15(b)中，隨著距離增大，峰值合速度先增大后減小。主要是當間距為0 m時，夾層底部巖體主要受到沖擊壓作用，所以其振動速度較大。隨著距離的增大，夾層底部在沖擊波的作用下，會發生反射，從而提高了振動速度。當距離較大后，夾層底部受到沖擊波和反射波都較小，故而其振動速度又發生降低。從H點和I點可以看出，觀察單元距離藥柱的絕對距離較大，其受到的沖擊波和反射波均較小，藥柱下端部與夾層的間距對其影響就相對較小，反映出來的速度變化就較小。

圖15 藥柱下端部與夾層底部間距對爆破振動的影響Fig.15 The influence of the distance between the bottom of the grain and the bottom of the interlayer on blasting vibration

2.4 爆破效應綜合分析

爆破主要受兩個因素控制，一方面，需要防止鹽腔圍巖損傷過大，影響鹽腔的密閉性；另一方面，需要實現預期的爆破效果，獲得理想的泥巖夾層破壞半徑。由于巖石類材料抗拉強度遠低于抗壓強度，以抗拉強度為準則，根據簡單的平面波估算鹽腔圍巖的振動速度承受范圍，依此設定裝藥量和裝藥位置，驗算夾層破壞是否符合預期效果。

(2)

式(2)可改寫為

(3)

式中：V為巖石的振動速度；[V]為臨界振動速度；c為應力波在巖石長傳播的速度；[σt]為極限抗拉強度；E為彈性模量。

超聲波在巖石中傳播波速與圍壓的關系。考慮到應力波在巖石中傳播波速與聲波相近，此處引用超聲波波速代替應力波波速，近似取4 350 m/s。

2.4.1 鹽腔破壞控制

根據李二兵的研究，鹽巖的抗拉強度為1.04 MPa。由于巖石的動態抗拉強度的應變率效應小，近似取動態抗拉強度等于準靜態抗拉強度。根據獲得的鹽巖抗拉強度(直接拉伸強度)和鹽巖的彈性模量，近似獲取臨界振動速度為16 cm/s。選取圖9中各單元為鹽腔破壞控制點。從圖15(a)可以看出，A點受到藥柱端部影響最大，且當藥柱底部距離夾層底部距離大于1.0 m時，距離的影響變小。因此，藥柱端部距離夾層底部1.0 m處為合理的安放位置。分析圖14(a)可知，在藥柱下端部距離夾層下部1 m時，裝藥量應控制不超過49.2 kg(為四分之一藥柱的重量)。

2.4.2 夾層破壞控制

依據夾層破壞效應，為了充分反映夾層的破壞范圍(包含爆破產生的直接效應和間接效應)，以使夾層產生裂縫為破壞標準，即夾層拉應力達到泥巖的抗拉強度。李二兵的研究給出泥巖抗拉強度2.19 MPa。為了更好的對計算結果定量分析，在模型圖中選擇高度4.5 m處水平面為觀察面，如圖16所示。在LS-DYNA中，以受壓為負，受拉為正，故選擇單元的最大主應力為觀察項目，最大主應力達到極限抗拉強度，即認為夾層受拉開裂。

圖16 破壞范圍觀察示意圖Fig.16 Observation map of damage range

圖17為破壞效果與裝藥量的關系。從圖中可以看出，破壞效果與裝藥量有直接的關系，隨著裝藥量的增加，破壞范圍增大，但隨著裝藥量增加，破壞半徑增大幅度減小。按照四分之一裝藥量考慮，19.2 kg時破壞半徑為8.66 m，49.2 kg時破壞半徑為11.3 m，可見裝量從19.2 kg增加1.56倍到49.2 kg，破壞半徑僅從8.66 m擴大了0.3倍到11.3 m。難溶夾層厚度8 m，一般8 m左右的破壞半徑足以滿足難溶夾層可控垮塌要求。因此，難溶夾層爆破四分之一裝藥當量可選在19.2 kg左右，整體裝藥量約76.8 kg。

圖17 破壞范圍與裝藥量的關系Fig.17 The relationship between the range of damage and the amount of charge

3 結 論

本文通過SHPB試驗獲取難溶夾層材料參數。使用塑性隨動硬化材料模型擬合，并使用LS-DYNA驗證材料模型的正確性和模型參數的準確性。對鹽穴實例建立數值模型進行研究分析，具體結論如下：

(1)在LS-DYNA中采用塑形隨動強化模型可以較好地擬合泥巖動力特性，與25 MPa圍壓下SHPB試驗結果具有較好的一致性。

(2)作用于鹽腔圍巖應力波通過鹵水和夾層巖體兩種介質傳遞，應力波在巖體沖傳播的速度明顯大于在水中傳播的速度。

(3)通過對本鹽穴實例進行數值模擬分析，為滿足鹽腔穩定性和密閉性的要求，若選取藥柱半徑7.5 cm，高度6 m，爆破炸藥當量宜選在76.8 kg左右，并布置在藥柱下端面與夾層板底部距離1 m的位置可以達到促使夾層垮塌的目的，同時保證鹽巖溶腔圍巖的氣密性。