一種沖擊波作用下結構毀傷算法研究*

周 朗，徐春光

（中山大學航空航天學院，廣東 深圳 518107）

對于常規彈藥而言，爆炸沖擊波是造成結構破壞的主要因素，研究爆炸沖擊波傳播與演化規律，對準確評估彈藥的毀傷效能具有重要意義。針對自由大氣內爆炸沖擊波的傳播規律，研究人員構建了多個經驗公式[1-4]，并應用于工程實際中，是研究爆炸沖擊波傳播與演化的重要手段，但對于建筑物等復雜空間，經驗公式不能反映沖擊波的演化過程，與實際物理過程存在比較明顯的差異，且無法給出沖擊波演化的細節，也不能反映沖擊波與結構的相互作用過程及影響。

隨著數值模擬技術和計算機技術的發展，研究人員利用數值模擬的方法對沖擊波與障礙物存在相互作用的問題進行研究。劉偉等[5]在考慮地面因素的影響條件下，研究了不同當量TNT 裝藥空爆下沖擊波超壓規律，獲得了有意義的結果。寧建國等[6]基于JWL 狀態方程和瞬時爆轟模型，研究了存在防護擋墻情況下，爆炸流場的傳播及演化規律，所得超壓計算結果比已有經驗公式計算結果偏高。王飛等[7]采用ALE 方法，對沖擊波繞過障礙物后的超壓分布情況進行了數值模擬，數值模擬結果能較好的反映超壓變化的規律。張曉偉等[8]對爆炸沖擊波作用于隔爆墻的繞射進行實驗和數值仿真，得到墻后沖擊波超壓分布規律，并給出了墻后近地面超壓峰值的工程計算公式。文獻[9-12]對沖擊波在坑道、隧道內的傳播規律進行了研究。上述工作對沖擊波與結構的相互作用進行了較為詳細的研究，但研究中均假設結構或障礙物為剛體，沒有考慮結構的變形和運動，難以描述結構在沖擊波作用下的變形及破壞。

對于考慮沖擊波作用下結構變形、運動及其對沖擊波傳播影響的問題，多采用國外軟件進行數值模擬，譬如LS-DYNA、AUTODYN、ABAQUS 等軟件。鄧榮兵等[13]利用LS-DYNA 對爆炸沖擊波與幕墻結構的作用過程進行了研究，分析了幕墻在沖擊波作用下的動態響應特性和內外層玻璃的破壞情況，與實驗結果符合較好。董湘乾[14]利用ABAQUS 對爆炸沖擊波與混凝土框架結構的作用過程進行了模擬，分析了多層框架結構在沖擊波作用下的動態響應。張秀華等[15]利用LS-DYNA 對爆炸沖擊波與框架結構的作用過程進行了模擬，分析了爆炸沖擊波的傳播特性,結果表明合理的設置泄壓口能一定程度的減小沖擊波對結構的破壞。游俊[16]利用ABAQUS 對沖擊波與破片聯合作用下砌墻體進行了數值模擬，對結構動力響應進行分析，得到了沖擊波傳播以及與結構相互作用導致結構毀壞的損傷機理。LS-DYNA 等商業軟件是目前此類問題數值模擬的主要手段，國內未見具有類似功能的算法軟件，同時軟件詳細處理模型、方法也缺乏完整的介紹。

基于上述問題，本文發展一種沖擊波作用下結構毀傷算法，并對爆炸沖擊波與建筑物結構的相互作用過程進行模擬，研究沖擊波作用下建筑物的毀傷范圍、建筑物破壞以后沖擊波的傳播特性等，驗證該算法的適用性。

1 沖擊波與結構相互作用機理分析

構建如圖1 所示的建筑物墻體簡化模型，采用沖擊波沖量[17]作為結構是否破壞的指標，分析沖擊波與結構相互作用機理。圖1 中將包含建筑物墻體等結構的整個計算空間全部劃分為網格單元，共有1～4 個網格單元，計算每個單元受到的比沖量，判斷該單元是否被破壞：

⑴ 如果載荷不足以破壞結構，建筑物單元按照剛性壁面假設處理，沖擊波發生反射，形成比入射波高很多的超壓；

⑵ 如果沖量超過破壞載荷則該單元毀壞，那么它就變成可以在沖擊波作用下向前運動的碎塊。

