蘇怡然，唐文勇，鄭紹文，李德聰

(1.上海交通大學(xué) 海洋工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，上海 200240;2.中國(guó)艦船研究設(shè)計(jì)中心，武漢 430064)

在實(shí)際工程中，氣泡脈動(dòng)引發(fā)沖擊作用的實(shí)例有很多:如水下爆炸氣態(tài)產(chǎn)物對(duì)船體結(jié)構(gòu)造成沖擊、螺旋槳工作時(shí)產(chǎn)生的毫米量級(jí)的空化氣泡對(duì)螺旋槳的剝蝕，以及超聲輻照產(chǎn)生的微米量級(jí)氣泡對(duì)微生物的沖擊作用等。當(dāng)氣泡距離結(jié)構(gòu)足夠近時(shí)，其脈動(dòng)過(guò)程往往伴隨著高速水射流［1］。水射流產(chǎn)生后會(huì)穿透氣泡并對(duì)附近結(jié)構(gòu)造成沖擊，此過(guò)程為氣泡破壞效應(yīng)的主要機(jī)理之一。

關(guān)于水下氣泡的運(yùn)動(dòng)特性，早期研究［1－2］主要基于一些理想假設(shè)，如氣泡呈球狀變化等。近年來(lái)，隨著計(jì)算機(jī)技術(shù)的發(fā)展，氣泡非球狀演化過(guò)程的數(shù)值模擬已成為可能。邊界積分方法假設(shè)流場(chǎng)無(wú)粘、無(wú)旋且不可壓縮，使用低維模型模擬高維空間，有效減小了計(jì)算量，在模擬水下爆炸氣泡射流現(xiàn)象中已有許多成功的應(yīng)用［3－7］。但上述假設(shè)限制了邊界積分法的應(yīng)用范圍:如螺旋槳空泡或超聲激發(fā)的氣泡，氣泡尺寸往往較小，此時(shí)粘性對(duì)氣泡脈動(dòng)的影響顯著;再如水下爆炸氣泡回彈階段，激波引起的能量耗散將對(duì)氣泡演化規(guī)律造成影響，在該過(guò)程的數(shù)值模擬中流體可壓縮性將不可忽略;此外，在射流穿透氣泡后，氣泡表面速度勢(shì)變?yōu)槎嘀岛瘮?shù)，計(jì)算將被終止，需引入環(huán)狀氣泡模型才能對(duì)該氣泡進(jìn)行模擬。因此，使用更具一般性的Navier-Stokes方程(以下簡(jiǎn)稱NS方程)模擬氣泡脈動(dòng)具有實(shí)際意義。

針對(duì)水下氣泡的非球狀脈動(dòng)現(xiàn)象，Popinet等［8］使用MAC法成功模擬了毫米量級(jí)氣泡的射流過(guò)程，Hao等［9］將該法拓展到三維直角坐標(biāo)系中，Johnsen、McKee等［10－11］實(shí)現(xiàn)了可壓縮流體中射流現(xiàn)象的數(shù)值模擬，張凌新等［14］使用VOF法模擬了水下爆炸氣泡的潰滅現(xiàn)象。這些方法均需要將流體邊界取在遠(yuǎn)離氣泡的位置，導(dǎo)致計(jì)算效率低下，尤其對(duì)于一些高強(qiáng)度氣泡，巨大的計(jì)算耗時(shí)將使數(shù)值模擬缺乏實(shí)用性。

本文在MAC法［13］基礎(chǔ)上提出一種適用于模擬氣泡脈動(dòng)問(wèn)題的近場(chǎng)壓力邊界條件，通過(guò)合理減小計(jì)算區(qū)域，提高了MAC法的計(jì)算效率，解決了高強(qiáng)度氣泡數(shù)值模擬中計(jì)算耗時(shí)較大的問(wèn)題。通過(guò)引進(jìn)Popinet［8］和Hirt［12］等對(duì)自由液面的處理方法，保證了自由液面邊界條件的計(jì)算精度和數(shù)值穩(wěn)定。

1 數(shù)學(xué)模型

1.1 控制方程

對(duì)于水下氣泡首次周期內(nèi)的脈動(dòng)問(wèn)題，使用軸對(duì)稱不可壓縮NS方程進(jìn)行描述:

