I February 1999
PHYSICS
LETTERS
A
Physics Letters A 25 I t 1999) 303-310
ELSEVIER
Laser-induced bubble formation and collapse in a nonlinear Rayleigh-Taylor unstable interface in a thin layer of molten metal S. Lugomer,
A, MaksimoviC
Rudjer Bo&~viC Instifute, Bijenih
c. 54. 10001 Zqveb,
Cnuria
Received 6 July 1998; revised manuscript received 28 September 1998; accepted for publication 11 November I998 Communicated by A.P. Fordy
Abstract Bubble formation and bubble collapse zones in a laser-induced nonlinear Rayleigh-Taylor (R-T) instability were shown to follow the contours of a series of R-T shock-fronts generated on an irregular planar (terrace-like) target surface. Bubble formation occurs in the regime of target planar vaporization. Spatial variation of the bubble density distribution, the bubble size, and the bubble-bubble distance, as a function of distance from the shock front envelope, were determined. Bubble collapse occurs in the regime of planar-to-volume boiling transition and proceeds by the so-called “chain reaction” collapse mechanism inside a 2D bubble array. The contours of bubble generation and of bubble collapse were simulated by using the analytical model of Ott with variable phase and variable amplitude of R-T modes. @ 1999 Elsevier Science B.V.
1. Introduction
There are two limits interesting for the study of laser-induced nonlinear Rayleigh-Taylor instability. The first is the limit Plight/heavy + 1 (p = fluid density ), when the Kelvin-Helmholtz (K-M) instability evolves, that results in the formation of micron-scale vortex filaments, and their self-organization [ 1,2]. The second is the limit Ptighr/neavy --+ 0, when the bubbling process takes place, which is the subject of this paper. This paper deals with the formation and collapse of bubbles appearing as a result of development of the Rayleigh-Taylor instability of the interface induced by laser-metal interaction. The main aim of the work is to reconstruct processes of the formation of bubble arrays and their collapse from the experimental micrographs of the laser-treated surface
by means of scanning electron micrography and by using numerical integration of the Ott model which is applicable for thin fluid layers. When the central light fluid is generated by the heavy peripheral fluid, the R-T instability takes place. The small amplitude, or the linear behavior of the RT instability is understood by the dispersion relation between the growth rate and the wavelength of the unstable modes. In the nonlinear and chaotic stage, the unstable modes develop into bubbles of light fluid and spikes of heavy fluid each penetrating into the other phase. The interaction among modes causes the chaotic mixing of the two fluids [ 3-51. Consequently, the late phase (stage) of R-T instability is characterized by bubble formation, bubble merging, spike break-up, and turbulent mixing [ 31, as the well-known phenomena in laser-matter interactions on planar and
03759601/09/$ - see front matter @ 1999 Elsevier Science B.V. All rights reserved. PI/ SO3759601(98)00862-7
304
S. Lugomer, A. MakFimoviC/Ph.wics
spherical targets [ 451. This paper contains two important results: (i) The first one concerns the bubble formation that occurs by the propagation of the R-T shock wave over the surface in the regime of planar vaporization. (ii) The second one concerns the bubble collapse that occurs by the propagation of the R-T shock wave over the surface in the regime of volume vaporization. Bubble generation and collapse were studied by using a Q-switched Nd:YAP laser of E = 0.3 J with a pulse duration of IO ns to irradiate the target. The pulse power density Qp was varied below the volume boiling threshold ( 106-lo7 W/cm2) by the variation of the spot size, until the bubbles start to appear on the surface. For the target material tantalum plates of 1 x 1 x 0.05 cm in size have been used. The target was prepared in such a way that surface non-planarity of N lo-20 ,um was introduced. In order to demonstrate that the process of bubble generation/distribution and of bubble collapse is associated with R-T instability, we have generated a series of shocks (or the oscillatory shock wave) on an irregular planar (terrace-like) target surface. The growth rates of the instability driven by the nonuniform shock were shown to depend on the phase of the oscillatory shock at the time when the shock hits the contact surface. We have shown that the R-T pattern and therefore the bubble generation and organization depend on the phase of the R-T instability modes. Ultrafast cooling after pulse termination causes that all surface dynamical structures stay permanently frozen, thus enabling a posteriori analysis. Micrographs of the laser-treated surface were obtained by scanning electron micrography (SEM) , and then were numerically filtered.
2. Nonlinear R-T instability in the regime of planar vaporization 2.1. Shock-induced generation of the bubble array When the laser power density Qr, is adjusted to planar vaporization, the shock wave generated in the center of a Gaussian spot spreads radially over the surface causing the acceleration of a light central fluid. Acceleration of a hot, light central fluid into a periphery of dense fluid under condition Piighi/heavy --+ 0 leads to the formation of bubbles in the unstable R-T inter-
Letters A 251 (1999) 303-310
face [ 71. Some authors have studied the development of bubbles of unique size (single mode) which form a periodic array that excluded the bubble-bubble interaction in a fixed system, i.e. steady state [8-IO]. It was found that in such an ideal situation the bubble has a constant velocity v, = Ci G, where Ci is a constant, A is the Atwood number. In addition to the Atwood number A, dimensionless compressibility parameter M* = Au/C: (where Ct, is the sound speed in the heavy fluid, and A is the characteristic length of the bubbles) was introduced [ 451. Two important dynamics were found: bubble motion and bubble merging. The dominant mechanism of bubble interaction is merging, whereby bubbles of different radii propagate with different velocities and the leading bubbles grow in size at the expense of their neighbours [4]. A series of independent shock fronts (the shock formed unstable R-T interfaces) was generated by a single laser pulse on the irregular planar surface. Every shock front is the unstable interface or the turbulent mixing zone in which the bubbling has taken place. The results are shown in Fig. la and the schematic reconstruction in Fig. lb. This micrograph shows the situation of bubble dynamics reached in the moment of pulse termination, i.e. it gives a time integrated picture. Schematic reconstruction of the bubbling process which forms a 2D array of bubbles inside the shock width zone is given in Fig. 2a. The comparison is done with the bubbling process obtained by numerical simulation for the compressible case with random interface initialization by Glimm et al. [4], Fig. 2b. 2.2. Characteristics of bubbling along direction I
(in Fig. lb) Spatial variation of the bubble density p, of the bubble diameter 2R and of the bubble-bubble distance d, along the distance L counted from the shock outer envelope, is given in Fig. 3a. The bubble density distribution, p = f(L), has a peak which coincides with the shock front position. The peak maximum is N 75 ,um behind the shock front line; it decreases to a value that stays constant in the isothermal region behind the shock. The bubble diameter distribution 2R = f(L) , has a broad maximum which is shifted to the back side of the bubble density peak (see Fig. 3a), and positioned at N 150 pm behind the shock front. The bubble di-
S. Lugomer, A. MakrimoviC/Physics
Letter.? A 251 (1999) 303-310
TIP SPLlTING HURBLE
COLUMN
ENvF:LOPE INNER
a)
Ib)
ENVEI,OPF:
so(AI
,
=
oy, 4
I
Fig. 2. (a) Development of the bubble array in the unstable R-T interface (schematically). (b) Plot of the interface for R-T instability at mesh 200 x 400. The results show that the interface becomes more complicated as the mesh becomes finer. However. the large-scale structure of the fluid interface is unchanged after mesh refinement. In these simulations A = 3 and M2 = 0.2. Only the portion of the computational domain near the interface is shown here. From Ref. [ 41.
ameter distribution in the peak is 2R 5 7.5 ,um, which indicates that they have grown by a merging process of 2 (or three) smaller bubbles. Small bubbles of diameter 2R < 4 ,cLmcan be seen on the outer envelope of the shock. The bubble-bubble distance distribution, d = f(L) has a maximum which coincides with the maximum of the bubble size distribution 2R = f(L), and decreases toward the front- and toward the back side of the shock width lines in Fig. 3a. The width of the shock front corresponds to the region between the outer and to the inner envelopes of chaotic mixing, i.e. to the region defined by the closest and furthest penetrating bubbles [ 51. 2.3. Characteristics (in Fig. lb) Fig. 1. Bubbling along the shock fronts of the R-T instability. (a) A series of bubbling fronts obtained on a Ta surface by a r = 10 ns laser pulse, on the irregular planar surface, M N 250x. (b) Reconstruction of a series of shock fronts generated by a single pulse on the nonplanar (terrace-like) surface.
of bubbling along direction 2
Spatial variation of the p, 2R and d distribution which shows a series of shocks were obtained along direction 2 on the micrograph of Fig. 1a,b. Two of the peaks in the p = f(L) distribution variation which correspond to the two strongest shock waves are
306
S. Lugorner; A. MahirnoviC/Physics
SHOCK
I?”
FRONT
_.
15
._
10
E
3
I
-0
d N
5
I 100
200
300
‘loo
500
600
700
L(vm)
2. SHOCK
-
nario of Gardner et al. [ 51, initially the bubble mergers are randomly distributed. As time progresses the size of the bubbles increases, and their radial distribution in space changes with time. After an initial period of slow increase, the velocity accelerates almost uniformly as a consequence of the merging process, so that the bubble density decreases monotonically and their radius increases [ 3-51. Interestingly, the spatial variation of the bubble density distribution strongly reminds of the bubble density distribution on the shocked foam targets observed by Hoarthy et al. [ 111. Bubbling that is observed ahead of the shock front (Fig. la) may be interpreted as a consequence of preheating of material by radiation from the shock. By measurement of the shocked foam material, Hoarthy et al. [ 111 have found that temperature in the shocked region rises sharply, then levels off before increasing rapidly close to the ablation surface, while the region behind the shocks is isothermal.
60
“E $ a
Letters A 251 (1999) 303-310
3. Nonlinear R-T instability in the regime of volume vaporization
‘lo zo
3.1. Shock-induced collapse of a 20 bubble array I
I
I
I
100
200
300
400
500
603
L(km)
Fig. 3. Spatial variation of the bubbling characteristics: p = f(L), 2R = f(L) and d = f(L), (a) along direction I in Fig. lb, (b) along direction 2.
clearly seen on the positions of 180 pm and 350 pm, counted from the outer envelope of the first shock; two shoulders in the density distribution at 450 pm and 500 pm correspond to the two smaller shocks that follow behind two large shocks. The bubble size distribution 2R = f(L) shows a wavy-like growth from the shock front toward its back side in Fig. 3b. They appear as local maxima and minima of distribution and again (as above) do not coincide with maxima of the p = f(L) distribution. Distribution of the bubble-bubble distance d = f(L) shows the increase (with local maxima and minima) toward the front side and toward the back side of the shock. The above bubble size distribution is characteristic for a “resonance” behavior, or a multimodal nature, clustered about the values nyo, for the bubbles formed from the merger of n initial bubbles. Assuming the sce-
When the laser power density Qp is slightly increased, the planar vaporization is transferred into a volume one. Since the bubbles spontaneously form on the surface and form a 2D bubble array, the accelerated fluid will spread from the center to the periphery through the existing bubble array. The prominences of accelerated fluid in the 2D R-T instability will interact with the bubble arrays. Such a case is realized when the onset of surface bubbling (the transition from planar to volume boiling) occurs simultaneously with the R-T instability. In that case the R-T fronts (the series of shocks) pass over the bubbling surface. Passage of a series of shocks through a 2D bubble array results in the bubble collapse that occurs along the shock fronts, Fig. 4. The bubble collapse under transient pressure pulses is caused by the low specific heat of bubbles which results in high shock temperatures in the bubbles [ 111. Larger magnification reveals that craters which form the curved channels (generated by the collapse of bubbles) follow the contours of the R-T instability which are the contours of a high pressure and temperature. The width of the channels (the width of
S. Lugotne,: A. MaksimoviC/Ph_ysics
Fig. 4. A series of the bubble collapse parts that coincide with a series of R-T instabilities (shock fronts) are the nonplanar (terrace-like) surfaces. (a) Micrograph obtained on a Ta surface of the micrograph (M ,-- 250x ). (b) Schematic reconstruction (a).
the shock indicating of bubbles lapse of a
Lefters A 251 (19991 303-310
lapse of rectangular, horizontal, vertical and triangular array of bubbles were studied by many authors [ 12161. Dear and Field [ 131 have studied the dynamics of the collapse of 2D arrays of bubbles. They found that with the bubble clusters, the collapse proceeded step by step with pressure waves from one collapsed row then collapsing the next row of bubbles.. In particular, the shock (horizontal column) that starts with the first bubble causes its collapse and the formation of micro-jets. Once the bubble totally collapses, the rebound shock wave is formed. During this time of collapse the second bubble has been shielded from the initial shock wave and has only experienced a slight lateral compression, but when the rebound shock wave from the first bubble reaches its lower surface, it too starts to collapse to produce a jet. The third bubble in the line is collapsed in a similar way by the collapse and rebound of the second bubble. A chain reaction along a line of bubbles is then feasible given the right conditions of shock strength, cavity diameter and spacing [ 131. Applied on the rectangular array of bubbles the above process indicates the interaction scenario, in which each row directly behind the next row would allow little of the main wave front to pass the first row [ 131. Hence, the collapse of a second row of bubbles should be mostly due to shock waves radiating from the collapse of the first row [ 13 1. Therefore. there appear to be two types of interaction between bubbles in well-ordered arrays. The first occurs between adjacent bubbles in a row collapsed at the same time by a shock wave. The second is the chain reaction in a column of bubbles where one bubble causes the collapse of another. For a row of bubbles parallel to the incident shock wave the jet direction and its form are affected by the reflected tensile waves from neighbouring bubbles. These reflected waves from adjacent bubbles affect the flow of liquid into the jet and the overall collapse of the bubble. For a row of bubbles, the regions of liquid between the bubbles where the tensile reflected waves overlap experience a drop in pressure below the ambient value [ 131. 3.2. Simulation of contours of bubble generation of bubble collapse
fronts) ranges from 3-4 bubble diameters, the existence of a bubble array of 3-4 rows in the shock front. The shock-induced colsingle bubble, and the shock-induced col-
?07
and
3.2. I. Model The model of dynamics capable of describing the bubbling contours as well as of the bubble collapse
308
S. Lugome< A. MaksimoviC/Physics
contours as the R-T instability, can be obtained in the analytic way. We start with rough and simple assumptions to avoid complex theoretical considerations and numerical solutions, usual for 3D and 2D cylindrical geometry (Sakagami and Nishihara [ 171) . Therefore, we follow Ott [ 181 and consider the R-T instability of a thin fluid layer of infinitesimal layer thickness, that is stable against the inertial forces caused by the fluid pressure. The analytical model of Ott, that we used in Ref. [ 1 ] for simulation of hypocycloidal contours of the vortex filament selforganization, is also the basis for simulation of the bubbling contours and of the bubble collapse contours. Namely, the contours in both cases represent the acceleration fronts of the surface fluid in the R-T instability. We assume that the fluid particle in the beginning has the coordinates x(O) = &J and y(0) = r]a, where ~7 and 5 are the Lagrangian coordinates of the moving particle. Thus, in t = 0 the fluid particle at the boundary has the coordinates x(~,~, 0) = 60, ~(5, 7, 0) = ~0. The neighbour particle has in t = 0 coordinates x(t+dt, v+dv, 0) = to + d.$e, y(dc, dv, 0) = v. + d?;lo. In the case of the “plain” boundary we have y(dc, dr), 0) = 70, where the arguments in x, y coordinates dt, dv represent the neighbour particle: (5 + dr). The element of the boundary, ds (the arch element), is given by [ 1,181 d.s=dm,
(1)
with ds = dl (for the plane element). The mass of this element is given by dm = uo ds, where a0 is the mass per unit of area. The force of the arch element is given [ l] by dmY=
-admj-Apdr
x k,
(2)
where dr = r’ - r, r’ is the position of the neighbour particle in t, and r is the position of another particle; a is acceleration (horizontal) of the surface layer. The i, j, k are the unit vectors in x, y and z direction, respectively. Since dr = dx i + dy j, one has (dxi-tdyj)
x k=-dxj-tdyi.
(3)
Using dx = xc d[ + y,, dv and dy = yc de + y? dq, the equation of motion becomes f=
-$$@‘+y,,dq),
Letters A 251 (1999) 303-310
j;= -a+
g(xld[+x,dg),
which is equal to (by using definition
of arc length)
K = -ay, ,
(6)
j;= -afax,.
(7)
The solution is the hypocycloide, ten [ 1,181 x( s, t) = s - c
y(s, t> =
1
AkeG
Bkefi
which may be writ-
sin( ks + Bx,k) ,
(8)
cos( ks + 6$k) ,
k
where Ak and Bk are the constant coefficients and 0k is the phase shift of the kth mode which is the solution of Ott [ 181. The evolution starts from the sinusoidal perturbation in t = 0, and leads to realistic hypocycloidal structures if the model with a larger number of modes is used. For reasons of simplicity we have taken only a short modal series k;. A simulation performed with only two modes with kl = 1, k:! variable, and with the phase shifts 8,,i = 0, 8 .v.I = 0, ox,2 = 0, 8,,2 variable, has given the R-T instabilities which for t = 0.1 (nondimensional time) are shown in Fig. 5a and b. In the case studied, we have a series of shocks generated by the nonplanar, terrace-like target surface instead of a single shock (characteristic for an ideally planar target). Thus, a series of independent R-T contours intersecting between themselves should be generated. Their number corresponds to the number of terraces in the interaction space; the three R-T contours correspond to the surface with three terraces. For reasons of simplicity we have assumed that every hypocycloidal contour starts in the same time (t = 0) from the common center, with the same amplitude. Comparing the contours in Fig. 5a and b (the upper left quadrant) with the bubbling contours on the micrograph (Fig. la) and its schematic reconstruction (Fig. lb), it can be seen that the system of contours in Fig. 5b (obtained for kl = 1, k2 = 3 (curve 1) ; kl = 1, k2 = 5 (curve 2); kl = 1, k2 = 8 (curve 3) and 8,~ = $,I = OX,2= 8,,2 = 0 in all cases), is the most similar one. Simulation of the bubble collapse contours seen in Fig. 4a, was performed under different assumptions,
S. Lugomec A. MaksimoviC/Physics LettersA 251 (1999) 303-310
-1500
1
’
’
’
-000-600-400-200
’
0
’
’
’
’
,
,
,
(
200 400 600 000
X
15001
,
,
,
1 I I1 -1500 -800-600-400-200
,
I‘ 0
I
200 400 600 800
-800
-400
X Fig. 5. Hypocycloidal (9)
bubbling contours obtained from Eqs.
as a solution of a nonlinear R-T
0
800
400
X
(8),
instability. (a) Three curves
Fig. 6. Hypocycloidal Eqs. (8).
(9)
bubble collapse contours obtained
in the set are obtained under following conditions: k, = 1. kZ = 3
bubble array (beginning of the volume boiling).
(curvr I): kl = I.kz =S (curve2);
in the set were obtained for: kl = I, k2 = 5 and /I,,,
phase shifts of the modes
I
kl = I,kz =8 (curve 3). The
and 2 were as follows: IS),,, = HY,I = 0,
0r.z = 0, $.2 = n. (b) Three curves in the set are obtained under the following conditions: kl = I, k? = 3 (curve 1); kl = I,kZ = 5 (curve
2); kl = I. k2 = 8 (curve 3); while all phase shifts are
equal lo zero.
19~.z= 0. B,,z = x. The amplitude of the modes increased from curve Atc”r”e 2) _ _ A,c”Wr II Afc”r”e 3) = ,:CUCvC?) Thee
from
for the case when series of shocks pass over the 2D
I
(a) Three curves
I
= 8,,,
= 0.
and 2 has been
to curve 2 and to curve 3 as follows:
curves i; the set
t9,,, as given above. The amplitudes were increased from curve to curve in the same way.
and with different values of k and 8 parameters. The set of hypocycloidal contours was generated under the assumption that amplitudes (the coefficients Ak and Bk in Eqs. (8), (9) ) are different for every shock generated on various surface terraces. The results shown in Fig. 6a are compared to the micrograph (Fig. 4a) and show some similarity. However, it may be said that the system of contours (only 3 contours for reasons of simplicity) shown in Fig. 6b (obtained for: k, = I. k? = 4, @,,I = OX,, = OS,2 = 0Y,2 = n-) is the most similar to Fig. 4. Although the system of contours (hypocycloidal) obtained in Figs. Sb and 6b are not exactly the contours in Figs. lb and 4a, respectively, they are the most realistic ones.
Generated on the nonplanar (terrace-like) target surface, a series of R-T shocks spreads in radial direction. Two cases of R-T instability and of bubble behavior were studied: (i) The R-T instability in the regime of planar surface vaporization. (ii) The R-T instability in the regime of planar-tovolume vaporization transition. In case (i) the prominences of R-T instability move on the vaporizing surface while in case (ii), the prominences move on the surface which has already started bubbling so that they move through a 2D array of bubbles. In case (i) the generation of bubbles, while in
S. Lupmer;
310
A. Makritnovid/Physics
case (ii) the bubble collapse follows the high pressure contour of the R-T prominences. Simulation performed on the basis of the multimodal model of Ott with variation of the mode numbers, of the phase shifts and of the mode amplitudes, has satisfactorily reproduced the bubbling and of the bubble collapse contours.
Acknowledgement The authors are thankful to Professor Edward Ott, Department of Electrical Engineering, University of Maryland for critically reading the manuscript.
References [ I J S. Lugomer,
A. MaksimoviC, J. Appl. Phys.. to be published.
Letters A 251 (1999) 303-310 [21 S. Lugomer, Phys. Lett. A 242 (1998) 319. [31 D.H. Sharp, Physica D 12 (1984) 3. 141 J. Glimm, X.L. Li, R. Menikoff, D.H. Sharp, Q. Zhang, Phys. Fluids A 2 ( 1990) 2046. 151 CL. Gardner, J. Glimm, 0. McBrain, R. Menikoff, D.H. Sharp, Phys. Fluids 31 ( 1988) 447. 161 R. Ishizaki, K. Nishihara, H. Sakagami. J. Ueshida, Phys. Rev. E 53 (1996) 5592. [71 J.A. Zufiria, Phys. Fluids 31 (1988) 440. [8] G.I. Taylor. Proc. R. Sot. London A 201 (1950) 192. [91 G. Birkhoff. D. Carter, J. Math. Mech. 6 ( 1957) 6. [ 101 P.R. Garabedian, Proc. R. Sot. London A 241 (1957) 423. 11I I D. Hoarty. A. Iwase, C. Meyer, J. Edwards, 0. Willi, Phys. Rev. Lett. 78 ( 1997) 3322. 121 J.M. Vanden-Broeck, J.B. Keller, Phys. Fluids 23 (1980) 1491. 131 J.P. Dear. J.E. Field, J. Fluid. Mech. 190 (1988) 409. 141 H. Takabe, A. Jamamoto, Phys. Rev. A 44 ( 1991) Sl42. IS I N.K. Boume. J.E. Field, J. Fluid. Mech. 244 (1992) 225. 161 A.P. Szumowski, K. Falkowski, Phys. Fluids 7 ( 1995) 2529. 171 H. Sakagami. K. Nishihara. Phys. Fluids B 2 (1990) 2715. 181 E. Ott, Phys. Rev. Lett. 29 (1972) 1429.