田曉建，姚安林,2，徐濤龍，蔣宏業(yè)，李又綠

(1.西南石油大學(xué) 石油與天然氣工程學(xué)院，四川 成都 610500； 2.油氣消防四川省重點實驗室，四川 成都 610500)

0 引言

油氣管道的建設(shè)過程中，由于受各種影響因素的制約，多條管道并行敷設(shè)逐漸增多[1]。在一些地區(qū)因規(guī)劃或集約化用地而往往要求盡量減小并行管道的間距[2]。加之城鄉(xiāng)建設(shè)發(fā)展較快，管道多處占壓，地區(qū)等級升級的情況越來越多，對在役管道的安全提出了挑戰(zhàn)[3]，管道一旦發(fā)生爆炸或燃燒等事故，在沖擊波或熱輻射的作用下必將影響相鄰管道的安全。因此，開展管道在爆炸荷載下的動力響應(yīng)方面的研究顯得尤為重要。

國內(nèi)外學(xué)者開展了相關(guān)的研究工作并取得了一定成果。Parviz等[4]采用數(shù)值模擬的方法研究了爆炸荷載下土質(zhì)、管材、TNT當量等對管道動力響應(yīng)的影響，發(fā)現(xiàn)低密度土壤中的爆炸對管道的破壞較??；Adibi等[5]研究了爆破強度、管材厚度、埋深、混凝土保護層對管材變形的影響，結(jié)果表明埋深的增加使管道的變形減?。籚ivek等[6]利用激波管產(chǎn)生沖擊波進行室內(nèi)模擬球形爆炸對管道的影響，結(jié)果表明管道受到?jīng)_擊波的沖擊立刻發(fā)生變形，壓力阻尼隨深度增加而增大；Mokhtari等[7]建立了基于歐拉-拉格朗日耦合法的三維有限元模型，分析了爆炸荷載下埋地管道的擾動(變形)行為；文霞等[8]基于ALE方法建立了隧道并行輸氣管道爆炸模型，分析不同工況下鄰管對爆炸沖擊的動力學(xué)響應(yīng)；王德國[9]基于數(shù)值模擬的方法計算不同并行間距下架空天然氣管道受并行管線爆炸沖擊的超壓及變形量。

已有研究均假定管道本體完好，缺乏對于含內(nèi)腐蝕缺陷管道的動力響應(yīng)研究。鑒于此，研究假定并行管道中的一條發(fā)生物理爆破，基于爆能相等原則確定TNT當量，在驗證模型可行性的基礎(chǔ)上建立基于SPH-FEM的管土-炸藥耦合模型，重點分析了爆炸荷載下缺陷參數(shù)和爆心距對管道動力響應(yīng)的影響，基于彈塑性失效準則求得管道安全平穩(wěn)運行的臨界深度和極限距離。研究可為埋地輸氣管道極端災(zāi)害下的風(fēng)險評估提供技術(shù)支撐，為并行管道可能的抗爆隔爆設(shè)計提供模擬數(shù)據(jù)支持。

1 爆炸能量與TNT當量換算

1.1 物理性爆炸能量估算

壓力容器因本構(gòu)缺陷或壓縮氣體超壓引發(fā)物理爆炸時，氣體在容器破裂之前不會產(chǎn)生物態(tài)變化，只有簡單的降溫膨脹過程。容器內(nèi)氣體絕熱膨脹所做的功轉(zhuǎn)變?yōu)楦邷貧怏w的爆炸能量，假定容器內(nèi)氣體的起始狀態(tài)參數(shù)為(P0，V0)，經(jīng)過絕熱膨脹后其狀態(tài)參數(shù)為(Pa,Va)，根據(jù)絕熱方程PVk=C得到P0V0k=PaVak,則氣體膨脹所做的功為:

(1)

若氣體絕熱膨脹所做的功全部轉(zhuǎn)化為爆破釋放的能量，則容器內(nèi)壓縮氣體的爆破能量可以表示為：

(2)

式中：E為容器內(nèi)壓縮氣體的爆破能量，kJ；P為容器內(nèi)氣體的絕對壓力，MPa；V為容器的體積，m3；k為氣體的絕熱指數(shù),取值范圍為1.2～1.4，燃氣一般取為1.29[10]。

1.2 等效TNT當量換算

當前TNT當量的換算方法主要有基于超壓、基于爆熱的等效計算方法和蒸汽云爆炸的當量計算方法。研究采用基于超壓模型的等效TNT當量方法，計算公式[11]：

(3)

式中：WTNT為TNT當量，kg；E為超壓爆炸能量，kJ；QTNT為TNT的爆炸熱，取值范圍為4 230～4 836 kJ/kg,通常取4.52×103kJ/kg。

