水下爆炸氣泡運動對水面形態影響的數值模擬

劉籌資，查恩堯，歐陽和平，程永舟，郭錦澤，鄭長青

(1.長沙理工大學水利工程學院，長沙 410114；2.中交第四航務工程勘察設計院有限公司，廣州 510000；3.湖南水利水電職業技術學院，長沙 410100；4.珠海爆破新技術開發公司，廣東 珠海 519099)

水下爆炸技術作為一種快速有效的水下施工手段，在航道疏浚[1]、水下巖石開挖[2]、水下結構拆除[3-4]等工程中應用廣泛。水下爆炸在不同條件下會形成多階段的物理現象，當水下爆炸發生在水面附近時，除沖擊波傳播、氣泡脈動現象外，氣泡在脈動過程中還會與水面相互作用在水面附近形成形態不同的水冢[5]。在重力作用下，水冢沖擊水面形成的波浪，會對周圍結構特別是艦船等浮體結構產生較大影響。因此開展水下爆炸氣泡運動對水面形態變化影響的研究，對深入了解近水面水下爆炸現象，指導水下爆炸實施時對周邊結構物的防護具有重要意義。

在爆炸興波方面，沈國光等[6-7]基于MAC算法將爆炸氣泡與自由液面的相互作用過程簡化，這一方法雖然能夠簡化水下爆炸的過程，但賦予氣泡的內壓是以與試驗水面現象的異同作為標準的，因此只能作為特定裝藥量的計算，此外由于模擬中忽略了氣泡脈動、射流的過程，對于興波機理的討論是不充分的。Safiyari等[8-9]基于SPH方法，將興波過程簡化為空腔的擴張、坍縮，進而引起水面波動的過程，這一研究同樣忽略了氣泡脈動、射流現象在水面形態形成中的關鍵作用。王志凱[10]基于ALE算法研究了水下爆炸氣泡破碎興波對浮動沖擊平臺的影響，該研究中雖然考慮了氣泡脈動，但其計算工況中，水深距離參數達到了2.7，依據文獻[11-13]，該條件下氣泡與自由液面相互作用形成的水冢對自由液面的擾動很小，能夠形成的興波波浪也較小，并且由于模型設置了浮體結構，水冢的發展以及興波波浪的形成也會被抑制。劉云龍等[12]基于BEM方法也討論了氣泡破碎興波對浮體結構的影響，但僅限于氣泡脈動成簡化的錐形水冢沖擊水面形成的波浪與浮體結構的相互作用。李汪諱等[14]、Wang等[15]基于BEM方法對水下爆炸興波現象進行了模擬，但文獻[14]未討論各種水冢形態對水面興波的影響，文獻[15]簡化了水冢形態，忽略不連續水體，將氣泡運動在水面附近形成的水冢分別簡化為水下空腔和水面水柱兩種形態。前人在興波方面的研究由于采用的BEM方法對不連續邊界的模擬存在困難，基本都建立在一定的水冢形態簡化條件下進行，這與實際的水面現象不相符，對于水面波動發展過程的討論也存在一定的誤差。

基于此，本數值模擬研究采用LS-DYNA非線性有限元軟件，基于ALE算法，結合水冢形成過程及其特征幾何參數、氣泡運動特性與水面波動過程，就不同水深距離參數下水下爆炸氣泡運動對水面形態的影響機理開展研究。

1 氣泡運動數值模型的建立

1.1 算法選擇與幾何模型的建立

由于炸藥爆炸的氣體產物在水下形成的氣泡存在著膨脹收縮過程，在氣泡運動的影響下水面會存在復雜的變形，同時為了捕捉氣泡的運動特性，對于氣泡邊界的捕捉也具有一定的要求。因此，本研究采用ALE算法進行三維數值模擬。模型中涉及到的炸藥、水和空氣均采用ALE算法。

