水面艦船附近氣泡運(yùn)動(dòng)數(shù)值模擬

王詩(shī)平 張阿漫 程瀟歐 姚熊亮 陳海龍

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

0 引言

當(dāng)在自由面與水中結(jié)構(gòu)物附近產(chǎn)生氣泡時(shí)，需要同時(shí)計(jì)入兩者對(duì)氣泡脈動(dòng)的影響［1-2］。水下爆炸被廣泛應(yīng)用于國(guó)防和工業(yè)等領(lǐng)域，其生成氣泡的運(yùn)動(dòng)特性由于浮力和邊界的特性而非常復(fù)雜。當(dāng)氣泡在剛性壁面附近運(yùn)動(dòng)時(shí)，剛性壁面會(huì)對(duì)氣泡產(chǎn)生吸引而誘導(dǎo)其產(chǎn)生朝向剛性壁面方向的射流；當(dāng)氣泡在自由面附近運(yùn)動(dòng)時(shí)，自由面會(huì)對(duì)氣泡產(chǎn)生排斥而使射流方向遠(yuǎn)離自由面，浮力的存在會(huì)使氣泡的運(yùn)動(dòng)形式更加復(fù)雜［3］。在數(shù)值模擬中描述在壁面、自由面和浮力共同作用下氣泡運(yùn)動(dòng)特性的研究較少，為此，本文將以計(jì)入自由面效應(yīng)的格林函數(shù)替代基本格林函數(shù)，采用邊界元法研究水面艦船附近氣泡的運(yùn)動(dòng)特性，同時(shí)計(jì)入浮力的影響。

在對(duì)自由面進(jìn)行建模時(shí)，需要將自由面區(qū)域設(shè)置得較大，或者將自由面建成封閉的四面體，這就導(dǎo)致自由面網(wǎng)格劃分較多［4］，嚴(yán)重影響計(jì)算速度。因此，本文將采用修正格林函數(shù)，在不對(duì)自由面進(jìn)行建模的情況下考慮自由面對(duì)氣泡脈動(dòng)的影響。

在考慮水下爆炸氣泡與水面艦船的相互作用時(shí)，由于藥包在水下爆炸生成的氣泡與自由面較近，因此需要考慮自由面對(duì)氣泡運(yùn)動(dòng)的影響，此時(shí)，在建立氣泡運(yùn)動(dòng)模型時(shí)，需要同時(shí)考慮氣泡表面、自由面和艦船水線面以下的部分。如此，就會(huì)產(chǎn)生一系列的問(wèn)題：首先，自由面需要大量的三維網(wǎng)格，這樣會(huì)增加數(shù)值計(jì)算量；其次是自由面與結(jié)構(gòu)的交界，這不僅會(huì)導(dǎo)致數(shù)值的不穩(wěn)定，而且在處理結(jié)構(gòu)與自由面的交界處時(shí)會(huì)帶來(lái)復(fù)雜的數(shù)值步驟。因此，本文將采用一種新的方法計(jì)算自由面與水面艦船附近氣泡的運(yùn)動(dòng)，通過(guò)在氣泡和水面艦船鏡像位置分布點(diǎn)匯，使得在計(jì)算過(guò)程中避免了對(duì)自由面建模的必要，從而提高計(jì)算速率。

1 數(shù)值模型及有效性驗(yàn)證

假定水下爆炸氣泡周圍的流場(chǎng)為無(wú)粘、無(wú)旋且不可壓縮的理想流體，滿足拉普拉斯方程［5-7］：

下面，建立此鏡像法的數(shù)學(xué)模型。按照此方法，只需建立水面艦船的水下部分，艦船水面以上的部分則可忽略。圖1所示為水面艦船（水下部分）、爆炸氣泡及其鏡像，其中水面艦船的船長(zhǎng)為L(zhǎng)，船寬為B，吃水為H。在數(shù)值計(jì)算中，將水面艦船及其鏡像劃分為3個(gè)區(qū)域：水面艦船自由面以下部分（S＋）、水線面結(jié)點(diǎn)（I）以及水面艦船自由面以下部分的鏡像（S-）。氣泡表面及氣泡表面的鏡像分別表示為B＋和B-。在此種情況下，可將拉普拉斯方程表示成如下邊界積分方程：

式中，rp為固定點(diǎn)；rq為積分點(diǎn)；λ為流場(chǎng)中的立體角。方程（2）聯(lián)系速度勢(shì)Φ和法向速度

式中，r′q表示rq關(guān)于水平面z=0的鏡像點(diǎn)。則此時(shí)邊界積分方程可表示為：

從方程（4）中可以看出，采用考慮自由面效應(yīng)的格林函數(shù)取代基本格林函數(shù)后，邊界積分方程得到了簡(jiǎn)化，僅需在船體水線以下部分和氣泡表面分布源和偶極即可得到計(jì)入自由面效應(yīng)的邊界積分方程。

