Full text
PHYSICAL REVIEW FLUIDS 8, 073606 (2023) Theory of the jets ejected after the inertial collapse of cavities with applications to bubble bursting jets José M. Gordillo *and Francisco J. Blanco–Rodríguez Área de Mecánica de Fluidos, Departamento de Ingeniería Aeroespacial y Mecánica de Fluidos, Universidad de Sevilla, 41092 Sevilla, Spain (Received 27 February 2023; accepted 13 July 2023; published 26 July 2023) The dynamics of the axisymmetric jets originated from the bursting of bubbles in a liquid of density ρ, viscosity μ, and interfacial tension coefficient σcan be rationalized as a two-stage process in which, initially, the pressure jump ∼σ/Rbaccelerates the liquid towards the axis of symmetry, inducing a far-field flow rate per unit length Q∞∝VcRb, with Rband Vc=√σ/(ρRb) indicating the radius of the bubble and the capillary velocity, respectively. The second stage, during which a fast jet of radius Rjet (T)Rband velocity Vjet(T)Vcis ejected, is driven by the far-field radial velocity field established initially, which forces the collapse of the cavity walls while keeping Q∞practically constant in time because liquid inertia and mass conservation prevent appreciable changes of this quantity during the very short timescale characterizing the ejection of the jet. Our theoretical predictions for Rjet(T)andVjet(T) reproduce fairly well the time evolution of the jet width and of the jet velocity for over three decades in time, obtaining good agreement with numerical simulations from the instant of jet inception until Rjet ∼Rb.The analytical expressions for the jet width and for the jet velocity provided here constitute the initial conditions for the explicit solution of the ballistic equations deduced in Gekle and Gordillo [J. Fluid Mech. 663, 293 (2010)], which, hence, can be straightforwardly used in order to quantify the size and velocity of the first drop ejected and the fluxes of mass, momentum, and energy transferred from the ocean into the atmosphere. In addition, motivated by the results obtained for the particular case of bubble bursting jets, we also present here a unified theoretical framework aimed at quantifying the dynamics of the type of generic jets produced by the collapse of axisymmetric gas cavities of arbitrary shape when their implosion is forced by the radial velocity induced by a far-field boundary condition expressing that the dimensionless liquid flow rate per unit length directed towards the axis of symmetry, q∞, remains constant in time. Making use of theory and of full numerical simulations, we first analyze the case of the collapse of a conical bubble with a half-opening angle βfinding that, when the value of q∞is fixed to a constant, this type of axisymmetric jets converge towards a purely inertial β-dependent self-similar solution of the inviscid Navier-Stokes equations, described here for the first time, which is characterized by the fact that the jet width and velocity are respectively given, in the limit β1, by rjet ≈2.25 tan β√q∞τand vjet ≈3q∞/(2 tan β√q∞τ), with τindicating the dimensionless time after the jet is ejected. For the case of parabolic cavities with a dimensionless radius of curvature at the plane of symmetry rc, our theory predicts that rjet ∝(2rc)−1/2(q∞τ)3/4and vjet ∝q∞(2rc)1/2(q∞τ)−3/4, a result which is also in good *[email protected] Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article’s title, journal citation, and DOI. 2469-990X/2023/8(7)/073606(29) 073606-1 Published by the American Physical Society
GORDILLO AND BLANCO–RODRÍGUEZ agreement with full numerical simulations. The present results might also find applications in the description of the very fast jets, with velocities reaching up to 1000 m s−1, produced after a bubble cavitates very close to a wall and in the quantification of the so-called bazooka effect. DOI: 10.1103/PhysRevFluids.8.073606 I. INTRODUCTION It is known that the type of high speed liquid jets produced after the implosion of a bubble, widely used in applications related with the cleaning of surfaces [1], is partly responsible for the structural damage in cavitating flows [2] and plays a key role in the dispersion of contaminants [3,4]. These high speed jets, which share similarities with those formed by the so-called bazooka effect [5] also possess applications in medicine, where they could be used in needle-free drug injection systems [6,7] and are commonly observed after the impact of either a drop [8–13]orasolid[14,15] against a gas-liquid interface and after the collapse of bubbles [16–19]. To our knowledge, and in spite of the number of recent contributions on the subject, there does not exist a commonly accepted framework capable of predicting the dynamics of the myriad of types of fast and thin inertial liquid jets produced after the implosion of a cavity. Hence, it will be our main purpose here to present a theory to describe the dynamics of the jets ejected after the collapse of slender axisymmetric cavities whose predictions will be compared with the results of full numerical simulations carried out using GERRIS [20,21]. The velocities of the jets produced after the implosion of a cavity depend very much on the initial geometry of the interface from which they are ejected, which justifies classifying the jets in, at least, two different categories: (i) the jets produced at the axis of an initially spherical surface and (ii) the jets produced at the base of a truncated conical surface. For instance, the type of jets reported Refs. [7,22], which can be analyzed as a function of the initial near-axis nonuniform radial velocity field at the interface [23] and are employed, for instance, in ink-jet printing applications [24,25], belong to the type of jets emerging from the bottom of a spherical cap. The jets produced by the cavitation of bubbles near boundaries [26,27] also belong to the same type of jets when the so-called standoff parameter, which expresses the ratio between the initial distance of the bubble center to the boundary and the maximum bubble radius, is of order unity or larger. Indeed, when the value of the standoff parameter is O(0.1), Ref. [2] presented numerical results retaining compressibility effects in the simulations revealing that the jets produced by cavitation bubbles near a rigid boundary do not emerge from an initially smooth, locally spherical cap, but from the base of a truncated conical interface. This numerical result was later confirmed experimentally in Ref. [28], in which velocities of the order of ∼100 m/s for the case of jets produced by cavitation bubbles emanating from a spherical cap were measured, i.e., for values of the standoff parameter of order unity, and much larger velocities, of the order of ∼1000 m/s, for the case of jets emerging from the base of a truncated conical bubble, a result in agreement with the original finding reported in Ref. [2]. Moreover, Ref. [28] realized that the jets emanating from the base of the locally conical bubble are the result of a boundary-parallel cylindrical flow converging towards the axis of symmetry and, motivated by this fact, Reuter and Ohl [28] also reported the results of simplified numerical simulations neglecting compressibility effects, obtaining results in agreement with their own experimental measurements. It will be the main purpose in this contribution to describe both numerically and theoretically the “conical jets,” namely the jets emerging from the base of a truncated conical interface. The conical jets, which are comparatively faster than those emerging from a spherical cap, are found not only in the collapse of cavitation bubbles for small values of the standoff parameter but also in other types of natural flows, like in the bursting of bubbles at an interface [12,29–31] or after the impact of a drop on a liquid pool [9–12]. Indeed, in these two latter cases, before the jet is ejected, capillary 073606-2
THEORY OF THE JETS EJECTED AFTER THE INERTIAL … waves traveling along the interface towards the bottom of the cavity deform the initially spherical bubble into a truncated cone [32]. Moreover, our analysis on the collapse of conical cavities is also extended here to the case of cavities with a generic shape, paying special attention to parabolic cavities, because this is the type of void geometry found when a solid impacts a free surface [14,15]. The paper is structured as follows: Section II is devoted to present the results of numerical simulations of bubble bursting jets, which motivates the theory in Sec. III, where we analyze the ejection of jets from slender axisymmetric cavities whose collapse is forced by a constant value of the flow rate per unit length. In this section, we also compare the predictions of our theory with the results of the numerical simulations presented in Sec. II and with the numerical results in Appendix A, where we also deduce a new self-similar axisymmetric solution of the inviscid Navier-Stokes equations that, up to our knowledge, had not been reported before. Conclusions are presented in Sec. IV. II. BUBBLE BURSTING JETS: NUMERICAL RESULTS Due to their widespread presence in nature and in technological applications, an increasing number of numerical [14,29,30,33], experimental [10,11,13,19,34], and theoretical [35,36] studies have been devoted to characterize the width and speed of the jets ejected after the collapse of a cavity or a bubble, a process which is ubiquitously present in our daily life experience. A type of jets exhaustively studied in recent times because of their crucial role played in the generation of nanometric sea spray aerosol are the jets produced after the bursting of bubbles see, e.g., Refs. [4,31,37–45], which break into a number of droplets that inject mass, heat, and momentum into the surrounding atmosphere, hence, influencing climate [33,46,47]. The only result in the literature suggesting a prediction for the time evolution of the width Rjet(T) and velocity Vjet(T) of bubble bursting jets; see Fig. 1, which is the necessary intermediate step to properly quantify the different fluxes transported by these jets, is the study of Lai, Eggers, and Deike [36], where the authors suggested that the inertiocapillary self-similar solution first reported in Ref. [35] governs the dynamics of bubble bursting jets. Therefore, the solution proposed in Ref. [36] is characterized by the fact that inertial terms in the momentum equation for a liquid of density ρ, viscosity μ, and interfacial tension coefficient σare supposed to be in balance with the capillary pressure terms, a fact implying that ρV2 jet(T)∝σ/Rjet(T) and, consequently, Wel=ρV2 jet(T)Rjet(T)/σ should remain constant in time [35]: In this contribution, we show that this is not the case and deduce, for the first time, the correct equations for Rjet(T) and Vjet(T). In order to first test whether Welremains constant in time, we have performed simulations of the type depicted in Fig. 1using GERRIS [20,21] because this numerical code has been proven to accurately reproduce experimental results related with the bursting of bubbles [30,31]. In the following, the bubble radius Rband the capillary velocity Vc=√σ/(ρRb) are used to build dimensionless variables, written in lowercase to differentiate them from their dimensional counterparts, written in capitals, i.e., rjet(τ)=Rjet(τ)/Rb,vjet(τ)=Vjet(τ)/Vc, and Wel(τ)=v2 jetrjet with t=TVc/Rband τ=t−tbubble referring to the dimensionless time after the jet is ejected; see Figs. 1and 2.Herewe consider that the Bond number verifies the condition Bo =ρgR2 b/σ 1, with gindicating gravity and, therefore, the jet ejection process is controlled, for fixed values of the density and viscosity ratios see Fig. 1, by just one dimensionless parameter, namely the Ohnesorge number, defined as Oh =μ/√ρRbσ, related with the Laplace number as La =Oh−2, a parameter which is varied here within the range 625 ⩽La ⩽17 000, with each value of La identified using the color code of Fig. 2(a). The jets depicted in Fig. 1are emitted once the capillary wave with a dimensionless wavelength λ∗∝Oh1/2, propagating with the Oh-independent dimensionless speed ≈5, reaches the bottom of the cavity at the instant tbubble(La) depicted in Fig. 2(a) [32]. The capillary wave deforms the initial spherical bubble into a truncated conical surface with a half-opening angle β(La) such that β≃45◦ for La ⩾2500; see the green dashed line in the third panel of Fig. 1,aswellasFigs.3and Fig. 1in the Supplemental Material [48]. For La ⩾2500, the minimum radius of the truncated cone is given 073606-3
GORDILLO AND BLANCO–RODRÍGUEZ FIG. 1. Bubble collapse and subsequent jet ejection processes for La =Oh−2=7200 with rand zindicating the dimensionless radial and axial cylindrical coordinates. As we did in Ref. [49], the numerical simulations have been carried out for a value of the Bond number Bo =0.01 and for fixed values of the density and viscosity ratios: ρg/ρ =1.2×10−3and μg/μ =1.8×10−2with ρgand μgindicating the gas density and viscosity. The last panel shows the base of the jet, which is the point at the interface r=rs(z,τ) of coordinates r=rjet(τ), z=zjet(τ)where∂rs/∂z→∞. The jet velocity is defined as vjet (τ)=vz(r=0,z=zjet ). The origin of times, t=0, is set at the instant when the interface of the initially spherical bubble starts deforming and tbubble(La) denotes the instant of time when the capillary wave reaches the bottom of the cavity; see Fig. 2(a) by [32,49]: rjet0(La) =0.22151−Oh 0.0305,(1) where the dependence with Oh1/2comes from the fact that λ∗∝Oh1/2[32]; note also that, for our subsequent purposes, it proves convenient to define here r∗ jet0 =rjet0(La =2500) ≃0.05. It was shown in Ref. [49] and in Fig. 3that the topology of the interface changes for La <2500 since, in these cases, a tiny bubble is entrapped beneath the cavity before the jet is ejected; see Fig. 3. For La <2500, the jet is also issued from the base of a truncated conical surface with a radius rjet0(La <2500) <r∗ jet0 and with an opening semiangle which verifies β(La <2500) <45◦;see Fig. 2(b) and Fig. 3[51]. 073606-4
THEORY OF THE JETS EJECTED AFTER THE INERTIAL … (a) (b) FIG. 2. (a) The values of tbubble, with tbubble the instant of time at which the capillary waves with a wavelength λ∗∝Oh1/2reach the bottom of the cavity are almost independent of La because capillary waves propagate with a velocity which is independent of La [32]. (b) The value of the opening semiangle β depends on La: For La ⩾2500 a bubble is not entrapped beneath the bottom of the cavity and β≈45◦;for La <2500, a satellite bubble is entrapped and β<45◦;seeFig.3and the results in Fig. 1of the Supplemental Material [48]. FIG. 3. Local shape of the cavity at the instant tbubble(La) the capillary waves reach the bottom of the bubble. Notice that, for La ⩾2500, a tiny satellite bubble is not entrapped beneath the cavity and the semiangle of the truncated conical surface from which the jet is ejected is very close to β=45◦, with this angle being indicated in the figure using red dashed lines. In contrast, for 625 La <2500, a satellite bubble is entrapped and β≃37◦for La ≈1000; see also Fig. 2(b). 073606-5
GORDILLO AND BLANCO–RODRÍGUEZ (a) (b) FIG. 4. (a) Wel=rjet v2 jet 1, with Wel∝τ−1/2for τ10−3. Inset: The values of the radial flow rate per unit length calculated using GERRIS as q=−rvrfor La =2400 and different values of r/rjet(τ)andz/zjet(τ) remain nearly constant along two decades in time. (b) Streamlines calculated at τ=10−5and cavity shapes at three instants of time for La =1000. The streamlines in the liquid are mostly horizontal at the side of the cavity walls, a fact indicating that the cavity collapse is driven by a radial velocity field. In contrast, the streamlines also reveal that the flow in the gas cavity is mostly directed along the axial direction. The velocity vectors plotted at two radial positions r≈0.5atτ=10−5(blue), τ=2.18 ×10−3(red), and τ=4.18 ×10−3(black) do not change with time. The results in (b), together with the ones depicted in the inset of (a), indicate that the velocity field driving the jet ejection process cannot be approximated by the velocity field generated by a threedimensional sink with a time-varying flow rate plus additional higher-order terms, a type of approximation which, however, can be used to model the velocity field before the jet is ejected; see Ref. [50] for details. In contrast, our results in Ref. [12] reveal that the axial velocity field at the bottom of the jets generated either after either a bubble bursts at a free interface or after a drop impacts a liquid pool, can be very well approximated by the flow generated by a line of sinks located at the axis of symmetry. In Ref. [12] we also show that the good agreement between predictions and the numerical results is kept along the whole time evolution of the jet. The result for Wel(τ)inFig.4(a) show that (i) Wel=v2 jetrjet 1 for all values of τ>0, a fact already suggesting that the jet dynamics cannot be driven by capillary forces and (ii) Welis far from being a constant since Wel∝τ−1/2for τ10−3. These results indicate that bubble bursting jets cannot be described using the inertiocapillary solution of Refs. [35,36]. It could be argued that the result in Fig. 4(a) refers to magnitudes evaluated right at the base of the jet: However, our previous results in Refs. [15,49] reveal that the time evolution of the jet radius and and of the jet velocity, and even the diameter and the velocity of the drops ejected, can be calculated in terms of the functions rjet(τ), zjet(τ), and vjet(τ). Indeed, in Ref. [49], we coupled the ballistic equations deduced in Ref. [15], valid in the limit Wel1 of interest here [see Fig. 4(a)] with the mass and the momentum balances at the top drop, with this latter balance incorporating the relative flux of momentum, the capillary retraction term, and the drag force exerted by the gas, finding an excellent agreement between the numerical results and the predictions. Then, since rjet(τ), zjet(τ), and vjet(τ) do not follow the inertiocapillary scaling in Refs. [35,36], as revealed by the results shown in Fig. 4(a), we can conclude that the jets emerging from the bursting of bubbles are not driven by balance between inertia and capillarity, as is claimed in Ref. [36]. In addition, Fig. 5 reveals that the scaled shapes of bubble-bursting jets superimpose onto the purely inertial, selfsimilar solution, found in Appendix A, which clearly differs from the inertiocapillary one reported in Ref. [36] because our new self-similar solution verifies the condition that the value of the farfield flow rate per unit length remains constant in time, a fact constituting an additional evidence 073606-6
THEORY OF THE JETS EJECTED AFTER THE INERTIAL … (a) (b) FIG. 5. The time-evolving jet shapes in the top row corresponding, from left to right, to La =17 000, La =7200, La =4000, and La =2500, overlap onto the universal self-similar and purely inertial solution depicted in Fig. 16, represented here in dashed lines, for over a decade in time. Here χ=r/rjet(τ)andη= (z−0)/(zjet(τ)−0), with z−0representing the distance measured from the vertex of the cone with an opening semiangle β≃45◦;seeFigs.1and 3and Appendix A. supporting our description of the velocity field in terms of a line of sinks located at the axis of symmetry [12,32], as it will become even more clear in Sec. III. In summary, the analytical results in Refs. [15,49] permit us to explicitly calculate the spatiotemporal evolutions of the jet shape and velocity, and even the sizes and the velocities of the drops ejected as a function of rjet(τ), zjet(τ), and vjet(τ), this being the reason why the main purpose in this contribution is to deduce the equations governing the three time-dependent variables illustrated in Fig. 1. In order to explain the results in Fig. 4(a), where it is shown that Wel∝τ−1/2,wefirstreview our previous results in Ref. [12], where we calculated the initial value of the jet velocity, vjet0(La), in terms of rjet0(La) given in Eq. (1) once the velocity field is expressed as the irrotational flow induced by a line of sinks with a flow rate per unit length q∞≈const. In Ref. [12] we showed that, when the flow is induced by a line of sinks located at the axis of symmetry, the jet width and the jet velocity verify the equation rjet(τ,La)vjet(τ,La) ∝q∞(2) 073606-7
GORDILLO AND BLANCO–RODRÍGUEZ with the value of q∞fixed during the initial cavity collapse process, namely before the jet is ejected [52]. Indeed, the capillary suction at the liquid side of the cavity induces a velocity V directed towards the axis of symmetry resulting from the balance ρV2∼σ/Rb⇒V∝Vcand, hence, the value of the far-field radial flow rate per unit length driving the collapse of the cavity is Q∞∝VcRb→q∞≈const. It was explained in Refs. [12,15] that, once the radial inflow is set during the initial stage of the collapse, the value of q∞remains approximately constant in time for τ>0 because liquid inertia and mass conservation prevent appreciable changes in this quantity during the very short instants of time following the emergence of the jet. Our prediction for the velocity field in terms of a line of sinks with a constant value of q∞agrees very well with the numerical velocity field calculated at the base of the jet, as shown in Ref. [12] and, moreover, our description is further supported by (i) the result depicted in the inset of Fig. 4(a), where it is shown that the flow rate per unit length q=−rvrremains approximately constant in time for different values of r/rjet(τ) and z/zjet(τ), with vrdenoting the radial velocity, and by (ii) the results shown in Fig. 4(b), where it is seen that (a) the streamlines at the side of the cavity walls are horizontal and (b) the velocity does not appreciably changes with time for a fixed value of r. Before proceeding to deduce in the next section the equations for the time-dependent functions rjet(τ), zjet(τ), and vjet(τ) up to prefactors, in the following we briefly explain the reason behind the result in Fig. 4(a), which shows that Wel∝τ−1/2for τ3×10−3. In order to calculate both rjet(τ,La) and vjet(τ,La), we need to add to Eq. (2) a relationship for zjet(τ), defined in Fig. 1, which we combine with the equation vjet =2dzjet/dτ: Indeed, since Wel1, the liquid pressure remains approximately constant at the portion of the interface surrounding the base of the jet. Hence, by virtue of the Euler-Bernoulli equation expressed in a frame of reference moving vertically with a velocity dzjet/dτ, the modulus of the velocity remains constant and equal to dzjet/dτalong the interface and, since the direction of the velocity is reversed in the moving frame of reference at the jet base we conclude that, in the laboratory frame of reference, vjet(τ)=2dzjet/dτ; see also Sec. III. The remaining equation for zjet can be understood in simple terms taking into account that zjet(τ) is located, for each value of τ, at the vertical position where the interface r=rs(z,τ) meets the axis of symmetry, r=0. Indeed, the integration of the kinematic boundary condition yields to drs dτ∝−q∞ rs⇒r2 s(z,τ)−(ztan β)2∝−2q∞τ, (3) where we have taken into account that vr≃−q∞/rand that rs(z,τ =0) =ztan β; see Fig. 1. Then zjet(τ) is deduced making rs(zjet(τ),τ)=0inEq.(3), from which it follows that zjet ∝ √q∞τ/tan β⇒vjet(τ)∝dzjet/dτ∝√q∞/τ/ tan β. Hence, by virtue of the mass balance, Eq. (2), rjet(τ)∝tan β√q∞τand, consequently, Wel=rjetv2 jet ∝τ−1/2, explaining the result depicted in Fig. 4(a). The scalings in the paragraph above are properly quantified in the next section. III. A THEORY ON THE EJECTION OF JETS FROM COLLAPSING CAVITIES The results in Sec. II reveal that the bursting of a bubble with a radius Rb√σ/(ρg), with gindicating the gravitational acceleration in a liquid of density ρ, viscosity μ, and interfacial tension coefficient σ, can be described in terms of a two-stage process characterized by an initial acceleration stage where the pressure difference p∝σ/Rbinduces a mostly radial velocity field characterized by a value of the flow rate per unit length given by Q∞∝VcRbwith ρV2 c∼p⇒Q∞∝Rbp ρ1/2 .(4) The flow rate Q∞in Eq. (4) remains mostly constant in time along the second stage, during which the jet is issued, because liquid inertia and mass conservation prevent appreciable changes of this quantity during the very short timescale characterizing the jet ejection process [12,15]. Clearly, 073606-8
THEORY OF THE JETS EJECTED AFTER THE INERTIAL … the previous description of the flow can also be applied to other similar physical situations where the initial acceleration stage is induced by either the hydrostatic pressure difference p∼ρgRb and, hence, Q∞∝Rb√gRb[19] or by the jump between the liquid and vapor pressures, namely p∼pa−pv, from which it follows that Q∞∝Rb[(pa−pv)/ρ]1/2[2,28]. Then, in order to describe the ejection of jets produced by the implosion of cavities triggered by the radial flow of a liquid which has been previously accelerated towards the axis of symmetry by the pressure difference p, here we will analyze the collapse of an axisymmetric bubble or cavity like the ones sketched in Figs. 6(a) and 6(b) under the incompressible approach. The flow rate per unit length Q∞which is established during the initial acceleration process, remains constant in time and acts as the far-field boundary condition driving the ejection of the jet because it forces the inward motion of the cavity walls and, therefore, the jet is issued as a mere consequence of mass conservation. We then consider the axisymmetric bubble or cavity sketched in Figs. 6(a) and 6(b) of characteristic length Lcfilled with a gas of density ρgand viscosity μgwhose collapse is driven by a far-field radial flow characterized by the velocity Vc.UsingLc,Lc/Vc, and ρV2 cas the characteristic scales of length, time, and pressure, the physical situation under study here and sketched in Figs. 6(a) and 6(b) can be described in terms of the following dimensionless parameters: We =ρV2 cLc σ,Re =ρVcLc μ,m=μg μ,=ρg ρ.(5) As in Sec. II, dimensionless variables will be written using lowercase letters to differentiate them from their dimensional counterparts, written in capitals, rand zwill respectively indicate the dimensionless radial and axial spatial coordinates in a cylindrical coordinate system, and τ will indicate the dimensionless time after the jet is ejected. In Sec. II,Lc=Rb, with Rbindicating the initial radius of the spherical bubble and Vc=√σ/(ρRb) is the so-called capillary velocity. Therefore, for the case of bubble bursting jets considered in Sec. II,We=1 and Re =Oh−1, with Oh =μ/√ρRbσthe so-called Ohnesorge number, which must be such that Oh 1forajetto be ejected from the bottom of the collapsing bubble [49]. As was also pointed out above, for the case of jets ejected as a consequence of the gravitational collapse of a bubble of length H,Lc=H and Vc=√gH [19], and, for the case of the cavitation of a bubble of radius Rbconsidered in Ref. [28], Lc=Rband Vc=[(pa−pv)/ρ]1/2. The results in Sec. II show the relevance played by the two time-dependent functions zjet(τ) and rjet(τ) characterizing the position of the interface of equation r−rs(z,τ)=0 in the description of the jet dynamics; see Figs. 6(c) and 6(d). Then, it will be our main purpose here to find analytical expressions for rjet(τ), vjet(τ), and zjet(τ)inthe limit τ1. The results in Sec. II indicate that rjet ∝√q∞τand vjet ∝dzjet/dτ∝√q∞/τ when the collapse of a conical cavity of semiangle β(see Figs. 6and 13) is driven by a far-field boundary condition expressing that q∞remains constant in time. Hence, the value of the local Reynolds number characterizing the flow at the base of the jet, Rel=Oh−1vjetrjet ∝Oh−1, will remain constant in time and the value of the local Weber number will diverge as Wel=v2 jetrjet ∝τ−1/21inthe limit τ1 of interest here; see Sec. II. Then, provided that the Reynolds number at the scale of the bubble is large, i.e., provided that Oh−11 for the case considered in Sec. II, both capillary and viscous effects can be neglected, an approximation which is further supported by the numerical results presented in Appendix A. Moreover, since we will not consider the effects of the gas on the jet ejection process, the solutions to be deduced next will not depend on any of the dimensionless parameters defined in Eq. (5); however, as will be shown below, our results will strongly depend on the value of the dimensionless flow rate per unit length q∞and on the initial geometry of the cavity, i.e., on the function rs(z,τ =0). Consequently, since in view of the discussion in the paragraphs above, viscous and capillary effects can be neglected, we consider here the case of an incompressible flow in which vorticity is zero, enabling us to express the velocity field in terms of a velocity potential φ, i.e., v=∇φ, with 073606-9
GORDILLO AND BLANCO–RODRÍGUEZ (a) (b) (c) (d) FIG. 10. Comparison between the numerical results in Appendix Acorresponding to the case of conical cavities with different values of the opening semiangle βwith the predictions in Eqs. (16), (23), and (25). Here vjet(τ)=vz[z=zjet(τ)+αrjet (τ)] with α=1 for the case of slender cavities, namely the conical cavities with β=5◦,β=10◦,andβ=20◦, whereas α=0 for the cases of nonslender conical cavities, namely those with β=30◦and β=40◦. The reason for the different values of αis that the liquid pressure relaxes to that of the gas at a different vertical locations within the jet depending on the value of the opening semiangle; see Fig. 9(c) for the case β=20◦. Here 0(β=5◦)=0.0080, 0(β=10◦)=0.0045, 0(β=20◦)=0.0028, 0(β=30◦)=0.0022, and 0(β=40◦)=0.0012. where we have taken into account that (i) the rate of volume increase in the bulk of the liquid is πr2 jetdzjet/dτ; (ii) by virtue of the Euler-Bernoulli equation, the jet velocity in the moving frame of reference is dzjet/dτ; and (iii) the function q, represented in Fig. 10 is such that −q(¯x→0) → 0.7 and, consequently, zjet+rjet zjet |q|dz ≃|q(¯x=0)|rjet ≈0.7rjet.(24) In order to validate the results in Eqs. (15), (16), (19), and (23) we make use of the results of the incompressible numerical simulations described in Appendix A. Indeed, in Appendix Awe analyze the results of numerical simulations using the GERRIS [20,21]scriptusedinRef.[28] to simulate the emergence of very fast jets, with velocities up to 1000 m s−1, which are emitted when a bubble cavitates in very close proximity to a solid wall once the value of the flow rate per unit length is prescribed as a far-field boundary condition. At this point, notice that our own analysis of the governing equations and of the numerical results reported in Appendix Apermits us to conclude that both the jet shapes and the velocity field converge towards a β-dependent self-similar solution of the 073606-16
THEORY OF THE JETS EJECTED AFTER THE INERTIAL … (a) (b) (c) (d) FIG. 11. Comparison between the numerical results of pmax(τ)−pgas,zjet (τ), vjet(τ)=vz(r=0,z= zjet +rjet ), and rjet(τ) and the respective values predicted by Eqs. (16), (22), and (23). The reason for the small deviation between the predicted and the calculated jet velocity is because the liquid pressure relaxes to that of the gas at a vertical position within the jet which is slightly above z=zjet(τ)+rjet(τ). Here 0(rc=1) =0.015. Euler-Bernoulli and Laplace equations such that lengths and velocities are respectively proportional to √q∞τand √q∞/τ for arbitrary values of βa result which, to our knowledge, had not been reported before. Figure 10 compares the predictions given in Eqs. (16), (20), (21), and (23) with the numerical values of rjet(τ), vjet(τ), pmax(τ), and of rjetdzjet/dτ(τ) obtained from the numerical simulations detailed in Appendix Afor the case of conical cavities with different values of the opening semiangle βtaking q∞=1. The results depicted in Fig. 10 confirm our predictions in the slender limit β1 and permit us to extend our results to generic values of βby just writing zjet(β,τ)=K(β)zjet(β1,τ)=K(β)1.5 tan β√q∞τwith K(β⩽5◦)=1 and (25) K(β=10◦)=1.13, K(β=20◦)=1.33, K(β=30◦)=1.47, and K(β=40◦)=1.63. Moreover, Fig. 11 compares the predictions in Eqs. (16), (22), and (23) with the numerical values of rjet(τ), vjet(τ), pmax(τ), and rjetdzjet/dτ(τ) calculated using the same type of simulations as those detailed in Appendix Abut replacing the equation for the initial shape of the interface by a parabola of equation rs(z,τ =0) =z2/2. The results shown in Fig. 11, which do not include any type of adjustable constant, provide further support to our theory. In order to compare our predictions with the numerical results corresponding to the case of bubble bursting jets, we make use of the equations for rjet and vjet corresponding to the case of 073606-17
GORDILLO AND BLANCO–RODRÍGUEZ conical cavities; see Eqs. (18), (23), and (25): vjet =1.5K(β) tan βq∞(rjet) τ,rjetvjet ≃3.4q∞.(26) with K(β≈40◦)≈1.6; see Fig. 2(b). The results in Figs. 12(a) and 12(b) reveal that the predictions of Eqs. (26)forWe l(τ), rjet(τ), and zjet(τ) closely follow the numerical results and that, in contrast, the analogous scalings deduced using the inertiocapillary ansatz [35,36], rjet ∝τ2/3,vjet ∝τ−1/3, and rjetdzjet/dτ∝τ1/3, clearly deviate from the data. The results in Fig. 12 have been obtained for the slowly varying values of the forcing flow rate per unit length q∞given by the following: If rjet ⩾r∗ jet0 =0.05 : q∞=0.82 If rjet <r∗ jet0 =0.05 : q∞(rjet)=−r0˙r0(0) [r0/(2rc)(0)]×exp ⎛ ⎝−ln r0(0) 2rc(0)2 −ln rjet r0(0)2⎞ ⎠, (27) see Fig. 3in the Supplemental Material [48], where we find that r0(0) =0.2, r0˙r0(0) =−1.5, and r0(0)/[2 rc(0)] =0.16 for the whole range of La <2500 considered here. The reason for the different expressions of the far-field forcing q∞in Eq. (27) can be explained as follows: For rjet >r∗ jet0 ≃0.05, the radius of the jet is sufficiently large for the flow rate per unit length q∞≈1 driving the jet ejection process to be fixed by the radial flow induced during the initial capillary collapse process of the cavity. However, for the cases corresponding to La <2500, rjet(τ1) 1 and, in addition, in these cases a tiny satellite bubble is entrapped beneath the bursting bubble. The satellite bubbles depicted in Fig. 3(see also Fig. 1in the Supplemental Material [48]) are produced as a result of the purely inertial collapse of a cylindrical gas thread and, consequently, when rjet 0.05, the flow rate per unit length q∞<0.82 forcing the ejection of the jet is fixed by the value of the local flow rate per unit length driving the purely inertial pinch-off of the bubble at the instant when the collapsing bubble radius coincides with that of the jet. Hence, the result in Eq. (27) for the case rjet <0.05, which has been deduced making use of the theory in Ref. [17], simply reflects the well-known fact that the flow rate driving the pinch-off of bubbles [18,58–62], see also the Supplemental Material [48] for details, is a decreasing function of the minimum radius of the entrapped bubble. Next we deduce the initial values of the jet width and velocity, namely the minimum value of the jet radius and the maximum value of the jet velocity for a given value of La. For La ⩾2500, the initial jet radius, rjet0, is given by Eq. (1) and, from Eq. (26), rjet0(La >2500) =0.22151−Oh 0.0305,vjet0(La >2500) ≃2.5 rjet0(La >2500).(28) For the cases in which a bubble is entrapped, namely La <2500, a jet will be ejected following the solution of Eqs. (26) and (27) if viscous effects are negligible, namely if the local Reynolds number Rel=Oh−1rjet vjet O(1) ⇒q∞⩾KOh [32], with Ka constant. Therefore, using the result in Eq. (27), the equations for the initial jet width and velocity deduced taking q∞=KOh when a bubble is entrapped, read: rjet0(La <2500) =r0(0) ×exp −1 2ln −KOh r0(0) (2rcr0˙r0)(0) 2 −ln r0(0) 2rc(0)2 vjet0(La <2500) =3.4q∞(rjet0) rjet0 =3.4KOh rjet0(La <2500),(29) where we have made use of the second of the equations in (26). The predictions in Eqs. (29) taking K=10 are pretty close to the numerical values depicted in Figs. 12(c) and 12(d). 073606-18
THEORY OF THE JETS EJECTED AFTER THE INERTIAL … (a) (b) (c) (d) FIG. 12. (a) Wel=rjetv2 jet 1, with Wel∝τ−1/2for τ10−3. The green line indicates the predicted values of Wel[rjet(τ)] =v2 jet rjet ={3.4q∞[rjet(τ)]}2/rjet(τ), with rjet and vjet calculated using Eqs. (26)and(27). (b) The values of vjet(τ), rjetdzjet/dτ,andrjet(τ) calculated using Eqs. (26)and(27), in green, reproduce the numerical results. In contrast, the inertiocapillary predictions [35,36], in pink, vjet ∝τ−1/3,rjetdzjet/dτ∝ τ1/3,andrjet ∝τ2/3, substantially deviate from the data. The numerical values in the figure correspond to τ>τ ∗(La), with τ∗(La) >0 the instant of time such that (ztip −zjet )/rjet ⩾1.25; see Fig. 1and Fig. 2in SM. [(c) and (d)] Comparison of the numerical values of rjet0(La) and vjet0(La) with the predictions (blue continuous line) in Eqs. (28)and(29). The numerical results in Ref. [57]forvjet0 are also included in Fig. 12(d).The differences existing between their results and ours for La <2500 are attributable to the differences existing in the instant of time at which the numerical values are calculated: In our case, the numerical values corresponds to instants of time τ∼10−5. Notice that the diameter and the velocity of the first drop ejected can be calculated using the ballistic equations in Refs. [15,49] using the results in this figure. Notice, however, that, as was pointed out in Ref. [12], for La 2500, the tip of the jet is hardly decelerated by capillary forces and by the viscous drag exerted by the gas, a fact implying that Eqs. (28) can be used to approximate the diameter and the initial velocity of the first drop ejected. 073606-19
GORDILLO AND BLANCO–RODRÍGUEZ As pointed out in Sec. II, the values of rjet0(La) and vjet0(La) deduced above when a bubble is entrapped are not the ones of the initial drop radius and velocity which, however, can be expressed as a function of rjet(τ) and of vjet(τ) following the results in Ref. [49], where the ballistic equations deduced in Ref. [15] are coupled with the mass and the momentum balances at the top drop [63]. As a final remark, notice that both the numerical simulations and the theoretical results have been deduced under the incompressible approximation and under the assumption that both the cavity collapse process and the jet produced, are axisymmetric. In real applications, however, asymmetries are always present but Ref. [64] showed that the amplitude of these asymmetries does not grow during the cavity collapse process. Then, when the minimum radius of the cavity equals the amplitude of the asymmetry, the liquid is no longer accelerated inwards, the cavity breaks and a jet is issued [15]. This fact possesses two important implications for the dynamics of the jet: On the one hand, the ejected jet is not axisymmetric, as it can be clearly appreciated in the experiments reported in Ref. [15] and in the experiments with cavitation bubbles reported in Ref. [28] and, on the other hand, provided that the local Reynolds number at the instant of jet ejection is large, the equation expressing that the minimum radius of the imploding cavity equals the amplitude of the asymmetry, sets the value of the minimum jet radius rjet0 and also the maximum value for the jet velocity, vjet0, through the mass conservation equation derived above, rjet0vjet0 ≃3.4q∞; see Eq. (23). Consequently, the viscous cut-off values of rjet0(La) and vjet0(La) deduced here will describe a real experiment provided that the amplitude of the asymmetries is smaller than rjet0(La). IV. CONCLUSIONS To conclude, here we have analyzed the high-Reynolds number implosion of cavities under the incompressible approach when the collapse is driven by a far-field boundary condition expressing that the flow rate per unit length remains constant in time. We have presented a self-consistent theory which describes the collapse of slender axisymmetric cavities with a generic geometry and provides with algebraic equations, with no adjustable constants, for the time evolution of the speed, width, and of the vertical position of the jets. When particularized to the case of conical bubbles, our algebraic equations reproduce the self-similar results obtained numerically in the limit of small values of the opening semiangle and, when particularized to the case of parabolic cavities, we also find that our predictions for the jet radius, the jet velocity, and the maximum liquid pressure are in good agreement with simulations. We have also shown here that the jets produced right after a bubble bursts at an interface can be described in terms of our purely inertial theory, which differs from the common belief that the dynamics of the jets produced after the inertiocapillary collapse of cavities can be described in terms of an inertiocapillary balance; see, e.g., Ref. [36]. For the case of bubble bursting jets, we have also deduced the equations describing the time-varying jet radius rjet(τ,La) ⩾rjet0(La) and the velocity vjet(τ,La) ⩽vjet0(La), as well as the initial values of the jet radius and of the jet velocity, i.e., rjet0(La) and vjet0(La), respectively. Our predictions agree well with the numerical results, which could reproduce a real experiment provided that the amplitude of the unavoidable asymmetries do not disturb the interface before the jet is ejected. As shown in Ref. [49], the time-dependent jet shape and the time-dependent velocity field within the jet can be expressed in closed form in terms of the time-dependent functions rjet(La,τ), zjet(La,τ), and vjet(La,τ), which constitute the initial conditions of the ballistic equations deduced in Ref. [15]. Hence, here we provide with the necessary information for the correct quantification of the different fluxes transported by bubble bursting jets from the ocean into the atmosphere. In spite of the results presented here resort on the incompressible approximation, they might also be applied to describe the type of very high speed jets generated when a conical bubble implodes very close to a wall, a phenomenon which could play a role in the degradation of the material produced by liquid cavitation. 073606-20
THEORY OF THE JETS EJECTED AFTER THE INERTIAL … ACKNOWLEDGMENTS This research work has been partially supported by the Grant No. PID2020-115655GB-C21 financed by the Spanish MCIN/AEI/10.13039/501100011033. J.M.G. designed the research and the theory, analyzed the data, and wrote the paper, and F.J.B.–R. performed the numerical simulations and analyzed the data. All authors reviewed the results and approved the final version of the paper. APPENDIX A: SELF-SIMILAR CONICAL JETS In this Appendix we present numerical results describing the high-Reynolds-number implosion of conical voids such as those found in the collapsing depressions of standing waves [35], during the collapse of bursting bubbles [29,32], or as a consequence of the collapse of the cavity formed after a drop impacts a liquid pool [9,10,12]. The numerical results have been obtained using the open-source package GERRIS [20,21] and therefore the effects of compressibility on the collapse and jet ejection processes are not retained in the simulations. In fact, the results in this section have been obtained using the GERRIS script provided in Ref. [28] in their study of the dynamic of jets produced in cavitation bubbles near a rigid boundary for small values of the standoff parameter. The numerical setup in Fig. 13(a) illustrates the radial and axial coordinates, Rand Z, respectively; the initial shape of a conical cavity with a half-opening angle β; and Lc, namely the radius of the cylindrical surface where the value of the far-field radial velocity, Vc, is imposed. In the following, ρ,μ, and σwill respectively denote the liquid density, viscosity, and interfacial tension coefficient, whereas ρgand μgwill be used to refer to the gas density and viscosity. The numerical code is solved fixing the value of the flow rate per unit length, LcVc,atR=Lcand imposing outflow boundary conditions at Z=±Lc/2, i.e., the Zderivatives of the velocity and pressure fields are zero at the top and bottom boundaries. The minimum grid size is 2−12Lcand the jet radius at the region where it meets the conical surface [see Fig. 13(b)] will be resolved with at least six numerical cells. Using Lc,Lc/Vc, and ρV2 cas the characteristic scales of length, time, and pressure, the flow sketched in Fig. 13(a) can be described in terms of the following dimensionless parameters: We =ρV2 cLc σ,Re =ρVcLc μ,m=μg μ,=ρg ρ,and β, (A1) but, since here we fix the values of the density and viscosity ratios to those characterizing the airwater system, =1.2×10−3,m=1.8×10−2, the numerical solutions will only depend on We, Re, and β. Notice also that, from now on, dimensionless variables will be written using lowercase letters to differentiate them from their dimensional counterparts, written in capitals, and τand n will indicate, respectively, the dimensionless time after the jet is ejected and the unit normal vector to the interface of equation r=rs(z,τ)⇒F=r−rs(z,τ)=0.(A2) Figure 14 shows that the numerical results are virtually independent of both We and Re for the range of values of these two parameters considered in this study, namely Re ⩾50 and We ⩾0.5. While the results regarding the effect of the Reynolds number could have been anticipated in view of the fact that the minimum value of Re is already large, Re =50, the independence of the solution with We indicates that the mechanism of jet ejection is not driven by capillarity. The influence of the opening semiangle βon the temporal evolution of the jet shapes is analyzed in Fig. 15, where it can be clearly appreciated that the jet velocity increases when βdecreases. It is our purpose now to explore whether there exists a self-similar solution for the high-Reynolds number jets depicted in Fig. 15 for instants of time such that τ1. The analysis starts by considering that the vorticity production at the interface can be neglected in the description of the jet dynamics in the limit Re 1 and then the velocity field vcan be expressed in terms of the velocity 073606-21
GORDILLO AND BLANCO–RODRÍGUEZ FIG. 13. (a) Sketch of the numerical setup indicating the far-field boundary condition, the initial shape of the conical interface, the material properties of the liquid and the gas, as well as the radial and axial coordinates, Rand Z, respectively. (b) Numerical results corresponding to Re =100, We =1, and β=40◦. [(c) and (d)] Numerical values of rvrcorresponding to the simulations depicted in (b) at two different instants of time: τ=4×10−4(c) and τ=4×10−3(d). Notice that the type of numerical simulations carried out using the boundary conditions sketched in (a) impose, as a far-field boundary condition, a constant value of the flow rate per unit length, rvr→−1. potential φas v=∇φ. (A3) Therefore, in the incompressible limit of interest here, continuity demands that ∇·v=0⇒∇2φ=0,(A4) where we have made use of Eq. (A3). The Laplace equation (A4) must be solved subjected to the value of the velocity potential calculated particularizing the Euler-Bernoulli equation ∂φ ∂τ +|∇φ|2 2−We−1∇·n=C(τ)(A5) at the interface of equation F=0, with Fthe function defined in (A2) satisfying the kinematic boundary condition, DF Dt =0⇒∂φ ∂r=∂rs ∂τ +∂φ ∂z ∂rs ∂z,with rs(z→∞,τ)→ztan β(A6) 073606-22
THEORY OF THE JETS EJECTED AFTER THE INERTIAL … (a) (b) FIG. 14. Analysis of the influence of Re and We on the jet ejection process for a fixed value of β=30◦ and two different instants of time τ=4×10−6(a) and τ=4×10−5(b). The results in the figure have been obtained for the following values of the control parameters: Re =100 and We =1 (blue continuous line); Re =50 and We =1(reddashedline);Re=200 and We =1 (red continuous line); Re =100 and We =5 (black continuous line); and Re =100 and We =0.5 (black dashed line). Notice that the differences between the jet shapes are practically indistinguishable and are only appreciable in the zoomed images corresponding to the jet tip, a fact indicating that the spatiotemporal evolution of the jet is virtually independent of Re and We. and with D/Dt indicating the material derivative. In Ref. [35], the self-similar structure of the flow field was deduced, introducing the scaled variables φ=(τ−τ0)δhz−0 (τ−τ0)ε,r (τ−τ0)εand rs(z,τ)=(τ−τ0)εgz−0 (τ−τ0)ε(A7) (a) (b) (c) FIG. 15. This figure analyzes the influence of the opening semiangle on the spatiotemporal evolution of the jet: (a) β=20◦,(b)β=30◦,and(c)β=40◦. In the three cases considered, Re =100 and We =1. Notice that the jet speed increases when βdecreases. 073606-23
GORDILLO AND BLANCO–RODRÍGUEZ into the system of Eqs. (A4)–(A6), finding that both Eqs. (A4) and (A6) can be written in terms of the self-similar variables defined in (A7)ifδ=2ε−1, with this condition also ensuring that the terms ∂φ/∂τ and |∇φ|2/2 in the Euler-Bernoulli equation (A5) possess an identical functional dependence with time. In Ref. [35], the value of the exponent ε=2/3 characterizing the length scale in their inertiocapillary self-similar solutions was deduced imposing the additional condition that the time dependence of the three terms at the left-hand side of Eq. (A5) must be the same. However, the results depicted in Fig. 14 reveal that the time evolution of the jet is independent of the Weber number, a fact indicating that the type of jets considered here are not forced by capillarity. Consequently, in the present case, the only restriction imposed by Eqs. (A4)–(A6) for self-similar solutions to exist is that δ=2ε−1 and then the value of the exponent εin Eq. (A7) should be fixed by the far-field boundary condition of the Laplace equation (A4) which, in our case, reads; see Fig. 13: r∂φ ∂r→−1forr→∞.(A8) The far-field boundary condition (A8) differs from the one corresponding to the inertiocapillary collapse of a cavity, which happens to be analogous to that driving the self-similar wave created after the impact of a disk over a liquid pool. Indeed, in this latter case, Refs. [65,66] reported that the far-field boundary condition for the velocity potential is the well-known two-dimensional irrotational flow around a wedge of angle 2π, φ(r→∞)∝r1/2w(θ),(A9) with rand θindicating the polar coordinates. The different boundary conditions (A8) and (A9) explain the differences between the exponent ε=2/3 characterizing the self-similar solutions reported in Refs. [35,66] and the corresponding value of εwhich is imposed by Eq. (A8). Indeed, here ε=1/2 because, as pointed out above, δ=2ε−1 and, in addition, Eq. (A8) can only be written in terms of the self-similar variables defined in (A7)ifδ=0. Since ε=1/2, in the case the numerical solution converged to a self-similar solution of the system of equations and boundary conditions for sufficiently large values of τ>τ 0, the characteristic length scale of the jets depicted in Fig. 14 should depend on time as (τ−τ0)1/2, whereas velocities should depend on time as (τ−τ0)−1/2, with these predictions also implying that, in the limit (τ−τ0)1 of interest here, both Cand the capillary term in Eq. (A5) are subdominant with respect to the inertial terms ∂φ/∂τ and |∇φ|2/2. In view of the previous analysis, here explore whether the different jet shapes depicted in Fig. 15 can be collapsed onto a single curve. For this purpose, we define the time-dependent functions zjet(τ,β) and rjet(τ,β) illustrated in Fig. 6. Figure 16 shows that the time-dependent jet shapes depicted in Fig. 15 do indeed collapse for over two decades in time onto a time-independent and almost β-independent function when plotted in terms of the scaled radial and vertical coordinates suggested by our previous analysis and defined as χ=r rjet(τ,β),and η=z−0(β) zjet(τ,β)−0(β),(A10) with the values of 0(β) given in Fig. 10, where it is depicted that, indeed, zjet(τ,β)−0(β)∝τ1/2, dz jet/dτ∝τ−1/2, and rjet(τ,β)∝τ1/2. The numerical results corresponding the case of a parabolic cavity shown in Fig. 11, have been obtained using the same numerical code and the same boundary conditions as those detailed in this Appendix, being the only difference that the equation for the initial shape of the interface is, in this case, rs(z,τ =0) =z2/2. 073606-24
THEORY OF THE JETS EJECTED AFTER THE INERTIAL … FIG. 16. The spatiotemporal evolution of the jet shapes of the type shown in Fig. 15 collapse, along two decades in time and for four different values of β, namely β=10◦,20◦,30◦, and 40◦, onto the τ-independent and almost β-independent function depicted in the figure, with χand ηdefined in Eq. (A10). This function for the self-similar shapes of the jets is compared in Fig. 5with the numerical results obtained for the case of bubble bursting jets emerging from the base of a truncated conical surface. APPENDIX B: SOLUTION OF THE INTEGRAL EQUATION FOR q Equation (11) is solved integrating by parts, once we notice that the integrands in Eq. (11) can be written, for ¯z>¯z0,as: df/d¯z0 (¯z−¯z0)2=d d¯z0f (¯z−¯z0)2−2f (¯z−¯z0)3(B1) and for ¯z<¯z0as: −df/d¯z0 (¯z−¯z0)2=− d d¯z0f (¯z−¯z0)2−2f (¯z−¯z0)3.(B2) Notice that lim γ→0¯z−γ 1 d d¯z0f (¯z−¯z0)2d¯z0=f(¯z) γ2,(B3) where we have taken into account that f(¯z=1) =0. Moreover, notice that lim γ→0∞ ¯z+γ−d d¯z0f (¯z−¯z0)2d¯z0=f(¯z) γ2.(B4) We now discretize f, which is assumed to be constant and of value f(¯zi)=fialong Npanels of constant width h, centered at ¯zi=1+3h/4+(i−1)h. Then, for i>j: ¯zj+h/2 ¯zj−h/2−2fj (¯zi−¯z0)3d¯z0=−fj1 (¯zi−¯zj−h/2)2−1 (¯zi−¯zj+h/2)2,(B5) 073606-25