其中:Φ=p/ρ代表壓力，ρ為液體密度，g為重力加速度，γ為運(yùn)動(dòng)粘性系數(shù)，剪應(yīng)力分量計(jì)算方法如下:

軸對(duì)稱模型如圖1所示，方程需滿足的邊界條件為:

對(duì)歐拉域外圍邊界Sp使用壓力邊界條件，邊界上的壓力值由遠(yuǎn)端流場(chǎng)確定。在求解遠(yuǎn)端壓力分布時(shí)，本文將脈動(dòng)氣泡假想為空間源匯，使用勢(shì)流理論計(jì)算遠(yuǎn)場(chǎng)壓力，詳見(jiàn)本文第2章。

圖1 水下氣泡數(shù)值模型Fig.1 Numerical model for the simulation of underwater bubbles

1.2 離散方法

使用交錯(cuò)網(wǎng)格對(duì)NS方程進(jìn)行離散，壓力節(jié)點(diǎn)和速度節(jié)點(diǎn)分別位于網(wǎng)格中心和網(wǎng)格邊長(zhǎng)中點(diǎn)。NS方程在空間上采用中心差分離散格式，在時(shí)間上使用一階隱式歐拉法進(jìn)行離散。

使用預(yù)測(cè)－校正法解決不可壓縮流動(dòng)的壓力－速度耦合問(wèn)題［15－16］:在每個(gè)時(shí)間步內(nèi)，首先忽略壓力梯度項(xiàng)，通過(guò)動(dòng)量方程求出預(yù)測(cè)速度，進(jìn)而得到壓力泊松方程，使用逐次超松弛迭代法對(duì)該方程進(jìn)行求解，最后用壓力梯度對(duì)預(yù)測(cè)速度進(jìn)行校正。

在處理自由液面問(wèn)題時(shí)，本文采用改進(jìn)MAC法的思想，通過(guò)標(biāo)記點(diǎn)追蹤自由液面運(yùn)動(dòng)，實(shí)現(xiàn)了對(duì)自由液面的精確捕捉。首先在自由液面上布置一串具有先后順序的無(wú)質(zhì)量標(biāo)記點(diǎn)，并在每個(gè)時(shí)步內(nèi)使標(biāo)記點(diǎn)跟隨流場(chǎng)運(yùn)動(dòng)，從而實(shí)現(xiàn)自由液面形狀和位置的追蹤。標(biāo)記點(diǎn)的運(yùn)動(dòng)速度由四個(gè)相鄰節(jié)點(diǎn)速度進(jìn)行雙線性插值得出。在計(jì)算過(guò)程中，需要不斷進(jìn)行標(biāo)記點(diǎn)位置的重分配，以保證標(biāo)記點(diǎn)沿液面均勻分布，從而降低計(jì)算誤差。自由液面將歐拉網(wǎng)格劃分為流體網(wǎng)格、液面網(wǎng)格和空網(wǎng)格三類，其中空網(wǎng)格不參與計(jì)算。

液面法向邊界條件見(jiàn)式(11)，式中σ為表面張力系數(shù)，κ為液面曲率，ΦInt為自由液面附近流場(chǎng)壓力，ΦB為氣泡內(nèi)壓力。傳統(tǒng)處理方法簡(jiǎn)單地認(rèn)為邊界單元壓力等于液面壓力，這樣處理在邊界法向壓力梯度較小時(shí)可滿足精度要求，但對(duì)于本文研究的高壓力梯度問(wèn)題，會(huì)造成計(jì)算誤差與數(shù)值波動(dòng)。為使計(jì)算誤差得到控制，采用Hirt等［14］提出的邊界單元壓力計(jì)算方法:認(rèn)為壓力在邊界附近呈線性分布，邊界單元壓力由相鄰流體單元壓力與液面壓力插值得到。

液面切向邊界條件見(jiàn)式(12)。根據(jù)Popinet等［8］提出的方法實(shí)現(xiàn)該條件，即認(rèn)為自由液面運(yùn)動(dòng)速度是由周邊流體網(wǎng)格處流速和空網(wǎng)格虛擬速度插值得到，通過(guò)人為指定虛擬速度的大小，可使流場(chǎng)滿足切向邊界條件。于是，該問(wèn)題轉(zhuǎn)化為條件極值問(wèn)題:

2 近場(chǎng)壓力邊界

