Effect of a rigid structure on the dynamics of a bubble beneath the free surface

Shi-Min Li, A-Man Zhang , Nian-Nian Liu

College of Shipbuilding Engineering, Harbin Engineering University, Harbin 150 0 01, China

Keywords: Underwater explosion Bubble dynamics Boundary integral method

ABSTRACT The coupling effect between a pulsating bubble and a free surface near a rigid structure is a compli- cated physical process.In this study, the evolution of an underwater explosion bubble and the free sur- face near a rigid structure is modeled by the boundary integral method.An approach of “double-node method”is used to maintain the stability of fluid-structure junction in the simulations, and meshes on the free surface and the structure are transformed to an open domain to ensure the calculation accu- racy and efficiency.Validations are conducted against an underwater explosion experiment near a rigid structure.As a result, the simulations trace the jetting behavior of the bubble and the rise of the free sur- face.Finally, the bubble migration and the height of the free surface for different structure draughts are analyzed.

Bubble dynamics near the free surface have always been a pop- ular topic due to the wide applications in underwater explosions [1–4] , jet printing [ 5 , 6 ], and ultrasonic cleaning [ 7 , 8 ].Once a struc- ture is near the free surface, the three-phase coupling effect be- tween the bubble, the structure, and the free surface becomes very complicated.The most important flow feature is the migra- tion of the bubble and the formation of the water column on the free surface.The migration direction of the bubble determines the impact region of the high-speed liquid jet generated during the asymmetric collapse and the characteristics of the shear flow be- tween the bubble and the structure, which in turn changes the ef- fects of underwater explosions and surface cleaning.The formation of the water column is a valuable tool in jet printing and some military affairs.An in-depth understanding of the interaction be- tween the bubble and free surface near typical structures is essen- tial for rational use of the characteristics of bubbles and the free surface.

Contributions have been made to many numerical and experi- mental studies on the interaction between bubbles and free sur- faces [9–17] or other various boundaries [18–23] .Many interest- ing phenomena, including the secondary cavitation, toroidal bub- ble splitting, crown-shaped water doom, and feather broken sur- face, are discovered.Especially for numerical simulation, it is one of the most effective and applicable methods for systematic studies and is not limited by experimental conditions.However, most nu- merical studies are all based on axisymmetric configurations.Once there is a structure near the free surface, the entire fluid domain becomes an asymmetric three-dimensional flow field, which sig- nificantly increases the difficulty and time consuming of the nu- merical simulation.Therefore, the boundary integral method (BIM) is a good choice in studying three-dimensional bubble dynamics due to its high efficiency and accuracy [24–28] .Some literature is already studying the dynamic behavior of bubbles near structures using three-dimensional BIM in recent years [29–32] , but they fo- cus more on the infinite boundaries or specific engineering appli- cations.Study on systematic descriptions of the characteristics of the bubble and free surface is relatively rare.

In this study, a numerical model for calculating the interaction between the bubble and free surface near a rigid boundary with a certain draught is established firstly using BIM.A simple and ba- sic structure, a right-angled corner, is regarded as the rigid struc- ture.An underwater explosion experiment near a cabin model is conducted to verify the reliability of our numerical method.Qual- itative analyses are performed for the behavior of the bubble and free surface for two draughts of the structure.Finally, the effects of the structure draught on bubble migration, jet direction, and free surface height are discussed.

The physical scenario in this paper can be described as the interaction between the bubble and free surface when two rigid perpendicular semi-infinite rigid walls (a vertical and a horizon- tal wall) exist, as shown in Fig.1 .The generation position of the bubble is regarded as the origin to establish a Cartesian coordinate system.Thez-axis points opposite to gravitational acceleration, and they-axis is horizontal to the right.The directions of thex,y, andzaxes meet the right-hand rule.The vertical wall crosses the free surface, and the structure is symmetrical about the bubble along thex-axis with the draught ofD.The structural walls are entirely rigid, without considering the flow of fluids on the external do- main.

The boundary integral method is based on the potential flow theory; that is, the fluid is regarded to be inviscid, irrotational, and incompressible.These assumptions can be reasonably adopted dur- ing most of the bubble’s lifetime [33–35] .Therefore, the Laplace equation is used to control the dynamics of fluids:

whereφis the fluid velocity potential.

