復(fù)雜結(jié)構(gòu)內(nèi)爆炸沖擊波傳播的大規(guī)模計算方法研究

王成，王雪桃，李濤，張志剛，孫偉福

（1.北京理工大學(xué) 爆炸科學(xué)與技術(shù)國家重點實驗室，北京 100081；2.中煤科工集團重慶研究院有限公司，重慶 400039；3.煤礦災(zāi)害防控全國重點實驗室，重慶 400039）

爆炸沖擊波是礦井等復(fù)雜密閉建筑結(jié)構(gòu)發(fā)生爆炸事故產(chǎn)生的主要災(zāi)害效應(yīng)之一，沖擊波不僅會造成人員傷亡，還會對周圍結(jié)構(gòu)及建筑物造成大規(guī)模破壞.空中自由場等開放空間發(fā)生爆炸事故時，由于爆炸沖擊波未受任何約束，通常以指數(shù)形式快速衰減.然而，礦井巷道等密閉空間內(nèi)爆炸產(chǎn)生的沖擊波由于受到壁面等結(jié)構(gòu)的約束作用，在其擴散過程中會發(fā)生多次的反射、折射等效應(yīng)，并隨之發(fā)生多次疊加作用，使得沖擊波超壓峰值較多且難以衰減至常壓，在很大程度上提高了爆炸壓力和能量.因此，爆炸沖擊波在復(fù)雜地下密閉空間內(nèi)的傳播規(guī)律是礦井、巷道等制定防爆及抗爆安全的主要科學(xué)依據(jù)，在爆炸事故的快速應(yīng)急響應(yīng)及防護救援等過程中具有極其重要的理論意義和實際應(yīng)用價值.

關(guān)于類似巷道、礦井等復(fù)雜密閉結(jié)構(gòu)內(nèi)的爆炸問題，陳海天等[1]、龐偉賓等[2-3]及其他研究者[4-6]通過實驗對其進行了理論研究，并基于自由場爆炸沖擊波壓強預(yù)測公式，得到了直邊墻拱形隧道、方形坑道、T 型坑道以及不同分叉角度的坑道內(nèi)爆炸沖擊波的衰減規(guī)律.但由于巷道在彎曲處和交叉處的沖擊波流場分布極其復(fù)雜，理論研究無法對其傳播過程進行直觀的反映，所以僅通過實驗對巷道內(nèi)爆問題進行深入細致的研究存在許多實際困難.因此，近些年來數(shù)值模擬技術(shù)在結(jié)構(gòu)內(nèi)爆問題中得到了廣泛應(yīng)用.張奇[7]、屈康康等[8]及其他研究者們[9-18]利用Autodyn、LS-DYNA 等有限元軟件對長直坑道、L 形坑道、S 形坑道、分叉管道、變徑管道等不同形狀巷道內(nèi)的沖擊波傳播規(guī)律進行了研究.樊保龍等[19]、王成等[20]對密閉巷道內(nèi)瓦斯的爆炸過程進行了數(shù)值分析.胡洋等[21]、宋詩謙等[22]、張曉偉等[23]及吳文溢等[24]對密閉空間內(nèi)TNT 爆炸效應(yīng)進行了數(shù)值模擬分析.張學(xué)博等[25]提出了沖擊波傳播的分段接力數(shù)值模擬法，并實現(xiàn)了礦井大尺度爆炸問題的二維模型模擬沖擊波傳播.楊科之等[26]通過自主計算程序?qū)唵蜗锏纼?nèi)爆流場進行了數(shù)值計算，并得出沖擊波在巷道內(nèi)衰減的計算公式，但是Autodyn 等商業(yè)軟件無法對炸藥放置于地面或近壁面的爆炸過程進行高效計算，也無法對近、中場的爆炸問題進行準確計算[27].另外，目前現(xiàn)有的自主編制爆炸流場的程序更多地是對簡單、小尺寸的爆炸問題進行計算，若對大尺寸的復(fù)雜密閉空間結(jié)構(gòu)進行數(shù)值模擬，則存在一定的難度和局限性.對于類似于地下礦井等大尺度的復(fù)雜密閉空間結(jié)構(gòu)，若在計算時整體計算域采用較粗的網(wǎng)格，則會導(dǎo)致計算結(jié)果不準確；若全流場采用較為精細的網(wǎng)格尺寸，則會導(dǎo)致由于網(wǎng)格數(shù)目過多而無法正常計算，即使利用高性能的計算設(shè)備和大規(guī)模的并行處理器，也會浪費大量的內(nèi)存和計算資源.

