基于AUTODYN的氣泡與固定壁面相互作用數(shù)值模擬

張 偉 岳永威 張阿漫 孫龍泉

1 中國艦船研究設(shè)計中心，湖北武漢 430064

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

0 引 言

水下爆炸氣泡引起的結(jié)構(gòu)破壞可分為3種：爆炸氣泡脈動激發(fā)船體梁總體振動，造成整體失穩(wěn)甚至斷裂失效；遠(yuǎn)場爆炸時，氣泡脈動引起艦船上較敏感設(shè)備的共振，造成設(shè)備破壞；當(dāng)炸藥近場爆炸時，氣泡受艦船結(jié)構(gòu)邊界的影響，形成沖擊射流，造成艦艇結(jié)構(gòu)局部損傷。第3種情況屬氣泡近壁面運(yùn)動規(guī)律問題，進(jìn)行理論研究的依據(jù)主要是以勢流理論建立的水平及垂直剛性面附近在浮力作用下運(yùn)動的氣泡理論模型。該模型基本能反映水下爆炸氣泡和周圍流體介質(zhì)的運(yùn)動規(guī)律，但其忽略了邊界對氣泡形狀的影響，較適于遠(yuǎn)場氣泡脈動分析。在試驗研究方面，關(guān)于水下爆炸氣泡對結(jié)構(gòu)的毀傷作用試驗研究多采用規(guī)則結(jié)構(gòu)或縮比模型，鮮有實船試驗。

近年來，由水下爆炸引起的氣泡動力學(xué)問題成為海軍艦船生命力技術(shù)領(lǐng)域關(guān)注的重點。然而，水下爆炸氣泡從形成、膨脹到最終潰滅是一個復(fù)雜的物理演化過程，尤其是氣泡在運(yùn)動過程中與周圍結(jié)構(gòu)的作用受許多因素的影響，研究難度較大。目前，我國學(xué)者主要是以高速攝像的方法對電火花生成的氣泡進(jìn)行觀測，進(jìn)而對氣泡的運(yùn)動規(guī)律進(jìn)行研究［1-2］，而基于數(shù)值平臺的仿真模擬則開展得相對較少。

AUTODYN［3］是一種顯式有限元分析程序，主要用于解決固體、流體、氣體及其相互作用的高度非線性動力學(xué)問題。目前，我國關(guān)于AUTODYN軟件在水下爆炸領(lǐng)域的應(yīng)用與其他商業(yè)軟件比相對較少。就水下爆炸［4-6］而言，由于爆轟產(chǎn)物和水都屬于流體，所以相對于Lagrange描述方法，在AUTODYN中采用Euler方法描述更為方便、有效。數(shù)值模擬的過程就是對守恒方程，包括連續(xù)性方程、動量方程及能量方程的離散，單純的守恒方程無法求解動力學(xué)問題，必須與狀態(tài)方程（EOS）聯(lián)立，才能構(gòu)成封閉方程組，繼而對流體動力學(xué)問題進(jìn)行求解。

1 有效性驗證

為了驗證AUTODYN軟件的有效性，使用球?qū)ΨQ計算模塊分析0.229 kg的TNT在178.6 m水深處爆炸時的相關(guān)數(shù)據(jù)，并將計算值與文獻(xiàn)［7］中的實驗值進(jìn)行了對比。在深水爆炸過程中，靜水壓力梯度可以忽略不計，取氣泡周圍的靜水壓力一定，水的計算域取為50 m。表1和表2所示為流場中距藥包中心0.5 m和1 m處沖擊波和氣泡壓力峰值、氣泡最大半徑及脈動周期的計算值與實驗值的對比。

表1 沖擊波及氣泡壓力峰值的AUTODYN計算值與實驗值對比Tab.1 Comparison between experimental data and the calculated results about shock and bubble pressure peak

表2 氣泡最大半徑及脈動周期的計算值與實驗值對比Tab.2 Comparison between experimental data and the calculated results about the maximum bubble radius and pulse cycle

從表中可以看出，數(shù)值模擬的氣泡脈動最大半徑和脈動周期與實驗值間的誤差約為10%，誤差在可接受的范圍內(nèi)。一維球?qū)ΨQ求解器可有效模擬水下爆炸氣泡的脈動。