對(duì)無(wú)窮流域問(wèn)題進(jìn)行數(shù)值模擬時(shí)，普遍需要通過(guò)較大的模型范圍來(lái)逼近無(wú)窮遠(yuǎn)邊界，計(jì)算效率較低。尤其對(duì)于高強(qiáng)度氣泡，脈動(dòng)過(guò)程中最大半徑可達(dá)初始半徑的30余倍，巨大的計(jì)算量將使數(shù)值模擬失去實(shí)用性。本文根據(jù)Best等［17］得到的規(guī)律，提出一種由近場(chǎng)壓力邊界描述無(wú)窮流域的數(shù)值算法，有效縮小了模型尺寸，提高了計(jì)算效率。

2.1 邊界壓力算法

Best等［17］提出，距氣泡一定距離之外的流場(chǎng)可近似看作發(fā)散狀。因此，將氣泡假想為空間點(diǎn)源(匯)，可認(rèn)為遠(yuǎn)端流場(chǎng)由該點(diǎn)源(匯)激發(fā)。設(shè)V為氣泡體積，Q為點(diǎn)源(匯)強(qiáng)度，r為空間點(diǎn)到氣泡中心的距離，遠(yuǎn)端流場(chǎng)速度勢(shì)為:

當(dāng)氣泡靠近壁面時(shí)，流場(chǎng)會(huì)受到鏡面效應(yīng)影響。此時(shí)，需在氣泡鏡像點(diǎn)處施加一個(gè)等強(qiáng)度的點(diǎn)源(匯)。設(shè)r1和r2分別表示空間點(diǎn)至氣泡中心和鏡像位置的距離，遠(yuǎn)端流場(chǎng)速度勢(shì)為:

在忽略遠(yuǎn)端粘性影響后，流場(chǎng)壓力可由伯努利方程求出。設(shè)Φ∞為無(wú)窮遠(yuǎn)處壓力(靜水壓)，近場(chǎng)邊界壓力ΦBC可由式(16)求出:

將ΦBC作為近場(chǎng)邊界的壓力值，即可減小數(shù)值模擬所需的模型尺寸。

2.2 穩(wěn)定性證明

首先，對(duì)源強(qiáng)度施加小量誤差δQ，該誤差引起的近場(chǎng)邊界壓力誤差為:

由于時(shí)間步長(zhǎng)Δt足夠小，式中第二項(xiàng)可以忽略。設(shè)R與·R·分別為氣泡等效半徑R對(duì)時(shí)間的一階和二階導(dǎo)數(shù)，則有 δ(·R·)? δ (R·)? δ(R)。根據(jù)Rayleigh-Plesset方程的積分形式:

忽略無(wú)窮小量后可得:

忽略無(wú)窮小量后可得:

可見(jiàn)，當(dāng)近場(chǎng)邊界距氣泡中心距離大于兩倍氣泡最大半徑時(shí)，近場(chǎng)壓力邊界條件已具有穩(wěn)定性。

2.3 邊界合理性判斷與數(shù)值震蕩抑制

當(dāng)標(biāo)記點(diǎn)移動(dòng)至新網(wǎng)格時(shí)，其速度將由新的一組節(jié)點(diǎn)速度插值得到，這會(huì)引起標(biāo)記點(diǎn)速度的微幅波動(dòng)。由于時(shí)間步長(zhǎng)很短，該波動(dòng)對(duì)于氣泡尺寸的影響可以忽略，但卻會(huì)使產(chǎn)生較大誤差，從而直接影響近場(chǎng)邊界壓力的準(zhǔn)確性。若不加控制，還易引發(fā)數(shù)值震蕩，導(dǎo)致計(jì)算失敗。

3 可靠性驗(yàn)證

氣泡的球狀演化規(guī)律由 Rayleigh-Plesset方程給出:

其中:R為氣泡等效半徑，γ為運(yùn)動(dòng)粘性系數(shù)，S為表面張力系數(shù)，p∞為靜水壓力，pB為氣泡壓力。氣體狀態(tài)方程為:

式中:k為氣體絕熱指數(shù)，取1.4。