2 數(shù)值模擬

2.1 物理模型建立及有限元網(wǎng)格劃分

運行管線某段發(fā)生物理爆炸時，非爆破段管輸介質(zhì)對爆破段介質(zhì)的膨脹做功有一定的干涉效應(yīng)，忽略干涉效應(yīng)僅考慮爆破段介質(zhì)的做功可能與實際不符，但是目前對于干涉效應(yīng)的研究尚缺乏相應(yīng)的實驗驗證。所以，為了開展后續(xù)工作，筆者采用已有研究[12-14]所作的簡化處理，取φ1 016 mm×26.2 mm，輸送壓力為10 MPa的爆管長度15.6 m，由式(2)求得爆破能量為280 646.387 4 kJ，由式(3)求得等效TNT當量為62.089 9 kg。由于結(jié)構(gòu)的對稱性，利用對稱條件，建立如圖 1所示的1/4物理模型，為方便建模炸藥量取16 kg，藥包尺寸為0.17 m × 0.34 m ×0.17 m，土體尺寸為8 m×6 m×6 m。炸藥及炸藥周圍大變形區(qū)域采用SPH粒子，其余部分采用8節(jié)點六面體實體單元，為提高計算精度，在炸藥、管道和炸藥周圍大變形區(qū)域細化網(wǎng)格邊長至50 mm，缺陷網(wǎng)格為5 mm，遠離爆炸近區(qū)單元邊長過渡至500 mm。在模型對稱面施加對稱約束，底面和其余兩側(cè)面施加無反射邊界約束，頂面為自由表面。

圖1 物理模型及網(wǎng)格劃分Fig.1 The physical model and grid partition

圖2為含內(nèi)點蝕缺陷管道示意圖，A點為迎爆面正對爆心缺陷中心處，固定ai，bi(ai=bi)變化ci來改變?nèi)毕莸纳疃?；固定ci變化ai，bi(ai=bi)來改變?nèi)毕莸谋砻娉叽?；固定ai，bi和ci，移動缺陷到B，C，D和E點來改變?nèi)毕莸奈恢茫还潭╝i，bi(ai=bi)和ci變化圖 1 中的R來改變爆心距離。由此開展多工況下的算例分析。

圖2 缺陷管道示意Fig.2 Schematic diagram of defective pipe

2.2 管道本構(gòu)模型

材料本構(gòu)采用MAT_PLASTIC_KINEMATIC模型。該模型遵從Von-Mises屈服準則，對于爆炸沖擊荷載作用下的金屬材料較為適用，利用其得出的結(jié)果與實驗數(shù)據(jù)吻合較好[14]，因而被廣泛應(yīng)用。其本構(gòu)方程為：

(4)

表1 管道的材料參數(shù)值Table 1 Material parameter value of pipe

2.3 土壤本構(gòu)模型

土體本構(gòu)選用帶有失效準則的MAT_SOIL_AND_FOAM砂土模型，其理想塑性屈服函數(shù)為：

(5)

式中：Sij為偏應(yīng)力分量；p為壓力，MPa；a0，a1和a2為剪切屈服面常數(shù)。主要參數(shù)值：密度為1.8 g/cm3，剪切模量為63.8 MPa，體積卸載模量為30 GPa，拉伸截止值為-6.9 kPa，剪切屈服面參數(shù)a0=3.4×109Pa2，a1=7.033×104MPa，a2=0.3。其他材料參數(shù)[17]見表2。

表2 砂土狀態(tài)方程參數(shù)Table 2 Parameters of state equation of sand

2.4 炸藥本構(gòu)模型

炸藥采用高能引爆燃燒材料模型MAT_HIGH_EXPLOSIVE_BURN和JWL狀態(tài)方程：

(6)

式中：p為壓力，MPa；E0為單位體積爆轟產(chǎn)物的內(nèi)能，GJ/m3；V為爆轟產(chǎn)物的相對體積；A，B，R1，R2和ω為定值。炸藥材料參數(shù)[18-19]取值見表3。

表3 炸藥材料參數(shù)Table 3 Parameters of TNT charge

3 含缺陷管道動力響應(yīng)結(jié)果分析

3.1 模型可靠性驗證

