水下爆炸圓柱殼沖擊波繞射特性分析

姚熊亮， 屈子悅， 張樂， 王志凱

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

在水下爆炸作用下，由于高能炸藥的爆炸、沖擊波的傳播以及復(fù)雜的流體-結(jié)構(gòu)相互作用現(xiàn)象的存在，使得結(jié)構(gòu)周圍的水動力場變得十分復(fù)雜。潛艇作為現(xiàn)代海戰(zhàn)的重要組成部分，由于其圓柱形結(jié)構(gòu)的特殊性，在分析水下爆炸沖擊波波動效應(yīng)時，除了要考慮反射、透射沖擊波外，在圓柱殼背爆面出現(xiàn)的繞射沖擊波也不容忽視。

姚熊亮等[1]使用ABAQUS有限元軟件，利用聲固耦合法對水下爆炸載荷作用下水面艦艇和潛艇結(jié)構(gòu)的動態(tài)響應(yīng)進(jìn)行了數(shù)值模擬。朱正洋[2]對近地空爆沖擊波進(jìn)行了數(shù)值模擬，得出了沖擊波在空氣介質(zhì)中的傳播過程和作用在結(jié)構(gòu)上的繞射載荷效應(yīng)。Ofengeim等[3]通過數(shù)值方法研究了平面激波在剛性圓柱體上的繞射。年鑫哲等[4]分析了空氣沖擊波作用于柔性防爆墻的透射和繞射效應(yīng)，比較了壓力波形的變化，得到了墻后壓力場的變化規(guī)律。邵宗戰(zhàn)[5]為研究水下爆炸對魚雷、水雷等水下目標(biāo)的破壞機(jī)理開展了一次水下爆炸試驗(yàn)，測量了圓柱殼正面、側(cè)面和背面的壓力場，得到了沖擊波壓力時程曲線。

圓柱殼體在水下爆炸沖擊波作用下的壓力空間分布情況的認(rèn)識尚不明確。本文基于空氣動力學(xué)理論，從膨脹波與激波的形成、波在膨脹波與激波中流動過程出發(fā)，來解釋沖擊波繞射的形成，并利用ALE流固耦合算法，探究圓柱殼結(jié)構(gòu)周圍壓力分布特點(diǎn)與繞射波特性。

1 數(shù)值分析模型

1.1 物理模型

本文所采用的物理模型以圓柱殼中心為坐標(biāo)原點(diǎn)建立模型，圓柱殼半徑為R、厚度為h的薄圓柱殼，殼體內(nèi)部為空氣，殼體外部為水，水流場中有TNT當(dāng)量為W的藥包，藥包中心與殼體間的距離為r，如圖1所示。

圖1 物理模型Fig.1 Physical model

本文采用LS-DYNA中的任意拉格朗日歐拉(ALE)算法實(shí)現(xiàn)沖擊波與圓柱殼的相互作用模擬。依據(jù)上述物理模型并結(jié)合文獻(xiàn)[6]中的實(shí)驗(yàn)?zāi)Ｐ统叽纾㈤L1.35 m，直徑1 m，板厚為10 mm的圓柱殼體。模型橫剖面圖如圖2所示，殼體內(nèi)部為空氣，外部為水，流場網(wǎng)格大小為0.04 m，結(jié)構(gòu)網(wǎng)格大小為0.02 m。流場邊界采用無反射邊界條件。基于計(jì)算模型的對稱性，在殼體上設(shè)置A、B、C、D、E5個測點(diǎn)。

1.2 模型參數(shù)設(shè)置

1)結(jié)構(gòu)材料特性。

本文所采用高強(qiáng)度低碳鋼，考慮材料應(yīng)變率敏感性的本構(gòu)方程有很多，本次圓柱殼體采用與實(shí)驗(yàn)數(shù)據(jù)吻合較好的Cowper-Symonds模型來描述材料的應(yīng)變率效應(yīng)，動態(tài)屈服強(qiáng)度：

(1)

式中：σd為動態(tài)屈服強(qiáng)度；σ0為屈服極限；E為楊氏模量；ET為切線硬化模量；εP為有效塑性應(yīng)變；εd為等效塑性應(yīng)變；D、δ為應(yīng)變率參數(shù)，可由試驗(yàn)獲得，文獻(xiàn)[7-10]均提到了低碳鋼有實(shí)驗(yàn)得到的與應(yīng)變率相關(guān)的參數(shù)：D=40.4 s-1，δ=5。材料相關(guān)參數(shù)見表1。其中，ρ1為鋼材密度；μ為泊松比。

圖2 計(jì)算模型Fig.2 Computation model

表1 模型參數(shù)Table 1 Model parameters

2)空氣與水狀態(tài)方程。