針對大尺寸復(fù)雜結(jié)構(gòu)內(nèi)爆數(shù)值計算浪費大量計算資源的問題，在此提出了一種基于數(shù)據(jù)流的分塊續(xù)算方法，并通過文獻實驗結(jié)果以及標準模型算例驗證了該計算方法的準確性和有效性，實現(xiàn)了復(fù)雜地下密閉空間內(nèi)爆問題的大規(guī)模、高效率、高精度數(shù)值模擬.同時，以復(fù)雜地下密閉空間結(jié)構(gòu)中的爆炸問題為研究對象，采用基于數(shù)據(jù)流的分塊續(xù)算方法對一個包含長方體形工房、圓柱形泄壓洞的復(fù)雜結(jié)構(gòu)內(nèi)的爆炸沖擊波傳播過程進行了數(shù)值模擬，并給出了沖擊波在巷道內(nèi)的傳播規(guī)律.

1 數(shù)學(xué)模型及算法

1.1 數(shù)學(xué)模型

由于炸藥在爆轟傳播過程中，沖擊波的傳播速度很快，因此分子擴散、熱傳導(dǎo)及黏性效應(yīng)在模型計算中可以忽略.控制方程組為由質(zhì)量守恒、動量守恒和能量守恒組成的非定常可壓縮三維Euler 方程組：

式中： ρ為密度；p為壓強；u、v、w為速度矢量分量；E為單位質(zhì)量的總能量，是單位質(zhì)量的比內(nèi)能e與動能的總和，滿足

為了使上述控制方程組封閉，需引入描述壓強、密度和比內(nèi)能之間關(guān)系的狀態(tài)方程.爆轟產(chǎn)物采用JWL 狀態(tài)方程，空氣介質(zhì)采用理想氣體狀態(tài)方程.

JWL 狀態(tài)方程為：

式中：p為壓強； ρ為密度；e為單位質(zhì)量的比內(nèi)能；A、B、R1、R2、 ω 和 ρ0為爆轟產(chǎn)物的材料常數(shù).

理想氣體狀態(tài)方程為：

式中：p為壓強； ρ為密度；e為比內(nèi)能； γ為多方氣體指數(shù)，對于空氣 γ一般取為1.4.

1.2 數(shù)值計算方法

在復(fù)雜結(jié)構(gòu)內(nèi)爆炸沖擊波流場的數(shù)值模擬過程中，為確保沖擊波的追蹤精度，采用5 階WENO（weighted essentially non-oscillatory）高精度有限差分格式[28]進行空間離散，采用3 階精度的TVD Runge-Kutta 方法進行時間離散.針對多介質(zhì)爆炸流場問題，采用RGFM（real ghost fluid method）[29]和 Level-set[30-31]相結(jié)合的物質(zhì)界面處理方法，通過Level-set 方法對界面進行隱式追蹤，通過求解法向黎曼問題確定物質(zhì)界面附近流場相互接觸作用之后的流動狀態(tài)，并利用RGFM 方法對界面處流場進行處理，實現(xiàn)了復(fù)雜地下密閉空間結(jié)構(gòu)內(nèi)爆轟產(chǎn)物及爆炸沖擊波傳播過程的高精度數(shù)值模擬.另外，通過采用一種基于WENO格式的自適應(yīng)擴大計算域高精度算法[32]，提高了多介質(zhì)界面追蹤與處理的精度和分辨率，加快了計算效率.在此基礎(chǔ)上，引入MPI（message passing interface）并行計算方法，對自適應(yīng)擴大計算域進行并行計算，即可實現(xiàn)高精度大規(guī)模爆炸問題的數(shù)值模擬.