氣泡脈動(dòng)過(guò)程滿足如下伯努利方程：

式中，δ為無(wú)量綱浮力參數(shù)，伯努利方程的無(wú)量綱化參考。本文將均采用無(wú)量綱變量，上述方程的無(wú)量綱化過(guò)程參考文獻(xiàn)［9］。

圖1 艦船模型和水下部分及其鏡像Fig.1 Ship model’s underwater part and its image

采用此格林函數(shù)可以計(jì)及自由面對(duì)氣泡動(dòng)態(tài)特性的影響，下面將檢驗(yàn)該方法在不建立自由面網(wǎng)格的基礎(chǔ)上考慮自由面效應(yīng)的有效性。采用本文的格林函數(shù)運(yùn)用三維模型進(jìn)行計(jì)算，并將計(jì)算結(jié)果與自由面附近軸對(duì)稱氣泡運(yùn)動(dòng)特性的數(shù)值結(jié)果進(jìn)行了比較。本文選取初始無(wú)量綱距離d＝1.4，強(qiáng)度參數(shù) ε＝200，浮力參數(shù) δ＝0，比熱比 γ＝1.25，對(duì)應(yīng)氣泡初始半徑 R0＝0.116。將氣泡表面劃分成642個(gè)結(jié)點(diǎn)，1 280個(gè)單元。兩種結(jié)果的對(duì)比如圖2所示。圖2（a）～圖2（c）是采用考慮了自由面效應(yīng)的三維氣泡脈動(dòng)結(jié)果，在計(jì)算過(guò)程中沒(méi)有對(duì)自由面進(jìn)行建模，圖中，z＝0的位置即為自由面的位置。圖 2（e）～圖 2（f）是在相同參數(shù)下采用軸對(duì)稱程序的計(jì)算結(jié)果。對(duì)于三維數(shù)值模擬，氣泡初始為球形，然后迅速膨脹，在t＝0.75時(shí)刻膨脹至最大體積（圖2（a）），隨后，氣泡開始坍塌，并在 t＝1.53時(shí)刻形成向下的射流（圖2（b）），在t＝1.64時(shí)刻，形成扁平的寬射流（圖2（c））。在整個(gè)計(jì)算過(guò)程中，氣泡的整體位置是向下移動(dòng)的，其計(jì)算結(jié)果與圖2（d）～圖2（f）中列出的軸對(duì)稱結(jié)果相似。由圖中可看到，采用考慮了自由面效應(yīng)的格林函數(shù)的結(jié)果與采用基本格林函數(shù)的二維軸對(duì)稱氣泡的結(jié)果基本相同，氣泡都在自由面的Bjerknes力作用下產(chǎn)生了遠(yuǎn)離自由面的射流。在基于不同方法的氣泡脈動(dòng)過(guò)程中，無(wú)量綱體積隨無(wú)量綱時(shí)間的變化曲線如圖3所示。由圖3可看出，采用考慮了自由面效應(yīng)的格林函數(shù)計(jì)算的體積與軸對(duì)稱模型吻合較好，驗(yàn)證了本文數(shù)值方法的有效性。

在處理水面艦船問(wèn)題時(shí)，水面艦船被自由面分割成水上部分和水下部分，忽略水上部分后，邊界不再封閉，因此，在計(jì)算奇異積分時(shí)，Hii=4π-Hij的性質(zhì)不再成立。邊界上立體角的求解可以通過(guò)如下2種方法求解。

圖2 自由面附近單個(gè)氣泡三維結(jié)果與軸對(duì)稱結(jié)果對(duì)比Fig.2 Comparison of three dimensional and axisymmetrical single bubble near free surface

圖3 三維與軸對(duì)稱氣泡無(wú)量綱體積隨無(wú)量綱時(shí)間的變化曲線Fig.3 Time history of bubble dimensionless volume of 3D and axisymmetrical result

1）每個(gè)小的球三角形（例如，球三角形NBC）的面積可以通過(guò)下面的積分求得［10］：

其中φ0和θ的關(guān)系可以從圖4中找出。圖中，點(diǎn)O′為邊界上的結(jié)點(diǎn)，ON是點(diǎn)O′的單位法向量，弧AB，BC，CD，DE和EA分別是以 O′點(diǎn)為原點(diǎn)的單位圓與線性單元（延長(zhǎng)線）的交界線。這樣，點(diǎn)O′的立體角λ便可以表示為弧ANB，BNC等球三角形的面積之和。

2）邊界立體角可通過(guò)如下方法求解［11］：

式中，Nb為結(jié)點(diǎn) p周圍三角形的個(gè)數(shù)；δi為兩個(gè)相鄰三角形單元的角度。

2 數(shù)值結(jié)果

