Computer Modeling in Engineering & Sciences

DOI: 10.32604/cmes.2021.015259

ARTICLE

The Nonlinear Coupling of Oscillating Bubble and Floating Body with Circular Hole

Ming He, Yunlong Liu*, Shaofei Ren and Wentao Liu

College of Shipbuilding Engineering, Harbin Engineering University, Harbin, Heilongjiang, China
*Corresponding Author: Yunlong Liu. Email: yunlong_liu@hrbeu.edu.cn
Received: 05 December 2020; Accepted: 08 February 2021

1  Introduction

Since the last century, bubble dynamics [1–4] has been receiving increasing attention from researchers, and it has been widely used in resource detection with air-gun array [5], underwater explosion [2,6–8], and biomedical fields [9]. With the rapid development of weapons and equipment technology, the damage caused by the underwater explosion of ships and marine structures has become more severe, which seriously restricts the vitality of vessels and marine systems. When the underwater explosion occurs near the marine structures, the generated shockwave will cause the marine structures to rupture. After the hole is created, the air inside the structures will bring about the appearance of the free surface. Simultaneously, generated bubbles will further affect the vitality of structures. The structural deformation induced by pulsating pressure, the further impact of the bubble jet load on the structures, and the water hump which causes water to flow into the cabin are all we should consider. Essentially, they are nonlinear coupling problems among the incomplete moving boundary, free surface and bubble. Therefore, to master the mechanics of the entire underwater explosion process, it is imperative to research the fluid-structure interaction between bubbles and broken structures. In recent years, many researchers have conducted extensive researches on related issues [10], but there are relatively few references considering the fluid-structure interaction effect.

Theoretical analysis [11], model testing [1,2,12–14], and numerical analysis [13,15–20] are three main methods for studying the underwater explosion problems. Extensive researches have been accomplished on bubble dynamics. In the early stage of research, theoretical research was the first to appear. The prediction model of the underwater explosion bubble was established by Geers et al. [11], who aim to solve spherical bubble problems. For nonlinear coupling problems under complex boundaries, theoretical analysis is arduous. With the development of the boundary element method [19–21], many researchers began to study non-spherical bubbles. However, the effect of the boundary element method to solve complicated engineering problems is not ideal. Because it requires many specific numerical processing skills. Therefore, researchers began to find more general numerical methods to study shock waves and bubbles, including the finite volume method [22], discontinuous Galerkin method [23], and smoothed particle hydrodynamics method [24]. Li et al. [25] used the finite volume method to research the bursting behaviors of bubble near free surface. Wang et al. [18] utilized the discontinuous Galerkin method to research the characteristic of pressure produced by the jet penetration, considering the compressibility of the fluid. Various numerical methods have deepened people’s understanding towards the underwater explosion problems. At present, relatively few articles have refined and full-cycle simulation, considering the movement of the structures under the ultra-close boundary conditions. Simultaneously, experimental research is also fundamental and can verify numerical algorithms. Cui et al. [1,14] did some basic rational experiments in this field.

Nowadays, numerical methods are continually evolving. The Eulerian finite element method is good at solving transient problems related to the large deformation of structures. This method was first established by Benson. Then he gave the main framework of the Eulerian finite element method, and proposed a brand-new advection approach for the node-centered parameters like momentum [26–29]. Liu et al. [30–32] applied it to the research of bubble dynamics, and obtained good numerical results. Besides, the volume of fluid method (VOF) is a mature interface capture method [22,33–36]. It can simulate the fracture of the free surface in solving the hydrodynamic problems. This method has been embedded in many commercial software, and has achieved good results.

The immersed boundary method was initially proposed by researchers to simulate the flow of blood in the heart [37,38], and this method is widely used over the years [39–41]. Fai et al. [42] proposed the lubrication format of the immersed boundary method. Yuan et al. [41] applied the immersed boundary method to deal with the compressible viscous fluid flow. Afra et al. [43] combined the immersed boundary method with the lattice Boltzmann method to establish a reliable numerical solution model. Kim et al. applied the penalty immersed boundary method for the research of a rigid body in fluid, which proved the penalty immersed boundary was quite significant. By comparison, the penalty immersed boundary method is a relatively accurate and stable processing method for fluid-structure interaction problems [44,45].

