樊凱旋 彭 旭 任家帆 饒國寧

南京理工大學化工學院(江蘇南京，210094)

引言

近年來，隨著我國航天航空技術的不斷快速發(fā)展，對于航天航空產品作業(yè)過程中的安全要求越來越高。 作為航天航空產品的動力裝置，固體火箭發(fā)動機是不利用周圍介質，而是利用自身攜帶的固體推進劑來生成工質的直接反作用發(fā)動機。 固體推進劑是化學推進劑的一種，它利用化學反應所釋放出的熱能作為推進系統的熱源[1]。 隨著現代固體推進劑能量的提高和火炸藥技術的融合，固體推進劑中會含有高能炸藥成分，在外界刺激作用下一旦發(fā)生爆炸，將會導致災難性事故的發(fā)生[2-4]。

從20 世紀60 年代起，國外對固體火箭發(fā)動機的意外爆炸及其后果開展了深入的研究[5-7]，不僅進行了大量的模擬仿真和理論分析工作，還進行過多次的固體發(fā)動機爆炸試驗，積累了大量的試驗數據，并且依此制定了安全防范措施。目前，國內的研究工作主要集中在幾個方面：1)針對固體發(fā)動機在外來沖擊下是否安全[8]，并未涉及到爆炸所產生的危害；2)固體火箭發(fā)動機自毀的危險性研究[9-10]；3)高能固體火箭發(fā)動機爆炸沖擊波毀傷效應研究[11]，這里只是涉及在無遮蔽的情況下固體發(fā)動機爆炸沖擊波對人的殺傷作用，并沒有考慮到固體發(fā)動機在工房中爆炸的情況。 而針對航天企業(yè)裝配作業(yè)時工房內的固體發(fā)動機爆炸事故的模擬和后果評估方面鮮有研究報道。

本文中，以航天產品裝配工房作業(yè)為研究對象，采用顯式動力學軟件AUTODYN 模擬固體發(fā)動機作業(yè)過程中的意外爆炸，分析工房內的沖擊波流場分布規(guī)律以及對人員傷害的影響，研究結果可以為裝配作業(yè)的風險評估、安全對策與防范提供理論支撐。

1 有限元模型的建立

以某型號固體火箭發(fā)動機裝配作業(yè)為研究對象，只考慮沖擊波對人的傷害作用，所以可以將固體發(fā)動機等效為相應質量的TNT。 有限元模型如圖1所示。

1.1 材料方程及參數的確定

1)空氣。 空氣采用Ideal Gas 理想氣體狀態(tài)方程，其表達式為

式中：p為壓力；ρ為空氣密度，ρ ＝1.225 kg/m3；e為氣體單位質量內能，e ＝206.8 kJ/g；γ為絕熱指數，γ ＝1.4。

2)炸藥。 根據固體發(fā)動機的TNT 當量，采用JWL 狀態(tài)方程來描述其爆轟過程，具體表達式為

式中：p為爆壓；V為爆轟產物的相對體積；E為初始比內能；A1、B1、R1、R2和ω為材料常數，其具體參數取值如表1 所示。 表1 中，ρ為TNT 的密度；D為炸藥的爆速[12-13]。

3)水泥混凝土。 CONC-35MPA 混凝土的狀態(tài)方程為p-alpha，其表達式為

表1 炸藥的JWL 狀態(tài)方程參數Tab.1 Parameters of JWL equation of state for explosives

式中：ρ0為材料初始密度；ρ為壓縮過程中材料的密度；e為初始內能；A1、A2、A3、B0、B1、T1、T2為材料常數，其具體參數取值如表2 所示[14]。

表2 CONC-35MPA 混凝土的材料方程參數Tab.2 Material equation parameters of CONC-35MPA concrete

強度模型則選擇適用于描述爆炸過程大變形、高應變的RHT Concrete 本構模型。 RHT 本構模型引入了3 個極限面， 即彈性極限面、失效面和殘余強度面，分別描述混凝土的初始屈服強度、失效強度及殘余強度的變化規(guī)律。 失效面的方程為

表3 RHT 本構模型強度參數Tab.3 Strength parameters of RHT constitutive model

1.2 數值計算模型

AUTODYN 里建立的模型為空氣模型、產品模型與工房模型。 氣域為長方體結構，將產品和工房包括在內，將三者用完全耦合的方式連接，空氣域的邊界條件為Flow-out (無限域、無反射邊界條件)。

1.2.1 空氣模型

空氣模型采用Euler，3D-Multi-Material 算法，空氣域為100 m×100 m×15 m 的長方體結構，采用六面體結構化網格，將網格設置為漸變網格，最小網格尺寸為50 mm。 網格劃分如圖2 所示。

1.2.2 產品模型

根據《軍工燃燒爆炸品工程安全技術規(guī)范征求意見稿2019》可以得出，此型號固體火箭發(fā)動機的TNT 當量為0.4[17]，該數值是通過實驗得出的。

此產品的TNT 當量比為0. 4，TNT 的密度為1.63 g/cm3，密度體積公式為