在結構單元已經破壞的情況下，實際物理過程為結構碎塊在沖擊波作用下飛散。對于建筑物類脆性材料，由于結構碎塊四周沒有固支，碎塊的變形很小，可以將結構碎塊作為剛體處理，則沖擊波作用下結構碎塊的飛散問題是動邊界繞流問題。對于此類問題，可以采用非結構動網格技術進行模擬，即將流體力學方程與6 自由度彈道方程耦合求解，求解流體力學方程獲得碎塊表面壓力，進而得到碎塊受到的氣動力，將其作為彈道方程的輸入，求解彈道方程獲得碎塊的新位置，并據此對計算網格進行變形移動，獲得新的計算網格，重復這一過程即可獲得氣動力作用下的碎塊軌跡。L?hner 等[18]和Baum 等[-19]對大規模破片飛散過程進行數值模擬，模擬破片數量達到上百枚；爆炸波作用下的鋼筋混凝土墻發生變形、斷裂、碎裂過程[20]也得到了較好模擬。王巍等[21-23]對非結構動網格技術開展了大量研究，并將其應用于多種場景下多體運動問題的模擬，如飛行器頭罩分離[21]，助推器分離[22]，子母彈拋殼[23]等，計算精度得到了確認。本文采用上述發展的非結構動網格技術及程序，對如圖2 所示單塊碎塊在沖擊波作用下的運動過程進行模擬，分析沖擊波與結構相互作用機理。

計算模型如圖2 所示，封閉空間一側墻體上設置一洞口，洞口內存在一個堵塊，堵塊與墻體之間留有縫隙。其中堵塊的尺寸為1.14 m×0.18 m×0.6 m，質量為50 kg；洞口尺寸為1.18 m×0.18 m×0.67 m；封閉空間尺寸為5.8 m×6.0 m×3.3 m。

圖3 為計算得到的堵塊運動位置和壓力等值線。從圖中可以清晰看出，在碎塊完全脫離墻體之前，碎塊運動的稀疏波引起超壓降低；在碎塊剛剛脫離墻體時，氣體透過碎塊與墻體之間的縫隙傳播，由于縫隙面積小于洞口入口處面積，縫隙處氣流速度達到聲速，形成壅塞，此時沖擊波既有繞射，也有反射，流動非常復雜；在碎塊脫離墻體較遠后，碎塊和墻體之間縫隙形成的氣流壅塞消失，爆炸點所在空間內部高壓氣體在洞口處形成壅塞，且由于碎塊運動速度小于氣體速度，大量氣體繞過碎塊，在繞過碎塊后規整成球面波。

根據以上分析，圖4 給出了沖擊波與結構碎塊作用機理，沖擊波與結構碎塊作用過程中主要存在如下物理效應：(1) 碎塊運動會產生稀疏波效應，如圖4(a)示意；(2) 碎塊剛脫離墻體后的壅塞效應，如圖4(b)示意；(3) 繞射波重新規整為平面波的效應，如圖4(c)示意。從建立數學模型的角度看，能夠刻畫以上機理的模型就可以模擬沖擊波與結構碎塊作用過程。

上述結果表明，在建筑物結構破壞以后，利用動網格技術能夠模擬沖擊波與結構碎塊的相互作用過程，能夠獲得較為真實的數值模擬結果，并且可以直接利用現有的動網格計算模塊進行計算。但在實際應用中，對于整棟樓的數值模擬，結構單元的破壞會產生成百上千，甚至上萬的結構碎塊。利用動網格技術處理，需要進行網格變形計算，并且由于每個碎塊獨立運動，使得每時每刻都可能會出現局部網格變差的情況，就需要將計算程序停止，進行局部重構和插值，計算量巨大，大大降低了計算效率。因此，盡管非結構動網格技術能夠很好地處理沖擊波推動結構碎塊飛散的問題，技術上也不存在不可克服的難點，但對計算機硬件要求太高，不利于工程應用，需要更高效率的計算方法。

2 模型建立

從上一節分析可知，沖擊波作用下的結構毀傷過程，存在結構未破壞形成的沖擊波剛性壁面反射、結構完全破壞后形成的沖擊波透射以及結構破壞過程中因碎塊運動形成的稀疏波和壅塞效應。對于未被破壞結構，按照剛性壁面假設處理，不是本文關注的重點。對于已破壞結構單元，沖擊波在碎塊下游的繞射波重新規整為平面波效應有比較好描述，不需要特別的處理。因此，碎塊引起的稀疏波和壅塞兩種效應是建模的關鍵。

2.1 運動稀疏波效應的模型

采用任意拉格朗日-歐拉方法（ALE）坐標描述三維無量綱可壓縮非定常流動的Euler 方程：

式中：S是控制體V的表面，n為其外法向單位矢量；Q為流體守恒變量；Fc(Q,xc) 表示相對流通矢量，xc是網格運動速度。

