Chinese Journal of Ship Research
近自由液面气泡与冲击波的相互作用
叶曦1，初文华2，3，4，陈林1，张阿漫5
1 200011中国船舶及海洋工程设计研究院，上海2 201306上海海洋大学 海洋科学学院，上海3 201306上海海洋大学 国家远洋渔业工程技术研究中心，上海4 201306上海海洋大学 大洋渔业资源可持续开发教育部重点实验室，上海5 150001哈尔滨工程大学 船舶工程学院，黑龙江 哈尔滨 摘 要：［目的］为了研究自由液面、气泡与冲击波三者之间的相互作用，［方法］基于间断迦辽金法，结合Level Set Real Ghost Fluid与 方法，分析复杂流场中冲击波传播特性及气泡运动特性，描述流场内各种波系的生成与发展过程。［结果］结果表明：在相互作用过程中，流场生成的复杂波系中包含多个稀疏波和冲击波。自由液面减缓了气泡的溃灭速度，而入射冲击波则加快了气泡的溃灭速度，并使自由液面的上拱运动增大。［结论］
所得结果可为水下爆炸对舰船结构的毁伤机理提供参考。关键词：气泡；自由液面；冲击波；间断迦辽金法中图分类号：U661.144 文献标志码：A DOI：10.3969/j.issn.1673-3185.2017.05.011
Interaction between bubble near free surface and shock wave
YE Xi1，CHU Wenhua2，3，4，CHEN Lin1，ZHANG Aman5 1 Marine Design and Research Institute of China，Shanghai 200011，China 2 College of Marine Sciences，Shanghai Ocean University，Shanghai 201306，China 3 National Engineering Research Center for Oceanic Fisheries，Shanghai Ocean University，Shanghai 201306，China 4 Key Laboratory of Sustainable Exploitation of Oceanic Fisheries Resources，Ministry of Education，Shanghai Ocean University，Shanghai 201306，China 5 School of Shipbuilding Engineering，Harbin Engineering University，Harbin 150001，China Abstract：［Objectives］This paper presents research into the interaction between a free surface，shock wave and bubble.［Methods］Based on the Discontinuous Galerkin method and combined with the Level Set method and Real Ghost Fluid method，the characteristics of shock wave propagation and bubble motion are analyzed，the generation and evolution of waves in fluid fields are described in detail.［Results］The results show the presence of complex waves in the fluid field after the interaction，including multiple rarefaction waves and shock waves. The collapse speed of the bubble is slowed down with the existence of the free surface，while the shock wave of the incident accelerates the collapse speed of the bubble and increases the up warp motion of the free surface.［Conclusions］The conclusions drawn from this paper can be used as reference points for further research into the damage mechanisms of ship structures subjected to underwater explosions. Key words：bubble；free surface；shock wave；Discontinuous Galerkin method
0引言
在船舶与海洋工程领域，经常会出现冲击波、气泡与自由液面间的相互作用，例如，水下爆炸或结构反射的冲击波与爆炸产生的气泡，以及船舶航行过程中在自由液面附近产生的空泡（船舶尾流中存在大量空泡）之间的相互作用。这种相互作用会改变冲击波的传播特性并影响气泡的运动特性，由此导致舰船结构载荷特性发生变化。因此，分析冲击波、自由液面以及气泡之间的相互作用对于深入研究水下爆炸工况的舰船结构毁伤机理具有重要意义。由水中气泡界面两侧介质差异造成的流场间Johnsen 等［1-2］断将给数值模拟带来一定的困难。采用体积分数法结合有限体积法进行加权本质无Weighted Essentially Non-Oscillatory，WE⁃振荡（ NO）格式重构，计算了水中气泡在强冲击波和弱冲击波作用下的演化历程，并3考虑了气泡下游存］采用水平集（Level在壁面的情况。王革和关奔［ Set）方法结合虚拟流体（Ghost Fluid）方法，并采用Euler有限体积法求解 方程组，分析了水中冲击波与气泡的相互作用。此外，自由拉格朗日法［4］、边界［5］ ［6-10］元法 、间断迦辽金法 等也适用于该物理现象的求解。然而，目前尚无自由液面、气泡与冲击波三者相互作用的研究成果。Level Set方法［11］、本文将基于间断迦辽金法，结合Real Ghost Fluid 方法［ 12 ］以及自适应网格加密（Adaptive Mesh Refinement，AMR ）技术［ 13 ］研究水下气泡、冲击波、自由液面之间的相互作用，分析自由液面对气泡运动及流场特性的影响，详细阐述新波系的产生与发展过程。
1 基本理论与方法
1如图 所示，计算模型采用 zor轴对称坐标系。在初始状态，气泡与自由液面相对静止，气泡半径为R，气泡中心与自由液面的初始距离为 d 。f强度为M的冲击波波前与自由液面平行，与气泡中心的距离为 d ，自下而上冲击气泡。本文所有s参数均进行无量纲化处理，参考量分别为气泡初始半径、标准大气压和空气密度。由于本文主要研究流场中的冲击波特性，故假设流场无粘且初Euler始无旋，并采用轴对称 方程描述流场： ¶E ¶F ¶G =S （1） + + ¶t ¶z ¶r
其中
E =[ ρ ρu ρv Σ]T
=[ F ρu ρu2 + p ρuv u ( Σ + p)]T =[ G ρv ρuv ρv 2 +p v ( Σ + p)]T [
S =- 1 ρu ρu2 ρuv u ( Σ + p)]T （2） z式（1）~式（2）中：t 为时间；ρ为流场密度；p为流场压力； u v 分别为 z r 轴方向的速度； ( 2) Σ = ρe + 0.5ρ u + v ，为总能量（内能与动能之2 in
和），其中 e 为单位质量的内能。in ［14］气泡内部采用理想气体状态方程描述： =( - 1) （3） p γ ρe g in =1.4，为气体比热比常数。式中：γ g Tait状态方程［6］进行描述：气泡外流场采用 （4） p = ρe ( N - 1)- NPw in =7.15 =式中： N 和 P 为介质常数， N ，P w w 3.31×108 Pa。Euler采用间断迦辽金法求解 方程，半离散格式为
pm dE å Ω Φ lΨ dΩ dtk = k k =0 pm ¶Φl pm ¶Φ å å Ω F E kΨk +G E kΨk l dΩ¶z ¶r k =0 k =0 pm pm
å å 5 ¶Ω F E kΨk nz e +G E kΨk nr e Φl dA （ ） k =0 k =0
式中：p 为多项式的最高次数；Ω 为单元面积； m Φl (l = 0 1 p ) 和 Ψk (k = 0 1 p )为试函m m数； E 为第k次多项式中Ψk 的系数； nz 和 nr为e e k单元边界的法向量；A为单元边界。Legendre选取正交多项式作为试函数，则流场中的变量E 为pm
å （6） E ( z r t )= E (t )Ψk ( z r) k k =0式（5）采用三阶Runge-Kutta显式格式进行时HLLC格式［15］计算数值通量间步进，采用F和G ，通过斜率限制器［16］抑制高精度离散格式带来的间断处非物理震荡。为了提高计算精度与效率，采
AMR 2），在流场中压力、用 技术构建多层网格（图密度等物理量出现间断的位置进行局部加密。
2 数值方法验证
3图 所示为水中冲击波与柱状气泡相互作用的数值计算结果对比。介质界面两侧拥有较大的密度与压力差，冲击波强度较大，对间断捕捉算法的计算稳定性及精度要求较高。计算区域取为24 mm×24 mm（此处仍采用带量纲的物理量，以便150×与文献结果［17-19 ］进行对比），底层网格数为150，采用3 AMR 3mm，层 网格。柱状气泡半径为气泡中心位于（12 mm，12 mm）处，冲击波马赫数1.72。气泡密度为1 kg/m3，压力为1 bar，气体常为 =1.4；气泡外部水的密度为1 000 kg/m3，压力数 γ g 1 bar，介质常数 N=4.4，P =6×108 Pa。冲击波为 w 54 mm。计算域的左边初始位置与气泡中心相距界为入口边界，右边界和上、下边界均为无反射边界。本文将计算结果与文献［17-19］的结果进行AMR了对比，虽然本文采用的底层网格数和 优化献［17］的结果，但仍能较好地模拟层数均小于文水中气泡与冲击波之间的相互作用。
3 数值计算分析
为了简化说明，设SWi为入射波； SWb为水下爆炸产生的冲击波；SWir为 SWi在自由液面的反射波；SWir*为 SWi在气泡上游表面的反射波；SWirr为在气泡下游表面的反射波；SWj SWir 为生成环状气泡时产生的冲击波；SWjr为 SWj在自由液面的反射波；SWt为气泡内部的透射冲击波。首先，考虑自由液面附近冲击波与初始内外压平衡的气泡之间的相互作用。气泡中心与自由=1.5，冲击波波前与气泡中液面的初始距离为 d f =2。计算域大小为3×6，采用 3层心的距离为 d s AMR 120×240，下边界为网格，最底层网格数量为入口边界条件，上边界为出口边界条件，左、右边界为壁面边界条件。各介质内的物理量记为［ ρ， ］，则气泡内部参数为［1，1，0，0，1.4］， p ，u ，v ，γ外部流场中水域参数为［1 000，1，0，0，7.15］，空气M=1.1，冲击域参数与气泡内部相同，冲击波强度4波与近自由面气泡的相互作用如图 所示。由于流场中未采用空化模型，因此压力存在负值。4（a）和 4（b）可由图 图 知，SWi传播至气泡表3面后分为 部分：一部分被反射形成稀疏波；另一部分透射至气泡内部并引起气泡内部气体的高速运动；剩下的一部分则绕过气泡继续向前传播。
4（c 4（d）可由图 ）和图 知，当 SWi到达自由液面后，一方面将从水中入射空气，由于水的特征阻抗高于空气，因此SWi以稀疏波的形式被反射回来，形成SWir；另一方面，反射回来的SWir作用于气泡表面，也从水中入射空气，因此稀疏波SWir将以冲击波的形式被反射回来，形成SWirr。以此类推，在气泡与自由液面之间，冲击波与稀疏波将不断地往复传播。此时，两侧自由液面在冲击波的作用下将出现上拱，而气泡上游表面将内凹形成射流，且在射流附近出现局部高压区域。4（e）和图4（f）可知，此时气泡在冲击波由图诱导下形成的射流已突破下游表面形成环状气泡，形成新的冲击波SWj，并向四周传播。由于SWj的强度远高于SWi，因此 SWj在自由液面处反射形成的稀疏波 SWjr所造成的负压也远高于SWir，此时，自由液面中心在SWj的作用下出现上拱。而在环状气泡的上游则出现了若干局部高压区域，这是由气泡形成射流时在上游产生的局部高压与向上游传播的SWj相交所产生。此外，当SWi和 SWj作用于自由液面时，有低强度冲击波透射入气相，将引起液面另一侧的气体运动。5 r=0图 所示为不同时刻 处的压力分布。5（a图 ）中单向箭头与各曲线相交的时间点从左t=0.402，0.702，0.772，0.816，0.852，至右分别为0.916，1.030，1.145；图 5（b）中 t1=0.702，t2=0.744， t3=0.775，t4=0.809，t5=0.863，t6=0.954，t7=0.982，t8= 1.016。SWi到达自由面后将以稀疏波SWir的形式5（b反射回来，并形成局部负压。由图 ）可知，当稀疏波 SWir作用在气泡下游表面后，将反射回冲击波SWirr，使附近流场逐渐改变为正压状态。而SWirr将再次到达自由面，并反射回稀疏波作用于气泡下游表面，最终再以冲击波的形式被反射。此时，气泡内部的冲击波SWt已接近气泡下游表5（b）中，t8面。图 时刻曲线即为气泡射流冲击下游表面后在冲击处附近产生的局部高压，并形成冲击波 SWj。当 SWj作用于自由液面时也会形成相应的稀疏波SWjr，并形成局部负压。由于SWj的强度远高于 SWi，则 SWjr所造成的负压也远高于
6图 所示为有、无自由液面时的气泡动态特性。自由液面反射的稀疏波作用于气泡表面时会使气泡下游界面发生朝向下游的运动，导致相应时刻的气泡体积增大，则形成最小体积的时间更长。由于反射波的作用时间较短，因此其产生的影响也较小。由于射流是由背向自由液面的气泡上游表面所形成，故自由液面的反射波对其影响较小，因此有、无自由液面时，气泡的射流顶点速度基本一致。水下多点爆炸比较常见，在产生气泡的同时 向四周传播冲击波，该冲击波与其他方向传播而来的冲击波相互耦合，进而对气泡运动及流场特性造成影响。冲击波（M=1.1）与初始内部高压的7气泡在自由液面附近的相互作用如图 所示，计8×12，采用 4 AMR算域大小为 层 网格，最底层网120×180。气泡内部介质参数为［ 1 630，格数为1 500，0，0，1.25］，其他参数与前文相同。 图7 冲击波（M=1.1）与初始内部高压的气泡在自由液面附近的相互作用（分图中左侧为压力分布，右侧为马赫场分布） Fig.7 Shock-bubble interaction with initial high inner pressure near free surface for M=1.1（In figures：left are pressure distribution，and right are Mach number field distribution）
7（a由图 ）可知，气泡释放的冲击波SWb与入SWi射波 在相汇处叠加，形成高压区。SWi与气泡SWir*，与先上游表面接触后发生反射形成稀疏波前的高压区叠加形成高压带状区并继续向四周传播。在气泡与自由液面之间的区域，SWb在自由液面处以稀疏波的形式被反射，随后在气泡下游表面以冲击波的形式被反射，如此往复，形成了间隔的冲击波区域与稀疏波区域。在初始时刻，气泡除了向四周释放冲击波SWb外，还在气泡内部形成内聚稀疏波，最终在气泡中心形成高压区，如图7（b）所示。8图 所示为不同初始条件下气泡（气泡内部1 500）与自由液面形态随时间的变化。初始压力与没有入射冲击波时相比，在SWi的作用下，同时刻的气泡体积减小，气泡膨胀的同时向自由液面偏移，自由液面两侧发生明显的上拱运动。当入射冲击波SWi强度较大时，自由液面的上拱运动更明显，气泡下游表面较快地形成明显的尖峰，并与自由液面发生融合。
4结论
Level Set本文采用间断迦辽金法，结合 和Real Ghost Fluid方法求解近自由液面气泡与冲击波间的相互作用，详细分析了气泡的运动特性、各类冲击波与稀疏波的生成和发展过程以及流场中压力、速度等物理量的分布特性。主要结论如下： 1）在气泡、冲击波与自由液面三者相互作用的过程中，流场内存在多个复杂的稀疏波系和冲击波系。2）当气泡初始内外压平衡时，气泡在入射冲击波的作用下溃灭形成环状气泡，在此过程中产生的冲击波强度远高于入射冲击波。在该冲击波作用下，自由液面中心形成向上凸起的水冢。当自由液面存在时，气泡的运动过程与自由场中相似，但是在自由液面反射的稀疏波作用下，气泡的溃灭速度减缓。3）当气泡内部存在一定的初始压力时，气泡运动开始时释放的冲击波与入射冲击波发生耦合，形成更为复杂的流场分布。在入射冲击波作用下，气泡并未形成环状气泡，且与无入射波的情况相比，气泡向自由液面偏移，体积减小，自由液面的上拱运动增大。
参考文献：
［ 1 ］ JOHNSEN E，COLONIUS T. Shock-induced collapse of a gas bubble in shockwave lithotripsy［J］. The Jour⁃ nal of the Acoustical Society of America，2008，124（4）： 2011-2020. ［2］ JOHNSEN E. Numerical simulation of non-spherical bubble collapse：with applications to shockwave litho⁃ tripsy［D］. California：California Institute of Technolo⁃ gy，2007. 3］ 奔. R22 究［J］. ［ 王革，关 激波作用下 气泡射流现象研力学学报，2013，45（5）：707-715. WANG G，GUAN B. A study on jet phenomenon of R22 gas cylinder under the impact of shock wave［J］. Chi⁃ nese Journal of Theoretical and Applied Mechanics， 2013，45（5）：707-715（in Chinese）. ［4］ JAMALUDDIN A R，BALL G J，TURANGAN C K，et al. The collapse of single bubbles and approximation of the far-field acoustic emissions for cavitation induced by shock wave lithotripsy［J］. Journal of Fluid Mechan⁃ ics，2011，677：305-341. ［5］ KLASEBOER E，FONG S W，TURANGAN C K，et al. Interaction of lithotripter shockwaves with single in⁃ ertial cavitation bubbles［J］. Journal of Fluid Mechan⁃ ics，2007，593：33-56.
［6］ PARK J. A coupled runge kutta discontinuous galer⁃ kin-direct ghost fluid （RKDG-DGF） method to near-field early-time underwater explosion（UNDEX） simulations［D］. Blacksburg：Virginia Polytechnic In⁃ stitute and State University，2008. ［7］ COCKBURN B， HOU S C ，SHU C W. The runge-kut⁃ ta local projection discontinuous Galerkin finite ele⁃ ment method for conservation laws IV：the multidimen⁃ sional case［J］. Mathematics of Computation，1990， 54（190）：545-581. ［8］ COCKBURN B，SHU C W. The runge-kutta discontin⁃ uous Galerkin method for conservation laws V：multidi⁃ mensional systems［J］. Journal of Computational Phys⁃ ics，1998，141（2）：199-224. ［9］ QIU J X LIU T G ，KHOO B C. Simulations of com⁃ ， pressible two-medium flow by runge-kutta discontinu⁃ ous Galerkin methods with ghost fluid method［J］. Com⁃ munications in Computational Physics，2008，3（2）： 479-504. 10］ 叶曦. ［ 水下气泡运动的噪声及近场载荷特性研究［D］.哈尔滨：哈尔滨工程大学，2014. YE X. Study on noise and near-field load characteris⁃ tics induced by bubble motion［D］. Harbin：Harbin Engineering University，2014（in Chinese）. ［11］ OSHER S，FEDKIW R. Level Set methods and dy⁃ namic implicit surfaces［M］. New York：Springer， 2002. ［12］ WANG C W， LIU T G ，KHOO B C. A real ghost flu⁃ 89 （上接第 页） ［ 5］ 邓聚龙. 灰理论基础［M］.武汉：华中科技大学出版社，2002.
6］ 中国船级社. 小水线面双体船指南［S］. ［ 北京：人民交通出版社，2005.
7］ 真. ［ 尹群，尹 小水线面双体船船体结构和强度初步研究［J］. 造船技术，1997（5）：22-24. 8］ 刘斌. 小水线面双体船极限强度研究［D］. ［ 武汉：武 id method for the simulation of multimedium com⁃ pressible flow［J］. SIAM Journal on Scientific Comput⁃ ing，2006，28（1）：278-302. ［13］ BERGER M J，COLELLA P. Local adaptive mesh re⁃ finement for shock hydrodynamics［J］. Journal of Computational Physics，1989，82（1）：64-84. ［14］ LANEY C B. Computational gasdynamics［M］. Ox⁃ ford：Cambridge University Press，1998. ［15］ HARTEN A， LAXPD ，LEER B V. On upstream dif⁃ ferencing and Godunov-type schemes for hyperbolic conservation laws［J］. SIAM Review，1983，25（1）： 35-61. ［16］ COCKBURN B， LIN S Y ，SHU C W. TVB run⁃ ge-kutta local projection discontinuous Galerkin fi⁃ nite element method for conservation laws III：one-di⁃ mensional systems［J］. Journal of Computational Phys⁃ ics，1989，84（1）：90-113. ［17］ NOURGALIEV R R，DINH T N，THEOFANOUS T G. Adaptive characteristics-based matching for com⁃ pressible multifluid dynamics［J］. Journal of Computa⁃ tional Physics，2006，213（2）：500-529. ［18］ BALL G J，HOWELL B P，LEIGHTON T G，et. al. Shock-induced collapse of a cylindrical air cavity in water：a free-lagrange simulation［J］. Shock Wave， 2000，10（4）：265-276. ［19］ BOURNE N K，FIELD J E. Shock-induced collapse of single cavities in liquids［J］. Journal of Fluid Me⁃ chanics，1992，244（1）：225-240. 汉理工大学，2009. ［9 黄雪忠. ］ 小水线面双体船的结构分析与优化设计［D］.上海：上海交通大学，2009. ［10］ 尹群，任鸿，朱安庆.小水线面双体船结构加强方式研究［J］. 造船技术，2003（3）：14-17，13. YIN Q，REN H，ZHU A Q. Research of SWATH ship structural strength methods［J］. Marine Technol⁃ ogy，2003（3）：14-17，13（in Chinese）.