為了兼顧計算效率，采用三維1/4模型，在對稱面上施加對稱邊界條件，在空氣域和水域外邊界處設置pressure outflow邊界單元層，以模擬無限邊界條件，模型的x、y方向和裝藥下方的水域深度應至少為氣泡脈動最大半徑的4倍，空氣層的高度則需滿足水冢中心水柱向上發展的要求。為了提高計算效率同時保證模擬結果的正確性，在氣泡最大半徑1.5倍范圍內采用與裝藥半徑尺寸相同的均勻網格，其他區域采用由均勻網格區域至流域外邊界逐漸漸變的網格，漸變網格最小處與均勻網格相同，最大處網格尺寸為均勻網格的10倍。幾何模型如圖1所示。

圖1 幾何模型Fig.1 Geometric model

1.2 材料模型的設置

模型主要涉及炸藥、水和空氣3種物質，其中炸藥采用JWL狀態方程進行描述，水采用Grüneisen狀態方程進行描述，空氣采用線性多項式狀態方程進行描述。各材料的具體參數如表1～表3所示[16-17]。

表1 TNT炸藥的材料方程和狀態方程參數取值

表2 水的材料方程和狀態方程參數取值

表3 空氣的材料方程和狀態方程參數取值

1.3 氣泡運動模型有效性驗證

為了驗證氣泡運動模型的有效性，設置與文獻[18]實驗相同的工況，使用5.2 g TNT的球形裝藥(與實驗中4 g PENT炸藥等效)，裝藥半徑為0.009 132 m，模型計算尺寸為1 m×1 m×4 m，漸變網格劃分尺寸最小為0.009 m，最大為0.1 m，均勻網格處網格尺寸為0.009 m。

由自由場條件下文獻[18]中實驗結果和數值模擬結果的對比(見圖2)可以發現，在自由場條件下數值模擬的氣泡運動趨勢與文獻中的各個階段保持一致；由在近自由液面條件下文獻[18]中實驗結果和數值模擬結果進行的對比(見圖3)可知，兩者氣泡運動趨勢大致一致。

圖2 自由場條件下氣泡運動趨勢對比Fig.2 Comparison of bubble movement trends under free field conditions

圖3 近自由液面條件下氣泡運動趨勢對比Fig.3 Comparison of bubble movement trends under the condition of near free liquid surface

對于水下爆炸氣泡脈動的第一周期及氣泡半徑計算，已經形成了能夠應用于工程實際的經驗公式[19]：

(1)

(2)

式中：R為氣泡第一次脈動的氣泡最大半徑，m；T為脈動周期，s；W為裝藥當量，kg；D為裝藥位置處水深，m；10.3為大氣壓強對應液柱高度，m；K4、K5為與炸藥相關的計算參數。

常見氣泡脈動參數取值如表4所示[19]。

表4 氣泡脈動經驗公式計算參數

另外也對自由場和近水面條件下的氣泡半徑時程進行了分析，數值模擬結果數據與文獻[18]中的工況數據進行對比分析發現兩者之間的誤差在10%以內(見表5)，為了體現研究的嚴謹性又對經驗公式[19]的結果進行了計算對比，分析發現在兩種條件下氣泡的半徑和周期變化誤差均在10%以內(見表5)，說明氣泡運動數值模型與實驗結果一致性良好。

表5 氣泡半徑與周期誤差計算結果

1.4 網格尺寸對計算精度的影響

為了討論網格尺寸大小對計算精度的影響，以5 kg TNT裝藥，在水面下1 m處爆炸為例，在理論氣泡脈動第一周期最大半徑的1.5倍范圍內分別劃分0.9、1.0、1.1、1.2倍裝藥半徑(0.09 m)的均勻網格，其余水域網格劃分由內至外逐漸增大，最大處為內部均勻網格尺寸的10倍。模型大小為20 m×20 m×30 m，水域高15 m，空氣域高15 m。

將計算所得的氣泡體積等效為規則圓的半徑并與經驗式(1)結果進行對比，得到了不同網格條件下的氣泡第一周期脈動最大半徑(見表6)，可以看出在幾種工況下，不同的網格尺寸的相對誤差均不超過5%。其中，當網格尺寸與裝藥半徑為0.9倍和1.0倍時，氣泡第一周期脈動最大半徑分別為2.629 m和2.633 m，兩者與經驗公式的理論半徑都較為接近，但采用0.9倍裝藥半徑與網格尺寸的計算時長相較于1.0倍裝藥半徑與網格尺寸約增加一半。在綜合考慮計算精度和計算效率的條件下，本文計算模型中心區域的均勻網格尺寸選取為1.0倍裝藥半徑大小。