To numerically solve Eq.(2) , the integral form of Laplace equa- tion (boundary integral equation [ 25 , 26 , 36 ]) is employed:

Based on the Bernoulli equation [ 37 , 38 ], the velocity potentials of the bubble and free surface are updated:

wherep∞ is the hydrostatic pressure at the initial position of the bubble, andpis the current pressure on the boundary surface.

To make the research universal, we organize all simulations on a non-dimensional system.Three fundamental dimensionless quantities are the maximum radius of the bubbleRm, the density of the liquidρ, and the hydrostatic pressure at the initial position of the bubblep∞.On this basis, the reference quantities of time, velocity and velocity potential can be easily obtained.Subsequent research contents are all based on the non-dimensional system un- less otherwise specified.Here, we give three dimensionless param- eters used in this study:

whereδis the buoyancy parameter used to characterize the bubble scale,dfanddware the distance between the initial bubble and the free surface and vertical wall, respectively.

In addition, the Kelvin impulse [ 33 , 39 ] is used here to quantita- tively measure the movement trend of bubbles:

wherenis the normal vector of nodes on the bubble surface, andsbis the area of the bubble wall.

The boundary conditions of the bubble surface and the free sur- face can be written as:

where the superscript indicates that the corresponding physical quantity is dimensionless; as ignoring the viscous effect and sur- face tension,pof nodes on the bubble surface is equal to the inter- nal bubble pressure calculated according to the adiabatic assump- tion (p′ =ε(V0 ′/V′)ς,εis the strength parameter,V0 ′ andV′ are the initial and current bubble volume, andςis polytropic coeffi- cient).

For nodes that belong to the free surface and rigid wall simul- taneously, they need to meet the boundary conditions of the free surface and cannot penetrate the rigid boundary, causing their ve- locities to be unstable and eventually collapse the simulation.In previous works [40] , the “weighted interpolation method”is used to deal with this problem, but it is only suitable for the slightly de- formed free surface.When the free surface dramatically deforms, a more general method is employed: the “double-node method”[ 27 , 41 , 42 ].Firstly, we define twon×nmatrices,GijandHij(i,j= 1, 2, …,n), representing the influence of the elements around nodejon nodei.GijandHijare calculated by ∫ ∫

sG(i,j)dsandrespectively.The detailed discrete pro- cess of the boundary integral method and the principle of the ’double-node method’ can refer to the works of the literature [ 24 , 43 ].Here, only the coefficient matrix used in the calculation is given:

whereκis the solid angle;sis the boundary surface composing the flow field; P and Q are discrete points on the boundaries (the source point P and the integration point Q);∂/∂nis the normal deviation of nodes on the boundaries.Once the boundaries is dis- cretized by a certain number of nodes and meshes, a calculation relationship between P and Q can be obtained, and thus a solv- ablen-element linear equation system (nis the number of nodes) is formed.After solving this linear equation system, the normal ve- locities of the bubble and the free surface and the velocity poten- tial of the wall are determined because the velocity potential of the bubble and free surface is known, and the normal velocities of rigid walls are always zero.

where the subscript ’f’ represents the nodes on the free surface, ’w’ represents that on the rigid walls, and ’c’ denotes that on the junction which is divided into two parts: the influence of the free surfacecf and the rigid boundarycw .

It can be seen that we need to calculate the effects of free sur- face elements and wall elements separately when calculating the matrixGij.Since the normal velocity of the node on the rigid walls is zero (∂φcw/∂n= 0), Eq.(8) can be simplified as:

Since∂φw/∂nis zero, the values ofGfw ,Gcfw andGww can also be set to zero; thus, only the influence of the free surface elements on the nodes needs to be taken into account for the matrixGij.

Another challenging problem in simulations is the scale of the boundaries.The boundary integral Eq.(2) requires the boundary surface to be closed, resulting in numerous nodes built outside the computational domain.Therefore, we need to artificially regard the nodes outside the computing domain to be on rigid walls to keep them stationary, which will undoubtedly cause errors and signifi- cantly reduce the calculation efficiency.Therefore, we use a vortex made by the free sides of boundaries to compensate for this de- ficiency.The basic idea is to replace the solid angle of the outer domain surface to nodes with the induced velocity potential of the vortex ring.The specific process is that all nodes on the free side of boundaries are considered on a three-dimensional vortex ring, and a local cylindrical coordinate system is established with the normal direction of the vortex ring as thez-axis.The detailed the- oretical foundation and calculation details of this approach can be found in Liu’s works[ 27,31 ] .The result is directly provided here:

