Delayed bubble entrapment during the drop impact on a bounded liquid bath

Numerical researches on the impact of a drop onto a bounded liquid bath are carried out using GERRIS code. A particular region of delayed bubble entrapment is presented. During a typical impact, a first crater, a liquid column and a second crater appear successively. It is found that gravity, restricted bath diameter and impact velocity play important roles in the whole process. With a proper bath diameter, the outward propagating waves at the liquid surface can be reflected back to the center in time and lead to a liquid column. The delayed bubble can be entrapped when the second crater is deep enough. When the bounded liquid bath is too narrow or too broad, no delayed bubble is formed, as the liquid column breaks before the second crater becomes deep enough. The mechanism of delayed bubble entrapment is explained through surface profiles and velocity fields.


I. INTRODUCTION
Various phenomena come out after a liquid droplet impacts onto a liquid surface, such as floating, bouncing, coalescence, jetting, splashing, and bubble entrapment. 1 These phenomena are associated with a lot of applications, such as metallurgical processes, 2 spray cooling, spray printing, and dissolved oxygen increment. 3,4 When a droplet impacts onto a broad deep bath, phenomena vary with the Weber number We and the Froude number Fr. They are defined as where ρL is the liquid density, D the drop diameter, V its velocity, σ the surface tension and g is the acceleration due to gravity. With a low impact velocity (Fr < 7), coalescing, floating and bouncing occur without bubble entrapment or jet. [5][6][7] Primary bubble occurs in higher velocity region (36.2Fr 0.186 < We < 48.3Fr 0.247 ), together with a high-speed jet. 2 Oguz and Prosperetti 8 explained the formation of bubble was determined by a delicate balance between the times at which the outward motion of the crater walls was reversed at different positions. Large bubble is another bubble phenomenon whose size is as large as the droplet. Thoraval et al. 9 used experiments and simulations to show that a concentrated vortex-ring controlled the formation process. At even higher velocities close to the terminal velocity, a crown rises high and close, forming a canopy bubble that floats at the surface. 10,11 Besides the impact velocity, restrictions of liquid zone, including bath depth in vertical direction and surrounding walls in horizontal direction, play important roles on the impact process. Based on the dimensionless ratio h * , defined as where h is the bath depth, the vertical restriction can be divided into three regions: deep baths (h * >>1), shallow baths (h * ≈1) and thin films (h * <<1). 6,12 Vander Wal et al. 13 investigated liquid baths with 0.1<h * <10, they found that thin films promoted splashing while shallow and deep baths inhibited it. Zou et al. 14 found that the success probability of drop bounce almost remained constant, as the depth exceeded five times the maximum crater depth. Nevertheless, restriction in horizontal direction is rarely studied, until Zou et al. did the pioneer works. 1,15 The behaviors of drops impacting on bounded liquid baths with various Dt * were studied using a high-speed video camera. The nondimensional bath diameter Dt * is defined as where Dt is the bath diameter. Different from the broad surface, the outwards propagating waves can be reflected back to the center and promote the flow processes. The outcomes, such as bouncing, large bubble entrapment and canopy, are dependent on the distance between the surrounding walls and the impact point. For the bouncing phenomenon, if the reflected waves reach the impact point before the drop rebounces from the liquid surface, the restitution coefficient will increase to about four times than that on broad surface. The canopy phenomenon can be observed frequently in small liquid baths, as the boundary walls stop the expansion of the crater and squeeze liquid into the crown. A unique phenomenon named delayed bubble was first observed by Zou et al. 1 It is quite similar with the primary bubble, [5][6][7] which is formed after the collapse of the first crater. However, the delayed bubble is formed after the crater-column-crater oscillation. Both the oscillation and the delayed bubble phenomenon only happen when a droplet impact onto a bounded liquid bath, and the mechanism is still unclear.
In this paper, a study of the dynamics of a drop impact onto a bounded liquid bath is presented. A numerical method is used to analyze the underlying mechanism. An approximate region of delayed bubble entrapment is carried out. Effect of gravity, diameter of the liquid bath, and impact velocity are investigated in detail. The liquid behaviors can be characterized by two nondimensional parameters: the position of the bath at the wall boundary H b * , and the height of the surface at the center Hc * . They are defined as where H b is the position of the bath at the wall boundary, and Hc the height of the surface at the center. At the initial status we have H b * =0, Hc * =1.05.