表6 不同網格條件下氣泡參數計算結果

2 氣泡運動對水面形態的影響

2.1 模型設置

由于驗證模型采用的是5.2 g的球形TNT炸藥，不便對氣泡運動對水面形態變化的影響進行研究，因此建立尺寸為10 m×10 m×30 m的模型，水域、空氣域各高15 m，采用5 kg球形TNT裝藥。

2.2 裝藥深度對氣泡運動特性的影響

水下爆炸氣泡在不同裝藥水深條件下將會出現不同的運動特性，為了研究裝藥深度對其運動特性的影響，首先定義裝藥水深H為水下裝藥位置距水面的距離，炸藥起爆位置在水面以下1～5 m的模型中心處，設置9種不同工況，各個工況之間間隔0.5 m。

由不同裝藥深度條件下氣泡垂直方向位移的變化情況(見圖4)可以看出，在H為1.0～3.0 m時，氣泡垂直方向位移呈現先上后下的變化趨勢，且移動大小隨裝藥水深的增大而減小，當H>3.0 m時，氣泡垂直方向位移逐漸向上移動。

注：a～i分別代表H為1.0、1.5、2.0、2.5、3.0、3.5、4.0、4.5、5.0 m。圖4 水深對氣泡垂直方向位移的影響Fig.4 Influence of water depth on the vertical displacement of bubble

由不同裝藥深度下氣泡半徑的變化情況(見圖5)可以發現，氣泡脈動第一周期時長隨著裝藥深度的增大逐漸延長，從以上分析可知在2.0～3.0 m裝藥水深范圍內存在一個臨界水深。當裝藥水深小于臨界水深時，氣泡脈動第一周期最大半徑與脈動第一周期時長隨裝藥水深的增大而增大，氣泡中心點的位移先向上后向下，且移動大小隨裝藥水深增大而減小；當裝藥水深大于臨界水深時，氣泡脈動第一周期最大半徑與脈動第一周期時長隨水深的增大而減小，氣泡中心點向上位移大小隨水深增大而增大。

注：a～i分別代表H為1.0、1.5、2.0、2.5、3.0、3.5、4.0、4.5、5.0 m。圖5 不同裝藥深度氣泡半徑變化關系Fig.5 Change relationship of bubble radius with different charge depth

2.3 水冢形成過程

氣泡運動的變化也會導致水面水冢形態發生變化，因此對水冢的形成過程進行研究非常有必要。為了研究水冢的形成過程，對近水面條件下5 kg球形TNT炸藥在改變不同水深距離參數時產生的不同類型的水冢形成過程進行研究。同時為了與水面中心水柱相區分，定義在水面附近形成的水體形態稱之為水冢。經前人研究發現，水冢現象的形成與水深距離參數γf=H/Rmax有關[5]，H為炸藥裝藥位置距水面的距離，m；Rmax為氣泡脈動過程的最大半徑，m。依據水下爆炸氣泡最大半徑公式，在水深距離參數γf范圍在0.000 1～1.1內，討論不同水深距離參數γf下水冢形態的特征及形成機理(空腔半徑為氣泡脈動第一周期所形成的實際最大半徑)。

當起爆位置距水面相近時，會產生零碎型水冢(見圖6)，此時水深距離參數γf為0.000 195，裝藥水深0.000 4 m，氣泡膨脹形成的水下空腔半徑為2.053 5 m。炸藥起爆后，形成的氣泡突出水面與空氣相連通而破碎(見圖6a、圖6b)；在t=354 ms氣泡膨脹至最大半徑，形成了水下空腔，水面水體形成了零碎型水冢；在重力作用下空腔四周水體向內匯聚并形成匯聚水柱(見圖6d、圖6e)；t=2 094 ms時刻匯聚水柱高度到達最大值；隨后受重力作用沖擊水面形成波浪傳播。

