水下爆炸載荷作用下變截面結構中應力波的傳播特性研究

孟子飛，徐思博，劉文韜，陳海龍

(哈爾濱工程大學 船舶工程學院，黑龍江 哈爾濱 150001)

0 引 言

近年來，國內外眾多學者對水下爆炸結構毀傷的力學機理進行了系統的研究[1-2]，但對于在水下爆炸載荷作用下應力波在艦船結構中的傳播特性研究甚少。目前，對于應力波在變截面結構中傳遞規律的研究主要以變截面桿為主，其中，王禮立等[3]運用一維桿中應力波的初等理論對應力波在錐形放大器中傳播特性進行了理論分析，指出了Leftheris 計算公式的錯誤之處并與實驗值進行對比，驗證了理論值的正確性，但計算得到的峰值應力比實測值偏大。周光泉、劉孝敏等[4-5]在此基礎上運用二維理論研究該問題并取得了較好的結果。但是，數值模擬與理論推導相比有著特殊的優勢，如果模型比較復雜并且載荷呈現強非線性時則難以利用理論推導出某處的應力狀態。另外，近年來隨著SHPB 實驗裝置研究領域的拓寬，其傳統的2 個基本假定（一維應力假定和均勻性假定）均受到挑戰[6-8]，而數值模擬對于優化實驗設計有著明顯的優勢，一方面數值模擬可以為實驗提供預測，如果數值模擬精度足夠的話完全可以代替實驗；另一方面在數值模擬中可以對應力波在不同材料、不同截面壓桿中的傳播特性進行系統的數值分析[9]。此外，在水下爆炸沖擊載荷作用下研究變截面桿以及舷側外板結構中應力波傳播的問題仍比較少見。因此，本文嘗試結合應力波理論，同時借助有限元的CEL 方法來探討在水下爆炸載荷作用下彈性波在變截面桿以及舷側外板結構中的傳播特性問題，并且利用沖擊響應譜分析方法探討舷側外板結構中的響應規律，為艦船遭遇水下爆炸后結構中應力波的傳播特性以及沖擊響應規律研究提供參考。

1 變截面桿中的應力波理論

本節在材料為線彈性的基礎上利用沖量定理、牛頓第三定律分析了彈性波在變截面桿中的反射和透射規律。

在如圖1 所示坐標系中研究一等截面均勻桿微元段在應力波下的縱向運動，其中A0為截面面積，為密度，F 為力。

圖 1 等截面均勻桿微元段示意圖Fig. 1 The sketch of a uniform cross-section bar

在推導控制方程之前，先作2 個基本假設：

1）桿在變形時橫截面保持為平面；

2）應力和應變之間遵循胡克定律。

規定應力和應變均以拉為正，質點速度以X 軸正向為正。

根據沖量定理，有

進而

將式（2）代入式（3）可得

聯立式（4）和式（5）得

2 水下爆炸流固耦合作用的CEL 方法

2.1 CEL 方法介紹

基于網格的數值方法有拉格朗日方法和歐拉方法，這2 種方法各有優缺點，拉格朗日方法能夠很好描述物質結構的運動，而歐拉方法能夠較好處理流體大變形的問題。由于拉格朗日方法和歐拉方法具有互補的性質，于是發展形成了拉格朗日方法和歐拉方法相互耦合的CEL 方法（耦合的歐拉-拉格朗日方法）。在CEL 算法中，歐拉網格和拉格朗日網格相互獨立。歐拉網格固定不動，材料在網格里流動，材料運動的同時它的各個物理量被映射到歐拉網格中，歐拉網格是獨立于空間位置的“監測點”實時追蹤流體的通量。而在拉格朗日網格中，網格節點固定在材料內部，網格隨著材料的變形而變形，質量、動量和能量隨著網格單元的運動而傳遞。在CEL計算時，結構利用拉格朗日方法求解，流體采用歐拉方法求解。在流體和結構相互耦合的界面上，拉格朗日單元受到歐拉材料給予的壓力載荷，歐拉單元接受拉格朗日單元施加的幾何約束，二者計算信息交換通過2 組網格之間的映射實現。基于以上原理，CEL 算法在解決物體的大位移，如碰撞、流體動力學及流體與固體之間的相互作用時，具有明顯優勢[11]。

2.2 控制方程

在CEL 方法中，拉格朗日數值方法的控制方程[12]如下：

根據質量守恒有

根據動量守恒有

根據能量守恒有

式中： ρ，v，σ，b，e分別為拉格朗日流體介質的密度，速度，應力，質量力和內能。

在CEL 方法中，歐拉數值方法的控制方程[12]如下：