在該大規(guī)模計算方法中，復(fù)雜結(jié)構(gòu)內(nèi)爆炸問題的數(shù)值模型一般包括兩部分區(qū)域：計算域和障礙物區(qū)域Ωobstacle，如圖1 所示.其中，在計算域中設(shè)置炸藥區(qū)域ΩTNT、空氣域Ωair兩種介質(zhì)以及剛性壁面區(qū)域，爆轟產(chǎn)物采用JWL 狀態(tài)方程計算，空氣采用理想氣體狀態(tài)方程.障礙物區(qū)域不參與計算，但由于礦井、巷道等地下空間的壁面一般為鋼筋混凝土墻，所以在障礙物與空氣/爆轟產(chǎn)物的交界處給定固壁無滑移邊界條件，用于模擬剛性壁面的反射過程.

2 分塊續(xù)算方法

由于地下空間或礦井等大規(guī)模復(fù)雜結(jié)構(gòu)爆炸問題的數(shù)值研究都是百米、千米量級，而炸藥爆炸形成的爆轟波陣面厚度通常在微米量級；另外，三維的爆炸問題存在時間和空間上的多尺度特性，若計算區(qū)域采用粗網(wǎng)格劃分會帶來較大的累積誤差，導(dǎo)致計算結(jié)果不準確甚至不收斂；若采用精細網(wǎng)格則只能對小尺寸的簡單結(jié)構(gòu)進行模擬，大規(guī)模復(fù)雜密閉空間結(jié)構(gòu)由于網(wǎng)格數(shù)目巨大，即使采用并行技術(shù)也很難滿足對其全尺寸快速計算的需求，對其進行全尺寸網(wǎng)格劃分更會浪費大量的CPU 計算資源.以圖2所示的二維結(jié)構(gòu)的爆炸問題為例，若對其進行全尺寸計算時需要72 個CPU 并行計算，且其中有36 個處于障礙物區(qū)域不參與計算（圖中陰影區(qū)域所示），極大地浪費了計算資源，如圖2（a）所示.為了減少大量計算資源的浪費，提高計算效率，實現(xiàn)對大規(guī)模復(fù)雜結(jié)構(gòu)內(nèi)爆炸問題的全流場數(shù)值模擬，提出了一種基于數(shù)據(jù)流的分塊續(xù)算方法，通過合理劃分模型并對每個部分進行單獨計算可對CPU 進行有效利用，如圖2（b）所示.對該二維結(jié)構(gòu)模型利用分塊續(xù)算方法計算時，共需要44 個CPU 但只有6 個處于障礙物區(qū)域，與全尺寸計算方法相比，節(jié)約了39%的CPU核數(shù)，有效利用了計算域，減少了計算資源的浪費.另外，若全尺寸計算需72n核時，則分塊續(xù)算共需22n核時，計算核時相比于全尺寸計算節(jié)約了69%.若以三維結(jié)構(gòu)來看，全尺寸計算浪費的CPU 計算資源更多，因此與全尺寸計算相比，數(shù)據(jù)流分塊續(xù)算方法極大地減少了計算資源的浪費，對復(fù)雜結(jié)構(gòu)進行合理劃分即可有效實現(xiàn)其內(nèi)爆流場的數(shù)值模擬.

以下簡要敘述數(shù)據(jù)流分塊續(xù)算方法的流程：

① 劃分模型.

在對礦井、巷道等復(fù)雜地下空間大尺度爆炸問題進行模擬時，可根據(jù)實際情況將計算模型進行簡化，并根據(jù)模型幾何形狀和爆炸流場的流體動力學(xué)特征將模型各部分進行合理劃分，即將整個計算模型劃分為若干部分（PART-i，i=1,2,3,···），使得在對每個部分進行數(shù)值計算時盡可能少地浪費計算資源.

② 數(shù)據(jù)儲存.

通常情況下，將放置炸藥的模型部分劃分為PART-1 對其先進行計算.首先根據(jù)具體的爆炸流場尺度以及計算資源，確定該部分的計算域尺寸、每個方向的CPU 數(shù)以及每個分區(qū)中的網(wǎng)格數(shù)，并給定計算域各CPU 的序號.另外，在PART-1 模型的出口處設(shè)置一個數(shù)據(jù)儲存區(qū)域，記為Ω，如圖2（b）所示.在計算若干時間步后，該部分模型計算域中的爆炸沖擊波向外傳播，當沖擊波傳播至PART-1 中Ω區(qū)域的80%～85%時，程序停止計算并將這一時刻Ω區(qū)域內(nèi)各CPU 所有計算網(wǎng)格點的物理量（包括密度、3個方向的速度、總能量、壓力等）保存至文件中，即數(shù)據(jù)儲存.