本次采用流體常用模型MAT _null模擬流場本構(gòu)關(guān)系，空氣狀態(tài)方程采用EOS_linear_polynormal方程，在LS-DYNA程序中，水的狀態(tài)方程采用Grüneisen狀態(tài)方程[11]，相關(guān)參數(shù)[12]見表2、表3。

表2 空氣狀態(tài)方程參數(shù)Table 2 Air flow field state equation parameters

表中ρ2代表空氣密度，E0代表初始內(nèi)能，V0代表初始相對體積，C0～C6為方程系數(shù)，線性多項(xiàng)式狀態(tài)方程：

P=C0+C1ζ+C2ζ2+C3ζ3+E′(C4+C5ζ+C6ζ2)

(2)

式中：ζ為當(dāng)前密度與初始密度之比；E′為當(dāng)前內(nèi)能。

表3 水狀態(tài)方程參數(shù)Table 3 Water flow field state equation parameters

線性多項(xiàng)式狀態(tài)方程見式(3):

(γ0+αμ)E1

(3)

式中：P為壓力；μ=ρ/ρ0-1；ρ3為水密度；E1為單位體積內(nèi)能；S1、S2、S3、γ0、α為常數(shù)。水的密度約為1 000 kg/m3，波在水中的傳播速度c約為1 500 m/s，水的運(yùn)動粘性系數(shù)約為ν=1.13×10-3N·s/m2。

3)炸藥狀態(tài)方程。

H. Hamashima等[13]通過水下爆炸實(shí)驗(yàn)擬合了JWL狀態(tài)方程，使其能夠更好地反應(yīng)炸藥膨脹過程中的二次反應(yīng)和約束條件。對于TNT炸藥，采用JWL[14]狀態(tài)方程：

(4)

式中：V當(dāng)前相對體積，A1、A2、r1、r2為常數(shù)，ω為比熱，各參數(shù)的取值表4。

表4 JWL狀態(tài)方程參數(shù)Table 4 JWL state equation parameters

1.3 壓力對比分析

1.3.1 自由場壓力對比分析

為驗(yàn)證爆炸沖擊波載荷的準(zhǔn)確性，對自由場中的水下爆炸沖擊波進(jìn)行數(shù)值模擬。計(jì)算中TNT當(dāng)量為1 kg，藥包初始半徑R0=0.053 m，采用中心點(diǎn)爆的形式。壓力測點(diǎn)與爆心的距離進(jìn)行無量綱化r→r/R0，測點(diǎn)距爆心的距離r=10，r=13，r=16，r=19，計(jì)算域?yàn)檫呴L為2.5 m的正方體區(qū)域，由于流場對稱，所以只取1/2流場進(jìn)行計(jì)算。流場邊界采用無反射邊界條件，流場剖面采用對稱邊界條件，流場網(wǎng)格大小為0.04 m。

圖3中給出了自由場中水下爆炸計(jì)算的數(shù)值模型和典型時刻自由場中的壓力云圖。圖4中給出了r=10，r=13，r=16，r=19處測點(diǎn)的沖擊波壓力與Zamyshlyayev[15]公式值的對比曲線，從圖中可得數(shù)值模擬結(jié)果衰減過程與經(jīng)驗(yàn)公式吻合良好，壓力曲線在峰值過后出現(xiàn)多次壓力震蕩，Cole[16]給出的解釋是由于裝藥中不同位置炸藥顆粒(對應(yīng)于數(shù)值計(jì)算模型中不同炸藥單元)爆炸后沖擊波疊加引起的[17]。

圖3 計(jì)算模型與典型時刻壓力云圖Fig.3 Computation model and pressure cloud chart

1.3.2 實(shí)驗(yàn)值與數(shù)值計(jì)算值的比較

根據(jù)文獻(xiàn)[6]中的實(shí)驗(yàn)，設(shè)置1 kg藥量，4 m爆距的工況下進(jìn)行數(shù)值模擬計(jì)算，得到圓柱殼迎爆面中心、上表面中心、背爆面中心的壓力峰值，與實(shí)驗(yàn)值比較，如表5所示。從表中可以看出，實(shí)驗(yàn)值與仿真值的誤差在10%左右。

表5 壓力峰值對比表Table 5 Peak pressure contrast

2 沖擊波繞射特性分析

2.1 繞射波形成原理

關(guān)于傳統(tǒng)繞射波現(xiàn)象的研究雖然比較普遍，且主要集中在空氣場中。有關(guān)水下爆炸沖擊波作用下結(jié)構(gòu)毀傷模式和特征的研究雖然已有不少，但水下強(qiáng)沖擊波與圓柱體耦合繞射波形成原理分析較少，也沒有形成相關(guān)理論體系。當(dāng)水下爆炸沖擊波在傳播過程中遇到壁面時，其傳播方向會發(fā)生改變，但傳播方向是如何改變的，其中的作用原理是什么？本文基于空氣動力學(xué)理論，從膨脹波與激波的形成、波在膨脹波與激波中流動過程出發(fā)，來解釋沖擊波繞射的形成。