whereris the vector from the node on the boundary to the center of the microelement of the vortex ring,rzandrrare the normal and radial components ofr, respectively;Cis the three- dimensional vortex ring formed by the free side of the free surface and walls, and dlis the microelement on vortexC.

An underwater explosion experiment was carried out to verify the feasibility of our numerical method.The experiment was con- ducted in a 4 m ×4 m ×4 m water tank.The layout of the exper- iment site can refer to the previous works [44] .5 g RDX explosives generates a bubble withRm of about 29.8 cm, with the buoyancy parameterδof 0.18.A cabin structure is welded above the water tank on a fixed steel frame and crosses through the free surface.Both the length (along thexaxis) and width (along theyaxis) of the structure are 80 cm.The distance from the bottom of the structure to the free surface is 40 cm.The depth of the explosive is 65 cm, and the relative position between the bubble at inception and the structure is shown in the first panel.The high-speed cam- era shoots the underwater explosion process at a speed of 50 0 0 frames per second, and the typical image is shown in Fig.1 .The high-pressure bubble expands rapidly after the explosive is deto- nated, remaining nearly spherical in the expansion stage.After the bubble volume reaches its maximum, the pressure of the flow field caused by the free surface and the wall is asymmetry, resulting in the obvious non-spherical collapsing.The upper part of the bubble is flattened due to a local high-pressure zone (3-4 frame) caused by the atmospheric boundary condition of the free surface [ 13 , 45 ].Under the combined action of the free surface and the rigid walls, the upper left part of the bubble collapses obliquely downwards towards the end of the first cycle.

At the same time, the numerical model simulates the behav- ior of the underwater explosion bubble.The identification of initial conditions can refer the works of Li [43] , and we do not describe it here in detail.In this and subsequent cases, the adiabatic parame- terςis 1.25 (for explosion bubbles [46]), and the initial radius and inner pressureεof the bubble are 0.149 and 100, respectively.The numerical simulation reproduces the bubble oscillation and jetting behavior well.Att= 51.284 ms, the bubble also exhibits the trend of collapsing obliquely downwards.The inclination angle is calcu- lated by Kelvin impulse to be about 56 °(arctan(Iz /Iy),Iz andIy are thez-component andy-component ofI, respectively).The direct measurement result in the image (that is, 58 °, the direction per- pendicular to the visible outline of the upper right part of the bub- ble, which matches the movement direction of the black flocculent impurities in frame 10) is very consistent with the calculation re- sult, verifying the feasibility of the numerical model on predicting the bubble migration.

In the numerical simulation, the left part of the bubble is more rounded and swollen than that in the experimental image (frame 9 and 10), which is attributed to the slight movement of the struc- ture.In fact, it is a challenging task to completely rigidify the structure as attacked by underwater explosions.In the experiment, the fluids between the bubble and the rigid walls flow freer since the rigid structure is slightly bounced during the bubble expan- sion stage, while the liquids are completely confined to one side of the walls in the simulation.So the bubble surface nearer the rigid structure shrinks slower in the experiment than in the simulation.Nevertheless, in terms of the overall oscillation process and migra- tion direction of the bubble, the weak motion of the structure does not affect the collapse process of the bubble to a large extent.

Figure 2 shows the evolution of the bubble and free surface as the right-angled structural draughtDis set to 0.5Rm when the bubble is located near the corner of the structure.In subsequent cases, the buoyancy parameterδis set to 0.4, equivalent to a bub- ble withRm of about 1.8 m.The bubble expansion causes the free surface to rise significantly, in which the part directly above the bubble rises fastest (we call it ’water doom’).The entire free sur- face is in the shape of an ’inclined umbrella’.At the maximum bubble volume (frame 2), the bubble surface close to the free sur- face is sucked in, and the part near the structure is repelled.Sub- sequently, an oblique downward liquid jet emerges as a result of the local high-pressure area upper the bubble [47] .The effect of the free surface drives the bubble to move downward, the effect of rigid walls attracts the bubble, and the buoyancy effect causes the bubble to migrate upwards.These effects are collectively re- ferred to as the ’Bjerkness effect’ [ 4 8 , 4 9 ].The oblique downward jet (3-5 frames) results from the combined influence of the free surface, rigid structure and buoyancy effect.Toward the end of the first cycle, the liquid jet pierced the opposite bubble wall, forming an annular bubble (frame 6).