③ 數(shù)據(jù)賦值.

在進行下一部分即PART-2 的計算時，同樣先確定該部分的計算域尺寸，并按PART-1 的網(wǎng)格大小對其進行網(wǎng)格劃分，使得相鄰PART 間能夠進行“映射”數(shù)據(jù)流的傳遞.在開始計算前，先更新時間步，讓兩部分程序的時間同步.另外，在PART-2 中設(shè)置一塊與PART-1 完全重疊的區(qū)域Ω，在計算該部分模型時，將上一部分模型PART-1 保存的區(qū)域Ω的數(shù)據(jù)作為該部分的初始條件，在開始計算時讀入此文件并重新啟動程序，將區(qū)域Ω內(nèi)各CPU 所有網(wǎng)格點的物理量直接映射在該部分區(qū)域Ω內(nèi)相對應(yīng)的計算網(wǎng)格點上，然后繼續(xù)求解歐拉方程進行計算，即數(shù)據(jù)賦值.若進行計算的整個大規(guī)模模型被劃分為很多部分，則根據(jù)流場中爆炸沖擊波的傳播過程對所有PART依次進行數(shù)據(jù)儲存、數(shù)據(jù)賦值的數(shù)據(jù)流傳遞，即可實現(xiàn)對整個模型的爆炸流場進行數(shù)值模擬.

分塊續(xù)算方法由于在相鄰PART 間分割處都會設(shè)置重疊區(qū)域，且各部分模型之間是通過數(shù)據(jù)儲存、數(shù)據(jù)賦值的方式來進行數(shù)據(jù)流傳遞的，因此，相對于全尺寸計算，這種計算方法增加了一定的計算量.但是，相鄰PART 間的這種直接映射的數(shù)據(jù)流傳遞方式也保證了數(shù)據(jù)的完整性，區(qū)域Ω在一定程度上可以代表上一個PART 計算得到的結(jié)果，使得接續(xù)計算時產(chǎn)生的誤差可以忽略（下文對此進行了驗證）.更為重要的是，分塊續(xù)算方法解決了大尺度復(fù)雜結(jié)構(gòu)由于網(wǎng)格數(shù)目過大而無法計算的難題，更解決了計算資源的大量浪費問題，實現(xiàn)了復(fù)雜地下密閉空間內(nèi)爆炸沖擊波傳播的大規(guī)模數(shù)值模擬，在礦井、巷道等復(fù)雜地下空間的防爆抗爆設(shè)計中具有廣闊的應(yīng)用前景.

3 驗證性算例

3.1 計算方法驗證

為了驗證該計算方法與其數(shù)值計算結(jié)果的準確性，結(jié)合文獻[6]對方形坑道內(nèi)爆炸沖擊波傳播的實驗結(jié)果，設(shè)計了兩端開口的TNT 方形坑道內(nèi)爆的計算模型，并利用該計算方法對其爆炸流場中的沖擊波衰減過程進行數(shù)值計算.該計算模型示意圖如圖3所示，方形坑道截面長×寬為1.2 m×1.8 m，坑道總長為41 m.圓柱形TNT 裝藥的密度為1.57 g/cm3，質(zhì)量為3.25 kg.柱狀裝藥距離地面0.9 m，且懸掛于距離坑道端面5 m 的位置.在坑道一側(cè)壁面的爆心所在截面上分別距離爆心6、10、16、20 、28、36 m 處設(shè)置超壓測點，其布置如圖3 所示.

圖3 方形坑道數(shù)值計算測點布置示意圖Fig.3 Diagram of measuring point layout for numerical calculation of square tunnel

表1 列出了文獻[5]所得實驗值以及數(shù)值模擬所得典型測點處的超壓峰值.從表中可以看出，20 m處沖擊波超壓的數(shù)值計算結(jié)果與實驗結(jié)果基本一致，而越靠近爆炸中心的測點，超壓峰值的偏差相對越大.因為在實驗中距離炸藥越近的地方?jīng)_擊波壓力越高，傳感器受震動會產(chǎn)生較多雜信號，而且靈敏度（分辨率）不同，會影響讀數(shù).但總體而言，利用該數(shù)值計算方法預(yù)測的峰值壓力與實驗所測得的壓力值基本一致，相對誤差在9.55%以內(nèi)，證明該計算方法能夠較為準確地預(yù)測結(jié)構(gòu)內(nèi)爆問題.