超聲速氣流沿外折微小角度的壁面流動時，產(chǎn)生弱膨脹波，氣流經(jīng)過此波后速度增大，壓強(qiáng)、密度、溫度降低。二維無粘性超聲速氣流繞折角為dθ的外凹角流動，如圖5所示。在這點(diǎn)O之前，氣流平行于壁面AO，流至點(diǎn)O，壁面向外轉(zhuǎn)折dθ，流動空間擴(kuò)大，氣流膨脹，并伴隨著壓強(qiáng)、溫度減小。壓強(qiáng)變化的擾動從點(diǎn)O向外傳播，形成一道弱膨脹波，它和波前氣流方向之間的夾角μ1只取決于馬赫數(shù)Ma1，即

(5)

圖4 爆炸沖擊波數(shù)值壓力與經(jīng)驗(yàn)公式值對比Fig.4 Comparisons numerical value and empirical formulas

圖5 微小外折角產(chǎn)生的膨脹波Fig.5 Expansion wave produced by small external bend angle

二維無粘性超聲速繞外凹曲面的流動如圖6所示。曲面的彎曲部分O1On可看成由無限多個連續(xù)轉(zhuǎn)折無限小折角構(gòu)成。從彎曲部分上的每一點(diǎn)都發(fā)出一條弱膨脹波，這無限多條弱膨脹波構(gòu)成“連續(xù)膨脹波”，氣流經(jīng)過這些膨脹波后，馬赫數(shù)由Ma1增大到Ma2，方向由平行于壁面AO1折轉(zhuǎn)到平行于壁面OnB，這一過程是漸變等熵過程。

根據(jù)極限的概念，一個連續(xù)的圓柱形拱殼可以看成是由無數(shù)段外折的微元折壁所組成的，如圖7所示。當(dāng)爆炸沖擊波繞壁面?zhèn)鞑r，圓柱殼壁面上的每一個點(diǎn)都相當(dāng)于一個折點(diǎn)，因此，每一點(diǎn)都將發(fā)出一條弱膨脹波，沖擊波每經(jīng)過一道膨脹波，參數(shù)值有一個微小的變化，通過一段膨脹波后，參數(shù)值及折角發(fā)生一個有限量的變化，從而形成了傳播方向發(fā)生改變的繞射波。

圖6 超聲速繞外凸曲面的流動Fig.6 Supersonic flow around convex surfaces

圖7 繞射波的形成Fig.7 Formation of diffracted waves

2.2 爆炸沖擊波的傳播過程

采用ALE流固耦合法對模型進(jìn)行數(shù)值仿真計(jì)算。計(jì)算得到的流場壓力云圖如圖8所示。裝藥起爆后，當(dāng)沖擊波作用在殼體后，由于殼體內(nèi)部介質(zhì)為空氣，大部分能量以沖擊波的形式返回到流場中，形成反射沖擊波，如圖8(a)所示；反射沖擊波向爆心方向傳播，當(dāng)傳播至氣泡表面后，由于氣泡內(nèi)部介質(zhì)的阻抗低于氣泡外部介質(zhì)的阻抗，又有一部分能量以稀疏波的形式反射回到流場中，另一部分則透射到氣泡內(nèi)部，氣泡反射回的稀疏波與殼體反射的稀疏波相遇后會相互疊加，使壓力迅速降低，如圖(b)；初始沖擊波繼續(xù)向前傳播，當(dāng)沖擊波繞過圓柱殼頂部和底部時，形成膨脹波，沖擊波的傳播方向發(fā)生改變，繞過臨界位置向背爆面?zhèn)鞑ィ鐖D8(c)、(d)所示。

圖8 流場壓力云圖Fig.8 Pressure distribution map of flow field

2.3 近壁面壓力分析

2.3.1 近壁面壓力空間分布

取圖2中所示的5個近壁面測點(diǎn)為研究對象。工況1～4的殼體厚度為10 mm，對爆距r進(jìn)行無量綱化處理，r→r/R(R為殼體半徑)，爆距分別為0.5、1、2、3；工況5～8在的爆距為2，殼體厚度分別為5、10、15、20 mm。

