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

張 偉 岳永威 張阿漫 孫龍泉

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

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

0 引 言

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

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

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

1 有效性驗證

為了驗證AUTODYN軟件的有效性，使用球對稱計算模塊分析0.229 kg的TNT在178.6 m水深處爆炸時的相關數據，并將計算值與文獻［7］中的實驗值進行了對比。在深水爆炸過程中，靜水壓力梯度可以忽略不計，取氣泡周圍的靜水壓力一定，水的計算域取為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

從表中可以看出，數值模擬的氣泡脈動最大半徑和脈動周期與實驗值間的誤差約為10%，誤差在可接受的范圍內。一維球對稱求解器可有效模擬水下爆炸氣泡的脈動。

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

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

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

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

2.1 概 述

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

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

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

2.1.1 空氣的狀態方程

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

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

2.1.2 炸藥的狀態方程

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

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

2.1.3 水的狀態方程

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

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

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

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

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

圖5～圖6所示為γf=2.0，3.0時的氣泡形狀圖。可以看出，距離壁面較遠時，壁面的影響減弱，γ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

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

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

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

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

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

近距離（γf=0.6，0.8）的氣泡由于其形狀不再保持為球形，因而其最大半徑為最大體積的等效半徑。由表3可見，距離壁面越近，氣泡的最大半徑越小，距離壁面越遠，其最大半徑就和自由場的相差越小。在γf=3.0的情況下，氣泡最大半徑已和自由場的相差無幾?？傊诿娴拇嬖趯馀莸淖畲蟀霃接绊懖淮?，在γf=0.6時僅與自由場相差（2.842-2.821）/2.842=0.7%。

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

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

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

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

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

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

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

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

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

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

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

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

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

圖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可見，當水射流躍擊到壁板上時，速度降至零，同時在水中產生壓縮波，作用在平板上的壓力即為水錘壓力。水錘壓力為：

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

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

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

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

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

3 結 論

本文運用AUTODYN軟件分別模擬了球對稱氣泡、重力場中氣泡以及剛性和彈性壁面與氣泡之間的相互作用，計算結果與實驗數據的對比和分析表明：

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

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

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

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

［1］張阿漫，王詩平，白兆宏，等.不同環境下氣泡脈動特性實驗研究［J］.力學學報，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］張阿漫，姚熊亮.近自由面水下爆炸氣泡的運動規律研究［J］.物理學報，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］張阿漫.水下爆炸氣泡三維動態特性研究［D］.哈爾濱：哈爾濱工程大學，2006.