提出如下模型：如圖5 所示，假設結構單元2、3 已經被破壞，在計算網格位置不變情況下，根據流固界面耦合算法、流場物理量和結構單元質量等參數，可以計算得到結構單元2、3 的網格運動速度 →xc，將其反映在流體控制方程中，計算流場參數，即可模擬結構單元2、3 運動引起的泄壓過程。在該計算過程中，盡管結構單元2、3 的網格位置沒有變化，但其單元運動速度真實的反映在流體控制方程中，從而在一定程度上模擬了結構單元運動引起的泄壓作用。

2.2 碎塊壅塞效應的計算模型

利用文獻[24]提出的虛擬網格透氣技術解決壅塞效應這一問題，圖6 為虛擬網格透氣技術原理圖，圖中單元1～3 為描述建筑物墻體結構的網格單元，單元4 為沖擊波一側流場中的氣體網格單元。完全氣體假設下的控制方程：

式中：Φ 的形式可簡寫為：

對于計算區域中單元4 物理量的計算可以表述為如下半離散形式：

式中：下標k用以表示不同表面。

對于圖6 所示情況，令單元4 與2 的邊界面為Sb，計算中采用如下形式：

式中：Φb為邊界面上的對流項，Sb表示單元4 和2 邊界投影面積。由于單元2 為結構單元，通常計算中式(5)中的Sb只能表示為固壁邊界，完全不通氣，不能模擬結構單元2 破壞后，沖擊波透過破壞單元傳播的過程，因而也就不能模擬碎塊破裂后的壅塞效應。

為了實現對結構單元破壞后邊界面由不透氣逐漸變為透氣，并最終達到完全透氣問題的模擬，文獻[24]引入如下計算模型：

式中：Sf、Sw分別是邊界面Sb通氣部分和不通氣部分的面積，即將邊界面分為通氣和不通氣兩部分，并計算對流項Φb：

式中：Φf是作為內部單元面處理時的對流項，可按一般式處理；Φw是作為固壁處理時的對流項。其對表面的通量積分為：

式中：U為流體相對于網格的速度；at為網格運動的法向速度；nx、ny、nz是n的三個分量;Sw為面積矢量。

對于具有虛擬擋板表面的單元，應該針對該表面通氣和壁面兩種狀態分別求解變量梯度，計算Φf和Φw時分別采用相應的梯度進行基本變量的重建。

將式(7)代入式(5)，則有如下形式：

分析式(9)，可以看出當Sf=0，即圖5 中單元4 與2 之間邊界面完全為固壁條件時，式(9)退化為式(5)，邊界面不透氣；當Sw=0 ，即單元4 與2 之間邊界面完全為流場內部面時，式(9)退化為式(4)，邊界面完全通氣；當處于邊界面狀態處于兩者之間時，單元4 與2 之間的邊界面可以實現按比例通氣。

根據第1 節中圖4 的分析結果，定義系數к為如下形式：

式中：S1為碎塊與墻體之間縫隙環帶的面積，可由碎塊運動速度、碎塊周長及運動時間等計算獲得；S2為碎塊的在來流方向的法向投影面積，固定不變，如圖7 所示。將к代入式(9)，S2代替其中Sb，則有如下形式：

當碎塊處于圖4(a)狀態時к=0，單元4 與2 之間邊界面完全為固壁條件，完全不通氣；當碎塊處于圖4(b)狀態時，к表示縫隙環帶面積與結構在來流方向的法向投影面積比值，模擬了縫隙對通氣量的限制；當碎塊處于圖4(c)狀態時к=1，單元4 與2 之間邊界面變為完全通氣狀體，此時結構碎塊對繞射后形成的第一道沖擊波強度影響很小。式(11)即為網格位置不變條件下，結構對沖擊波壅塞效應計算模型。

3 結果與分析

3.1 與非結構動網格模擬結果對比確認

對本文所提方法的計算過程，采用如圖6 所示的計算模型，模型的形狀、參數與圖2 中的模型一致，但洞口處沒有堵塊，不存在動網格，而是設置一個“虛擬擋板[24]”。計算過程中，當沖擊波到達洞口后，采用本文提出的計算模型，對沖擊波的傳播過程進行模擬，監測與圖2 中模型相同位置處超壓變化，并進行對比，分析驗證計算模型的準確性。其中，具體監測點位置A、B如圖8 所示，在墻壁中心位置；給定高壓氣團球心位置時，坐標定義如圖9 所示，坐標原點在堵塊中心。