表1 典型測點處的文獻實驗及數(shù)值模擬超壓峰值Tab.1 Literature experiments and numerical simulations of overpressure peaks at typical measurement points

3.2 分塊續(xù)算方法驗證

由于數(shù)據(jù)流分塊續(xù)算方法在相鄰PART 間用于“數(shù)據(jù)儲存、數(shù)據(jù)賦值”的Ω區(qū)域只保存了PART-i中爆炸沖擊波的主要流場數(shù)據(jù)，實際上Ω區(qū)域以外的爆炸流場對后續(xù)沖擊波的傳播也有一定的影響.為了研究Ω區(qū)域保存的流場數(shù)據(jù)能否作為下一部分PART-（i+1）的初始條件，實現(xiàn)PART 之間的沖擊波續(xù)算過程，設(shè)計了如圖4 所示的標準模型算例對分塊續(xù)算方法的準確性進行了驗證.

圖4 標準模型算例示意圖Fig.4 Schematic diagram of standard model calculation example

該計算模型由一個4 m×4 m×3 m 的長方體和半徑為1 m、長10 m 的管道組成，長方體內(nèi)放置一塊藥量為50 kg 的TNT 正方體狀裝藥，位于距離地面高0.5 m 處.整個計算域有兩種介質(zhì)：空氣介質(zhì)和TNT炸藥，空氣介質(zhì)采用理想氣體狀態(tài)方程，TNT 炸藥采用JWL 狀態(tài)方程.另外，在計算域長管道中心軸線方向每隔1 m 設(shè)置一個沖擊波超壓峰值測點，用于記錄沖擊波超壓在管道內(nèi)隨距離的衰減過程.

為驗證數(shù)據(jù)流分塊續(xù)算方法的準確性，對該計算模型內(nèi)的爆炸問題進行兩次數(shù)值計算即全尺寸計算和分塊續(xù)算.全尺寸計算時，計算模型不進行拆分直接對整個計算域劃分網(wǎng)格進行計算，如圖4 所示，網(wǎng)格尺寸dx=dy=dz=0.01 m，網(wǎng)格總數(shù)約為2.94 億，采用600 個CPU 并行計算，總的計算時長為19 h；分塊續(xù)算時，將整個計算模型劃分為PART-1 和PART-2兩部分，并在PART-1 和PART-2 中分別設(shè)置一個長為2 m 的區(qū)域Ω，用于接續(xù)計算時數(shù)據(jù)流的傳遞，如圖4（b）所示.續(xù)算時，網(wǎng)格尺寸不變，網(wǎng)格總數(shù)約為1.17 億，共采用240 個CPU 并行計算，總的計算時長為8.5 h 左右.與全尺寸計算方法相比，分步續(xù)算方法節(jié)約了60.00%的CPU 核數(shù)，極大地減少了計算資源的浪費；另外，計算時長相比于全尺寸計算降低了55.26%左右，即分步續(xù)算方法在給定算例下相比于普通全尺寸計算，節(jié)約了88.68%的計算核時，極大地提高了計算效率.

表2 列出了全尺寸計算和分塊續(xù)算得到的管道內(nèi)典型測點處沖擊波的超壓峰值.可見分塊續(xù)算得到的超壓峰值與全尺寸計算時的基本一致，相對誤差在0.33%以內(nèi).表明分塊續(xù)算方法可實現(xiàn)基于數(shù)據(jù)流的數(shù)據(jù)傳遞，并能夠保證數(shù)據(jù)在傳遞過程中的準確性.另外，對于復(fù)雜地下密閉空間等大規(guī)模結(jié)構(gòu)，分塊續(xù)算方法可以有效利用計算資源，節(jié)約計算時長，提高計算效率，實現(xiàn)爆炸流場中沖擊波傳播的大規(guī)模、高效率、高精度數(shù)值模擬.

表2 全尺寸計算和分塊續(xù)算時典型測點處的超壓峰值Tab.2 Peak overpressure at typical measuring points during full-scale calculation and block continuation

4 爆炸沖擊波在復(fù)雜結(jié)構(gòu)中的數(shù)值模擬