下面將以重力場為例進(jìn)行對比分析。在模擬重力場中水下爆炸氣泡的運(yùn)動時，采用AUTODYN軟件中獨(dú)有的映射技術(shù)，將一維球?qū)ΨQ計算結(jié)果映射至二維軸對稱求解器，從而解決了網(wǎng)格尺寸過小、計算時間過長的問題。設(shè)置長、寬、高分別為18 m，18 m和7 m的流場，以35 g藥量在流場中心下3.5 m處引爆，觀察該工況下氣泡在重力場中的運(yùn)動規(guī)律。在流場中預(yù)設(shè)A，B，C等3個測點以便測量氣泡在運(yùn)動過程中的流場壓力，它們分別位于爆心水平方向0.7 m處；爆心下方水平方向0.7 m、垂向0.71 m處；爆心上方水平方向0.7 m、垂向1.095 m處。圖1所示為測定A，B，C的數(shù)值模擬壓力時歷曲線，通過與文獻(xiàn)［7］中的相似工況及測點壓力曲線進(jìn)行對比，發(fā)現(xiàn)曲線的時間發(fā)展趨勢以及壓力峰值基本吻合，進(jìn)一步驗證了AUTODYN在模擬氣泡在重力場中運(yùn)動的精確性。

圖1 流場中測點的壓力數(shù)值模擬曲線Fig.1 Bubble pressure of flow field numerical simulation curve of three measuring points

圖1給出了80～83 ms時間段的壓力值，此時，氣泡收縮至最小體積，流場中輻射氣泡二次壓力波。由圖可見，氣泡二次壓力波的峰值和持續(xù)時間模擬值與實驗值吻合良好，平均誤差在10%以內(nèi)，但壓力峰值發(fā)生的時間略有提前。究其原因，可能是數(shù)值上雖然采用了邊界處理，但有限的邊界仍會對氣泡脈動產(chǎn)生影響，使氣泡的周期較實驗值略小。

2 近固壁面的水下爆炸氣泡射流

2.1 概 述

氣泡與壁面的相互作用一直是研究人員關(guān)注的問題。處于固壁面附近的氣泡在受到壁面Bjerknes力［8］的同時還受重力的作用，為此，設(shè)置不同的無量綱距離參數(shù)［9］（爆心距壁面的距離與氣泡最大半徑的比值），將炸藥置于剛固平板的下方，如圖2所示。

圖2 氣泡在固壁面下的計算模型Fig.2 Model of bubble under a solid boundary

水下爆炸氣泡的模擬涉及到的流體材料包括空氣、炸藥和水，下面將分別介紹這3種材料的狀態(tài)方程。

2.1.1 空氣的狀態(tài)方程

計算中可能會涉及到自由面，因此，需要對空氣進(jìn)行模擬。對空氣進(jìn)行模擬采用理想氣體狀態(tài)方程（Ideal Gas）：

式中，ρ=1.225×10-3g/cm3；γ為絕熱指數(shù)，γ=1.4；e為比內(nèi)能，e=2.06785e5kJ/kg。

2.1.2 炸藥的狀態(tài)方程

此處，選取TNT炸藥材料模型作為水下爆炸氣泡形成的條件。采用JWL狀態(tài)方程描述炸藥的爆轟過程：

需要說明的是，當(dāng)炸藥膨脹到相當(dāng)?shù)捏w積時，JWL方程右端的前兩項可以忽略，此時，可以用理想氣體的狀態(tài)方程模擬炸藥的行為：p=ρ(γ-1)e，其中γ=ω+1。

2.1.3 水的狀態(tài)方程

AUTODYN的材料庫中，水的狀態(tài)方程有2種，多項式（Polynomial）狀態(tài)方程和沖擊（shock）狀態(tài)方程。由于需要考慮靜水壓力，因此，本文選用多項式狀態(tài)方程進(jìn)行計算：

2.2 不同距離固壁面下的氣泡形狀

為考察固壁面下不同距離的氣泡形狀，計算50 kg TNT在75 m水深處爆炸時的相關(guān)數(shù)據(jù)，使用STEEL材料模擬固壁，計算中采用映射技術(shù)。設(shè)置無量綱距離γf=0.6，0.8，1.0，1.5，2.0，2.5，3.0，5.0。圖3～圖6所示為距離參數(shù)γf=0.6，1.0，2.0，3.0時氣泡形狀的演變。

圖3所示為γf=0.6時的氣泡形狀圖。在膨脹階段，氣泡受到壁面的阻礙，在103.5 ms時氣泡體積達(dá)到最大，其上端成扁平狀；在收縮階段，氣泡上部仍被吸附在固壁表面，在183.5 ms時下端產(chǎn)生指向壁面的射流，196.5 ms時射流沖破氣泡，213 ms時體積達(dá)到最小。

圖4所示為γf=1.0時的氣泡形狀圖。值得注意的是，在收縮階段，氣泡上部距固壁表面的距離變近，這主要是受浮力的影響。