根據式(7)可以得到TNT 的尺寸為0.55 m ×0.40 m ×0.40 m。 TNT 采用Euler 算法，單發(fā)產品炸點為炸藥的中心，雙發(fā)產品的炸點有兩個，分別布置在兩個炸藥的中心，起爆時間為同時起爆，雙發(fā)產品間距為1 m，炸藥嵌在空氣域里。 炸藥采用的是六面體結構化網格，網格尺寸為50 mm，單發(fā)和雙發(fā)產品網格劃分如圖3 和圖4 所示。

1.2.3 工房模型

工房為長方體狀，尺寸為30 m×15 m×9 m，工房采用六面體結構化網格，網格尺寸為200 mm。 單發(fā)產品放置在工房地面的中心處，雙發(fā)產品在工房地面的中軸線上并排放置。 網格劃分如圖5 所示。

為了定性和定量地分析工房內部壓力分布規(guī)律，在空氣域內設置了一系列監(jiān)測點，監(jiān)測點的位置均為距工房地面1 m 高。 為了便于比較，設置單發(fā)產品和雙發(fā)產品爆炸時監(jiān)測點的位置相同，如圖6所示。

2 數值模擬結果分析

產品在工房地面爆炸時，形成的爆炸沖擊波向四周傳播，當遇到壁面時反射形成反射沖擊波，這些反射波及入射波相互疊加，使空間內的壓力增大。TNT 在密閉空間爆轟后，能量的釋放主要分為兩個過程：爆轟和后燃燒過程。 爆轟形成了沖擊波壓力，后燃燒形成了密閉空間中的準靜態(tài)壓力，高頻沖擊波壓力波與低頻準靜態(tài)壓力波疊加[18]。 這里只研究爆轟產生的沖擊波超壓的影響，以此來進行定性和定量分析。

2.1 不同時刻的壓力分布

單發(fā)產品爆炸過程中，不同時刻工房的壓力云如圖7 所示。

當t＝0 ms 時，由于爆炸還沒有開始，工房處在空氣域內，所以此時工房內的壓力應該與環(huán)境壓力相同，壓力為101 kPa。

當t為10 ms時，窗口處的壓力大于工房內其他位置的壓力，窗口處的最大壓力為1 018 kPa。這是由于：炸藥爆炸沖擊波以球形的方式向外傳播，炸藥距離窗口比距離墻面更近，所以沖擊波到達得更早，窗口的壓力要比其余位置大；泄壓時，工房內部的壓力從窗口泄放，就會導致窗口附近的壓力急劇升高。

當t為20 ms 時，工房頂部的最大壓力可以達到2 800 kPa，工房內大部分區(qū)域的壓力為負值。

當t為30 ms 時，工房內的壓力值維持在450 kPa 左右，工房頂部的壓力開始下降。

當t為40 ms 時，工房內大部分區(qū)域壓力在620 kPa 左右，最大壓力為5 600 kPa。

t從50 ms 到70 ms，窗口處的壓力為工房內的最大壓力，最大壓力從5 700 kPa 上升到8 000 kPa，工房內出現負壓區(qū)域。

雙發(fā)產品爆炸過程中，不同時刻工房的壓力云如圖8 所示。

模擬未開始時，整個工房的壓力與環(huán)境壓力相同；當t為10 ms 時，工房頂部和窗口處的壓力率先到達最大值，最大壓力為2 600 kPa；當t為20 ms時，工房頂部的最大壓力可以達到3 380 kPa，工房內大部分區(qū)域的壓力為負值；當t為30 ms 時，工房內的壓力為580 kPa，工房頂部的壓力開始下降；當t為40 ms 時，窗口處的壓力依然是工房內的最大壓力，最大壓力為9 100 kPa；t從50 ms 到70 ms，窗口處的壓力為工房內的最大壓力，最大壓力從10 730 kPa 上升到13 200 kPa。

2.2 沖擊波超壓與傷害半徑

在整個空氣域內要對沖擊波超壓進行定量分析，就必須在空氣域的不同位置添加監(jiān)測點，得到實時壓力變化情況，根據壓力變化曲線來判斷沖擊波超壓的大小，再根據超壓對人的傷害作用，從而求出傷害半徑。 考慮窗口對壓力的變化會有影響，所以分兩種情況來對單發(fā)產品和雙發(fā)產品爆炸的傷害范圍進行分析。 其中，X方向與窗口平行，Y方向與窗垂直。 查閱文獻得到沖擊波超壓對人的傷害作用如表4 所示[19]。

2.2.1X方向上的壓力分析

表4 沖擊波超壓對人的傷害作用Tab.4 Damage effect of shock wave overpressure on humans

X方向上監(jiān)測點的分布為：沿著工房的中心向X軸分布，所有的點均在與X軸平行的直線上。 測點的位置如表5 所示。

表5 X 方向上測點的位置Tab.5 Location of points on the X direction

各測點的壓力曲線見圖9、圖10。

可以看出，單發(fā)產品爆炸時，測點1?！?＃的壓力均從101 kPa開始變化，不同位置的沖擊波超壓到達時間也不相同，到炸點的距離越遠，超壓到達的時間也越大。 并且隨著到炸點距離的增大，沖擊波超壓也越來越小。 這也符合沖擊波超壓作用的原則。