4.1 計算模型設(shè)置

針對復(fù)雜地下密閉空間結(jié)構(gòu)中的爆炸問題，采用上文提出的基于數(shù)據(jù)流的分塊續(xù)算方法，設(shè)計了一個包含長方體型工房及圓柱形泄壓洞的大規(guī)模復(fù)雜結(jié)構(gòu)對爆炸流場中的沖擊波傳播及衰減過程進行數(shù)值模擬，其計算模型示意圖如圖5 所示.

圖5 大規(guī)模復(fù)雜結(jié)構(gòu)計算模型示意圖Fig.5 Diagram of large-scale complex structure calculation model

整個計算模型由工房、泄壓洞A、B、C及泄壓洞B和C出口的自由場區(qū)域組成，工房內(nèi)放置一塊質(zhì)量為392 kg 的方形炸藥，TNT 裝藥密度為1.58 g/cm3.其中，工房長×寬×高為14 m×5 m×5 m，泄壓洞A和B的直徑均為2 m，泄壓洞C直徑為1 m.泄壓洞A長為10.5 m，泄壓洞B包括7 m 的水平直段部分和15 m 的向上傾斜45°部分，泄壓洞C長為25 m，整個模型的壁面都設(shè)置為固壁無滑移邊界條件.另外，在泄壓洞B和C的出口處分別設(shè)置一塊4 m×4 m×2 m及16 m×8 m×2 m 的自由場區(qū)域，該部分邊界條件設(shè)置為無反射出流.由于該計算模型較為復(fù)雜，為了減少計算資源的浪費，根據(jù)模型具體幾何形狀特征將其劃分為6 個部分（PART-1、PART-2、PART-3.1、PART-3.2、PART-4.1、PART-4.2），并利用基于數(shù)據(jù)流的分塊續(xù)算方法對其爆炸流場進行數(shù)值計算，如圖6 所示.在該計算模型中，用于數(shù)據(jù)儲存、數(shù)據(jù)賦值的Ω區(qū)域（即①～⑤區(qū)域）設(shè)置為1.6 m 長，整個模型各PART 間數(shù)據(jù)流的傳遞過程如圖6 所示.

圖6 數(shù)據(jù)流傳遞過程示意圖Fig.6 Diagram of data flow transmission process

4.2 計算結(jié)果及分析

4.2.1 工房內(nèi)部及工房出口處的沖擊波傳播規(guī)律

在工房內(nèi)前后左右每個側(cè)壁面的三等分點處分別設(shè)置2 個測點，其高度與炸藥中心平齊，監(jiān)測每個測點處的超壓-時間曲線.通過數(shù)值計算可以發(fā)現(xiàn)，距離炸藥最近的測點2、4 處的反射超壓高達281.1 MPa，主要是由于高壓沖擊波直接撞擊至壁面處形成反射波（如圖7（b）所示），而在撞擊壁面后沖擊波在向外傳播的同時壓力快速衰減，因此測點1、3 處的超壓峰值僅有29 MPa.

圖7 PART-1 典型時刻壓力云圖Fig.7 Pressure distributions at typical time of PART-1

在1.030 ms 時，沖擊波傳播至泄壓洞A內(nèi)，由于工房與泄壓洞A的交匯處存在T 字型拐角，致使此處產(chǎn)生極強的稀疏效應(yīng)，密度與壓強均接近于0，如圖7（c）所示；沖擊波在1.333 6 ms 時傳播至泄壓洞A右側(cè)壁面，并在1.650 ms 時刻從右側(cè)反射至管道中心軸線區(qū)域，因此使得泄壓洞A始段約0.8、1.0 m 處的壓力并不高，反而在泄壓洞A約1.2～1.6 m 處的壓力升高.在2.084、2.677 ms 時沖擊波傳播至工房右、左兩側(cè)壁面，并產(chǎn)生較高壓力的反射波，分別高達54、45 MPa.

4.2.2 泄壓洞A內(nèi)部沖擊波衰減規(guī)律