To further study the nonlinear coupling mechanism of the underwater explosion, the fluid-structure interaction model between the oscillating bubble and floating body with circular hole is established based on the Eulerian finite element method, the VOF method, the monotonic upwind scheme for conservation laws (MUSCL), and the penalty immersed boundary method. This paper focus on the basic research in engineering problems, so the physical model has been simplified to the axisymmetric model. A floating rigid body with a circular hole is initially placed at the free surface to represent the broken part of ships and marine structures, affected by the underwater explosion. Then, an oscillating bubble is beneath the hole and interacts with the floating structure above, as showed in Fig. 1. d is the distance from the free surface to the explosive location, Rh is the hole radius of the floating body, Rw is the radical width of the floating body, At is the axial thickness of the floating body, Rm represents the maximum radius of the bubble which can be used to calculate the initial conditions, and the origin of the axial coordinate is directly below the axis of symmetry.

Figure 1: The schematic diagram of the nonlinear coupling problem in this paper

This article is arranged as following. Above all, the numerical algorithm is presented, and verified by comparing the results of simulation with those of the electric discharge bubble experiment and explosion experiment. Secondly, this paper considers the influence of the hole size, the explosive location, the floating body density, and the buoyancy on the fluid-structure interaction characteristics. It shows from the results that the four factors play essential roles simultaneously. In addition, the free surface crushing, water hump, fluctuation, spike skirt, and jet penetration is carefully investigated. Finally, some conclusions are summarised.

2  Theoretical and Numerical Methods

2.1 The Eulerian Finite Element Method for Fluid Solver

The Eulerian finite element method can effectively deal with the explosion and impact problems [26,31,32], where the viscosity and thermal conduction can be ignored [12]. In this method, the governing equation for fluid is the following equation, i.e.,

∂φ∂t+∇⋅(φu)=S.(1)

where u={ur,uz} is the velocity vector, φ is the conserving variable, t is the time, and S is the source term. The above formula is based on a cylindrical coordinate system, where r and z are the radial and axial coordinates. Utilizing the operator split, the solution of the above equation can be divided into the Lagrangian phase and Eulerian phase, i.e.,

∂φ∂t=S,(2)

∂φ∂t=-∇⋅(φu).(3)

As for mass, momentum, and energy conservation, S and φ can be expressed as the following equation:

{S=0,φ=ρS=-∇p+f,φ=ρu,S=-p∇⋅u,φ=ρein(4)

where ρ is the material density, p is the pressure, ein is the internal energy per unit mass of fluid, and f is the body force. We should pay attention to that ∇⋅u=∂ur∂r+∂uz∂z+1rur.

The Lagrangian phase is to solve Eq. (2), and it is calculated by the explicit finite element method. Then, the Galerkin equation used to solve the momentum conservation can be written as

∬Ωρu̇ϕαϕβds=∬Ω(ρfϕα+p∇ϕα)ds-∫Γpϕαndl,(5)

where Ω is the discretized two-dimensional orthogonal computation domain, Γ is its boundary, n is the normal element vector pointing to the outside of the domain, ϕ is the computational shape function, and α, β are the number of nodes. Immediately, the fluid velocity is renewed by the explicit integral scheme, and the mesh advances with the fluid material, i.e.,

u(n+12)=u(n-12)+u̇(n)Δt,(6)

x(n+1)=x(n)+u(n+12)Δt,(7)

where the superscript n indicates the increment number, and x is the displacement, and Δt is the time increment. Immediately, the material density and specific internal energy are updated by

dρdt=-ρ∇⋅u(8)

and

d(ρein)dt=-(p+ρein)∇⋅u.(9)

Meanwhile, the Lagrangian phase is over, and the Eulerian phase starts. The core idea of the Eulerian phase is to achieve the transport of parameters by the VOF method and the MUSCL. The VOF method is a commonly used and effective method in computational fluid dynamics [33,34,46,47]. Its transport equation is [32]

dfdt=u⋅∇f+∂f∂t=0,(10)

where f is the volume fraction. The schematic diagram of the VOF method is explained in Fig. 2. Here, Δvi,y,j is the transport volume of element, i, j is the number of adjacent element situated at the side of the interface, and y is the number of the element boundary. In general, the interface between the two phases is difficult to determine. In this paper, the interface is solved by the least-squares method, i.e.,

∂∂kb Σj[fj-f(rj,zj)]2=0,(b=1,2,3…6),(11)

where f(r,z)=k1r2+k2z2+k3rz+k4r+k5z+k6 is the assumed distribution function of the corresponding element volume fraction as a function of position. Coefficients kb are used to obtain the characteristic parameters of the fluid interface, and they are calculated by solving Eq. (11). The normal vector of the fluid interface is required in numerical calculations, and it can be determined by k4 and k5, i.e., n=[k4,k5]k42+k52. After the transport volume is computed, the MUSCL [26] is used to transport the element-centered parameters. For a more detailed introduction, please refer to our previous article [32].

Figure 2: The schematic diagram of the VOF method

After the Eulerian phase, the fluid equation of state is used to complete the closed solution of equations. Referring to related literature, the equation of state of the fluid is given by the Tammann equation [31], i.e.,

p=ρein(χ-1)-χPw,(12)

where χ denotes the ratio of specific heats, and Pw is the reference pressure. The corresponding fluid parameters of water, air, and explosive production are noticed in Tab. 1. Here, W is the initial explosive mass, and V0 is the initial explosive volume.

Table 1: Fluid parameters of water, air, and explosive production

2.2 The Penalty Immersed Boundary Method

For the nonlinear fluid-structure interaction problems, the processing of the fluid-structure interaction interface is complicated. The penalty immersed boundary method is used in this paper, to perform the unbalanced constraint processing of the boundary position [44,45,48,49]. In numerical simulation, this method adopts a relatively simple Cartesian grid, which avoids the divergence caused by the body-fitted grid [29]. Submission of the penalty immersed boundary method changes the governing equations of the fluid, and introduces a new additional volumetric force. Now, this paper extends the governing equation from the two-dimensional format to axisymmetric format. And the body force is the resultant force of the additional volumetric force and gravity. Continuously, Eq. (13) is the momentum equation for the fluid where the solid boundary node is located, i.e.,

∂ρu∂t+∇⋅(ρuu)+∇p=fc+ρg.(13)

The interaction force fc can achieve the fluid-structure interaction condition, and it can be solved through continuous iteration of the immersed depth [29]. For the nonlinear problems, such as the explosion problems, the immersed depth formula is improved in this model. This paper adds a damping force term to increase the stability of the numerical simulation further, i.e.,

fRB=-κεn-Cε̇n,(14)

εn+1=εn+Δt(ub-∑m=14ϕmum)⋅n,(15)

ε̇=εn+1-εnΔt=(ub-∑m=14ϕmum)⋅n,(16)

fc=-ϕmfRB(m=1:4),(17)

where κ is the penalty stiffness coefficient, fRB is the total interface force, ε is the immersed depth, ub is the velocity of the contact boundary node, n is the normal vector, um is the node speed of the fluid element in which the solid boundary node is located, ϕm is the shape function of node m , and C is the damping coefficient. The schematic diagram of the penalty immersed boundary method is explained in Fig. 3. Bringing Eqs. (14)–(17) into Eq. (13), the fluid-structure interaction problems are solved in this paper. In addition, the time step obtained should be less than the numerical oscillation period of fluid-structure boundary, to ensure the stability of the calculation results. This oscillation is caused by the introduction of the penalty stiffness coefficient.

Figure 3: The schematic diagram of the penalty immersed boundary method

2.3 Non-Dimensionalization

To extend the law obtained in this paper, the dimensionless processing is performed on all physical variables, which helps us understand the physical problems better. The corresponding rules are showed in Tab. 2. Here, the subscript l represents the water, At is selected as the reference length, and P∞ is the hydrostatic pressure at the depth of the initial bubble center.

Table 2: Variable scales for non-dimensionalization

In addition, Rh* is the dimensionless radius of the hole, ρ* is the dimensionless density of the floating body, t* is the dimensionless time, d* is the dimensionless distance from the free surface to the explosive location, R* is the dimensionless equivalent radius, u* is the dimensionless axial velocity of floating body, r* is the dimensionless coordinate of radical direction, and z* is the dimensionless coordinate of axial direction. In previous studies, characteristic variables of bubbles are generally used for dimensionless processing, including maximum equivalent bubble radius. But in this article, more attention is paid to the fluid-structure coupling characteristics of the structure. So the axial thickness of the floating body, as the characteristic variable of structure, is used for dimensionless processing. Furthermore, the buoyancy parameter is defined as δ=ρgRm/p∞. And the bubble equivalent radius can be calculated by

R=(32∬Ωferdrdz)13.(18)

Here, fe is the volume fraction of explosive production.

3  Numerical Verification and Convergence Test

3.1 Near-Free Surface Cracking Numerical Verification

In this paper, the numerical method is verified by our electric discharge experiment. The experiment was performed in a 50×50×50cm3 water tank. The bubble is generated by a 220 V high voltage discharge device. Meanwhile, the free surface and bubble characteristics are captured by a high-speed camera. For more detailed experimental introduction, please refer to the previous paper [14,25]. The phenomenon of the water hump and bubble jet penetration can be accurately recorded. The thickness of the polymethyl methacrylate plate in the test is 2.0 mm, and the plate is fixed on the water. In addition, there is a circular hole in the middle of the plate. And Rh is 15.0 mm, Rm is 13.6 mm, and d is 7.0 mm. In order to display the boundary of the fluid vividly, the image of experimental results is subjected to the boundary capture processing. The experimental results are as following showed in Fig. 4 (above).

Figure 4: The experimental results diagram (above) at t = 0.60 ms, 1.40 ms, 2.80 ms, and 3.40 ms, and the numerical results diagram (below) at t = 0.52 ms, 1.32 ms, 2.75 ms, and 3.14 ms for case where Rm = 13.6 mm, d = 7.0 mm, and Rh = 15.0 mm

The numerical results are showed in Fig. 4 (below), and the boundary condition of simulation is same as that of experiment described above. In simulation, the initial conditions are obtained by the equation of bubble motion proposed by predecessors [12]. The size of the calculation domain is 0.05m×0.15m, and the element size is 0.00025 m. The boundary of the computational domain is a numerical non-reflective boundary condition [31,50]. It can be seen from the Fig. 4 that the numerical results agree well with the experimental results. The first moment shows the pressure propagating outward process of high-pressure bubble, in which the hole and free surface hold an obvious effect on the shock wave. At the second moment, the free surface climbed along the hole boundary. When the bubble expands to a particular volume, caused by inertia and pressure difference, the pressure inside the bubble turns to be lower than the surrounding ambient pressure, and the bubble starts to collapse with a downward jet. The free surface near the hole shrinks at the same time, which is consistent with the experimental results. At the fourth moment, the bubble emerges tearing deformation, and the water hump is also well simulated carefully simultaneously. In summary, the accuracy of the numerical method in this paper is very reliable.

3.2 Convergence Test

In numerical simulation, the convergence of the results is a vital factor that researchers must consider. In general, the grid size has a large effect on the results of numerical simulation. When the grid size is small, the calculation efficiency is reduced. And when the grid size is large, the calculation accuracy is poor. In this paper, grid length L = 0.04At, 0.03At, and 0.02At are selected for comparison. For the case in this section, At = 1.0 m, Rw is 2.0At, Rm is 1.0At, Rh is 1.25At, d is 0.63At, and the density of the floating body is half that of water. While the grid size is 0.02At, the result has reached the corresponding convergence standard, as showed in Fig. 5. So the grid of this size is used for all cases. In calculation, the size of the calculation domain is 7.5At×12.5At, and the overall domain consists of 150000 elements. The boundary of the computational domain also is a numerical non-reflective boundary condition.

Figure 5: The time evolution of the equivalent radius curves, to study the convergence with different grid sizes

3.3 The Verification of Pressure

Flow field pressure is verified through the explosion experiment from literature [12]. The comparison diagram between the numerical and experimental results is showed in Fig. 6. The working condition is the explosion of 35 g Hexocire explosive with 3.5 m from the free surface. The pressure measuring point is located at 0.7 m from the center of the explosive, and specific experimental parameters can be obtained from the literature [12]. Standard TNT in the simulation can substitute for the Hexocirex explosive by the coefficient 1.23. In the numerical simulation, the initial TNT explosive is 43.05 g, the explosive density is 1630Kg/m3, the computational domain size is 3.0 m × 6.0 m, and the length of element is 0.005 m. Through the comparison, the shock wave pressure and bubble pulse pressure can be calculated better. It is proved that the pressure change of the flow field can be simulated well by the numerical model in this work.

Figure 6: The comparison diagram between the numerical results and experimental results [12]

4  Results Analysis

4.1 Basic Features of Interaction

In this work, the nonlinear coupling characteristics will be researched systematically. The differences in the hole size, floating body density, and explosive location can affect the corresponding physical process. So the effects of these factors are analyzed in detail and qualitatively. In the last section of this chapter, the buoyancy that account for the scale effect is also discussed. In order to analyze the different factors through the method of controlling variables, a reference case is selected firstly. And the parameters of the reference case have been described in Section 3.2.

Figure 7: The diagram of free surface and flow field pressure at t* = 0.016, 0.41, 0.83, 1.30, 1.67, and 2.10 for case where Rh *=1.25, ρ*=0.5, and d* = 0.63. The jet, jet penetration, fluctuation, and spike skirt are showed in above figure

It can be seen from Fig. 7 that the bubble gradually expands to the maximum bubble radius, from the initial time to t* = 0.41, which is the result of bubble internal pressure continuously propagates outward. When the bubble expands to a certain volume, a jet away from the free surface is visible at t* = 0.83, due to the induction of the free surface. And there is a small high-pressure zone above the jet. The jet continuously moves downward, until the bubble is broken at t* = 1.30. With the bubble breaking, a high pressure temporarily propagates outward once again. At both moments of t* = 1.67 and 2.10, the bubble exhibits compressible effects. The fluctuation and spike skirt phenomenon occur at the free surface. Throughout the physical process, the energy generated by the bubble acts on the floating body, which manifests itself as the oscillation of a floating body. The floating body movement is accompanied by the expansion, collapsing, and rebound of the bubble, as showed in Fig. 8. In addition, the results with or without damping are compared, and it can be seen that the damping term can effectively reduce the numerical oscillation.

Figure 8: The floating body axial velocity curve for case where Rh *=1.25, ρ*=0.5, and d* = 0.63

4.2 Influence of the Hole Size

The impact pressure of the transient shock wave will cause the damage of structures, which will emerge a few holes with different sizes. Different hole sizes will change the characteristics of bubble motion, so it is of great significance to study the influence of the hole size. In this paper, compared with the reference case in Section 4.1, the hole parameter Rh* is selected to be 1.0 and 1.5 for numerical simulation, respectively. The free surface curves at typical moments for cases with different hole sizes are showed in Fig. 9.

Figure 9: The free surface curves at typical moments for cases with different hole sizes. The induction of structure, cavity, and crushing of free surface is showed in above figure

The effect of hole size reflects the coupling between the bubble and solid boundary. It can be seen from Fig. 9 that the size of the hole has a significant influence on the changes of the free surface. Compared with the reference case, when the radius of the hole Rh less than the maximum radius of the oscillating bubble, the dynamic effect of the broken floating body on the bubble embodies significantly. At t* = 0.52, it can be seen that the structure can induce the deformation of bubble, resulting in the no-smooth appearance of the bubble. At t* = 1.28, the air above the system rushes into the sides of the hole, and the cavity appears when Rh*=1.0. It is caused by the local pressure imbalance of the flow field. At t* = 1.82, the bubble movement comes to the middle-late stage. The free surface crushing and violently motion are effectively transmitted for the case where Rh*=1.0, and this case is relatively extraordinary. The effective simulation of free surface crushing reflects the superiority of the numerical method in this paper, which cannot be simulated by the traditional numerical method, such as the boundary element method.

The influence of the hole size is effectively discussed in this work by Fig. 10. It can be seen from the bubble equivalent radius diagram that when Rh*=1.0, the bubble equivalent radius changes more obvious. And the corresponding overall trend of this case is different from those of other cases. It proved that the induction of the broken floating body affects the development of the bubble oscillation. When Rh*>1.25, the floating body has a small effect on the pulsation of the bubble. For the time of the jet penetration, it is gradually delayed as the hole radius decreases. It can be seen from the floating body velocity diagram that when Rh*=1.0, the movement amplitude of the floating body is more severe, and the speed peak of this case is twice that of case where Rh*=1.25. From the perspective of the coupling mechanism, the axial energy generated by the bubble acts more on the floating body when the hole radius is small, causing the intense movement of the floating body.

Figure 10: The curves of the bubble equivalent radius (left) and floating body velocity (right) with different hole radiuses

4.3 Influence of the Floating Body Density

In accordance with the different floating body density, the floating body has various circumstances on the free surface, including floating and suspension. The different circumstances of the floating body will form different boundary conditions, which will induce some new physical features. In this paper, the density parameter ρ* is selected to be 1.00, 0.75, and 0.25 for numerical simulation, respectively. The evolutions of pressure and free surface for the case where Rh*=1.25, ρ*=1.0, and d* = 0.63 are showed in Fig. 11.

Figure 11: The diagram of free surface and flow field pressure at t* = 0.24, 0.62, 0.96, 1.32, 1.61, and 2.06 for case where Rh *=1.25, ρ*=1.0, and d* = 0.63. The spike skirt, jet penetration, small bubble, and multiple spike skirt are showed in above figure

It can be seen from Fig. 11 that the floating body density owns a certain influence on the evolution of the free surface. At t* = 0.96, the climb of the free surface near the hole is more apparent. And the climbing of liquid forms a new spike skirt after t* = 1.32. The new spike skirt phenomenon is a fascinating feature for this case, which leads to multiple spike skirt in the late stage of bubble movement. The mechanism of this new spike skirt in Fig. 11 is different from that in Fig. 7. The previous spike skirt is a peculiar appearance caused by the coupling of the hole and free surface. And it occurs due to the free surface climbs to the upper surface of the floating body, during the expansion of the bubble. The last spike skirt is caused by the rebound of the bubble. Apart from this, the differences of free surface and flow field pressure between this case and the reference case are not obvious.

Through the Fig. 12, the influence of the floating body density is quantitatively analyzed in this paper. It can be seen from the bubble equivalent radius diagram that the floating body density has a small influence on the motion of the bubble. Between t* = 0.0 and 1.0, there is almost no difference in the bubble equivalent radius for the four cases. At the end of the first period of bubble motion, with the decrease of the floating body density, there is a slight advance for the time of bubble jet penetration. It can be seen from the velocity diagram of the floating body that the density of the floating body has a significant influence on this process. With the decrease of the floating body density, the amplitude of floating body movement increases obviously. When the two objects are subjected to the force with same magnitude, it is known from the kinematics formula, that the acceleration of the object with a small mass is larger. That is, the corresponding movement response is more severe. In a word, a reasonable prediction is made by the fluid-structure interaction model of this paper.

Figure 12: The curves of the bubble equivalent radius (left) and floating body velocity (right) with different densities

4.4 Influence of the Explosive Location

In actual naval battle, the damaging effect of the underwater explosion is not linear with the explosive location. Therefore, it is of great significance to choose a reasonable explosive location. The fluid-structure interaction mechanism with the changes of the explosive location is analyzed in this paper. Combining the reference case in Section 4.1, the explosive location parameter d* is selected to be 1.10, 1.50, and 1.90 for numerical simulation, respectively. Then, the corresponding law is summarized. The free surface curves at typical moments for cases with different explosive locations are showed in Fig. 13.

Figure 13: The free surface curves at typical moments for cases with different explosive locations

The effect of explosive location reflects the coupling between the bubble and free surface. It can be seen from Fig. 13 that the depth of the explosive location has a great influence on the changes of the fluid field. For the reference case showed in Fig. 7, the water hump phenomenon occurs at the late stage of the bubble motion. However, as the explosive location moves away from the free surface, this phenomenon does not occur, and the free surface is relatively stable, as showed in Fig. 13. As d increases, the energy of the bubble propagating to the free surface, cannot cause the massive motion of the surface fluid. In addition, there are other several significant differences in bubble characteristics. As showed at t* = 0.31, when d* > 1.1, the bubble shape is more similar to a sphere in the early stage of the expansion process. With the change of d*, three jet modes are effectively displayed at t* = 1.80, including the counter jet, downward jet, and upward jet. As demonstrated at t* = 2.20, when the explosive location moves away from the free surface, the secondary rebound process is more intense after the bubble jet penetration.

Through Fig. 14, the influence of the explosive location is vividly displayed in this work. It can be seen from the bubble equivalent radius diagram that the explosive location has a great influence on the movement of fluid. For the oscillating bubbles with different explosive locations, there are differences in the time at which the bubble reaches the maximum radius and the jet penetration. As the distance between the bubble and free surface increases, the time is correspondingly delayed. And the change of the bubble maximum size is nonlinear due to the induction of free surface and hydrostatic pressure. It can be seen from the velocity diagram of the floating body that the broken floating body’s movement response is different, when the bubble migrates towards, and move away from the free surface. When d* = 0.63, the movement amplitude of the floating body is relatively small, which is determined by the action mechanism of the bubble energy. The energy of the bubble acts more on the radial direction of the floating body with circular hole under this condition. Finally, while d* = 1.10, d* = 1.50, and d* = 1.90, the initial floating body motion is consistent, and the floating body motion is determined by the phase of the bubble pulsation after t* = 0.75.

Figure 14: The curves of the bubble equivalent radius (left) and floating body velocity (right) with different explosive locations

4.5 Influence of the Buoyancy Effect

For small-scale bubbles, the buoyancy has less effect on them. But for the underwater explosion bubbles studied in this paper, the bubble size is very huge, and the buoyancy will have an excellent influence on the motion characteristics of bubbles. So it is reasonable to analyze the buoyancy effect through a non-dimensionless parameter. The buoyancy effect is investigated by the buoyancy parameter δ, which has been introduced in the Section 2.3.

Figure 15: The diagram of non-dimensional bubble location when the jet penetration

Combining the reference case in Section 4.1, δ=0.32, δ=0.55, and δ=0.71 three cases are chosen to study the influence of buoyancy in this work. Through the diagram of non-dimensional bubble location when the jet penetration, the effect of buoyancy on bubble movement can be seen in Fig. 15. While the bubble is larger, the time of the jet penetration is advanced, and the closer the bubble center is to the free surface. That’s because when the bubble is more gigantic, the buoyancy of the bubble is greater. Under the effect of buoyancy, the bubble has an upward movement response as a whole. In addition, the free surface curves at the later stage for cases with different buoyancy parameters are showed in Fig. 16. It can be concluded that the buoyancy of fluid can intensify the evolution of the flow field, and produce a lifting effect on the floating body in this condition.

Figure 16: The free surface curves at the later stage for cases with different buoyancy parameters

5  Conclusion

In this paper, the fluid-structure interaction model of oscillating bubble and floating body with circular hole is established by combining the Eulerian finite element method, the penalty immersed boundary method, and the VOF method. This complex nonlinear problem is effectively analyzed by numerical simulation. In our work, the numerical method is verified by the electric discharge bubble experiment and explosion experiment. The effect of factors including the hole size, floating body density, explosive location, and buoyancy on the fluid-structure interaction characteristics is summarized. It is showed from the results that the four factors play essential roles simultaneously.

The effect of the hole size reflects the coupling between the bubble and solid boundary. When the radius of the hole Rh less than the oscillating bubble’s maximum radius the changes in the whole system are incredibly intense, and the free surface crushing will emerge. When Rh*>1.25, the floating body has a small effect on the pulsation of the bubble. The influence of the floating body density is relatively small. The free surface climbs to the upper surface of the floating body during the expansion of the bubble, when the floating body has the same density as the water. And it leads to the multiple spike skirt in the late stage of bubble movement. The effect of explosive location reflects the coupling between the bubble and free surface. The near free surface bubble induces the occurrence of the water hump phenomenon. And when the explosive location is different, the form of the bubble energy acting on the floating body is also significantly different. The energy of the bubble acts more on the radial direction of the floating body, when the explosive location parameter is small. Taking into account the role of buoyancy, while the bubble is more gigantic, the time of the jet penetration is advanced, and the closer the bubble center is to the free surface. And the buoyancy of fluid can produce a lifting effect on the floating body. In summary, the effective numerical simulation of the coupling problem among the incomplete moving boundary, free surface and bubble is carried out in this paper.

Funding Statement: The National Key R and D Program of China (2018YFC0308900), the National Natural Science Foundation of China (Grant No. 51879050), the Industrial Technology Development Program (Grant No. JCKY2018604C010), and the China Postdoctoral Science Foundation (Grant Nos. 2016M600244, 2018T110271).

Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.

This work is licensed under a Creative Commons Attribution 4.0 International License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.