為驗證SPH-FEM耦合算法的可靠性，分別建立裝藥0.2 kg，藥包埋深0.8，0.9和1 m的驗證模型，研究土體爆炸空腔的變化發(fā)展規(guī)律。圖3為爆腔模擬尺寸與穆朝民等[20]提出的爆腔預(yù)估公式計算值對比曲線圖，包含了不同藥包埋深下的最終爆腔形態(tài)，兩者吻合較好，爆腔水平半徑R2和裝藥上方半徑R3均隨藥包埋深的增加而減小，前者減速平緩后者減速急劇；驗證模型中地表面形成鼓包而無可見彈坑，隨藥包埋深的增加鼓包隆起度減小，爆腔內(nèi)壁出現(xiàn)貫穿裂紋，爆腔由“倒梨形”向“橢圓形”過度。

圖3 爆腔尺寸與藥包埋深的關(guān)系Fig.3 The relationship between the size of brust cavity and the depth of TNT charge

土壤質(zhì)點的峰值振速可由薩道夫斯基峰值振速經(jīng)驗公式[21]計算：

(7)

兩端取對數(shù)：

(8)

式中，V為峰值振速，m/s；K為場地系數(shù)；W為裝藥量，kg；R為爆心距，m；η為衰減系數(shù)。

在驗證模型中提取不同距離下土壤質(zhì)點的峰值振速，分別進行函數(shù)擬合，擬合曲線如圖4所示。公式中場地系數(shù)和衰減系數(shù)均接近，擬合曲線呈線性關(guān)系，與經(jīng)驗公式規(guī)律吻合，驗證了SPH-FEM耦合法的可靠性。

圖4 土壤質(zhì)點峰值振速擬合曲線Fig.4 Fitting curve of peak velocity of soil particle

3.2 缺陷深度

為研究缺陷深度對管道動力響應(yīng)的影響，在R=4 m時，設(shè)定圖2中缺陷尺寸ai=bi=4 cm，缺陷深度ci由0.1t等距增至0.7t共4種設(shè)計工況，進行算例分析。

圖5為不同工況下缺陷中心橫截面上峰值位移分布。由圖5可知，缺陷中心橫截面上的峰值位移隨ci的增大而增大，增速由緩慢到急劇，最大位移位于臨近管道的頂部，這是由于隨著時間的推移，爆炸空腔由“圓形”向“倒梨形”過度，爆腔形心上移，水平半徑持續(xù)增大，因管體擠壓管頂明顯上移。

圖5 缺陷中心管道截面上位移峰值分布Fig.5 Distribution of peak displacement on the cross section of a defect center pipe

圖6 缺陷中心管道截面峰值應(yīng)力分布Fig.6 Peak stress distribution diagram of the cross section of a defect center pipe

圖6為不同工況下缺陷中心橫截面上峰值應(yīng)力分布。由圖6可知，截面上各處應(yīng)力分布不均勻，缺陷處存在明顯的應(yīng)力集中現(xiàn)象，在缺陷邊界處應(yīng)力呈現(xiàn)先減小后增加的現(xiàn)象。隨著缺陷深度的增加，橫截面上各處應(yīng)力呈增大趨勢，缺陷中心峰值應(yīng)力增速最快、應(yīng)力最大，最大應(yīng)力分別為315.26，355.61，426.18和626.46 MPa。橫截面上各處應(yīng)力分布沿水平線非對稱，鄰近爆心水平線下方120°范圍內(nèi)管體應(yīng)力峰值大于水平線上方相同角度范圍內(nèi)的應(yīng)力值，且ci越大，這種現(xiàn)象越明顯。這是由于隨著時間的推移，爆轟產(chǎn)物急劇擴張迫使周圍土壤粒子急劇壓縮移動并形成爆腔的過程中，缺陷中心水平線以下爆炸能量全部用于壓縮土壤粒子，而水平線上方爆炸能量除部分用于壓縮土壤粒子外，還有部分轉(zhuǎn)化為土壤粒子的機械能，土壤粒子透過地表自由面被拋向空中并形成鼓包，導(dǎo)致部分能量耗散，因此呈現(xiàn)應(yīng)力分布沿水平線非對稱的現(xiàn)象。

圖7為不同工況下缺陷中心橫截面上峰值振速分布。由圖7可知，橫截面上峰值振速分布呈“草莓”狀，峰值振速分布沿水平線趨近對稱，隨缺陷深度的增加呈增加的趨勢，缺陷中心處峰值振速最大。

圖7 缺陷中心管道截面峰值振速分布Fig.7 Peak vibration velocity distribution map of the defect center pipe section

綜合圖6～7可知，缺陷中心橫截面上的位移、應(yīng)力、振速分布沿水平線呈非對稱關(guān)系；隨缺陷深度的增加，均呈現(xiàn)增加的趨勢，增速呈現(xiàn)由平緩到急劇的轉(zhuǎn)變；缺陷處為響應(yīng)的熱點區(qū)域，缺陷中心為危險點。為此，將缺陷中心單元在不同ci下的峰值應(yīng)力進行擬合，擬合曲線如圖8所示，應(yīng)力隨著ci的增加呈指數(shù)增加。擬合公式為：