圖5～圖6所示為γf=2.0，3.0時的氣泡形狀圖。可以看出，距離壁面較遠(yuǎn)時，壁面的影響減弱，γf=2.0與γf=3.0時的氣泡形狀并沒有明顯的差別。因γf=5.0時的氣泡形狀變化與自由場中的相似，因而在此不再給出。

圖3 γf=0.6時的氣泡形狀演變Fig.3 Evolution of bubble shape with the distance parameterγf=0.6

圖5 γf=2.0時的氣泡形狀演變Fig.5 Evolution of bubble shape with the distance parameterγf=2.0

圖6 γf=3.0時的氣泡形狀演變Fig.6 Evolution of bubble shape with the distance parameterγf=3.0

將自由場氣泡與壁面下方的氣泡進(jìn)行對比可見，在距離參數(shù)較大（γf=2.0，3.0，5.0）的工況下，氣泡形狀的演變過程與自由場相似；距離壁面越近，射流的寬度越大。此外，距離參數(shù)不同，氣泡的脈動周期及射流發(fā)生的時間也不同，下面，將對脈動周期及射流時間進(jìn)行分析。

2.3 固壁面對氣泡最大半徑、脈動周期及射流時間的影響

藥包距壁面的距離不同，受到壁面的影響也不同，最直觀的體現(xiàn)可參見氣泡第1次脈動周期及最大半徑，此外，射流沖擊時間也有所不同。第1次脈動周期指的是氣泡第1次收縮至最小體積，將要進(jìn)行第2次膨脹的時刻，射流沖擊時刻指的則是射流將氣泡的上下表面擊穿的時刻。

50 kg TNT在75 m水深處爆炸，按照經(jīng)驗公式計算，可得氣泡第1次脈動周期為182.6 ms，第1次膨脹的最大半徑為2.85 m。下面，將給出不同距離參數(shù)壁面下氣泡的第1次膨脹最大半徑、第1次脈動周期及射流沖擊時刻的數(shù)值模擬值，如表3所示。

表3 不同距離參數(shù)壁面下的氣泡數(shù)值模擬值Tab.3 Numerical simulation results with different distance parameters

近距離（γf=0.6，0.8）的氣泡由于其形狀不再保持為球形，因而其最大半徑為最大體積的等效半徑。由表3可見，距離壁面越近，氣泡的最大半徑越小，距離壁面越遠(yuǎn)，其最大半徑就和自由場的相差越小。在γf=3.0的情況下，氣泡最大半徑已和自由場的相差無幾。總之，壁面的存在對氣泡的最大半徑影響不大，在γf=0.6時僅與自由場相差（2.842-2.821）/2.842=0.7%。

氣泡第1次脈動周期受壁面的影響較大，距離參數(shù)越小，脈動周期越大。在γf=3.0時，氣泡脈動周期已與自由場的值相差不大，這是因為壁面在氣泡膨脹過程中對其造成了阻礙作用，而在收縮時，因氣泡依附在壁面上，因此距壁面越近，氣泡脈動周期便越大。

為研究壁面的存在對射流時間的影響，使用氣泡脈動周期作為參考量將射流沖擊時間進(jìn)行無量綱化：Tj′=Tj/T，繪制其隨無量綱距離的變化而發(fā)生的變化，如圖7所示。圖中，圓圈為數(shù)值解；實線為自由場中的無量綱射流，時間Tj′=1.083；虛線代表射流沖擊時刻與氣泡脈動周期相同，Tj′=1 。

圖7 無量綱射流沖擊時間隨距離參數(shù)的變化曲線Fig.7 Variation of dimensionless jet shock time with respect to the distance parameter

由圖7可見，射流沖擊時刻與氣泡第1次脈動周期并非總是一致，但兩者相差不大，即射流沖擊是在氣泡第1次體積達(dá)到最小時刻附近發(fā)生的，此時也是氣泡開始輻射2次壓力波的時刻；對于壁面下的氣泡，其無量綱射流沖擊時間是隨著距離參數(shù)的增加而增加，這是因為壁面效應(yīng)是隨距壁面的距離增加而減弱，當(dāng)距離參數(shù)約為2.0（γf≈2.0）時，射流沖擊時刻與氣泡第1次脈動周期基本相等；當(dāng)距離參數(shù)大到一定（γf≥3.0）時，無量綱射流沖擊時間與自由場的就已非常接近。

2.4 近固壁面氣泡射流速度及壓力