We then illustrate the evolution of the bubble and free surface as the structural draft is 2Rm in Fig.3 when the bubble is located on the side of the vertical wall.In the expansion stage, the left part of the bubble is flattened due to the obstruction of the wall, and the upper part is sucked into the free surface.As the bubble col- lapses, the liquid jet tends to point downward and hit the opposite bubble wall under the influence of the free surface.At the same time, the lower right part of the bubble also tends to be concave upwards diagonally under the action of the buoyancy.Before the jet impacts the bubble wall, the rightmost part of the bubble be- comes pointed, resulting from the more significant influence of the vertical wall.Compared with the previous case, the more oversized vertical wall hinders the oscillation of the left part of the bubble so that the right part shrinks more quickly relative to the other parts of the bubble wall.Therefore, the upper right part and the lower right part of the bubble surface are respectively affected by the free surface and the buoyancy effect at a relatively faster con- traction speed, resulting in a sharpened shape of the bubble wall.

Figure 5 shows the displacement of the bubble centroidu(Panel a) and the evolution of the height of the water doomh(Panel b) for different structure draughts.The bubble migration process in the vertical direction differs slightly for various draughts because the distance from the bubble to the free surface and the buoyancy parameters are fixed,γf=γw= 1,δ= 0.4 .During most of the first cycle, the vertical migration range of bubbles increases as the structural draft increases.However, the vertical migration range does not change monotonously with the structural draft by the end of the first cycle, mainly because of the non-spherical bub- ble shape as the descriptions in Figs.3 and 4 .In general, when the scale and position of the bubble are fixed, the influence of the structural draft on the longitudinal bubble migration in the first period is not apparent, and the bubble does not significantly mi- grate downward or upward.However, the jet directions all point downwards (not displaying all cases), which is consistent with the prediction results of Kelvin impulse theory [ 49 , 50 ] (asγfδ<0.442 , the free surface effect is stronger than the buoyancy effect).In- stead, the influence of the structural draft on the horizontal migra- tion of bubbles is significant.The bubble is slightly repelled by the walls in the early stage and is intensely attracted at the late stage.As the structure draught decreases, the bubble is less affected by the rigid walls, so the migration range of the bubble away from and towards the walls is smaller.Panel b shows the time histo- ries of the height of the water doom.The free surface rises rapidly as the bubble expands and changes relatively slowly after reaching the maximum volume.As the structural draft increases, the free surface is raised higher.A physical explanation is provided as fol- lows: as the structural draft increases, the flow of liquids between the bubble and the structure is more limited, so the pulsating of the bubble makes the fluids near the walls flow to the top of the bubble instead of the bottom of the structure, leading to the in- crease in the pressure between the bubble and the free surface.

In this study, the boundary integral method (BIM) is used to model the dynamic behavior of the bubble and the evolution of the free surface for different structure draughts.The ’double-node method’ and ’open-domain correction’ are adopted to ensure the stability and accuracy of the numerical simulations.By comparing with the underwater explosion experiment, the feasibility of our numerical model has been verified.Under the influence of the free surface, the rigid structure and the buoyancy effect, a downward liquid jet is generated and deflected towards the structure surface.Simultaneously, the free surface is obviously raised to form a wa- ter doom taking an ’inclined umbrella’ shape.When the initial po- sition of the bubble is fixed, the vertical movement of the bubble changes little with the structural draft in the first period, and the jet direction is consistent with the prediction of Kelvin impulse theory.The structural draught significantly affects the horizontal migration of the bubble.The larger structural draught causes the bubble to be repelled more strongly in the early oscillation stage and move toward the structure more vigorously in the late os- cillation stage.The existence of the structure also causes the free surface to rise higher, which can provide design references for the ’water curtain anti-missile’ near a ship.

Declaration of Competing Interest

There exists no conflict of interest statement in this works.

Acknowledgments

This work is supported by the National Key R&D Program of China (2018YFC0308900), the National Natural Science Founda- tion of China (11872158 , 52001085), the Postdoctoral Science Foun- dation (2019M661256), and the Heilongjiang Postdoctoral Fund, China (LBH-Z19135).