(9)

式中：y為缺陷中心處的峰值應(yīng)力，MPa；x為缺陷的深度，cm。為保證管道的安全運行，基于塑性失效準則，由公式(9)求得滿足條件的臨界腐蝕深度為1.659 96 cm，約為0.633 6倍的管道壁厚。

圖8 峰值應(yīng)力擬合曲線Fig.8 Peak stress fitting curve

3.3 不同缺陷表面積

為研究缺陷表面尺寸對管道的動力響應(yīng)規(guī)律，在爆心距R=4 m時，設(shè)定圖2中缺陷尺寸ci=0.3t，ai=bi由4 cm等距增至10 cm共4種設(shè)計工況進行分析。

圖9為不同工況下缺陷中心橫截面上位移峰值分布。由圖9可知，缺陷中心截面上的位移隨缺陷表面積的增加呈減小趨勢，背爆面上一定角度范圍內(nèi)呈現(xiàn)先減小后增加再減小的波動，位移分布總體差異較小。

圖9 缺陷中心管道截面上位移峰值分布Fig.9 Distribution of peak displacement on the cross section of a defect center pipe

圖10為不同工況下缺陷中心橫截面上峰值應(yīng)力分布。由圖10可知，隨著缺陷表面尺寸的增大，應(yīng)力集中范圍增大，缺陷中心處應(yīng)力有小幅度的增加，應(yīng)力分別為355.79，359.89，360.64和366.29 MPa，缺陷外應(yīng)力分布無顯著變化。此外，水平線兩側(cè)管體截面應(yīng)力分布非對稱，水平線下方應(yīng)力較上方大。

圖10 缺陷中心管道截面峰值應(yīng)力分布Fig.10 Peak stress distribution diagram of the cross section of a defect center pipe

圖11為不同工況下缺陷中心橫截面上峰值振速分布。由圖11可知，缺陷中心管道截面上的峰值振速分布曲線呈“草莓”狀，迎爆面上的振速值大于背爆面的振速,水平線下方峰值振速略大于水平線上方峰值振速，隨ai=bi的增加，各處峰值振速無明顯變化，缺陷中心處峰值振速最大，最大分別為2.437 2，2.489 7，2.482 6和2.487 5 m/s。

綜合圖10～11可知，缺陷表面尺寸對缺陷中心管體橫截面上的位移分布、應(yīng)力分布和峰值振速分布影響較小。

圖11 缺陷中心管道截面峰值振速分布Fig.11 Peak vibration velocity distribution map of the defect center pipe section

3.4 不同缺陷位置

管道內(nèi)腐蝕一般發(fā)生于管道下部壁面[22]，為研究缺陷位置對管道動力響應(yīng)的影響，研究僅針對缺陷位于管道下半部的情況。R=4 m時，設(shè)定圖2中缺陷尺寸ai=bi=40 mm，ci=0.3 t，在A，B，C，D，E中變化缺陷位置，分析管體的動力響應(yīng)特征。

圖12為不同缺陷位置處的位移峰值分布。由圖12可知，不同缺陷位置處的位移分布曲線存在微小差異，位于同一水平高度的缺陷，迎爆面上缺陷的曲線值高于背爆面的值，這是由于不同缺陷位置處的爆心距不同。

圖12 缺陷中心管道截面上位移峰值分布Fig.12 Distribution of peak displacement on the cross section of a defect center pipe

圖13為不同缺陷位置處的峰值應(yīng)力分布。由圖13可知，缺陷處存在明顯的應(yīng)力集中現(xiàn)象，迎爆面上的應(yīng)力集中程度較背爆面大，缺陷外各曲線基本重合，反映了缺陷位置對管道動力響應(yīng)的響應(yīng)影響較小。此外，應(yīng)力集中最明顯的位置并非位于迎爆面上A處，而是位于迎爆面上的B處，這是由于炸藥中心上方能量部分轉(zhuǎn)化土壤粒子的動能，粒子被拋向空中形成鼓包，而下方能量全部用于壓縮土體，力最終作用于管道，因而出現(xiàn)B處應(yīng)力集中程度大于A處的現(xiàn)象。

圖13 缺陷中心管道截面峰值應(yīng)力分布Fig.13 Peak stress distribution diagram of the cross section of a defect center pipe