II. EXPERIMENTAL RESULTS
The delayed bubble was first discovered by Zou et al. 1 They investigated a droplet impact upon a bounded liquid surface. Highly purified water was held by glass tubes of inner diameters Dt=7.9mm, 12.0mm, 17.0mm and 26.0mm, while droplets with uniform diameter D=2.64mm were produced with a stainless steel needle. A high speed video camera (FASTCAM-ultima APX, USA) and a Nikkor 60-mm micro lens are employed to capture the images. The speed of the high-speed camera is 2000fps (frames per second) and the resolution is 68μm/pixel. Flickerless backlighting is produced by a high-intensity LED lamp with a thin sheet of drafting paper as a diffuser. The environment pressure was kept at 101.3kPa. The viscosity and surface tension coefficient of the liquid were 0.8937 × 10 −3 Pa⋅s and 0.07197N/m at temperature of 25 ○ C, respectively. The static contact angle between the tube wall and the liquid was 60.8 ○ with uncertainty of 1 ○ , which was measured from the images. The impact velocity was measured from the images with an uncertainty of 0.04m/s. The maximum standard deviation of drop diameter is 0.4mm. The error of Weber numbers is estimated to be less than 0.2. Each experimental point was performed several times (3 times as a minimum) to ensure the repeatability. 1 A typical case with a bath diameter Dt * =4.54 (Dt=12.0mm) and We=73.3 is shown in Fig. 1. After the impact, a first crater, a liquid column and a second crater appear in turn, which can be described as crater-column-crater oscillation. Here we name the jet which is totally above the surface as liquid column. Then a necking process (from 65ms to 70ms) occurs during the collapse of the second crater, and finally entraps a delayed bubble at 71ms.

A. Method and settings
We perform numerical simulations to explore drop impact on a bounded liquid bath using the GERRIS code. 16,17 It is an open source software, and it has been validated by researchers on similar drop impact phenomena. 9,18 The governing equations of incompressible Newtonian fluids can be written 17 ∂tρ + ∇ ⋅ (ρu) = 0, with u the fluid velocity vector, p the pressure, ρ ≡ ρ(x, t) the density, μ ≡ μ(x, t) the dynamic viscosity, σ the surface tension coefficient, D the deformation tensor defined as Dij = (∂iuj + ∂jui)/2. Note that x is the position vector. The Dirac distribution function δs expresses the fact that surface tension is concentrated on the interface, κ and n the curvature and normal to the interface. For two-phase fluid, the density and dynamic viscosity can be defined as where ρ 1 , ρ 2 and μ 1 , μ 2 the densities and viscosities of the first and second fluids, respectively. c(x, t) is the liquid volume fraction of the first fluid (VOF method), defined as 1 the first fluid exits at x 0 the first fluid not exits at x, while 0 < c(x, t) < 1 represents the interface. Fieldc is either identical to c or is constructed by applying a smoothing spatial filter to c. The density advection equation can be replaced by an equivalent equation for the volume fraction For the original continuum-surface-force σκδsn, it can be calculated approximately by the following method A staggered in time discretization of the volume fraction/density and pressure leads to the following formally second-order accurate time discretization Space is discretized using a graded quadtree partitioning (octree in three dimensions). The detail numerical methods can be found in Stephane Popinet's works. 16,17 GERRIS code solves the Navier-Stokes equations with a volume-of-fluid method for two-phase flow, by a dynamic adaptive grid refinement. We use axisymmetric simulation method in the GERRIS code, the computational domain is shown in a symmetrical form in Fig. 2, only the right half is for calculation. The left boundary is the axisymmetric axis. The top boundary is set outflow (p = 0, ∂u ∂y ). The right side and bottom boundaries are set no-slip solid walls. Concerning the contact angle makes the problem too complicated, which cannot be completed with the GERRIS code. It should deserve a separate study and is therefore out of the scope of the present research. Consequently, we use the default constant 90 ○ contact angle in GERRIS code. A spherical droplet with  diameter D is initialized set above the liquid bath (0.05D), with a certain impact velocity V. We add several calculating boxes (each box is a square with a side length of 0.5Dt) in vertical direction to form a rectangle computational domain, whose width and vertical height are 0.5Dt and L 1 +L 2 , respectively. The height above liquid bath L 1 is large enough (>5D, higher than the maximum column height) to capture the whole process and the bath depth L 2 is large enough 12 (>10D) to avoid the influence of the bottom. Hobbs and Osheroff 12 studied the collision of a droplet onto liquid pools with different depths. They found that with depth of 25mm or larger, the sequence of results is the same as that of splashing on a deep pool. Considering the calculating time, we set the pool depth as L 2 >10, note that, 25mm/2.64mm=9.47. Two tracers are used in the simulations, a main tracer to identify the air-water interface, and a passive tracer to identify the liquid originating from the drop. The two tracers are independent of each other.