從本質(zhì)上講，Rayleigh-Plesset方程是NS方程在理想條件下的簡(jiǎn)化形式，因此，算法可靠性的最基本要求即為可逼近Rayleigh-Plesset經(jīng)典解。對(duì)忽略重力作用與壁面作用的氣泡脈動(dòng)過(guò)程進(jìn)行數(shù)值模擬，考慮不同強(qiáng)度的氣泡(設(shè)氣泡強(qiáng)度系數(shù)為其壓力峰值與環(huán)境壓力的比值，這里取 50，100，150，200)，模型參數(shù)見(jiàn)表 1，計(jì)算結(jié)果如圖2所示。在強(qiáng)度系數(shù)為50時(shí)，氣泡半徑誤差為0.2%;隨著強(qiáng)度增高，半徑誤差略微變大，這是因?yàn)樵谀Ｐ统叽绮蛔兊那闆r下，隨著氣泡體積變大，近場(chǎng)邊界距氣泡距離不斷減小。計(jì)算結(jié)果表明，各項(xiàng)誤差均在可接受范圍內(nèi)，其中半徑誤差小于1%，脈動(dòng)周期誤差約為0.5% ～0.6%。

表1 數(shù)值模擬主要參數(shù)Tab.1 Parameters of numerical model

圖2 近場(chǎng)邊界方法與R－P方程計(jì)算結(jié)果對(duì)比Fig.2 Result comparisons between the near-field boundary method and the R-P equation

無(wú)窮流域問(wèn)題數(shù)值模擬的傳統(tǒng)方法需使用較大的模型范圍，且模型邊界條件設(shè)為恒定靜水壓力。為驗(yàn)證本文所述邊界條件在提高計(jì)算效率、降低計(jì)算誤差中的效果，表2將兩種邊界條件的誤差進(jìn)行了對(duì)比。對(duì)于不同模型尺寸，恒壓邊界均會(huì)產(chǎn)生更大的誤差，尤其是脈動(dòng)周期誤差。隨著模型尺寸的增大，兩種算法的誤差均有所減小，其中近場(chǎng)壓力邊界的誤差維持在0.5%左右，而恒壓邊界的誤差則從36%快速下降，但即使模型范圍達(dá)到10倍氣泡半徑，恒壓邊界誤差依然大于近場(chǎng)壓力邊界誤差。這說(shuō)明本文所述邊界壓力算法在提高計(jì)算效率的同時(shí)也降低了計(jì)算誤差。

表2 近場(chǎng)邊界與恒壓邊界誤差比較Tab.1 Error comparison between near-field boundary and constant pressure boundary

然而，上述模型中自由液面均未發(fā)生變形。為驗(yàn)證算法處理液面變形問(wèn)題的可靠性，本文對(duì)重力作用下的氣泡脈動(dòng)過(guò)程進(jìn)行數(shù)值模擬，并與實(shí)驗(yàn)數(shù)據(jù)進(jìn)行了對(duì)比。Shima等［18］采用脈沖放電的方法產(chǎn)生水下氣泡，并提出由氣泡平衡半徑Re來(lái)確定氣泡壓力:

對(duì)于該實(shí)驗(yàn)涉及的氣泡尺寸，流體粘性效應(yīng)影響可以忽略，但由于粘性作用顯著的微米量級(jí)氣泡的脈動(dòng)實(shí)驗(yàn)難以進(jìn)行，所以這里僅將上述實(shí)驗(yàn)數(shù)據(jù)同數(shù)值仿真結(jié)果進(jìn)行對(duì)比，作為證明算法可靠性的依據(jù)，計(jì)算結(jié)果如圖3所示。在首次脈動(dòng)周期內(nèi)，數(shù)值模擬結(jié)果與實(shí)驗(yàn)數(shù)據(jù)吻合良好，氣泡半徑平均誤差為氣泡最大半徑的3.91%，脈動(dòng)周期誤差為2.68%;當(dāng)氣泡進(jìn)入回彈階段時(shí)，數(shù)值模擬得到的氣泡等效半徑明顯大于實(shí)驗(yàn)結(jié)果。這是因?yàn)樵诨貜楇A段，多種物理過(guò)程產(chǎn)生的能量耗散效應(yīng)逐漸顯著，如氣泡界面上熱交換、可壓縮流體激波等。