圖8 給出了泄壓洞A內(nèi)部沖擊波超壓峰值隨距離的衰減規(guī)律.在泄壓洞A入口0 m 處，管道內(nèi)中心軸線的超壓峰值為9.78 MPa，隨后爆炸沖擊波以波動衰減的趨勢逐漸減小，經(jīng)過約10 m 長的泄壓洞A后，沖擊波波動衰減至約1.2 MPa 左右，與入口處相比降低了約88%.因為沖擊波從工房傳播出來時并不是平面波，從而導(dǎo)致沖擊波在泄壓洞內(nèi)的壁面上來回反射，因此在泄壓洞A內(nèi)部的沖擊波并不像空中自由場爆炸一樣隨著爆距的增加而快速衰減，而是以波動衰減的趨勢逐漸降低，其波動衰減過程如圖9（a）～9（d）所示.

圖8 泄壓洞A 內(nèi)超壓峰值衰減過程Fig.8 Attenuation process of overpressure peak in pressure relief tunnel A

圖9 PART-2 典型時刻沖擊波壓力云圖Fig.9 Pressure distributions at typical time of PART-2

4.2.3 泄壓洞A、B、C交匯處T 型通道沖擊波傳播規(guī)律

由圖9 可以看出，沖擊波在約7.424 0 ms 時刻傳播至泄壓洞A側(cè)R15 圓弧處，在9.268 4 ms 時撞擊至泄壓洞B側(cè)R15 圓弧處，并產(chǎn)生較高反射壓力，此后，該反射沖擊波向泄壓洞A、B和泄壓洞C內(nèi)傳播，使得T 型區(qū)域的超壓峰值有一定程度的升高.

由于沖擊波在T 型區(qū)域處進行無規(guī)則反射，致使該區(qū)域超壓較不穩(wěn)定，為研究沖擊波在T 型通道處的分流情況，此處取距離三管道中心軸線交點2 m 處的測點進行沖擊波分流分析.通過讀取相關(guān)測點處的超壓峰值可知，泄壓洞A內(nèi)距離交點2 m處的超壓峰值為1.39 MPa，泄壓洞B的超壓峰值為0.93 MPa，泄壓洞C的超壓峰值為0.89 MPa.即在T型通道分流處，與泄壓洞A同等直徑的泄壓洞B的沖擊波超壓峰值約占分流前的67%，而呈直角且直徑只有泄壓洞A1/2 的泄壓洞C的超壓峰值約占分流前的64%.由此可見，分流時管道間采用直角交匯和縮小管道直徑在一定程度上有利于沖擊波超壓的衰減.

4.2.4 泄壓洞B和C內(nèi)部沖擊波衰減規(guī)律

從T 型區(qū)域開始，隨著沖擊波向外傳播，泄壓洞C內(nèi)的沖擊波超壓值出現(xiàn)了明顯減小，其衰減過程如圖10（b）所示.沖擊波超壓峰值從泄壓洞C內(nèi)0 m的1.22 MPa 衰減至24 m 處的0.19 MPa，降低了84%，且其衰減速度逐漸減慢.由于泄壓洞C前端部分位于T 型區(qū)域內(nèi)，因此沖擊波在0～3 m 范圍內(nèi)有一定的反射，超壓峰值有所升高，在2.6 m 處達到最高值1.8 MPa.距離16 m 處開始波動衰減消失，超壓峰值以緩慢的速度逐漸降低，說明在泄壓洞C內(nèi)16 m 左右, 沖擊波經(jīng)過多次壁面反射后已基本形成平面波.

圖10 泄壓洞B 和C 內(nèi)部超壓峰值衰減規(guī)律Fig.10 Attenuation law of overpressure peak inside relief tunnels B and C

沖擊波于T 型通道分流后，在泄壓洞B內(nèi)經(jīng)過7 m的水平直段及15 m 的向上傾斜45°部分的傳播，其在斜段部分的衰減過程如圖10（a）所示.爆炸沖擊波的超壓峰值從斜坡起始處的0.82 MPa 衰減至15 m處的0.20 MPa，降低了76%，而且其衰減速度逐漸減慢.同樣，沖擊波在斜坡管道內(nèi)傳播時，其超壓峰值也是一個波動衰減的過程.

4.2.5 泄壓洞B和C出口自由場區(qū)域沖擊波衰減規(guī)律