雙發(fā)產品爆炸時，測點1?！?＃的壓力也是從環(huán)境壓力開始變化的，壓力變化趨勢和單發(fā)產品爆炸時的壓力變化基本相同，但是超壓峰值的增長與爆炸產品數量的增長呈現非線性關系。 根據表4 沖擊波超壓對人的傷害作用，得到單發(fā)產品和雙發(fā)產品爆炸的X方向上的傷害范圍，如表6 所示。

表6 X 方向上的傷害半徑Tab.6 Damage radius in the X direction m

所以，單發(fā)產品爆炸在X方向上無傷害半徑為40.0 m 之外；輕傷半徑＞29.0 ～40.0 m；重傷半徑20.5 ～29.0 m；死亡半徑為20.5 m 之內。 雙發(fā)產品爆炸在X方向上無傷害半徑為45.0 m 之外；輕傷半徑＞32.0 ～45.0 m；重傷半徑24.0 ～32.0 m；死亡半徑為24.0 m 之內。

2.2.2Y方向上的壓力分析

Y方向上監(jiān)測點的分布為：沿著工房的中心向Y軸分布，所有測點均在與Y軸平行的直線上。 測點的位置如表7 所示。 各測點的壓力曲線如圖11、圖12 所示。

表7 Y 方向上點的位置Tab.7 Location of the point in the Y direction

可以看出，單發(fā)產品爆炸時，測點7?！?2＃的壓力也是從環(huán)境壓力開始變化。 為了更好地比較各個點的沖擊波超壓的差異以及壓力變化規(guī)律，這里將測點6＃的壓力曲線舍去，因為測點6＃的沖擊波超壓過大，如果放在同一個坐標系中，很難根據壓力曲線對其他點做出比較。

圖11 反映出隨著到炸點的距離變遠，沖擊波超壓到達的時間隨之增大，沖擊波超壓也在逐漸減小，同樣符合沖擊波超壓作用的原則。

雙發(fā)產品爆炸時，測點7?！?2＃壓力增加的幅度明顯要比測點1?！?＃的壓力增加幅度大，這是因為測點7＃～12＃分布于兩個炸藥連線的中垂線上，由于沖擊波是以球形的方式向外擴散，所以測點7＃～12＃受到兩股沖擊波同時作用。 根據表4 沖擊波超壓對人的傷害作用，就可以得出在Y方向上的傷害半徑，如表8 所示。

表8 Y 方向上的傷害半徑Tab.8 Damage radius in the Y direction m

所以單發(fā)產品爆炸在Y方向上無傷害半徑為43.0 m 之外，輕傷半徑＞31.0 ～43.0 m；重傷半徑22.0 ～31.0 m；死亡半徑為22.0 m 之內。 雙發(fā)產品爆炸在Y方向上無傷害半徑為55.0 m 之外，輕傷半徑＞42.5 ～55.0 m；重傷半徑32.0 ～42.5 m；死亡半徑在32.0 m 之內。

2.3 模擬結果與實驗結果的比較

根據陽建紅等[20]做的無遮蔽情況下高能推進劑爆炸實驗，可以將實驗結果與模擬結果作比較，從而印證上述模擬方法是正確的。

具體做法為：由實驗中推進劑的質量和推進劑TNT 當量可以計算出TNT 質量。 由于雙發(fā)產品之間的沖擊波會相互影響，這里只對單發(fā)產品進行分析。 將文章中有限元模型的工房部分刪除，產品質量與實驗中換算得出的TNT 質量相同，然后再增加測點，模擬測點位置與試驗測點位置相同，其余條件保持不變。 可以得出測點處沖擊波超壓的數據，如表9 所示。

從表9 可以看出，在距離爆源大于6 m 時，數值模擬結果與實驗所測數據誤差范圍在10%以內；距離爆源越近，爆炸場內沖擊波傳播越復雜，馬赫反射的影響越顯著，從而造成了模擬數據與實驗數據相差較大。 通過實驗結果與模擬結果相比較，說明文中所建立的計算模型與方法是正確的。

3 結論

通過數值模擬的手段研究了航天產品裝配工房內單發(fā)產品和雙發(fā)產品的爆炸過程，得出如下結論：

表9 實驗結果與模擬結果相比較Tab.9 Result comparison of experiment with simulation

1)根據不同時刻的壓力云圖可以得到，產品在工房內爆炸要比在空氣中爆炸復雜得多，超壓峰值的大小與工房的結構有很大的關系。

2)當單發(fā)產品爆炸時，沖擊波超壓在X方向上對人的殺傷距離為29.0 m，爆炸后人員的安全距離為40.0 m；在Y方向上對人的殺傷距離為31.0 m，爆炸后人員的安全距離為43.0 m。 當雙發(fā)產品爆炸時，沖擊波超壓在X方向上對人的殺傷距離為32.0 m，爆炸后人員的安全距離為45.0 m；在Y方向上對人的殺傷距離為42.5 m，爆炸后人員的安全距離為55.0 m。