水下爆炸載荷作用下變截面結(jié)構(gòu)中應(yīng)力波的傳播特性研究

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

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

0 引 言

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

1 變截面桿中的應(yīng)力波理論

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

在如圖1 所示坐標(biāo)系中研究一等截面均勻桿微元段在應(yīng)力波下的縱向運(yùn)動，其中A0為截面面積，為密度，F(xiàn) 為力。

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

在推導(dǎo)控制方程之前，先作2 個基本假設(shè)：

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

2）應(yīng)力和應(yīng)變之間遵循胡克定律。

規(guī)定應(yīng)力和應(yīng)變均以拉為正，質(zhì)點速度以X 軸正向為正。

根據(jù)沖量定理，有

進(jìn)而

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

聯(lián)立式（4）和式（5）得

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

2.1 CEL 方法介紹

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

2.2 控制方程

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

根據(jù)質(zhì)量守恒有

根據(jù)動量守恒有

根據(jù)能量守恒有

式中： ρ，v，σ，b，e分別為拉格朗日流體介質(zhì)的密度，速度，應(yīng)力，質(zhì)量力和內(nèi)能。

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

歐拉算法中，材料在固定的網(wǎng)格里流動，根據(jù)質(zhì)量守恒定律，有連續(xù)性方程

根據(jù)動量守恒有

根據(jù)能量守恒有

2.3 歐拉體積分?jǐn)?shù)法

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

圖 2 歐拉體積分?jǐn)?shù)示意圖Fig. 2 The sketch of eulerian volume fraction

2.4 罰函數(shù)接觸方法

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

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

2.5 CEL 方法驗證

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

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

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

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

表 1 JWL 方程的參數(shù)Tab. 1 The parameter of JWL equation

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

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

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

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

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

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

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

圖 7 沖擊波進(jìn)入桿左端典型位置圖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 水下爆炸載荷作用下舷側(cè)外板中應(yīng)力波的傳播特性

4.1 模型建立

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

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

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

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

4.2 應(yīng)力波從舷側(cè)外板到T 型材腹板傳播的特性及沖擊響應(yīng)譜分析

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

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

將以上測點的加速度時歷數(shù)據(jù)進(jìn)行轉(zhuǎn)換得到?jīng)_擊響應(yīng)譜曲線，如圖13 所示?？梢园l(fā)現(xiàn)，對于[500,2 000] Hz 較高頻率范圍內(nèi)的沖擊響應(yīng)經(jīng)過第1 個T 型材后呈大幅衰減趨勢。這說明舷側(cè)外板在受到水下爆炸載荷作用下[500, 2 000] Hz 內(nèi)的T 型材腹板的結(jié)構(gòu)響應(yīng)與外板結(jié)構(gòu)響應(yīng)相比大為減小。而[0, 500] Hz 內(nèi)的響應(yīng)相對來說變化則較為平緩?？傮w來看，從舷側(cè)外板到第1 個T 型材，低頻段的結(jié)構(gòu)響應(yīng)譜值兩者相差不大，而在高頻段外板的譜值遠(yuǎn)高于T 型材腹板中的譜值。

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

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

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

將以上測點的加速度時歷數(shù)據(jù)進(jìn)行轉(zhuǎn)換得到?jīng)_擊響應(yīng)譜曲線，如圖15 所示?？梢钥闯?，沖擊響應(yīng)規(guī)律與經(jīng)過第1 個T 型材后的規(guī)律不同，舷側(cè)外板在受到水下爆炸載荷作用下，只有[1 300, 1 800] Hz 內(nèi)的結(jié)構(gòu)響應(yīng)才大于腹板的響應(yīng)，在其他頻段兩者則呈波動交替趨勢，幅值相差不大。

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

結(jié)合上述描述可以得出，當(dāng)應(yīng)力波再次經(jīng)過從小到大截面形式的變化時，舷側(cè)外板產(chǎn)生的中高頻段沖擊響應(yīng)大于T 型材腹板中的響應(yīng)，而其他頻段的響應(yīng)相差不大。

4.3 應(yīng)力波在舷側(cè)外板中的傳播特性及沖擊響應(yīng)譜分析

本節(jié)采用上節(jié)相同的分析方法，進(jìn)一步探究應(yīng)力波在舷側(cè)外板結(jié)構(gòu)中的傳播特性。

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

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

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

從圖16（c）與圖16（d）可以發(fā)現(xiàn)，#5，#6 測點在[0, 500] Hz 內(nèi)的幅值比對應(yīng)的#3，#4 測點小，這說明在該頻率范圍內(nèi)的應(yīng)力波經(jīng)過第2 個T 型材后強(qiáng)度呈衰減趨勢，衰減幅度約50%，而#5，#6 測點在[500,2 000] Hz 內(nèi)的幅值相對于#3，#4 測點來說衰減幅度并不大?；诖耍@說明應(yīng)力波再次經(jīng)過相同截面形式的T 型材后會濾掉大部分低頻信號，這與4.2 節(jié)的結(jié)論相吻合。

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

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

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

5 結(jié) 語

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

1）當(dāng)應(yīng)力波由大截面向小截面?zhèn)鞑r，透射波擾動強(qiáng)度大于入射波擾動，這種結(jié)構(gòu)起到“聚波”作用。

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

3）在舷側(cè)外板到T 型材腹板這種從小到大截面形式變化的結(jié)構(gòu)中，激勵產(chǎn)生的低頻段結(jié)構(gòu)響應(yīng)譜值，兩者之間相差不大，而外板的高頻段譜值遠(yuǎn)高于T 型材腹板中的譜值。當(dāng)應(yīng)力波再次經(jīng)過這種截面形式的變化時，舷側(cè)外板產(chǎn)生的中高頻段沖擊響應(yīng)大于T 型材腹板中的響應(yīng)，而其他頻段的響應(yīng)相差不大。

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

以上結(jié)論對研究船舶與海洋工程領(lǐng)域的型材、縱骨、縱桁等變截面結(jié)構(gòu)中應(yīng)力波傳播特性以及后續(xù)研究復(fù)雜艦船結(jié)構(gòu)的沖擊響應(yīng)規(guī)律具有一定的參考意義。