圖14為不同缺陷位置處的峰值振速分布。由圖14可知，峰值振速分布曲線呈“草莓”狀，缺陷位于A點時，峰值振速最大，其最大值為2.506 7 m/s；缺陷位于C，D，E處時，缺陷中心處的峰值振速近似相等，分別為2.440 7，2.433 3和2.445 3 m/s；位于B處時，峰值振速最小，其最小值為2.388 9 m/s。當缺陷位于背爆面上的E處時，缺陷中心處有最大振速，其最大值為1.784 8 m/s，小于迎爆面上的最小振速2.388 9 m/s。

圖14 缺陷中心管道截面峰值振速分布Fig.14 Peak vibration velocity distribution map of the defect center pipe section

綜合圖12～14可知，缺陷位置對管道的動力響應(yīng)影響較小，迎爆面上的動力響應(yīng)程度明顯高于背爆面。

3.5 不同爆心距

為了定量研究管道的安全間距，基于后果最嚴重原則，研究ai=bi=40 mm，ci=0.7 t，R由4 m等距增至5.2 m共4種設(shè)計工況下管體的動力響應(yīng)特征。

圖15為不同工況下缺陷中心橫截面上位移峰值分布。由圖15可知，隨著R的增大，缺陷中心橫截面上的位移值減小，迎爆面上臨近管頂處峰值位移最大，4種工況下最大位移分別為6.583 4，3.943 2，3.693 1和3.365 9 cm。

圖16為不同工況下的峰值應(yīng)力分布。由圖16可知，缺陷中心橫截面上的峰值應(yīng)力分布呈“滴管”狀，隨著R的增大，應(yīng)力分布值減小。缺陷處存在應(yīng)力集中現(xiàn)象，且R越小應(yīng)力集中現(xiàn)象越明顯，且缺陷中心處應(yīng)力峰值最大，4種工況下的應(yīng)力峰值分別為473.27，504.54，516.58和627.32 MPa。截面上各點的應(yīng)力分布關(guān)于水平線非對稱，近爆點管體水平線下方120°范圍內(nèi)等效應(yīng)力大于上方相同范圍內(nèi)的應(yīng)力，背爆面水平線上方35°范圍內(nèi)應(yīng)力大于水平線下方同等范圍內(nèi)的應(yīng)力。

圖17為不同工況下的峰值振速分布。由圖17可知，缺陷中心橫截面上的峰值振速分布呈“草莓”狀，水平線下方曲線值略大于水平線上方曲線值，振速分布關(guān)于水平線近似對稱。隨著R的增大，峰值振速減小，減速由急劇到平緩。缺陷中心處的峰值振速最大，4種工況下最大振速分別為1.174 6，1.478 0，1.869 8和3.534 7 m/s。

圖17 缺陷中心管道截面峰值振速分布Fig.17 Peak vibration velocity distribution map of the defect center pipe section

由上述可知，爆心距R是影響管道動力響應(yīng)的關(guān)鍵因素之一，缺陷中心為風(fēng)險點。將不同R下缺陷中心處的峰值應(yīng)力進行擬合，擬合曲線如圖18所示，擬合公式為:

y=33 408.62-208.478 67R+
0.440 3R2-3.099 74×10-4R3

(10)

圖18 不同爆心距下的峰值應(yīng)力擬合曲線Fig.18 Peak stress fitting curve under different blasting distances

式中：y為缺陷中心處的峰值應(yīng)力，MPa；R為爆心距，cm。為了保證管道安全平穩(wěn)運行，利用式(10)求得基于彈性失效準則的管道安全距離為5.16 m。

4 結(jié)論

1)采用SPH-FEM耦合法計算得到不同藥包埋深時的爆腔形態(tài)與爆腔預(yù)估尺寸較吻合，且峰值振速與比例距離成線性關(guān)系，與振速經(jīng)驗公式規(guī)律相吻合，驗證了SPH-FEM算法的可靠性，為進一步開展埋地結(jié)構(gòu)近場爆炸沖擊動力響應(yīng)提供新思路。

2)爆炸沖擊下，腐蝕管道動力響應(yīng)的熱點區(qū)域位于管體迎爆面上的缺陷處，且最大動力響應(yīng)參數(shù)值位于缺陷中心點處，隨著缺陷深度的增加或管間距的減小參數(shù)值增速由平緩到急劇。

3)相比缺陷深度和爆心距對管道的擾動程度，缺陷表面尺寸和缺陷位置對管道的動力響應(yīng)影響較小。通過改善腐蝕防護工作和增大管間距，可有效保證管道安全平穩(wěn)運行。

4)在本研究的條件下，為防止管道的物理爆破對并行缺陷鄰管的沖擊破壞，建議管道的最大腐蝕深度不宜大于0.633 6倍的管道壁厚，并行間距不應(yīng)小于5.16 m。