B. Grid dependence
GERRIS code have the dynamic adaptive grid refinement function, by setting the coarsest grid level and finest grid level. The case Dt * =4.54, We=73.3 is calculated using four kinds of grid to test the grid dependence, which are shown in Table I. The grid cell size can be calculated as 0.5Dt/2 n , n is the refinement level, 0.5Dt is the side length of the calculating box. The comparison of grid (I) and (III) is shown in Fig. 3, the left side is the axisymmetric boundary and the right side is the wall boundary. The grid is refined near the air/water interface, the passive tracer interface, and in the regions of high vorticity. It is observed a coarse grid could not capture the physics of the process very well, the delayed bubble phenomenon is not observed as shown in Fig. 4 for grid (I). However, with fine grid (III) (shown in Fig. 5) and (IV), the delayed bubble can be captured well. As a matter of fact, it is the grid cell size (0.5Dt/2 n ) that influences the simulation. 19,20 Thus we choose different refinement levels n for different Dt to keep the grid cell sizes approximately equal to that in grid (III), which gives consideration to both precision and computation time, shown in Table II. The finest grid cell size is 11μm approximately the same as 15μm in the work of Deka et al. 20

C. Validation
A typical time evolution of delayed bubble entrapment (Dt * =4.54, We=73.3) using grid (III), is shown in Fig. 5. The numerical results are in good agreement with the experimental results. A thin air sheet is formed and become a little bubble (shown in Fig. 5 at 5ms) during the coalescence of the drop and the bath, which is called the Thoroddsen bubble. 21 The sequence of crater-columncrater oscillation is identical to that in experiment, and finally a delayed bubble is entrapped. The delayed bubble is entrapped at 71ms for experiment and 77ms for numerical simulation. This difference may be caused by the neglect of contact angle, which affect the liquid structure close to the boundary walls.