首先，研究氣泡射流速度隨距離參數(shù)的變化規(guī)律。射流速度在射流形成過程中不斷提高。當(dāng)射流沖擊后，射流頂端的速度會有一定程度的損失，而氣泡射流沖擊后的射流速度是人們最為關(guān)心的。因此，定義射流沖擊完成后的射流頂端速度為Vj。圖8給出了γf=0.6，1.0，2.0，3.0，5.0時的射流頂端速度時歷曲線及重力場中的射流頂端速度的時歷曲線。

圖8 不同距離參數(shù)下的射流頂端速度時歷曲線Fig.8 Time history of the jet top velocity with different distance parameters

由圖8可見，在射流沖擊時刻，射流頂端速度急速增加，在達(dá)到最大值后開始下降。除γf=0.6的情況外，其余工況的射流頂端速度均趨于20 m/s。同時還可看出，距離參數(shù)越小，射流沖擊持續(xù)的時間越短，這是因為壁面的存在使得氣泡的射流在更短的時間內(nèi)完成。為研究距離參數(shù)對射流頂端速度的影響，繪制了射流頂端速度最大值隨距離參數(shù)的變化曲線，如圖9所示。

圖9 射流頂端速度峰值隨距離參數(shù)的變化曲線Fig.9 Variation of max top velocity with respect to the distance parameter

由圖9可見，射流頂端速度峰值隨距離參數(shù)的變化規(guī)律較復(fù)雜。在距離參數(shù)較小時，射流頂端速度峰值隨距離參數(shù)的增加而增加，在γf=2.0附近達(dá)到最大值，隨后又隨距離參數(shù)的增加而減小。當(dāng)距離參數(shù)較大時，射流頂端速度峰值與自由場的相差不大。據(jù)分析，距離參數(shù)不同，壁面對氣泡的作用也不同，氣泡射流沖擊發(fā)生的時間以及射流的形狀也不同。因此，當(dāng)距離參數(shù)很小時，雖然壁面的作用非常強(qiáng)烈，但射流沖擊發(fā)生時間相對較早，而且氣泡緊貼壁面，射流無法完全發(fā)展，形成的射流頂端速度峰值并不是最大的。隨著距離參數(shù)的略微增大，射流時間稍微推遲，因氣泡離壁面有一定的距離，此時射流能更進(jìn)一步發(fā)展，因此，射流頂端速度峰值不斷增加。當(dāng)γf=2.0時，射流時間與氣泡第1次脈動時間基本一致，此時射流發(fā)展最為完全，射流頂端速度峰值達(dá)到最大值。隨著距離參數(shù)的進(jìn)一步增大，壁面的影響削弱，在氣泡第1個脈動周期內(nèi)不再發(fā)生射流，因而射流頂端速度峰值也就不斷減小。

氣泡射流沖擊后，射流頂端速度經(jīng)過氣泡頂端與壁面間流體的衰減作用，即為作用在壁面上的射流速度Vjp。圖10所示為壁面附近0.1 m處射流速度峰值Vjp(max)隨距離參數(shù)的變化曲線。

圖10 壁面上射流速度峰值隨距離參數(shù)的變化曲線Fig.10 Variation of the peak velocity of jet near a solid boundary with respect to the distance parameter

由圖10可見，經(jīng)過流體的衰減作用，壁面射流速度峰值與氣泡頂端射流速度峰值隨距離參數(shù)的變化規(guī)律又有所不同。在小距離參數(shù)（γf=0.6，0.8）情況下，無流體的衰減作用，因此射流頂端速度即為壁面射流速度。但距離參數(shù)越大，流體的衰減作用便也越大，最后作用在壁面的射流速度就越小。當(dāng)γf＞3時，便可以認(rèn)為作用在壁面的已不再是射流，而是氣泡第2次膨脹引起的流體運(yùn)動。可見，壁面上的氣泡射流載荷是非常復(fù)雜的，氣泡攻擊壁面的威力也與多種因素有關(guān)。

圖11所示為γf=0.6時壁面附近0.1 m處的壓力時歷曲線。

圖11 γf=0.6時壁面附近0.1 m處的壓力時歷曲線Fig.11 Time history of pressure near a solid boundary about 0.1 m with the distance parameterγf=0.6

由圖11可見，當(dāng)水射流躍擊到壁板上時，速度降至零，同時在水中產(chǎn)生壓縮波，作用在平板上的壓力即為水錘壓力。水錘壓力為：

式中，Ph為水錘壓力；ρ為水的密度；C為水的聲速；Vj為作用在平板上的射流速度。按式（4）計算，γf=0.6時壁面的水錘壓力應(yīng)為145 MPa，數(shù)值模擬值為170 MPa，考慮到此測點并非壁面上的點，數(shù)值模擬值是在誤差允許的范圍內(nèi)。同時注意到，在射流壓力之后，程序可以捕捉到氣泡的二次壓力波，其峰值為20.56 MPa。