Klaseboer等［4］在實(shí)驗(yàn)水池中進(jìn)行了不同裝藥量的爆炸實(shí)驗(yàn)，并進(jìn)行高速攝影，圖4給出了爆炸氣泡脈動(dòng)過(guò)程與數(shù)值模擬結(jié)果的對(duì)比，兩者在趨勢(shì)上吻合較好。從數(shù)值模擬結(jié)果可以看出，在潰滅階段，受重力的作用，射流產(chǎn)生于氣泡底端并逐漸發(fā)展至穿透氣泡，射流方向與重力方向相反。實(shí)驗(yàn)中監(jiān)測(cè)了爆炸氣泡附近壓力場(chǎng)隨時(shí)間的變化，由近到遠(yuǎn)三個(gè)監(jiān)測(cè)點(diǎn)分別記為P1、P3、P5。圖5對(duì)比了三個(gè)監(jiān)測(cè)點(diǎn)壓力的實(shí)驗(yàn)數(shù)據(jù)和數(shù)值計(jì)算結(jié)果。可以看出，兩者趨勢(shì)一致，誤差保持在15%以內(nèi)。產(chǎn)生該誤差的主要原因是試驗(yàn)水池邊界和數(shù)值模型邊界之間的差異。數(shù)值模擬中沒(méi)有考慮池壁、水面等邊界影響。此外，由于算法沒(méi)有考慮流體壓縮性，所以不同監(jiān)測(cè)點(diǎn)同時(shí)達(dá)到壓力峰值，而非試驗(yàn)中先后達(dá)到峰值。

圖3 重力作用下氣泡等效半徑隨時(shí)間變化Fig.3 Time history of bubble equivalent radius under the influence of gravity

圖4 水下爆炸氣泡脈動(dòng)過(guò)程(虛線為數(shù)值模擬結(jié)果)Fig.4 Pulsation of underwater explosive bubbles(The dash line represents numerical results)

圖5 爆炸氣泡附近壓力隨時(shí)間變化Fig.5 Time evolution of the pressure near explosive bubble

4 近壁氣泡脈動(dòng)過(guò)程分析

近壁氣泡脈動(dòng)過(guò)程往往伴隨著射流現(xiàn)象，由于國(guó)內(nèi)外公開(kāi)的相關(guān)實(shí)驗(yàn)結(jié)果多為高速攝影，定量分析較為匱乏，故本文僅對(duì)典型近壁氣泡脈動(dòng)過(guò)程進(jìn)行定性分析。Tomita等［19］使用脈沖放電的方法制造出近壁氣泡，并進(jìn)行高速攝影，如圖6(a)所示，圖6(b)為該過(guò)程的數(shù)值模擬計(jì)算結(jié)果。數(shù)值模型中，氣泡位于7.5 m水深，強(qiáng)度系數(shù)取50，初始半徑為0.15 m，氣泡距壁0.7 m。

在膨脹階段，氣泡始終保持近似球形，并被固壁略微排斥，膨脹階段停止于約0.5倍脈動(dòng)周期;隨后進(jìn)入潰滅階段，該階段內(nèi)氣泡受固壁的強(qiáng)烈吸引而產(chǎn)生大幅位移和變形，在潰滅階段尾期，射流形成于氣泡遠(yuǎn)離固壁端;在124 ms時(shí)刻射流擊穿氣泡，使氣泡變?yōu)榄h(huán)狀。

射流形成時(shí)刻流線圖與壓力云圖的計(jì)算結(jié)果如圖7所示。此時(shí)，氣泡壓力大于靜水壓力，氣泡處于減速收縮階段，且射流后方存在局部高壓區(qū)域;流場(chǎng)邊界處流動(dòng)近似發(fā)散狀，這與Best［17］提出的流速規(guī)律是相吻合的。

圖6 近壁氣泡脈動(dòng)過(guò)程Fig.6 Bubble pulsation in the vicinity of a rigid boundary

圖7 射流形成時(shí)刻流線圖與壓力云圖Fig.7 Streamlines and pressure profiles at jet formation

5 結(jié)論

本文基于NS方程，提出一種近場(chǎng)壓力邊界算法，有效提高了計(jì)算效率。該方法適用于水下爆炸氣泡脈動(dòng)、螺旋槳空泡潰滅等多種過(guò)程的數(shù)值模擬和相應(yīng)射流沖擊載荷的預(yù)報(bào)。本文主要結(jié)論如下:

(1)本文所述近場(chǎng)壓力邊界允許通過(guò)較小區(qū)域的流體模型模擬無(wú)窮流域問(wèn)題，大幅提高了計(jì)算效率;