D. Parameter space
We explore the effect of Dt * and velocity on delayed bubble entrapment. The We-Dt * parameter space is shown in Fig. 6. Numerical results are generally in accordance with the experimental results. The dash oval represents an approximate region of delayed bubble. The delayed bubble forms in the range from Dt * =4.0 to Dt * =5.0, with proper impact velocities.
As can be seen in Fig. 6, the liquid column breaks when the Weber number is too high. For instance, when Dt * =4.54 the critical Weber number is 94.8 in the experiment (shown in Fig. 7(b)) and 105.6 in the simulation (shown in Fig. 8(b)). Note that the bubble in Fig. 7(b) is the primary bubble. When the Weber number is above this critical We, the column-break phenomenon always occurs. The delayed bubble can also be entrapped with columnbreak phenomenon, 1 as shown in Fig. 7(c). However, in our simulation, the delayed bubble disappears when the critical We is exceeded, as shown in Fig. 8(c). This may be due to the neglect of contact angle which causes capillary forces at the boundary. In experiments, the upward capillary forces restrain the rising of the liquid column and help to break it. In simulations, there are no capillary forces due to the 90 ○ contact angle. The critical Weber number increases and the column-break time delays (We=94.8 at 48ms in the experiment, and We=105.6 at 63ms in the simulation). Comparing Fig. 7(c) at 43ms to Fig. 8(c) at 50ms, it can be found that the column is much thicker in the simulation, which leads to quite different crater dynamics. The outside liquid keeps rising ( Fig. 8(c) at 68ms) but not descending and converging to the center ( Fig. 8(a) at 74ms), no necking process is formed. This indicates that our simulations are not competent column-break phenomenon.
Consequently, this research focuses on the region below the critical We.

A. Role of gravity
The role of gravity in the entrapment process is studied and presented as follow. It can be simply accomplished in numerical simulations by setting gravity to 0. Typical sequences are shown in Fig. 9 and time evolutions of center and boundary are shown in Fig. 10. Although both cases have the similar process, delayed bubble is not formed without gravity. The Froude number of droplet (D=2.64mm, We=73.3), Fr = V/ √ gD = 9 indicates that gravity can be neglected in the first few milliseconds, which is similar to large bubble entrainment. 9 This can be evidenced by the coherence of boundary and center movements before 12ms, shown in Fig. 10. However, as time goes on, gravity may play an important role in the later process. After 12ms, the curves of the center and boundary begin to bifurcate. Without gravity, only surface tension dominates the oscillation, which causes a longer oscillation period and a larger amplitude. The second crater is in V-shape instead of U-shape. Finally the crater bottom rebounds before side walls retract, consequently no delayed bubble is entrapped.