歐拉算法中，材料在固定的網格里流動，根據質量守恒定律，有連續性方程

根據動量守恒有

根據能量守恒有

2.3 歐拉體積分數法

在CEL 中，材料在每一個增量步都會計算每個單元材料的分布以此來描述流體的變形狀態，而材料的初始分布情況在Abaqus 中利用“體積分數工具”進行計算。該工具對歐拉體與參考體之間進行布爾運算來確定材料的分布。圖2 中利用“體積分數工具”將T 型材作為參考體映射到歐拉體后計算得到的歐拉體積分數示意圖，圖中灰色部分代表在歐拉網格中填充“水”材料，數字代表“水”在歐拉單元中所占的體積分數，白色部分代表T 型材結構其材料的敷設在拉格朗日體中進行。后處理中通過EVF 選項卡來觀察流體的運動情況。

圖 2 歐拉體積分數示意圖Fig. 2 The sketch of eulerian volume fraction

2.4 罰函數接觸方法

在CEL 方法中，歐拉域與拉格朗日域之間的接觸是通過罰函數接觸方法實現的。拉格朗日單元邊和面上建立節點，歐拉材料表面建立錨點，二者通過“硬”接觸-壓力過盈行為建立耦合關系[13]。這種方法允許歐拉材料和拉格朗日體之間存在一些微小的滲透。節點與錨點之間的接觸力Fp正比于二者的滲透距離d[14]。

式中：k 為罰函數剛度系數，取決于拉格朗日和歐拉材料的屬性。

2.5 CEL 方法驗證

本節采用自由場水下爆炸模型，以仿真值與經驗公式值對比的方式對CEL 方法進行驗證，具體的幾何模型如圖3 所示。

球形T N T 的藥量為0.0 5 5 g，裝藥密度1 600 kg/m3，爆速 。爆轟產物的壓力與爆轟產物相對比容和比內能之間的關系釆用JWL 狀態方程表示[13]。方程中參數具體值見表1。

水的密度ρ=1 000 kg/m3，粘性系數為0.001，采用US-UP形式的EOS 狀態方程，其中，水的聲速C0=1 480 m/s。

圖 3 自由場水下爆炸幾何模型示意圖Fig. 3 The geometry model for underwater explosion of free field

表 1 JWL 方程的參數Tab. 1 The parameter of JWL equation

為減小邊界對沖擊波影響，計算歐拉域的所有邊界都設定為流出無反射邊界，在無邊界效應影響的水域進行爆炸，其產生沖擊波的波陣面壓力峰值pm可由Cole 和Zamyshlyayev 給出的經驗公式[15-16]表示。圖4為在不同距離處CEL 計算壓力峰值與經驗公式計算壓力峰值的對比。從圖中可以發現，CEL 方法計算出來的結果與Cole 和Zamyshlyayev 給出的經驗公式計算結果相似，準確度較高。因此，采用CEL 方法可以充分模擬水下爆炸沖擊載荷。

圖 4 不同距離處的壓力峰值對比Fig. 4 Peak pressure contrast at different distances

3 水下爆炸載荷作用下變截面桿中應力波的傳播特性

下面采用CEL 方法模擬出水下爆炸沖擊載荷，然后通過有限元方法模擬沖擊波在變截面桿中的傳播問題。此問題涉及液固相互作用并且水下爆炸載荷呈現非線性難以利用理論公式得到解析解，本節運用CEL方法來處理流固耦合問題以期達到研究變截面桿在受到水下爆炸載荷下彈性波傳播特性的目的。

本節建立1 根兩段等長L 的變截面桿，L 為200 mm，桿左端直徑D 為4 mm，桿的截面比為2∶1，桿的材料為普通鋼（不考慮塑性），桿的兩端均為自由端，在距桿左端50 mm 處建立和第2.5 節一樣的TNT 炸藥模型，最后建立一足夠大的Euler 域以保證TNT 炸藥具有足夠的傳播空間，Euler 域四周均設置無反射邊界條件，具體模型示意圖如圖5 和圖6 所示。

圖 5 變截面桿幾何示意圖Fig. 5 The geometry model of a variable section bar