圖9、10為圓柱殼表面測點(diǎn)A～E的壓力在不同工況下的時程曲線，測點(diǎn)位置如圖2所示。當(dāng)初始爆炸沖擊波傳播到A點(diǎn)后，A點(diǎn)出現(xiàn)第1次壓力峰值，隨后壓力迅速衰減，由于殼體中應(yīng)力波的反射形成了第2個壓力峰值；B點(diǎn)處壓力峰值降低，壓力變化規(guī)律與A點(diǎn)相似；測點(diǎn)C位于臨界點(diǎn)，壓力峰值降低，脈寬增大；D、E2點(diǎn)壓力變化規(guī)律與迎爆面測點(diǎn)壓力變化規(guī)律差異較大，繞射沖擊波的到達(dá)形成的壓力值較低，脈寬增加。通過圖9和圖10的對比可以發(fā)現(xiàn)，板厚主要影響應(yīng)力波的透射規(guī)律，各點(diǎn)的載荷時歷特征受爆距參數(shù)影響更加明顯。

圖9 不同爆距下各測點(diǎn)壓力時程曲線Fig.9 Time-history curve of pressure at different explosion distances

圖10 不同板厚下各測點(diǎn)壓力時程曲線Fig.10 Time-history curve of pressure at different thicknesses

2.3.2 近壁面壓力峰值特征

為探究沖擊波近壁面壓力的峰值在圓柱殼周圍的分布特征，定義無量綱相對壓力峰值λPi=Pmi/PmA，其中Pmi表示第i位置處的壓力峰值。λP變化趨勢如圖11所示。

圖11 λP值折線圖Fig.11 The broken-line graph of λP values

從圖11中可以看出，λPB的值最大，臨界點(diǎn)和背爆面測點(diǎn)λPC、λPD、λPE的值依次減小；臨界點(diǎn)C處λPC與迎爆面處的λPB值相差較大，與背爆面處的λPD、λPE值相差較小；隨著爆距與板厚的增加，λPB-λPE值也隨之增大。不管爆距與結(jié)構(gòu)形式如何改變，背爆面λPD、λPE的值始終非常接近。

2.3.3 脈寬變化特征

壓力峰值所代表的是瞬時的壓力，而脈寬代表了沖擊波作用的時間，同樣是一個重要的參量，定義脈寬是指沖擊波從峰值衰減到峰值的1/e所需的時間，定義無量綱相對脈寬λτi=τi/τA，其中τi表示第i位置處的脈寬。圖12給出了各工況下每個測點(diǎn)的λτi=τi/τA。背爆面測點(diǎn)處的沖擊波脈寬明顯比迎爆面測點(diǎn)處的脈寬大。無論爆距與板厚如何改變，背爆面兩測點(diǎn)的λτD-E始終非常接近；并且隨著爆距與板厚的增加，λτD-E與界點(diǎn)處λτC的值趨于一致。迎爆測點(diǎn)λτB和臨界點(diǎn)λτC幾乎不發(fā)生變化。

圖12 λτ值折線圖Fig.12 The broken-line graph of λτ values

2.3.4 背爆面近壁面相鄰位置壓力分布特征

經(jīng)上述分析發(fā)現(xiàn)背爆面與迎爆面處的壓力分布特征存在較大差異，本節(jié)分析背爆面相鄰測點(diǎn)處的壓力分布特征。將D、E2處的峰值之比和脈寬之比在不同爆距與板厚下的變化趨勢繪于圖13中。

圖13 PmD/PmE與τD/τE值折線圖Fig.13 The broken-line graph of PmD/PmE and τD/τE values

可以發(fā)現(xiàn)，隨著爆距的增大，PmD/PmE從1.16降低至0.99，在爆距為3時，測點(diǎn)E的壓力峰值超過了測點(diǎn)D的壓力峰值，脈寬之比τD/τE從0.79降至0.97。

3 結(jié)論

1)水下裝藥爆炸后釋放大量能量，以沖擊波的形式向外傳播，當(dāng)沖擊波傳播到結(jié)構(gòu)物表面時，出現(xiàn)壓力峰值；一部分沖擊波會反射回流場形成反射沖擊波，入射波傳播到臨界點(diǎn)時會改變其傳播方向，形成繞射波到達(dá)平板背爆面。

2)迎爆面測點(diǎn)處壓力迅速達(dá)到峰值，壓力峰值大但脈寬較小；圓柱殼頂部及背爆面處測點(diǎn)壓力上升緩慢，壓力值小但脈寬較大。迎爆面的壓力峰值相比與背爆面大致相差一個數(shù)量級，而背爆面的沖擊波脈寬時間能達(dá)到迎爆面的數(shù)倍。

3)迎爆面測點(diǎn)的壓力分布特征與背爆面測點(diǎn)的壓力分布特征存在明顯差異，臨界點(diǎn)的壓力特征與背爆面壓力特征更為相似。無論爆距與板厚如何改變，背爆面各處的壓力峰值與脈寬不隨空間位置的變化而變化。

4)隨著爆距和板厚的增加，背爆面的沖擊波脈寬λτD-E與臨界點(diǎn)處的沖擊波脈寬λτC的值趨于一致。