(2)MAC方法適于追蹤具有較大曲率的自由液面，允許人為的指定自由液面拓?fù)湫螤畹母淖儯谒職馀莘乔驙蠲}動(dòng)過(guò)程的模擬中表現(xiàn)出色;

(3)本文旨在闡明近場(chǎng)壓力邊界，所以文中主要模擬了氣泡首次脈動(dòng)周期內(nèi)的近似不可壓過(guò)程，關(guān)于氣泡回彈時(shí)流體壓縮效應(yīng)和液面熱交換作用的影響有待進(jìn)一步研究。

［1］Blake J R，Gibson D C.Cavitation bubbles near boundaries［J］.Annual Review of Fluid Mechanics，1987，19(1):99－123.

［2］Plesset M S，Prosperetti A.Bubble dynamics and cavitation［J］.Annual Review of Fluid Mechanics，1977，9(1):145－185.

［3］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.

［4］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.

［5］ Zhang A M，Yao X L，Yu X B.The dynamics of threedimensional underwater explosion bubble［J］.Journal of Sound and Vibration，2008，311(4):1196 －1212.

［6］張阿漫，姚熊亮.水下爆炸氣泡與海底相互作用研究［J］.振動(dòng)與沖擊，2008，27(3):92 －98.ZHANG A-man， YAO Xiong-liang. Interaction between underwater explosion bubble and seabed［J］.Journal of Vibration and Shock，2008，27(3):92－98.

［7］宗 智，何 亮，張恩國(guó).水中結(jié)構(gòu)物附近三維爆炸氣泡的數(shù)值模擬［J］.水動(dòng)力學(xué)研究與進(jìn)展 A輯，2007，22(5):592－602.ZONG Zhi， HE Liang， ZHANG En-guo. Numerical simulation of a three-dimensional underwater explosion bubble near a structure［J］.Journal of Hydrodynamics Ser A，2007，22(5):592－602.

［8］ Popinet S，Zaleski S.Bubble collapse near a solid boundary:a numerical study of the influence of viscosity［J］.Journal of Fluid Mechanics，2002，464:137 －163.

［9］ Hao Y，Prosperetti A.A numerical method for threedimensional gas liquid flow computations［J］.Journal of Computational Physics，2004，196(1):126 －144.

［10］ Johnsen E， Colonius T. Numerical simulations of nonspherical bubble collapse［J］.Journal of Fluid Mechanics，2009，629:231 －262.

［11］ McKee S，Tome M F，F(xiàn)erreira V G，et al.The MACmethod［J］.Computers＆ Fluids，2008，37(8):907－930.

［12］ Nichols B D，Hirt C W.Improved free surface boundary conditions for numerical incompressible-flow calculations［J］.Journal of Computational Physics，1971，8(3):434 －448.

［13］Hawker N A，Ventikos Y.Interaction of a strong shockwave with a gas bubble in a liquid medium:A numerical study［J］.Journal of Fluid Mechanics，2012，701:59 －97.

［14］張凌新，尹 琴，邵雪明.水中氣泡潰滅的理論與數(shù)值研究［J］.水動(dòng)力學(xué)研究與進(jìn)展 A輯，2012，27(1):68－73.ZHANG Ling-xin，YIN Qin，SHAO Xue-ming.Theoretical and numerical studies on the bubble collapse in water［J］.Chinese Journal of Hydrodynamics Ser.A，2012，27(1):68－73.

［15］ Griffith B E.An accurate and efficient method for the incompressible Navier Stokes equations using the projection method as a preconditioner［J］.Journal of Computational Physics，2009，228:7565－7595.

［16］ Brown D L，Cortez R，Minion M L.Accurate projection methods for the incompressible Navier Stokes equations［J］.Journal of Computational Physics，2001，168(2):464－499.

［17］ Best J P.The dynamics of underwater explosions［D］.PhD dissertation. The University of Wollongong， New South Wales，1991.

［18］ Shima A，Tomita Y，Sendai.The behavior of a spherical bubble near a solid wall in a compressible liquid［J］.Ingenieur-Archiv，1981，51:243－255.

［19］Tomita Y，Shima A.High-speed photographic observations of laser-induced cavitation bubbles in water［J］.Acustica，1990，71(3):161－171.