B. Propagation of energy
After a drop impacts onto a bounded liquid bath, the impact kinetic energy E k is converted into gravitational potential energy Epg and surface energy Eps, also portion of the energy dissipates. As we use the axisymmetric simulation method, the velocity can be written as u=(u, v), u and v are the radial and axial velocities, respectively. The liquid energy can be defined in the cylindrical coordinates as where E k is the kinetic energy, Epg the gravitational potential energy, Eps the surface energy, Em the mechanical energy, c the volume fraction of liquid phase. x and y are the radial and axial coordinates, respectively. The initial gravitational potential energy E pg0 can be written as The kinetic energy E k and the gravitational potential energy Epg are calculated by the "sum function" in the GERRIS code. For the surface energy Eps, S is the liquid-air interface area and S 0 is the initial interface area. The interface profile line y = f (x) is first extracted, then the area S and S 0 can be simply calculated by If the profile line overturns, then it needs to be separated into several parts, each part will be calculated as above. All energy is nondimensionalized by the initial kinetic energy Em * =Em/E k0 , E k * =E k /E k0 , Epg * = Epg/E k0 , Eps * = Eps/E k0 . At the initial status we have Eps * =Epg * =0, Em * =E k * =1. According to the Navier-Stokes equations, the energy dissipation equation for incompressible fluid without heat conduction can be derived as where ε is the internal energy. It can be found that, the mechanism energy dissipates and finally becomes internal energy. Time evolutions of energy are shown in Fig. 11, after drop impact with We=73.3 onto a liquid bath (Dt * =4.54). Generally, the mechanical energy Em decreases constantly, while the kinetic energy, gravitational potential energy and surface energy exhibit oscillatory behaviors. Gravitational potential energy and surface energy are in phase, while in phase opposition with kinetic energy. The oscillatory behaviors correspond to the crater-columncrater oscillation. This oscillation does not occur when it is an infinite surface. The time sequence and the center movement when We=73.3 and Dt * =20.0 are shown in Fig. 12. This later value is large enough to be regarded as infinite surface. A high-speed thin jet is formed at 15ms. Under the action of gravity and surface tension, the thin jet gradually slows down and falls back to the bath. The maximum height is reached at 27ms, marked as t minf , shown in Fig. 12 (right). No second crater is formed and no delayed bubble is entrapped. Without boundary walls, the impact energy is propagated outside by waves without reflection after the impact. Outlines of first few milliseconds after the impact of Dt * =20.0, We=73.3, is shown in Fig. 13. Two main kinds of waves can be found, 2 small amplitude front wave (small arrows in Fig. 13) and large amplitude swell wave (large arrows in Fig. 13). The front wave is a disturbance formed during the initial coalescence of the droplet and the liquid bath, 2 and the swell wave is formed during the expanding of the first crater. A tongue is created above the cavity at 3ms, it erects up and moves outward under the action of surface tension, which forms the swell wave. In certain conditions, this tongue can induce a large bubble. 9 The horizontal crest motions of the front wave and the swell wave are shown in Fig. 14, whose velocity turns a constant after 3ms and 7ms, respectively. Using the least square method to calculate the slope, we have the average velocity V f ≈ 0.61m/s (front wave) and Vs ≈ 0.33m/s (swell wave), both higher than the velocity of capillarygravity wave (0.23m/s). As the waves can still be found when the gravity is set to 0 in the simulation, the waves are predominantly of capillary nature according to the linear solution of the water waves. 22 Thus we have the corresponding wavelength λ f = 2πσ ρV 1 2 ≈ 1.2mm and λs = 2πσ ρV 2 2 ≈ 4.2mm. λ f is about the radius of the droplet, different from Leng's works. 2 It shows that the wavelength of outward capillary wave is about the diameter of the droplet. This difference may be due to the front wave we study here is unmeasurable in experiments as it dissipates quickly. Thus Leng's outward wave should be the latter one. λs is about the diameter of the first crater, which indicates the relationship between the swell wave and the first crater.

C. Crater-column-crater oscillation
When the liquid bath is bounded, the waves will reflect back to the center with the energy. As shown in Fig. 13, the propagating energy is mainly carried by the swell wave, and it could promote the upward motion of center. According to V f and Vs, we can predict the time t cf and tcs when the front wave and swell wave reach the center. As shown in Fig. 15 (left), the reflected energy reaches between the blue line (front wave t cf ) and the red line (swell wave tcs), it delays when Dt * increases. For Dt * ≤6.0, the reflected energy converges to the center before t minf , which gives extra energy to the thin jet and finally forms a liquid column. This liquid column can be found in Fig. 16(a) at 43ms and Fig. 16(b) at 52ms. However, the thin jet fails to transform to liquid column when Dt * >6.0, since the reflected energy reaches the center too late. Note that, when Dt * =7.0, very few energy reaches before t minf , as the main energy are carried by the swell wave. The extra energy delays the column peak time from t minf to tm, and increase the column peak height from H cmaxinf * to Hcmax * . With increasing Dt * , both the peak time and the peak height tends to the value of infinite Dt * .
When Dt * ⩽6.0, the liquid column latter collapses and causes a second crater, as the liquid bath is confined in boundary walls. This second crater is a prerequisite for delayed bubble entrapment, as shown in Fig. 16(a) at 70ms and Fig. 16(b) at 82ms. When Dt * >6.0 no second crater occurs as no liquid column is formed, as shown in Fig. 16(c) at 50ms and Fig. 16(d) at 50ms.