針對上述兩個模型，模型1（堵塊模型），模型2（虛擬透氣計算模型），對不同參數設置的算例進行計算，算例情況見表1。

表1 算例參數設置Table 1 Example parameter setting

圖10～圖13 為不同算例的監測點A、B壓力變化，可以看出兩種模型情況下得到的監測點壓力峰值及規律均符合較好，表明本文提出的計算模型具有較高的準確度，能夠用于此類問題的模擬。

3.2 典型建筑物作用過程模擬結果分析

利用本文構建的爆炸沖擊波與建筑物結構作用計算模型，對沖擊波作用下建筑物的毀傷及沖擊波傳播過程進行了模擬，獲得了建筑物毀傷效果和沖擊波超壓分布。

3.2.1 計算模型

如圖14 所示，構建了一個大樓模型。為了提高計算效率，橫向寬度只包含3 個房間。其中，樓房左側的3、4、5 層建立結構計算單元，即墻體、樓板等結構部件內部建立網格單元，計算中根據結構單元毀傷載荷判據計算結構單元是否破壞，若破壞則變為透氣單元，樓房的其他部分不參與結構計算。計算中裝藥當量、結構毀傷載荷根據相關設計參數確定。

3.2.2 計算結果與分析

圖15 為建筑物結構毀傷計算結果，給出了結構單元破壞的過程。從計算結果來看，406 房間左右隔墻、上下樓板和走廊隔墻完全崩塌；405、407 房間外側墻體大部分崩塌；4 層走廊樓板大部分崩塌；306 房間下樓板和506 房間上樓板局部崩塌，較為真實地反映了試驗中建筑物的毀傷效果（其中，中間的房間是406，左邊405、右邊407，上邊506、下邊306）。計算中對每個單元獨立判斷其是否破壞，只與超壓分布有關，由于房間沖擊波反射、相交影響，使得超壓分布非常不均勻，導致出現墻體破壞也非常不均勻，且計算模型只考慮結構是否破壞透氣，未考慮結構間的支撐及運動情況，出現了四周均破壞而中間單元未破壞的情況，但計算結果仍具有較高的參考價值。

圖16 為爆炸流場演化云圖，與圖12 對照，可以看出沖擊波傳播與結構破壞特征一致，并且由于“虛擬擋板”對通氣量的限制，剛開始透過破壞單元傳播的沖擊波強度較弱，表明本文計算模型較好的模擬了結構碎塊對沖擊波的阻礙效應。

圖17 為406 房間監測點沖擊波超壓計算結果，由于沒有毀傷效應試驗超壓測量數據，將計算結果與無限空氣中爆炸沖擊波超壓經驗公式計算結果進行對比，其中TNT 爆炸當量為120 kg。表2 為經驗公式計算結果，經驗公式采用Baker 公式[25]，具體如下：

表2 監測點超壓經驗公式計算結果Table 2 The numerical calculation results of shock wave overpressure under different monitoring point

正反射超壓公式如下：

式中：w為TNT 藥量，kg；r為距離，m；Δp為超壓，kg/cm2。

比較圖14 與表2 計算結果，可以看出在爆炸房間內，沖擊波反射超壓數值計算結果與經驗公式結算結果符合較好，表明本文采用的爆炸流場計算方法具有較好的計算精度。

4 結 論

本文針對爆炸沖擊波與建筑物結構作用過程，在分析沖擊波與結構碎塊作用機理的基礎上，構建了一種計算模型，能夠在計算網格不動的條件下模擬沖擊波與結構的作用過程，主要結論如下：

(1)根據非結構動網格技術計算沖擊波推動結構碎塊飛散模擬中，非定常流場的變化規律，分析得出沖擊波與碎塊之間相互作用機理主要分為碎塊運動會產生稀疏波效應、碎塊剛脫離墻體的壅塞效應、繞射波重新規整為平面波效應。

(2)在機理清晰認識的基礎上，本文采用“虛擬運動”模擬稀疏波效應、“虛擬網格透氣技術[1]”處理碎塊和墻體之間的壅塞效應，發展了一種新的模擬爆炸沖擊波與建筑物作用過程的模型。通過算例對比，該模型計算結果與非結構動網格技術計算結果符合較好，表明該模型能夠用于沖擊波作用下結構破壞以及沖擊波傳播過程的模擬，具有良好的工程應用價值。

(3)建立建筑物模型，采用本文構建的爆炸沖擊波與建筑物結構作用計算模型 ，模擬了對建筑物的毀傷及沖擊波傳播過程，獲得了建筑物毀傷效果和沖擊波超壓分布。