圖6 零碎型水冢發展過程Fig.6 Development process of fragmented water mound

當水深距離參數γf增大到0.082 9，此時裝藥水深下降到0.215 m，氣泡膨脹形成的水下空腔半徑為2.589 m時會形成飛濺型水冢(見圖7)。炸藥起爆后，在t=6、54 ms時刻氣泡突出水面發生破碎，并在t=174 ms時刻達到最大半徑，形成空腔；t=294、594 ms時刻，由于氣泡的帶動，使得水體在空腔上部聚集，產生背向射流，形成飛濺型水冢。接著便在t=774、1 194 ms時刻產生類矩形的水面水柱，并在t=1 614 ms時刻達到最大值。

圖7 飛濺型水冢發展過程Fig.7 Development process of splashing water mound

接著當水深距離參數γf增加到0.549時，裝藥水深為1.5 m，氣泡脈動第一周期實際最大半徑為2.732 m會形成酒杯型水冢(見圖8)。在t=6 ms時刻，氣泡開始膨脹，并在t=118 ms時刻達到最大半徑；t=258 ms時刻，氣泡在收縮過程中，由于受到自由液面Bjerknes力的排斥作用，會形成向下的氣泡射流，進而帶動水流形成向上的射流，并逐漸發展為水面中心水柱；t=318 ms時刻，氣泡在向下射流的帶動下，分裂為兩部分，并在t=338 ms時刻達到第一周期的最小半徑；t=478 ms時刻氣泡再次膨脹，并且中部和上部氣泡相融合，同時在氣泡的推動和中心水柱的阻礙下，向上運動的水體在中心水柱周圍形成了環狀水體，環狀水體與中心水柱一起構成了酒杯型水冢。

圖8 酒杯型水冢發展過程Fig.8 Development process of wine glass water mound

最后當水深距離參數γf增大到0.998時，裝藥水深為2.692 m，此時氣泡脈動第一周期實際最大半徑為2.697 m，形成的水冢形態為皇冠型水冢(見圖9)。在t=6 ms時刻，氣泡開始膨脹，并在t=118 ms時刻達到最大半徑；t=258 ms時刻，氣泡在收縮過程中，由于受到自由液面Bjerknes力的排斥作用，會形成向下的氣泡射流，進而帶動水流形成向上的射流，并逐漸發展為水面中心水柱；t=318 ms時刻氣泡收縮至最小半徑，并在t=338 ms時刻開始二次脈動，同時向上的水射流也在不斷發展，并促進了中心水柱的發展，環狀水體也在氣泡膨脹的推動下開始突出水面向上發展，并與中心水柱一起構成了皇冠型水冢。

圖9 皇冠型水冢發展過程Fig.9 Development process of crown-shaped water mound

2.4 水冢氣泡運動特性

對上述4種不同類型的水冢氣泡運動特性進行分析，由于零碎型水冢和飛濺型水冢工況下的水下爆炸氣泡，炸藥裝藥水深較淺，不能形成穩定的氣泡脈動周期和射流，因此不對其運動特性進行討論。

對于酒杯型水冢和皇冠型水冢工況下的氣泡，裝藥水深能夠滿足氣泡脈動所需的條件，可以形成穩定的氣泡脈動與射流。通過兩者氣泡脈動半徑、周期的對比(見圖10a)可以看出，酒杯型水冢氣泡脈動第一周期的最大半徑為2.732 m，周期為336 ms；皇冠型水冢氣泡脈動第一周期的最大半徑為2.697 m，周期為360 ms。這表明氣泡最大半徑隨著水深的增大而減小，氣泡脈動周期則相反，這一現象與文獻[20]得到的結論相同。

從氣泡中心點在垂直方向的位移變化(見圖10b)可以看出，在膨脹階段，酒杯型水冢的氣泡中心點在垂直方向上的位移明顯高于皇冠型水冢，這是由于酒杯型水冢氣泡所處位置水深較淺，氣泡受到的外部靜水壓強也較小，因此膨脹階段的阻力也是偏小的。在收縮階段，由于受到Bjerknes力的作用形成向下射流，因此在射流和氣泡收縮趨勢的帶動下，氣泡中心位置不斷下移。