本文計(jì)算氣泡的代碼采用Fortran語(yǔ)言編寫，采用有限元軟件ABAQUS進(jìn)行后處理，以實(shí)現(xiàn)計(jì)算過(guò)程可視化。考慮初始質(zhì)量為N kg TNT的藥包安置在船體正下方0.4L的位置。氣泡脈動(dòng)過(guò)程如圖5所示。圖5（a）為初始?xì)馀荩跉馀輧?nèi)部充滿了高溫高壓氣體，隨后氣泡迅速膨脹，在船體表面中部壓力最大。在t＝0.51時(shí)刻，氣泡體積達(dá)到最大然后便開始坍塌，在t＝1.04時(shí)刻，氣泡開始產(chǎn)生向上的射流，在t＝1.08時(shí)刻，環(huán)狀氣泡形成，并在t＝1.11時(shí)刻體積達(dá)到最小。不同時(shí)刻艦船表面的壓力分布如圖6所示。由圖中可看出，在t＝0的初始時(shí)刻和環(huán)狀氣泡體積最小的時(shí)刻（t＝1.11）附近，短時(shí)間內(nèi)船體表面承受的垂向載荷數(shù)值很大，對(duì)艦船表面產(chǎn)生了類似“中拱”的巨大力矩，這時(shí)，結(jié)構(gòu)處于危險(xiǎn)狀態(tài)。在其余的絕大部分時(shí)間內(nèi)，艦船表面的垂向載荷均很小。

圖5 N kg TNT藥包距離水面0.4L爆炸生成氣泡的動(dòng)態(tài)特性Fig.5 Bubble dynamics generated by N kg TNT at 0.4L depth

圖6 艦船表面沿船長(zhǎng)方向的無(wú)量綱垂向載荷分布Fig.6 Dimensionless vertical load distributing on ship surface

艦船表面中點(diǎn)的氣泡載荷時(shí)歷曲線如圖7所示。艦船遭受氣泡作用的載荷與氣泡本身的運(yùn)動(dòng)息息相關(guān)，在氣泡脈動(dòng)的絕大部分時(shí)間內(nèi)，氣泡為近似球形，氣泡脈動(dòng)速度較慢，艦船表面始終維持著低于靜水壓力的壓強(qiáng)值；在氣泡初始時(shí)刻和環(huán)狀氣泡體積最小時(shí)刻附近，艦船遭受著危險(xiǎn)載荷，且射流在沖破氣泡壁的瞬間還會(huì)對(duì)艦船表面產(chǎn)生較大的射流沖擊載荷，計(jì)算值為14.83。射流沖擊主要是對(duì)結(jié)構(gòu)造成局部性損傷，并有可能在流場(chǎng)中產(chǎn)生空化，射流沖擊的特性還有待進(jìn)一步的探索。

圖7 艦船表面中點(diǎn)無(wú)量綱垂向載荷時(shí)歷曲線Fig.7 Dimensionless vertical load of bubble on middle ship surface

當(dāng)水下爆炸不是發(fā)生在艦船正下方，例如，發(fā)生在舷側(cè)與水面附近時(shí)（圖8），舷側(cè)的存在使得氣泡產(chǎn)生朝向舷側(cè)曲面的射流，而自由面的存在又使得氣泡產(chǎn)生朝向遠(yuǎn)離自由面的射流，在這兩股影響下，氣泡便產(chǎn)生斜向下的射流，具體實(shí)驗(yàn)過(guò)程參見文獻(xiàn)［12］。因此，當(dāng)藥包在自由面和舷側(cè)附近爆炸時(shí)，氣泡的攻擊效果會(huì)被削弱。下面，將模擬水下爆炸氣泡在舷側(cè)和自由面附近的情況，在計(jì)算過(guò)程中，采用鏡像方法考慮了自由面對(duì)氣泡脈動(dòng)的影響。

圖8 自由面附近氣泡對(duì)舷側(cè)的攻擊效果Fig.8 Attack effect of bubble near shipboard and free surface

初始質(zhì)量為N kg TNT的藥包安置在船體正下方0.3L，左側(cè)0.55B的位置處。氣泡脈動(dòng)過(guò)程如圖9所示，其中，圖中云圖表示壓力。此工況為靠近水面艦船和自由面的情況，艦船和自由面對(duì)氣泡脈動(dòng)的影響都很大。在此算例中，藥包距離船體較近，在氣泡膨脹至最大體積的過(guò)程中，艦船對(duì)氣泡有阻礙作用，而在氣泡坍塌的過(guò)程中，船體對(duì)其的吸引較劇烈，氣泡在靠近船體一側(cè)較大，形成了遠(yuǎn)離自由面和靠近艦船的斜射流。