圖7～圖9 給出了爆炸沖擊波從開始進入到桿的左端、傳遞到變截面附近、傳遞到桿右端附近3 個典型位置圖。當時間s 時，提取出入射波作用下桿左端到右端的應力隨距離傳遞曲線和速度隨距離傳遞曲線，計算距離從0～0.2 m 對應的應力峰值作為速度峰值作為 vI=0.28 m/s。當時間s 時，提取出透射波作用下桿左端到右端的應力隨距離傳遞曲線和速度隨距離傳遞曲線，計算距離從0.2～0.4 m 對應的應力峰值作為，速度峰值作為m/s。各曲線如圖10 所示，與理論解進行對比如表2 所示。可以發現變截面桿在受到水下爆炸載荷作用下彈性波規律十分接近理論解，并且可以看出當應力波由大截面向小截面傳播時，透射波擾動強度大于入射波擾動，這種結構起到“聚波”作用。由于在t=1.32×10-4s 時，桿左端到桿變截面處存在反射波與部分入射波的疊加，所以難以提取反射波作用下桿的應力傳遞曲線和速度傳遞曲線，在此不予處理。

圖 7 沖擊波進入桿左端典型位置圖Fig. 7 The shock wave enters the left of the bar

圖 8 沖擊波傳遞到變截面附近位置圖Fig. 8 The shock wave enters the center of the bar

圖 9 沖擊波傳遞到桿右端附近位置圖Fig. 9 The shock wave enters the right of the bar

4 水下爆炸載荷作用下舷側外板中應力波的傳播特性

4.1 模型建立

本節借助CEL 方法以期達到研究舷側外板在受到水下爆炸載荷下彈性波傳播特性的目的。此問題與上節變截面桿模型相比，不僅涉及液固相互作用而且模型更為復雜，無法利用理論公式得到解析解。

圖 10 桿的應力和速度隨距離變化曲線Fig. 10 Stress and velocity change with distance

表 2 仿真值與理論公式值對比Tab. 2 Comparison of simulation and theoretical values

圖 11 舷側外板遭遇水下爆炸仿真模型Fig. 11 The simulation model of side plating underunderwater explosion

4.2 應力波從舷側外板到T 型材腹板傳播的特性及沖擊響應譜分析

采用頻譜分析方法和沖擊響應譜分析方法，探究應力波從舷側外板結構往T 型材腹板中的傳播特性。分別提取出#1，#2，#7，#8 四個測點的加速度隨時間變化的數據，并進行快速傅里葉變換得到對應的頻譜曲線，如圖12 所示。從圖中可以看出，在[500,2 000] Hz 較高頻率范圍內的應力波從舷側外板經過第1 個T 型材后傳遞到T 型材腹板處的強度大為減小，而[0, 500] Hz 內的應力波強度則相差不大。這說明應力波在經過第1 個從小到大截面時濾掉了大部分高頻信號。

圖 12 各測點頻段分析Fig. 12 Analysis on frequency band at various gauging points

將以上測點的加速度時歷數據進行轉換得到沖擊響應譜曲線，如圖13 所示。可以發現，對于[500,2 000] Hz 較高頻率范圍內的沖擊響應經過第1 個T 型材后呈大幅衰減趨勢。這說明舷側外板在受到水下爆炸載荷作用下[500, 2 000] Hz 內的T 型材腹板的結構響應與外板結構響應相比大為減小。而[0, 500] Hz 內的響應相對來說變化則較為平緩。總體來看，從舷側外板到第1 個T 型材，低頻段的結構響應譜值兩者相差不大，而在高頻段外板的譜值遠高于T 型材腹板中的譜值。

圖 13 各測點處沖擊響應譜曲線Fig. 13 The curve of shock response spectrum at various gauging points

分別提取出#9，#10 測點的加速度隨時間變化的數據，并進行快速傅里葉變換，得到對應的頻譜曲線，并與#3，#4 進行對比，如圖14 所示。可以看出，在[0, 500] Hz 的應力波從舷側外板經過第2 個T 型材后傳遞到T 型材腹板處的強度大為減小，而在其他頻段中腹板的應力波強度與外板相比相差不大。這說明當應力波再次經過從小到大截面形式的變化時濾掉了大部分低頻信號。

圖 14 各測點頻段分析Fig. 14 Analysis on frequency band at various gauging points

將以上測點的加速度時歷數據進行轉換得到沖擊響應譜曲線，如圖15 所示。可以看出，沖擊響應規律與經過第1 個T 型材后的規律不同，舷側外板在受到水下爆炸載荷作用下，只有[1 300, 1 800] Hz 內的結構響應才大于腹板的響應，在其他頻段兩者則呈波動交替趨勢，幅值相差不大。

圖 15 各測點處沖擊響應譜曲線Fig. 15 Shock response spectrum at various gauging points

結合上述描述可以得出，當應力波再次經過從小到大截面形式的變化時，舷側外板產生的中高頻段沖擊響應大于T 型材腹板中的響應，而其他頻段的響應相差不大。