圖10 酒杯型水冢、皇冠型水冢氣泡運動特征對比Fig.10 Comparison of bubble movement characteristics between wine glass and crown-shaped water mound

2.5 不同類型水冢的特征參數對比

對比上述4種不同水冢的特征尺寸(見圖11)可知，幾種類型的水冢寬度基本處于4～5 m，與氣泡最大直徑相近，呈現出一定的相似性，結合上文對水冢發展過程進行詳細的研究，可以知道水冢的形成過程與氣泡的運動息息相關，氣泡的運動特性決定了水冢形態。

圖11 不同類型水冢特征尺寸對比Fig.11 Comparison of the characteristic sizes of different types of water mounds

零碎型和飛濺型水冢的最大寬度是氣泡在破碎前形成的空腔直徑。飛濺型水冢寬度大于零碎型水冢，這是由于零碎型水冢的起爆水深更淺，氣泡會更快速的膨脹到空氣中發生破碎。

酒杯型水冢和皇冠型水冢的水冢寬度，主要與氣泡進行第二周期脈動推動水面形成的射流位置相關，皇冠型水冢中心水柱的向上速度更低，形成的中心水柱底部寬度更大，使得氣泡第二周期脈動形成的環狀水體在更靠外側的位置形成，水冢寬度也因此更大。

對于零碎型水冢和飛濺型水冢的中心水柱高度，形成的主要原因是空腔四周水體在重力作用下向空腔底部匯聚形成的匯聚水柱。由圖6c及圖7c可知，飛濺型水冢空腔周圍水體范圍更大，且空腔高度更高，從能量角度分析，飛濺型水冢形成的匯聚水柱應該會比零碎型水冢匯聚水柱更高，但由于飛濺型水冢在空腔頂部形成的交匯向下射流與匯聚水流的運動方向相反，阻礙了匯聚水柱的形成，使得飛濺型水冢形成的中心水柱高度更低，這一現象也體現在飛濺型水冢的匯聚水柱的上下斷面寬度基本相同，而零碎型水冢的匯聚水柱為錐形。

對比酒杯型水冢和皇冠型水冢的水柱高度，由圖11可知，酒杯型水冢形成的中心水柱高度為35.213 m，皇冠型水冢形成的中心水柱高度僅為4.437 m，兩者相差較大。這是由于皇冠型水冢所處位置水深較大，向下氣泡射流的劇烈程度相較于酒杯型水冢是偏低的，因此形成的向上射流速度也較低。在向上射流形成的時刻，垂直方向距氣泡頂部1 m位置處，皇冠型水冢向上射流速度為9.27 m/s，而酒杯型水冢的向上射流速度為29.8 m/s，同時皇冠型水冢氣泡上部水體厚度更大，向上射流要穿透水體進入空氣中形成中心水柱所要消耗的能量也更大。綜上原因，在皇冠型水冢與酒杯型水冢的工況下雖然都能形成穩定的氣泡脈動、射流，但形成的中心水柱高度存在極大的差距。

2.6 不同類型水冢的水面波動特性

水下爆炸引起的水面波動在一定程度上會對水上結構物施加波浪壓力，對于類似核彈等超大裝藥量情況下的水下爆炸，形成的水面波動甚至可以等效為海嘯波。根據前文所形成的的水冢，分析不同類型水冢形成的水面波動特性。

距炸藥氣泡中心5 m處不同形態水冢形成的水面波動如圖12所示。對于零碎型水冢的水面波動，水冢形態形成的800 ms以內，由于氣泡破碎在100 ms左右已經破碎，自由液面只是在空腔膨脹過程中，水面出現了較低程度的上升。在1 000～2 200 ms左右的液面下降是空腔匯聚水柱開始形成，到2 650 ms匯聚水柱達到最高時刻之間，中心水柱不斷升高，四周水體向中心區域流動，在一定的時間延遲之后帶動5 m處自由液面下降。自2 650 ms水柱達到最高點之后，在重力作用下匯聚水柱開始沖擊水面形成水面波動，并持續向外傳播，2 200～4 000 ms水面上升，此時水面最大的上升高度為31.570 cm。