初始質(zhì)量為N kg TNT的藥包安置在船體正下方0.3L、左側(cè)0.8B的位置處。氣泡脈動(dòng)過(guò)程如圖10所示，其中，圖中云圖表示壓力。在此算例中，藥包距離船體較遠(yuǎn)，受船體的吸引較小，氣泡在膨脹至最大體積附近時(shí)受到了艦船輕微的排斥，而在坍塌過(guò)程中則受到了艦船的吸引。與上一算例相比，氣泡被船體吸引較小，氣泡射流向下。

如果水下爆炸氣泡出現(xiàn)在水中結(jié)構(gòu)物附近，通常會(huì)形成一股朝向水中結(jié)構(gòu)物的射流，這股射流非常大（速度將達(dá)到100 m/s，甚至更大），并有可能會(huì)對(duì)水中結(jié)構(gòu)物造成毀傷。根據(jù)氣泡與水中結(jié)構(gòu)物的相對(duì)位置，由于近自由面和浮力的影響，氣泡射流可能會(huì)作用在水中結(jié)構(gòu)物上，也有可能不會(huì)直接作用在結(jié)構(gòu)上。

3 結(jié) 語(yǔ)

對(duì)于水面艦船，自由面的存在不僅增加了網(wǎng)格的數(shù)量，而且在自由面與結(jié)構(gòu)的交界處也增加了數(shù)值處理的難度。本文采用考慮了自由面效應(yīng)的格林函數(shù)來(lái)取代基本格林函數(shù)，計(jì)入自由面的存在對(duì)氣泡運(yùn)動(dòng)的影響，并將計(jì)算結(jié)果與軸對(duì)稱數(shù)值模擬結(jié)果進(jìn)行了對(duì)比，數(shù)值結(jié)果吻合良好。

通過(guò)數(shù)值模擬和實(shí)驗(yàn)發(fā)現(xiàn)，當(dāng)藥包在舷側(cè)和自由面附近爆炸時(shí)，氣泡射流可能不會(huì)完全作用在艦船上，自由面效應(yīng)的存在削弱了氣泡的打擊能力。因此，在指導(dǎo)常規(guī)武器攻擊水面艦船時(shí)，為充分發(fā)揮氣泡射流的威力，應(yīng)盡量避免其在自由面附近爆炸。

［1］KLASEBOER E，KHOO B C，HUNG K C.Dynamics of an oscillating bubble near a floating structure［J］.Journal of Fluid and Structures，2005，21（4）：395-412.

［2］COLE R H.Underwater explosion［M］.New Jersey：Princeton University Press，1948：118-127.

［3］張阿漫，王詩(shī)平，汪玉，等.水下爆炸對(duì)艦船結(jié)構(gòu)損傷特性研究綜述［J］.中國(guó)艦船研究，2011，6（3）：1-7.ZHANG A M，WANG S P，WANG Y，et al.Advances in the research of characteristics of warship structural damage due to underwater explosion［J］.Chinese Journal of Ship Research，2011，6（3）：1-7.

［4］張阿漫，姚熊亮.近自由面水下爆炸氣泡的運(yùn)動(dòng)規(guī)律研究［J］.物理學(xué)報(bào)，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.

［5］WANG Q X，YEO K S，KHOO B C，et al.Unstructured MEL modelling of unsteady nonlinear ship waves［J］.Journal of Computational Physics，2005，210 （1）183-224.

［6］ZHANG Y L，YEO K S，KHOO B C，et al.3D jet impact and toroidal bubbles［J］.Journal of Computational Physics，2001，166（2）：336-360.

［7］宗智，何亮，孫龍泉.水下爆炸氣泡對(duì)水面艦船載荷的數(shù)值研究［J］.船舶力學(xué)，2008，12（5）：733-739.ZONG Z，HE L，SUN L Q.Numerical study of loading on the surface ship near an underwater explosion bubble［J］.Journal of Ship Mechanics，2008，12（5）：733-739.

［8］WANG C，KHOO B C.An indirect boundary element method for three-dimensional explosion bubbles［J］.Journal of Computational Physics，2004，194（2）：451-480.

［9］戴遺山，段文洋.船舶在波浪中運(yùn)動(dòng)的勢(shì)流理論［M］.北京：國(guó)防工業(yè)出版社，2008.

［10］ZHANG Y L，YEO K S，KHOO B C，et al.Three-dimensional computation of bubbles near a free surface［J］.Journal of Computational Physics，1998，146（1）：105-123.

［11］WANG Q X.Numerical modelling of violent bubble motion［J］.Physics of Fluids，2004，16（5）：1610-1619.

［12］張阿漫，王詩(shī)平，白兆宏，等.不同環(huán)境下氣泡脈動(dòng)特性實(shí)驗(yàn)研究［J］.力學(xué)學(xué)報(bào)，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］.Acta Mechanica Sinica，2011，43（1）：71-83.