王詩平 張阿漫 程瀟歐 姚熊亮 陳海龍

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

0 引言

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

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

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

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

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

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

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

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

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

氣泡脈動過程滿足如下伯努利方程：

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

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

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

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

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

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

1）每個小的球三角形（例如，球三角形NBC）的面積可以通過下面的積分求得［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)的單位圓與線性單元（延長線）的交界線。這樣，點(diǎn)O′的立體角λ便可以表示為弧ANB，BNC等球三角形的面積之和。

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

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

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

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

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

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

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

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

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

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

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

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

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

3 結(jié) 語

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

通過數(shù)值模擬和實(shí)驗(yàn)發(fā)現(xiàn)，當(dāng)藥包在舷側(cè)和自由面附近爆炸時，氣泡射流可能不會完全作用在艦船上，自由面效應(yīng)的存在削弱了氣泡的打擊能力。因此，在指導(dǎo)常規(guī)武器攻擊水面艦船時，為充分發(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］張阿漫，王詩平，汪玉，等.水下爆炸對艦船結(jié)構(gòu)損傷特性研究綜述［J］.中國艦船研究，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)動規(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.

［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］宗智，何亮，孫龍泉.水下爆炸氣泡對水面艦船載荷的數(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)動的勢流理論［M］.北京：國防工業(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］張阿漫，王詩平，白兆宏，等.不同環(huán)境下氣泡脈動特性實(shí)驗(yà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］.Acta Mechanica Sinica，2011，43（1）：71-83.