圖12 不同形態水冢形成的水面波動時程變化Fig.12 Time-history of water surface fluctuation formed by different forms of water mounds

對于飛濺型水冢的水面波動，在水冢形成階段，相較于零碎型水冢，飛濺型水冢對水面的影響更加明顯，這主要是由于其形成的空腔范圍更大，氣泡的膨脹帶動水體向上運動但并未破碎飛濺，形成了連續的水面隆起，帶動了監測點處的自由液面位置上升。在400～1 900 ms之間，監測點處自由液面下降，一方面是由于前一階段氣泡脈動破碎之后水面的隆起消失，更主要的是與零碎型水冢相同，空腔四周水體開始匯聚形成匯聚水柱，帶動周圍水體向中心匯集，使水面降低。在2 800 ms時刻匯聚水柱達到最大高度，此后開始在重力作用下沖擊水面，形成波浪向外傳播。監測點處的水面最大上升高度為64.060 cm，大致為零碎型水冢的2倍，雖然飛濺型水柱高度比零碎型水柱高度大0.325 m，兩者的體積分別為38.780 m3和15.450 m3，但是飛濺型水冢所具有的重力勢能更大，因此沖擊水面形成的興波波浪高度也更大。

對于酒杯型水冢的水面波動，在水冢形成初期，氣泡迅速膨脹頂起的水體在氣泡收縮階段的下落形成了初期的水面波動。在400 ms之后，隨著中心水柱和四周環形水體的上升，四周水體向中心流動，監測點的水面略微下降。隨后四周的環形水體和中心水柱根部附近的大量水體在重力作用下持續沖擊水面，形成了較大的水面波動，監測點最高的水面波動位移為83.085 cm。在最高水面波動通過監測點后，可以看到后續出現一個小的水面波動，以及一個持續的水面上升過程，則是殘余的部分環狀水體和部分下落的中心水柱沖擊水面形成的殘余水面波動。

對于皇冠型水冢的水面波動，與酒杯型水冢類似，在水冢形成初期，氣泡膨脹頂起的水面下墜形成了一定程度的水面波動。伴隨著水冢的形成，四周水流向中心流動，使得監測點處自由液面出現一定程度的下降。隨后在中心水柱和四周環狀水體形成的水冢下墜沖擊下，自由液面形成波動。皇冠型水冢在監測點處形成的最大水面波動位移為53.640 cm，相對酒杯型水冢較低，形成該現象的原因一方面是由于中心水柱的高度更低，僅為4.337 m。另一方面是由于氣泡多次脈動形成的自由液面破碎情況更加劇烈，后續的環狀水體與之前形成的環狀水體相互沖擊消耗了能量。從圖3中可以觀察到在最高水面波動通過監測點之后，還形成了兩個小的水面波動，前者是部分殘余的環狀水體沖擊水面形成的，后者是后續中心水柱上部水體沖擊水面形成的。

3 結論

1)隨著水深距離參數的增大，水下爆炸氣泡逐漸由不具備完整周期的破碎狀態變為具有脈動射流的完整周期狀態，同時由于Bjerknes力減小，氣泡向上運動趨勢逐漸增大。

2)隨水深距離參數的增大，水冢形態依次為零碎型、飛濺型、酒杯型、皇冠型；前兩者形成水冢形態與氣泡破碎時帶動的水體運動有關，水冢的特征尺寸受氣泡半徑大小和水冢形成機理影響；后兩者形成的水冢形態是氣泡脈動和射流現象疊加造成的，水冢的特征尺寸與氣泡大小及形成的射流速度相關。

3) 隨著水深距離參數的增大，零碎型水冢和飛濺型水冢的興波機理為：在重力作用下，氣泡破碎后空腔四周水體在空腔底部匯聚形成的匯聚水柱沖擊自由液面形成水面波動；酒杯型水冢和皇冠型水冢在重力作用下，由組成水冢的中心水柱和周圍環狀水體，沖擊自由液面形成水面波動。