4.3 應力波在舷側外板中的傳播特性及沖擊響應譜分析

本節采用上節相同的分析方法，進一步探究應力波在舷側外板結構中的傳播特性。

分別提取出#1～#6 六個測點的加速度隨時間變化的數據，并進行快速傅里葉變換，得到對應的頻譜曲線，如圖16 所示。

通過圖16（a）與圖16（b）兩組測點相互對比可以發現，在[500, 2 000] Hz 內#3，#4 測點的幅值比#1，#2 測點小，這說明在該頻率范圍內的應力波經過第1 個T 型材后強度呈大幅衰減趨勢，衰減幅度約60%，而#3，#4 測點在[0, 500] Hz 內的幅值比對應的#1，#2 測點大，這說明在該頻率范圍內的應力波經過第1 個T 型材后強度略微增強。結合這兩點，說明當應力波經過第1 個T 型材后濾掉了大部分高頻信號，而增強了低頻信號，前者是由于應力波在經過第1 個從小截面到大截面形式的變化時會濾掉高頻信號的特性所導致，后者是因為應力波從大截面進入小截面其強度會增強的原理所導致。

圖 16 各測點頻段分析Fig. 16 Analysis on frequency band at various gauging points

從圖16（c）與圖16（d）可以發現，#5，#6 測點在[0, 500] Hz 內的幅值比對應的#3，#4 測點小，這說明在該頻率范圍內的應力波經過第2 個T 型材后強度呈衰減趨勢，衰減幅度約50%，而#5，#6 測點在[500,2 000] Hz 內的幅值相對于#3，#4 測點來說衰減幅度并不大。基于此，這說明應力波再次經過相同截面形式的T 型材后會濾掉大部分低頻信號，這與4.2 節的結論相吻合。

將以上測點的加速度時歷數據進行轉換得到其對應的沖擊響應譜曲線，如圖17 所示。可以看出，#1，#3，#5 測點處沖擊響應譜曲線與#2，#4，#6 測點處沖擊響應譜曲線規律相似。為方便起見，下面以#2，#4，#6 測點處沖擊響應譜曲線規律進行闡述。在[0,2 000] Hz 內，#2，#4，#6 測點的譜峰值依次減小，說明舷側外板在受到水下爆炸載荷作用下，在[0,2 000] Hz 內的結構響應峰值在經過2 個T 型材后呈衰減趨勢。具體說來，在[0, 500] Hz 內#2 測點的譜峰值小于#4，則說明舷側外板在受到水下爆炸載荷作用下[0, 500] Hz 內的結構響應在經過第1 個T 型材后呈放大趨勢，而#6 測點的譜峰值小于#4 測點，這說明舷側外板在受到水下爆炸載荷作用下[0, 500] Hz 內的結構響應在經過第2 個T 型材后呈衰減趨勢。對于[500,2 000] Hz 較高頻率范圍內的沖擊響應經過2 個T 型材后整體呈衰減趨勢。

綜上所述，舷側外板中低頻段的沖擊響應譜峰值沿著從下往上的方向先升高再下降而高頻段的譜峰值則依次下降。

圖 17 各測點處沖擊響應譜曲線Fig. 17 The curve of shock response spectrum at various gauging points

5 結 語

本文基于應力波理論，結合CEL 方法，首先研究了在水下爆炸載荷作用下變截面桿結構中應力波的傳播特性，與理論公式對比驗證了數值精度，然后進一步探討了舷側外板結構在遭遇水下爆炸載荷后的應力波傳播特性，得到如下結論：

1）當應力波由大截面向小截面傳播時，透射波擾動強度大于入射波擾動，這種結構起到“聚波”作用。

2）應力波在經過第1 個從小到大截面形式的變化時濾掉了大部分高頻信號，當應力波再次經過從小到大截面形式的變化時會進一步濾掉響應信號。

3）在舷側外板到T 型材腹板這種從小到大截面形式變化的結構中，激勵產生的低頻段結構響應譜值，兩者之間相差不大，而外板的高頻段譜值遠高于T 型材腹板中的譜值。當應力波再次經過這種截面形式的變化時，舷側外板產生的中高頻段沖擊響應大于T 型材腹板中的響應，而其他頻段的響應相差不大。

4）在含有2 個T 型材的舷側外板中，低頻段的沖擊響應譜峰值沿著從下往上的方向先升高再下降而高頻段的譜峰值則依次下降。

以上結論對研究船舶與海洋工程領域的型材、縱骨、縱桁等變截面結構中應力波傳播特性以及后續研究復雜艦船結構的沖擊響應規律具有一定的參考意義。