D. Dynamics of necking process
After a drop impacts on a bounded liquid bath, a second crater could forms after the liquid column's collapse. However, as shown in Fig. 6, the delayed bubble can only be entrapped in the red dash region of the phase diagram. The question is therefore why not all the second craters entrap a delayed bubble?
Numerical results of drop impact on the liquid bath Dt * =4.54 with varying impact We are shown in Fig. 17 bubble occurs when We=73.3 and 84.1. As shown in the second line of Fig. 17, for We<105.6, the outside liquid moves downwards and converges to the center during the second crater's collapse. This converge causes the necking process, shown in Fig. 18. As shown in Fig. 17, the liquid column turns higher and thinner with increasing We, and finally breaks at We=105.6. When the liquid column breaks, the second crater becomes too shallow to form a necking process, and also the detached drop disturbs the collapse of the second crater. The collapse processes of the second crater are shown in Fig. 18 and the H cmin * -We map is shown in Fig. 19. Generally, the center of second crater descends and then rebounds, while the side walls of the crater keep converging to the center. For We<73.3, the second crater could go deeper with increasing velocity (Fig. 19), which forms a sharper bottom at the time just before it rebounds (Fig. 18, We=41.1, 69.8). The bottom rebounds before the side walls reach each other, with no entrapment of delayed bubble. For We≥73.3, the bottom of the second crater remains at the lowest position (Fig. 18, We=73.3) or still keeps moving down (Fig. 18, We= 84.1) when the side walls collide with each other. A nipple can be seen at the bottom of the crater, then a delayed bubble is pinched off. The maximum depth decreases slightly when We≥73.3 (Fig. 19), as the bottom's downward movement is disturbed by the pinch-off phenomenon. When the impact velocity is too high (Fig. 17, We=105.6), the column breaks and the second crater becomes quite shallow. The bottom rebounds very early, before the convergence of side walls starts. This relationship between the impact velocity, the maximum depth of second crater and entrapment of delayed bubble, can also be found in other Dt * . The H cmin * -We map of the second crater are presented in Fig. 20. For a constant Dt * , before the delayed bubble is entrapped, the maximum depth generally increases with velocity. When the second crater is deep enough, marked by black circles in Fig. 20, the delayed bubble is entrapped. With higher impact velocity, the liquid column goes higher and finally breaks, marked by black triangles in Fig. 20, which causes a much shallower second crater and no bubble entrapment. However, not all liquid baths can entrap the delayed bubble with a proper impact velocity. If the bath is too narrow (Dt * <4.0 in Fig. 6) or too broad (Dt * >5.0 in Fig. 6), no delayed bubble can be formed. If the surface is too narrow, the column breaks with a quite low Weber number (We=48.3 for Dt * =3.0). Thus the liquid column cannot stored enough energy, which later cannot cause an enough deep second crater. If the liquid bath is too broad, the reflected energy is not as effective as narrower surfaces. This causes a shallower second crater in a certain velocity and cannot entrap the delayed bubble. Consequently, it needs more initial impact energy to form a deeper second crater. However, with such high velocity, the liquid column is too thin and easy to be pulled off, which causes energy loss and no delayed bubble can be formed.

V. CONCLUSIONS
In this paper, numerical simulations are performed to investigate delayed bubble entrapment after a droplet impact onto a bounded liquid bath. An approximate region of delayed bubble entrapment is presented. Effects of gravity, surrounding walls and impact velocity are analyzed. With surrounding walls, a cratercolumn-crater oscillation can occur due to the timely reflected energy. Both the surface tension and gravity restrain the liquid motion during the oscillation, absence of gravity will cause a longer oscillation period and a larger amplitude. When 4.0≤Dt * ≤5.0, the delayed bubble can be entrapped with proper impact velocities. Lower or higher impact velocity will lead to a shallower crater and no bubble entrapment. When Dt * <4.0 or Dt * >5.0, the crater cannot be deep enough to entrap a bubble, as the liquid column breaks.