泄壓洞B和C外部自由場區(qū)域沖擊波超壓峰值測點的分布示意圖如圖11（a）和11（b）所示，其所對應(yīng)的超壓峰值分別如圖11（c）和11（d）所示，可以看出，與管道內(nèi)沖擊波衰減規(guī)律不同，沖擊波在自由場區(qū)域是持續(xù)衰減的，且其衰減速度比較快.爆炸沖擊波從泄壓內(nèi)傳出后，分別于1 m 和5 m 處“整合”為平面波.

圖11 泄壓洞出口測點分布及超壓衰減過程Fig.11 Distribution of measuring points and the process of overpressure attenuation

4.3 Ω 區(qū)域取值對數(shù)值計算結(jié)果的影響

上述采用基于數(shù)據(jù)流的分塊續(xù)算方法對爆炸沖擊波在復(fù)雜結(jié)構(gòu)內(nèi)的傳播過程進行了數(shù)值計算，為了在一定程度上能夠減少計算量，用于數(shù)據(jù)儲存、數(shù)據(jù)賦予的Ω區(qū)域不宜取得過大.上述計算模型中用于數(shù)據(jù)傳遞的Ω區(qū)域取值為1.6 m，在此通過對Ω區(qū)域進行適當?shù)臄U大，研究Ω區(qū)域的取值對分塊續(xù)算結(jié)果的影響大小，進而驗證數(shù)據(jù)流分塊續(xù)算方法的準確性.

選用PART-3.1 和PART-3.2 兩部分模型，對Ω區(qū)域的取值對計算結(jié)果的影響進行研究.用于數(shù)據(jù)傳遞的Ω區(qū)域擴大為原來的3 倍，即在計算PART-3.1 時，保存所有物理量的區(qū)域③由原來的1.6 m 擴大為4.8 m，在 PART-3.2 中也設(shè)置一塊與PART-3.1中區(qū)域③重疊的區(qū)域③.圖12 給出了Ω區(qū)域為1.6 m 和4.8 m 時，泄壓洞C出口區(qū)域超壓峰值衰減曲線比較圖.由圖中對比結(jié)果可以看出，沖擊波在軸線方向、左斜及右斜方向的超壓峰值衰減過程基本一致，因此，對上述大規(guī)模復(fù)雜結(jié)構(gòu)進行分塊續(xù)算時，進行數(shù)據(jù)傳遞的Ω區(qū)域取為1.6 m 是合理的，對整體爆炸流場的數(shù)值計算結(jié)果影響極小.

圖12 Ω 區(qū)域為1.6 m 和4.8 m 時泄壓洞C 出口區(qū)域超壓峰值衰減過程比較Fig.12 Comparison of peak overpressure at outlet of pressure relief tunnel C when Ω is 1.6 m and 4.8 m

5 結(jié) 論

①提出了一種基于數(shù)據(jù)流的可計算復(fù)雜地下密閉空間內(nèi)爆炸問題的分塊續(xù)算方法，解決了計算資源的大量浪費問題，實現(xiàn)了復(fù)雜地下密閉空間內(nèi)爆炸沖擊波傳播的大規(guī)模、高效率、高精度數(shù)值模擬，并利用典型實驗與標準模型算例驗證了該方法的準確性與有效性.

②利用該方法對一個包含長方體型工房和圓柱形泄壓洞的大規(guī)模復(fù)雜結(jié)構(gòu)內(nèi)的爆炸問題進行了數(shù)值模擬，結(jié)果表明：沖擊波在泄壓洞內(nèi)是以“波動”衰減的趨勢逐漸減小的，且其衰減速度比在自由場的衰減慢；T 型分流處管道間采用直角交匯或縮小管道直徑在一定程度上有利于沖擊波超壓的衰減.該方法在地下空間等大尺度復(fù)雜結(jié)構(gòu)的防爆抗爆設(shè)計中具有廣泛的應(yīng)用前景.

③在此提出的數(shù)據(jù)流分塊續(xù)算方法雖解決了計算資源的大量浪費問題，但是用于數(shù)據(jù)流傳遞的Ω區(qū)域只保存了爆炸沖擊波的主要流場數(shù)據(jù)，實際上Ω區(qū)域以外的爆炸流場對后續(xù)沖擊波的傳播也有一定的影響.然而以上算例中將Ω區(qū)域的邊界條件一律設(shè)置為固壁無滑移邊界條件進行計算，在下一步研究中，需要優(yōu)化對其邊界條件的處理，使其更接近于真實流場中爆炸沖擊波的傳播.