氣泡在發(fā)生射流之后，仍然會第2次膨脹，圖12所示為γf=0.6，1.5，3.0時爆心水平方向2.6 m處的氣泡脈動壓力時歷曲線，圖13所示為爆心水平方向2.6 m處氣泡脈動壓力峰值隨距離參數(shù)的變化曲線。

圖12 爆心水平方向2.6 m處氣泡脈動壓力時歷曲線Fig.12 Time history of bubble pulse pressure at 2.6 m from explosion center horizontally

圖13 爆心水平方向2.6 m處氣泡脈動壓力峰值隨距離參數(shù)變化曲線Fig.13 Variation of bubble pulse pressure peak with respect to distance parameter at 2.6 m from explosion center horizontally

由圖12和圖13可見，水平方向的氣泡脈動壓力峰值隨距離參數(shù)的變化規(guī)律與壁面附近射流速度的規(guī)律是相反的，在γf=0.8附近最小。這是因為氣泡射流消耗的能量越多，氣泡再次膨脹的能量就會越發(fā)縮小，相應(yīng)的脈動壓力也就越小。這也再次驗證了在γf=0.8時氣泡對壁面的射流作用最大。

3 結(jié) 論

本文運(yùn)用AUTODYN軟件分別模擬了球?qū)ΨQ氣泡、重力場中氣泡以及剛性和彈性壁面與氣泡之間的相互作用，計算結(jié)果與實驗數(shù)據(jù)的對比和分析表明：

1）運(yùn)用AUTODYN軟件計算得出的一維及重力場中水下爆炸氣泡的平均誤差在10%以內(nèi)，完全適用于氣泡動力學(xué)的研究。

2）剛性壁面的存在對氣泡的運(yùn)動特性具有較大影響。對于固壁面下方的水下爆炸氣泡，距離參數(shù)γf越小，壁面的影響就越大，氣泡射流的寬度越大，第1次脈動周期越大。

3）氣泡的無量綱射流沖擊時間隨距離參數(shù)γf的增加而增加，在γf≈2.0時，無量綱射流沖擊時間為1，氣泡射流發(fā)展最為完全，射流頂端速度峰值達(dá)到最大值。距離參數(shù)γf越小，氣泡射流沖擊持續(xù)的時間就越短。

4）壁面上的氣泡射流載荷非常復(fù)雜，與氣泡頂端射流速度、氣泡射流形狀以及射流頂端距離壁面的距離均有關(guān)。

［1］張阿漫，王詩平，白兆宏，等.不同環(huán)境下氣泡脈動特性實驗研究［J］.力學(xué)學(xué)報，2011，43（1）：71-83.ZHANG A M，WANG S P，BAI Z H，et al.Experimental study on bubble pulse features under different circumstances［J］.Chinese Journal of Theoretical and Applied Mechanics，2011，43（1）：71-83.

［2］張阿漫，姚熊亮.近自由面水下爆炸氣泡的運(yùn)動規(guī)律研究［J］.物理學(xué)報，2008，57（1）：339-353.ZHANG A M，YAO X L.The law of the underwater explosion bubble motion near free surface［J］.Acta Physica Sinica，2008，57（1）：339-353.

［3］VUYST T D，VIGNJEVIC R，CAMPBELL J C.Coupling between meshless and finite element methods［J］.International Journal of Impact Engineering，2005，31（8）：1054-1064.

［4］COLE R H.Underwater explosion［M］.Princeton：Princeton University Press，1948.

［5］GEERS T L，HUNTER K S.An integrated waveeffects model for an underwater explosion bubble［J］.The Journal of the Acoustical Society of America，2002，111（4）：1584-1601.

［6］BRUJAN E A，NAHEN K，SCHMIDT P，et al.Dynamics of laser-induced cavitation bubbles near an elastic boundary：influence of the elastic modulus［J］.Journal of Fluid Mechanics，2001，433：251-281.

［7］KLASEBOER E，HUNG K C，WANG C，et al.Experimental and numerical investigation of the dynamics of an underwater explosion bubble near a resilient/rigid structure［J］.Journal of Fluid Mechanics，2005，537：387-413.

［8］BJERKNES.Fields of force［M］.New York：Columbia University Press，1966.

［9］張阿漫.水下爆炸氣泡三維動態(tài)特性研究［D］.哈爾濱：哈爾濱工程大學(xué)，2006.