- Email: [email protected]
Lattice Boltzmann simulation of a drop impact on a moving wall with a liquid film
Computers and Mathematics with Applications (
)
–
Contents lists available at SciVerse ScienceDirect
Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa
Lattice Boltzmann simulation of a drop impact on a moving wall with a liquid film Cheng Ming ∗ , Lou Jing Institute of High Performance Computing, Agency for Science, Technology and Research (A*STAR), 1 Fusionopolis Way, #16-16 Connexis, Singapore 138632, Singapore
article
info
Keywords: Drop splashing Moving wall Two-phase flow Lattice Boltzmann method
abstract A liquid drop impact on a moving wall with a pre-existing thin film of the same liquid is simulated by using a two-phase flow lattice Boltzmann model. A perturbation to flow near the wall is introduced shortly after the moment of impact to reflect the instability of drop impact dynamics. The present simulations not only confirm the splash behavior of a drop impact on a stationary wall with the liquid film, but also reveal some interesting dynamics features of the drop impact on a moving wall with the liquid film. It is found that the moving wall can act to enhance or suppress the splash. The influence of a moving wall on the drop splashing is studied by examining the interface shapes, velocity and vorticity fields. The numerical results indicate that a critical threshold of wall velocity exists, which is determined by the Reynolds number, Weber number and the film thickness. If the wall velocity were larger than this critical value, only a part of the drop could splash on the wall. © 2013 Elsevier Ltd. All rights reserved.
1. Introduction Drop impact on a solid or liquid surface is an everyday occurrence. It comprises a rich variety of fluid mechanics facets, such as multiphase and interfacial flow. Drop–surface interactions are also of great importance in many industrial applications. These include spray cooling, welding in material processing, and fuel injection in internal combustion engines, ink-jet printing, soil erosion by rain, etc. The physical phenomena accompanying drop impact on a surface have fascinated researchers for a long time because of their intrinsic complexity and practical importance. Since the pioneering work of Worthington (1908) [1] who was generally considered as the first to investigate the drop impacts systematically, a lot of efforts have been made to understand the physics of drop–wall interactions and to comprehend the underlying mechanism governing flow development at different length scales. Many experimentalists [2–9] had shown that the outcome of droplet impact on thin fluid layers and dry surfaces can be dramatically altered, depending both on the physicochemical characteristics of the drop and on those of the surface itself. When a liquid drop hits a solid surface, it forms a rim and spreads radially until it reaches a certain maximum diameter. If the liquid does not wet the substrate, it then retracts again. During this phase, the drop can show different types of behavior, such as deposition, partial rebound, or splashing. When the liquid drop impacts a thin liquid layer at higher impact velocities, the collision also causes splashing. The drop can form the shape of crowns consisting of a thin liquid sheet with an unstable free rim at the top, from which numerous secondary droplets are ejected. Comprehensive reviews, covering most aspects of this subject in great detail, have been given by Peregrine (1981), Rein (1993), Prosperetti (1993) and Yarin (2006) [10–13]. Although some numerical simulations of the drop splashing on a thin liquid film have been conducted in the past [5,14–23], most of them are restricted to either the two-dimensional (2D) or axisymmetric flow on a stationary wall.
∗
Corresponding author. Tel.: +65 6419 1276; fax: +65 6467 4350. E-mail address: [email protected] (C. Ming).
0898-1221/$ – see front matter © 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.camwa.2013.07.003
2
C. Ming, L. Jing / Computers and Mathematics with Applications (
)
–
The limitation of the earlier numerical simulations could be due to the lack of efficient numerical method for this kind of multiphase flow. In fact, the multiphase flow modeling is one of the most challenging fields of research in applied mathematics and computational fluid dynamics. In the present work, a lattice Boltzmann method (LBM), which was previously presented by Lee and Lin (2005) [24] for incompressible two-phase flow at high density ratio, is used to simulate the three-dimensional (3D) drop splashing on a moving wall with a pre-existing liquid film. Unlike traditional numerical methods, which are based on the solution of the macroscopic variables such as velocity and density, the LBM is based on the kinetic theory for a single particle velocity distribution function [25,26]. The macroscopic variables can be obtained from the moment integrations of the particle distribution function. As such, the LBM has some computational advantages. For example, the lattice Boltzmann equation is linear, as opposed to the nonlinear Navier–Stokes equations. It is well known that the interface boundary between the two phases is mesoscopic by nature. In drop splashing flow, liquid structures of very different length scales develop, while the inertia and surface tension forces strongly act on the liquid at the same time. Hence, it is necessary to consider the effect of microscopic factors on large scale splashing behavior. In this regard, the LBM based on a mesoscopic kinetic equation is a natural multi-scale approach. The objective of the present work is to develop the LBM as an alternative to the experimental approach to investigate the effect of a moving wall on drop splashing. The paper is organized as follows. In Section 2, we provide a brief description of the lattice Boltzmann method for multiphase flows [24] and relevant initial and boundary conditions. We validate our code by comparing the LBM results with the existing experiment results for a drop impact on a thin film. In Section 3, we present some results on the impact of a drop on a moving wall with a thin liquid film. We illustrate the evolution of three-dimensional flow patterns that are observed for different wall velocities. We examine the velocity and vorticity dependence on the moving wall velocity. Conclusions are presented in Section 4. 2. Numerical method As mentioned earlier, the computational approach adopted in this study is based on the model proposed by Lee and Lin (2005) [24]. Since this approach has been discussed previously in Refs. [22,24], only the essential features, such as evolution equations, the initial and boundary conditions for the flow and the validation of the model are described in this section. 2.1. Lattice Boltzmann equation The basic idea of the model [24] is to use two particle distribution functions: one for the order parameter which tracks the interface between two different phases and other for the flow field of the two fluids. Since the primary roles of the two distribution functions differ, the stress and potential forms of the surface tension force are respectively adopted in accordance with the purpose of the distribution function. A stable discretization scheme is used for the intermolecular forcing terms. Based on the discrete Boltzmann equation, the evolution equations for the order parameter distribution function (fα ) and pressure distribution function (gα ) can be written in the following form: fα − fαeq Λf α |t + Λf α |t +δt fα − fαeq − + δt , fα (x + eα δ t , t + δ t ) − fα (x, t ) = − 2τ 2τ 2 (x,t ) (x+eα δ t ,t +δ t )
(1)
gα − gαeq gα − gαeq Λg α |t + Λg α |t +δt gα (x + eα δ t , t + δ t ) − gα (x, t ) = − − + δt , 2τ 2 τ 2 (x,t ) (x+eα δ t ,t +δ t )
(2)
and
Λf α = Λg α =
(eα − u) · [∇(ρ cs2 ) − ρ∇(ϕ − κ∇ 2 ρ)] cs2
(eα − u) · ∇(ρ cs2 ) cs2
gα = fα +
Γα (u) =
p cs2
fαeq
ρ
[Γα (u) − Γα (0)] +
Γα (u),
(3)
(eα − u) · (∇|∇ρ|2 − ∇ · ∇ρ∇ρ) cs2
[Γα (u)),
(4)
− ρ Γα (0),
(5)
,
(6)
where α is the discrete particle velocity direction, eα is the particle velocity, fαeq and gαeq are the equilibrium distribution functions, τ is the relaxation parameter, cs is the lattice speed of sound, ϕ is the chemical potential, and κ is the surface tension parameter. In this paper, we use a nineteen-velocity model in three dimensions (D3Q19). The discrete velocities of the model {eα } are
(0, 0, 0)c , eα = (±1, 0, 0)c , (0, ±1, 0)c , (0, 0, ±1)c , (±1, ±1, 0)c , (±1, 0, ±1)c , (0, ±1, ±1)c ,
α = 0, α = 1, . . . , 6, α = 7, . . . , 18,
(7)
C. Ming, L. Jing / Computers and Mathematics with Applications (
)
–
3
where c = δx /δt . For the athermal D3O19 model, the equilibrium distribution functions are
eα · u
fαeq = wα ρ 1 +
g α = wα eq
p cs2
cs2
+
+
ρ eα · u cs2
(eα · u)2
+
2cs4
−
ρ(eα · u)2 2cs4
u·u 2cs2
−
,
ρu · u
(8)
2cs2
,
(9)
where cs = 1/3 and wα are the corresponding integral weights: 1/3, wα = 1/18, 1/36,
α = 0, α = 1, . . . , 6, α = 7, . . . , 18.
(10)
The chemical potential ϕ is given by
ϕ(ρ) = 4β(ρ − ρ1 )(ρ − ρ2 )(ρ − 0.5(ρ1 + ρ2 )),
(11)
where ρ1 and ρ2 are the densities of the two fluids, respectively. The parameters β and κ depend on the surface tension (σ ) and the interface width (W ), as follows:
β= κ=
12σ
,
W (ρ2 − ρ1 )4
(12)
1.5W σ . (ρ2 − ρ1 )2
(13)
The density of fluid, ρ , the pressure, p, and the velocity, u, are calculated using
ρ=
18
fα ,
(14)
α=0
p=
18
gα cs2 +
α=0
ρu =
18
α=0
1 2
gα eα +
u · ∇(ρ cs2 )δ t , 1 2
(κ∇|∇ρ|2 − κ∇ · ∇ρ∇ρ)δ t .
(15)
(16)
The kinematic viscosity relates to the relaxation parameter by ν = τ cs2 δ t and is calculated by linear interpolation as
ν = ν1 +
ρ − ρ1 (ν2 − ν1 ), ρ2 − ρ1
(17)
where ν1 and ν2 are the kinematic viscosity of the two fluids, respectively. The corresponding macroscopic order parameter equation for Eq. (1) and macroscopic dynamical equations for Eq. (2) are
∂ρ + u · ∇ρ = ∇ · [τ (ρ∇ϕ − ∇ p)δ t ], ∂t 1 ∂p + ∇ · u = 0, ρ cs2 ∂ t ρ
∂u + (u · ∇)u = −∇ p + ∇ · Πν + ∇ · Πσ , ∂t
(18) (19) (20)
where
∇ · Πν = ρν(∇ u + ∇ uT ),
(21)
∇ · Πσ = κ(|∇ρ| I − ∇ρ∇ρ).
(22)
2
For low-Mach number flow, the time derivative of pressure in Eq. (19) is small and the incompressible condition is approximately satisfied.
4
C. Ming, L. Jing / Computers and Mathematics with Applications (
)
–
Fig. 1. Schematic sketch of a droplet impact on a moving plane wall with a thin liquid film.
Eqs. (1) and (2) are solved using the following three steps. Pre-streaming collision step:
eq δt ¯fα (x, t ) = fα (x, t ) − fα − fα + Λf α , 2τ (x,t ) 2 (x,t ) eq δt gα − gα . + g¯α (x, t ) = gα (x, t ) − Λ gα 2τ 2 (x,t ) (x,t )
(23)
(24)
Streaming step: fα (x + eα δ t , t + δ t ) = f¯α (x, t ),
(25)
gα (x + eα δ t , t + δ t ) = g¯α (x, t ).
(26)
Post-streaming collision step:
2τ + Λ , f α 2τ + 1 (x+eα δ t ,t +δ t ) 2τ + 1 (x+eα δ t ,t +δ t ) gα − gαeq 2τ g¯α (x + eα δ t , t + δ t ) = gα (x + eα δ t , t + δ t ) − + Λg α . 2τ + 1 (x+eα δ t ,t +δ t ) 2τ + 1 (x+eα δ t ,t +δ t ) f¯α (x + eα δ t , t + δ t ) = fα (x + eα δ t , t + δ t ) −
fα − fαeq
(27)
(28)
2.2. Flow configuration, boundary and initial conditions We consider that a liquid drop of diameter d0 , density ρ2 and dynamic viscosity µ2 normally impacts onto a wall with a thin liquid film of the same liquid, at an impact velocity u0 . The thickness of the film is h0 . The surrounding gas has density ρ1 and viscosity µ1 . The liquid–gas surface tension is σ . The drop is initially placed at a distance z0 = 0.5d0 + h0 + W from the wall; here, the interface width W = 5δx [24]. A schematic of a drop impingement onto a liquid film is shown in Fig. 1. We will limit ourselves to the flow stage of the crown spike formation. At that stage of drop splashing, gravity effects are typically not important [13,23]; therefore, the gravity is ignored in the present study. The computational domain is lx × ly × lz , where lx and ly are the lateral dimensions of the domain in the horizontal direction, and lz is the dimension in the vertical direction. The no-slip boundary condition is used at z = 0, while pressure boundary conditions are applied at z = lz . The periodic boundary conditions are used at the rest of the boundaries. The no-slip conditions in the LBM are realized with the bounce-back scheme in which all particles colliding with the moving wall not only reverse their momenta, but also gain the momentum imposed by the wall, gα¯ = gα − 6wα ρ2
eα · uw , c2
(29)
where gα¯ is the distribution function of the velocity eα¯ := −eα . A fixed pressure is imposed through the equilibria at z = lz . The initial velocity field is assigned by
−u0 k, u = uw i, 0,
(ρ = ρ2 , z > h0 ), (ρ = ρ2 , z ≤ h0 ), (ρ = ρ1 ),
(30)
and the initial condition of fα and gα can be set to its equilibrium value corresponding to the given velocity field: fα = fα(eq) (u) gα =
gα(eq) (u)
(t = 0),
(31)
(t = 0).
(32)
C. Ming, L. Jing / Computers and Mathematics with Applications (
)
–
5
Fig. 2. Comparison of time evolution of (a) the crown radius, and (b) the crown height at different mesh resolutions for Re = 2000, We = 300 and H = 0.1.
The dimensionless variables and parameters, position (R), well velocity (Uw ), time (T ), film thickness (H), density ratio (δ ), viscosity ratio (λ), Weber number (We) and Reynolds number (Re), are defined according to R=
r d0
,
Uw =
|uw | u0
,
T =
u0 t − W d0
,
H =
h0 d0
,
ρ2 ρ2 d0 u20 ρ2 d0 u0 µ2 , δ= , We = , Re = , µ1 ρ1 σ µ2 where r = xi + yj + zk is a spatial vector in the Cartesian coordinate system. λ=
(33)
2.3. Validation Numerical simulations of the 3D splash of a drop impact on a thin film liquid wall require massive computational resources. In this study, an LBM code with Message Passing Interface (MPI) is developed. In what follows, we will validate our LBM-MPI code by comparison with the existing results for a drop impact on a thin liquid film [5,27]. The method used in this paper has been previously verified and applied to 2D drop impact dynamics [22,24]. The influence of W on the results was studied by Lee and Lin [24]. The numerical result does not strongly depend on the interface thickness W in this particular case when W > 4. We find that the computational domain size of lx × ly × lz = 6d0 × 6d0 × 3d0 is sufficiently large for the present calculation (T ≤ 2). Therefore, we will use the computational domain size of lx × ly × lz = 6d0 × 6d0 × 3d0 unless otherwise stated. We conduct a grid convergence study using four different grid densities defined by the initial diameter of the drop, d0 = 120δx , 160δx , 200δx and 240δx , for a drop impact on a thin film (Uw = 0, H = 0.1, Re = 2000 and We = 300). In Fig. 2, we show a comparison of the time evolution of the crown radius Rc and the crown height Hc with different grid sizes. As the size increases, both the radius and height increase. When the size increases from d0 = 120δx to d0 = 160δx , the changes of the height with time are more obvious. When the size increases further from d0 = 160δx to d0 = 200δx , the values of Rc and Hc vary by less than 3% and 5%, respectively. It is found that the shapes obtained with d0 = 200δx and 240δx resolutions are slightly different but with the same flow features. However, it takes about 3 times the CPU time for the latter. Therefore, the size d0 = 200δx will be adopted for all subsequent calculations. We also check the LBM code performance using Np = 2i processors (i = 0, 1, . . . , 8). The execution time decreases with increasing number of processors. The speedup is linear up to Np = 128. Then the speedup performance slightly drops below the ideal linear curve when Np = 256. All simulations in this paper are carried out on the SGI Altix UV 1000 machine. A typical simulation on a 1200 × 1200 × 600 grid for a run of 104 time steps takes about 20 CPU hours using 256 processors. The present code performance in terms of million node updates per second (MNUPS) is about 1200 × 1200 × 600 × 104 /(20 × 60 × 60) = 120 MNUPS on 256 processors, or an equivalent performance of 0.486 MNUPS per processor. Here, it should be mentioned that if one is interested in the further evolution behavior of droplets beyond the present limitation (T ≤ 2), a relatively larger computational domain and a finer grid system will be necessary. To further validate our code, we simulate the case of drop splashing on a thin liquid film (Uw = 0, We = 300, Re = 2000 and H = 0.15). The flow is similar to, as close as possible, the experiment conducted by Lohse et al. (2008) [27]. In numerical calculations of such a flow, there are no destabilizing effects except for the numerical factors. Since the computational domain and the boundary conditions are symmetric, the calculation result is also symmetric. Fig. 3(a) shows the crown shapes at different times without any perturbation. It is found that as time increases the rim will remain axisymmetric. In real flows, once the drop impacts the surface, the instability takes place immediately following the moment of impact, which grows with time and eventually causes the crown spike formation [7]. Therefore an initial perturbation to flow is required
6
C. Ming, L. Jing / Computers and Mathematics with Applications (
)
–
Fig. 3. Crown shapes for a drop impact onto a stationary thin film (a) without and (b) with perturbation for Re = 2000, We = 300 and H = 0.15 at T = 0.3 and 0.6. The perturbation amplitude is ϵ = 0.02, the wave number n = 20.
to better reflect the dynamics of real drop impact. In our simulation, the perturbation is introduced into impacting velocity in the azimuthal direction by a cosine-wave perturbation for a short period of time. The perturbation is of the form u∗z = uz + ϵ uz exp
h0 − z d0
cos(nθ ).
(34)
The perturbation is specified in terms of its amplitude (ϵ) and wave number (n). We have checked the amplitude effects in the range of 0.01 ≤ ϵ ≤ 0.05. It is found that a higher ϵ does not significantly affect the overall characteristics of the crown spike, although it enhances the rim instability. The wave number determines the number of crown spikes, which can be estimated by [28]
n = 1.14 We .
(35)
Fig. 3(b) shows the crown shapes at different times for the drop impact with an initial perturbation of ϵ = 0.02 and n = 20. The perturbation is introduced into the fluid at T = 0.1 and withdrawn at T = 0.2. At the moment of introduction of the perturbation, only the equilibrium distributions of f and g are modified. We can see that the instability grows with time and induces the crown spike formation. A comparison of our numerical results and the experimental pictures of Lohse et al. (2008) [27] is shown in Fig. 4. Our simulation clearly captures the typical flow features of drop–film collision: a thin liquid sheet jets after the impact; the sheet grows into a crown and propagates radially from the drop; the rim grows at the edge of the crown becomes unstable, and followed by the spikes. The above perturbation parameters are adopted in all subsequent computation. 3. Results and discussion In this section we will study the normal impact of a liquid drop on a moving wall which is covered by a thin film of the same liquid. The ratio of the wall moving velocity to the drop impact velocity varies in the range of 0 ≤ Uw ≤ 1. The influence of the moving wall on the crown formation and evolution will be investigated systematically by varying Uw in the specified range. In order to isolate the effects of the film thickness, surface tension, density and viscosity the simulations are carried out for a set of fixed values: ϵ = 0.02, n = 20, H = 0.15, δ = 100, λ = 50, Re = 2000 and We = 300. 3.1. Flow patterns We firstly observe the flow patterns of drop splashing on the moving wall. Fig. 5 shows the crown formation and evolution of a drop impact on the moving wall at Uw = 0, 0.5 and 1. The azimuthal symmetry of the flow is the most notable feature in the case of Uw = 0, as shown in Fig. 5(a). When the wall moves at a constant velocity along the x direction, asymmetric splashing is observed in Fig. 5(b) and (c), and the splash is enhanced in the direction opposed to the motion and attenuated in the direction of motion. It is found that as Uw increases, the suppression effect becomes more obvious. Comparing Fig. 5(c) with Fig. 5(b), one can observe that when Uw = 1 the crown appears on both upward and downward sides, but the regular spikes do not occur on the downward side, and the crown height also falls. The present results suggest that when the wall velocity is sufficiently large, the splash can be suppressed completely on the downward side, which is similar to the experimental observation of Bird et al. (2009) for drop splashing on a moving dry surface [29]. 3.2. The influence of the moving wall on the ejecta formation We next study the effects of the moving wall on the ejecta which forms at the early stage of drop splashing. Based on the computational flow field data, we plot the interface (ρ = (ρ1 + ρ2 )/2) on the xz plane of symmetry for a given time
C. Ming, L. Jing / Computers and Mathematics with Applications (
)
–
7
Fig. 4. Evolution of the crown for a drop impact onto a thin liquid film at different times. (a) The experimental images of Lohse et al. (2008) [27] vs. (b) the LBM results for Re = 2000, We = 300 and H = 0.15.
Fig. 5. Formation and evolution of the crown for a drop impact onto a moving wall with a liquid film at different wall velocities. (a) Uw = 0, (b) Uw = 0.5, (c) Uw = 1.
T by use of the TECPLOT software, as shown in Fig. 6. Using the probe function of the software, we can obtain the x and z coordinates of all points on the interface curve and find out the local extreme value to determine the crown height (Hc ) by a point-by-point comparison. Fig. 6 shows the time evolution of the interface shapes near the neck of the drop impact at different Uw . For the case of Uw = 0, as shown in Fig. 6(a), once the drop impacts the wall, a layer of displaced liquid beneath it moves radially, and since it fails to accelerate the bulk of surrounding liquid, the inertial force squeezes the liquid layer radially and symmetrically from the neck so that the liquid layer is defected upwards and forms an ejecta. The ejecta sets in the middle of the neck between the oncoming drop and the liquid layer. Very shortly after initial impact, the ejecta experiences a very large acceleration in an approximately impulsive manner as can be seen from Fig. 6(a). When Uw = 0.5,
8
C. Ming, L. Jing / Computers and Mathematics with Applications (
)
–
Fig. 6. Interface shapes near the neck of a drop impacting onto a moving thin film at different times. (a) Uw = 0, (b) Uw = 0.5, (c) Uw = 1.
Fig. 7. Time evolution of the crown height on (a) the upward and (b) downward sides.
the evolutions of interface shapes near the neck are found to be basically the same as those in the case of Uw = 0. Although the ejecta is not uniform, the ejecta occurs on both sides (see Fig. 6(b)). However, when the wall velocity is increased to Uw = 1, as shown in Fig. 6(c), the ejecta only occurs at the upward side. The present results indicate that the moving wall influences the flux ejected into a splash jet on the upward and downward sides. The jetting velocity increases and decreases respectively on the upward and downward sides, resulting in a difference in crown height, as shown in Fig. 6(b) and (c). Fig. 7 plots the time evolution of the crown height for different Uw . We note that the height on the upward and downward sides obviously depends on the wall velocity. The larger the wall velocity, the higher the crown height on the upward side, while the lower the crown height on the downward side.
C. Ming, L. Jing / Computers and Mathematics with Applications (
)
–
9
Fig. 8. The velocity fields on the xz plane of symmetry for a drop impacting onto a moving thin film at different times. (a) Uw = 0 and (b) Uw = 1.
3.3. Velocity field around the ejecta Fig. 8 shows a comparison of the velocity fields on the xz plane of symmetry for the cases of the stationary wall and moving wall. When Uw = 0, at the time of impact and shortly thereafter, the inertial motion of the drop is transformed into the radial motion of liquid away from the axis; the radial movement of the liquid is indicated by the velocity vectors. When T = 0.2, the liquid in the film at radial distances greater than the position where the crown is formed has not started moving, as shown in Fig. 8(a). Initially, this discontinuity may aid in the formation of the crown. Subsequently, as liquid is fed into the crown, the velocity discontinuity persists. The velocity vectors in the crown point away from the wall in the radially outward direction as a result of its growth. As time increases, a characteristic feature of the flow appears near the tip of the crown; a vortex flow is formed near the tip of the crown due to the interaction between ejecta and surrounding fluid. We note that the velocity distribution around the drop is symmetric. We found that the maximum ejecta velocity can be several times higher in magnitude than the impact velocity, in agreement with the early observation of Toroddsen (1998), (2006) [30,31]. When Uw > 0, as shown in Fig. 8(b), at T = 0.2 the velocity discontinuity appears. As time increases, the velocity discontinuity gradually disappears on the downward side. This is because the moving wall causes the increase in relative velocity on the upward side and the decrease on the downward side. This rise and fall of the relative velocity explains the splash enhancement downward and attenuation upward. 3.4. The effect of the moving wall on the evolution of crown We now study the effect of the moving wall on the evolution of crown by examining the vorticity field. Fig. 9 shows the evolution of vorticity with time at Uw = 0 and 1 respectively, which allows us to quantitatively analyze the effect of moving wall on drop splashing. The vorticity distributions around the crown are directly correlated to the velocity field shown in Fig. 8. For the case of Uw = 0, Fig. 9(a) clearly illustrates that there is a pair of counter-rotating vortices around the crown due to the interaction of the ejecta with the surrounding fluid. The one above the ejecta is generated by the flow of fluid from the drop; another below the ejecta is generated by the flow of fluid from the initial layer. The strengths of vorticity are equal, which drive the liquid radially and symmetrically from the neck. The vortex expands radially rather rapidly in the initial stage of the impact during 0.2 ≤ T ≤ 0.4. The radial expansion slows down during T > 0.4. This vortex pair will assist in stripping the spike from the rim and lead to the secondary drop formation. For the cases of the moving wall, when Uw = 0.5 (not shown here), the patterns of vorticity contours’ evolution are found to be basically the same as those in the case of a normal impact, apart from the lack of the perfect rotational symmetry. A pair of counter-rotating vortices appears around the crown but their strengths are not uniform along the neck. When the wall velocity increases to Uw = 1, there is an obvious difference between the pair of counter-rotating vortices on the downward side, as shown in Fig. 9(b). As time increases, the viscous effect diffuses vorticity and significantly weakens the strength of the vorticity generated by the flow of the initial layer. When T > 0.2, the vortex on the downward side is quickly dissipated due to viscous. Thus, one vortex appears on the downward side. It means that the flow does not form the crown spikes on the downward side.
10
C. Ming, L. Jing / Computers and Mathematics with Applications (
)
–
Fig. 9. Evolution of the vorticity around the ejecta of a drop impacting onto a moving thin film at different times. (a) Uw = 0 and (b) Uw = 1.
4. Conclusions In this paper, a numerical study of a liquid drop impact on a moving wall with a pre-existing thin film of the same liquid is presented. A lattice Boltzmann model for incompressible two-phase flow at high density ratio is employed to investigate the dynamic behavior of drop splashing on the wall. In our simulation, an initial perturbation to flow is introduced to reflect the dynamics of drop evolution. The LBM-MPI code is validated by comparing the numerical results with the experimental observations for drop splashing. The qualitative features of the drop splashing on the moving wall at the early stage have been obtained by the present numerical simulation. We study the effects of the moving wall on the formation and evolution of ejecta and crown. It is found that an obvious effect is that the moving wall can act to enhance or suppress the splash. Further study will undoubtedly deepen our understanding about some important phenomena accompanying drop–wall interactions. Acknowledgments The authors gratefully acknowledge Professor Detlef Lohse for helpful discussions and his experimental photographs. They would also like to thank Dr. Z. Shang for his MPI-code development assistance. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17]
A.M. Worthington, A Study of Splashes, Green, London, 1908. H.E. Edgerton, J.R. Killian, Flash! Seeing the Unseen by Ultra-High-Speed Photography, Boston, Branford, 1954. G.E. Cossali, A. Coghe, M. Marengo, The impact of a single drop on a wetted solid surface, Exp. Fluids 22 (1997) 463–472. A.B. Wang, C.C. Chen, Splashing impact of a single drop onto very thin liquid films, Phys. Fluids 12 (2000) 2155–2158. C. Josserand, S. Zaleski, Droplet splashing on a thin liquid film, Phys. Fluids 16 (2003) 1650–1657. R. Rioboo, C. Bauthier, J. Conti, J. De Coninck, M. Voue, Experimental investigation of splash and crown formation during single drop impact on wetted surfaces, Exp. Fluids 35 (2003) 648–652. R.L. Vander-Wal, G.M. Berger, S.D. Mozes, Droplets splashing upon films of the same fluid of various depths, Exp. Fluids 40 (2006) 33–52. B. Ching, M.W. Golay, T.J. Johnson, Droplet impacts upon liquid surfaces, Science 226 (1984) 535–537. R. Krechetnikov, G.M. Homsy, Crown-forming instability phenomena in the drop splash problem, J. Colloid Interface Sci. 331 (2009) 555–559. D. Peregrine, The fascination of fluid mechanics, J. Fluid Mech. 106 (1981) 59–80. M. Rein, Phenomena of liquid drop impact on solid and liquid surfaces, Fluid Dyn. Res. 12 (1993) 61–93. A. Prosperetti, H. Oguz, The impact of drops on liquid surfaces and the underwater noise of rain, Annu. Rev. Fluid Mech. 25 (1993) 577–602. A.L. Yarin, Drop impact dynamics: splashing, spreading, receding, bouncing, Annu. Rev. Fluid Mech. 38 (2006) 159–192. J. Fukai, Z. Zhao, D. Poulikakos, C.M. Megaridis, O. Miyatake, Modeling of the deformation of a liquid droplet impinging upon a flat surface, Phys. Fluids A 5 (1993) 2588–2599. J. Fukai, Y. Shiiba, T. Yamamoto, O. Miyatake, D. Poulikakos, Wetting effects on the spreading of a liquid droplet colliding with a flat surface: experimant and modeling, Phys. Fluids 7 (1995) 236–247. M. Pasandideh-Fard, Y.M. Qiao, S. Chandra, J. Mostaghimi, Capillary effects during droplet impact on a solid surface, Phys. Fluids 8 (1996) 650–659. D.A. Weiss, A.L. Yarin, Single drop impact onto liquid films: neck distortion, jetting, tiny bubble entrainment, and crown formation, J. Fluid Mech. 385 (1999) 229–254.
C. Ming, L. Jing / Computers and Mathematics with Applications ( [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31]
)
–
11
M. Rieber, A. Frohn, A numerical study on the mechanism of splashing, Int. J. Heat Fluid Flow 20 (1999) 455–461. M. Bussmann, S. Chandra, J. Mostaghimi, Modeling the splash of a droplet impacting a solid surface, Phys. Fluids 12 (2000) 3121–3132. I.V. Roisman, C. Tropea, Impact of a drop onto a wetted wall: description of crown formation and propagation, J. Fluid Mech. 472 (2002) 373–397. S.D. Howison, J.R. Ockendon, J.M. Oliver, R. Purvis, F.T. Smith, Droplet impact on a thin fluid layer, J. Fluid Mech. 542 (2005) 1–23. S. Mukherjee, J. Abraham, Crown behavior in drop impact on wet walls, Phys. Fluids 19 (2007) 052103. G. Coppola, G. Rocco, L. de Luca, Insights on the impact of a plane drop on a thin liquid film, Phys. Fluids 23 (2011) 022105. T. Lee, C.L. Lin, A stable discretization of the lattice Boltzmann equation for simulation of incompressible two-phase flows at high density ratio, J. Comput. Phys. 206 (2005) 16–47. S. Chen, G. Doolen, Lattice Boltzmann method for fluid flows, Annu. Rev. Fluid Mech. 30 (1998) 329–364. L.-S. Luo, Theory of the lattice Boltzmann method: lattice Boltzmann models for nonideal gases, Phys. Rev. E 62 (2000) 4982–4996. W. van-Hoeve, T. Segers, H. Droes, D. Lohse, M. Versluis, A splash of red, in: The 61st Annual Meeting of the APS DFD, November 23–25, 2008 in San Antonio, Texas, 2008. N.Z. Mehdizadeh, S. Chandra, J. Mostaghimi, Formation of fingers around the edges of a drop hitting a metal plate with high velocity, J. Fluid Mech. 510 (2004) 353–373. J.C. Bird, S.S.H. Tsai, H.A. Stone, Inclined to splash: triggering and inhibiting a splash with tangential velocity, New J. Phys. 11 (2009) 063017. S.T. Thoroddsen, J. Sakakibara, Evolution of fingering pattern of an impacting drop, Phys. Fluids 10 (1998) 1359–1374. S.T. Thoroddsen, T.G. Etoh, K. Takehara, Crown breakup by marangoni instability, J. Fluid Mech. 557 (2006) 63–72.
)
–
Contents lists available at SciVerse ScienceDirect
Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa
Lattice Boltzmann simulation of a drop impact on a moving wall with a liquid film Cheng Ming ∗ , Lou Jing Institute of High Performance Computing, Agency for Science, Technology and Research (A*STAR), 1 Fusionopolis Way, #16-16 Connexis, Singapore 138632, Singapore
article
info
Keywords: Drop splashing Moving wall Two-phase flow Lattice Boltzmann method
abstract A liquid drop impact on a moving wall with a pre-existing thin film of the same liquid is simulated by using a two-phase flow lattice Boltzmann model. A perturbation to flow near the wall is introduced shortly after the moment of impact to reflect the instability of drop impact dynamics. The present simulations not only confirm the splash behavior of a drop impact on a stationary wall with the liquid film, but also reveal some interesting dynamics features of the drop impact on a moving wall with the liquid film. It is found that the moving wall can act to enhance or suppress the splash. The influence of a moving wall on the drop splashing is studied by examining the interface shapes, velocity and vorticity fields. The numerical results indicate that a critical threshold of wall velocity exists, which is determined by the Reynolds number, Weber number and the film thickness. If the wall velocity were larger than this critical value, only a part of the drop could splash on the wall. © 2013 Elsevier Ltd. All rights reserved.
1. Introduction Drop impact on a solid or liquid surface is an everyday occurrence. It comprises a rich variety of fluid mechanics facets, such as multiphase and interfacial flow. Drop–surface interactions are also of great importance in many industrial applications. These include spray cooling, welding in material processing, and fuel injection in internal combustion engines, ink-jet printing, soil erosion by rain, etc. The physical phenomena accompanying drop impact on a surface have fascinated researchers for a long time because of their intrinsic complexity and practical importance. Since the pioneering work of Worthington (1908) [1] who was generally considered as the first to investigate the drop impacts systematically, a lot of efforts have been made to understand the physics of drop–wall interactions and to comprehend the underlying mechanism governing flow development at different length scales. Many experimentalists [2–9] had shown that the outcome of droplet impact on thin fluid layers and dry surfaces can be dramatically altered, depending both on the physicochemical characteristics of the drop and on those of the surface itself. When a liquid drop hits a solid surface, it forms a rim and spreads radially until it reaches a certain maximum diameter. If the liquid does not wet the substrate, it then retracts again. During this phase, the drop can show different types of behavior, such as deposition, partial rebound, or splashing. When the liquid drop impacts a thin liquid layer at higher impact velocities, the collision also causes splashing. The drop can form the shape of crowns consisting of a thin liquid sheet with an unstable free rim at the top, from which numerous secondary droplets are ejected. Comprehensive reviews, covering most aspects of this subject in great detail, have been given by Peregrine (1981), Rein (1993), Prosperetti (1993) and Yarin (2006) [10–13]. Although some numerical simulations of the drop splashing on a thin liquid film have been conducted in the past [5,14–23], most of them are restricted to either the two-dimensional (2D) or axisymmetric flow on a stationary wall.
∗
Corresponding author. Tel.: +65 6419 1276; fax: +65 6467 4350. E-mail address: [email protected] (C. Ming).
0898-1221/$ – see front matter © 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.camwa.2013.07.003
2
C. Ming, L. Jing / Computers and Mathematics with Applications (
)
–
The limitation of the earlier numerical simulations could be due to the lack of efficient numerical method for this kind of multiphase flow. In fact, the multiphase flow modeling is one of the most challenging fields of research in applied mathematics and computational fluid dynamics. In the present work, a lattice Boltzmann method (LBM), which was previously presented by Lee and Lin (2005) [24] for incompressible two-phase flow at high density ratio, is used to simulate the three-dimensional (3D) drop splashing on a moving wall with a pre-existing liquid film. Unlike traditional numerical methods, which are based on the solution of the macroscopic variables such as velocity and density, the LBM is based on the kinetic theory for a single particle velocity distribution function [25,26]. The macroscopic variables can be obtained from the moment integrations of the particle distribution function. As such, the LBM has some computational advantages. For example, the lattice Boltzmann equation is linear, as opposed to the nonlinear Navier–Stokes equations. It is well known that the interface boundary between the two phases is mesoscopic by nature. In drop splashing flow, liquid structures of very different length scales develop, while the inertia and surface tension forces strongly act on the liquid at the same time. Hence, it is necessary to consider the effect of microscopic factors on large scale splashing behavior. In this regard, the LBM based on a mesoscopic kinetic equation is a natural multi-scale approach. The objective of the present work is to develop the LBM as an alternative to the experimental approach to investigate the effect of a moving wall on drop splashing. The paper is organized as follows. In Section 2, we provide a brief description of the lattice Boltzmann method for multiphase flows [24] and relevant initial and boundary conditions. We validate our code by comparing the LBM results with the existing experiment results for a drop impact on a thin film. In Section 3, we present some results on the impact of a drop on a moving wall with a thin liquid film. We illustrate the evolution of three-dimensional flow patterns that are observed for different wall velocities. We examine the velocity and vorticity dependence on the moving wall velocity. Conclusions are presented in Section 4. 2. Numerical method As mentioned earlier, the computational approach adopted in this study is based on the model proposed by Lee and Lin (2005) [24]. Since this approach has been discussed previously in Refs. [22,24], only the essential features, such as evolution equations, the initial and boundary conditions for the flow and the validation of the model are described in this section. 2.1. Lattice Boltzmann equation The basic idea of the model [24] is to use two particle distribution functions: one for the order parameter which tracks the interface between two different phases and other for the flow field of the two fluids. Since the primary roles of the two distribution functions differ, the stress and potential forms of the surface tension force are respectively adopted in accordance with the purpose of the distribution function. A stable discretization scheme is used for the intermolecular forcing terms. Based on the discrete Boltzmann equation, the evolution equations for the order parameter distribution function (fα ) and pressure distribution function (gα ) can be written in the following form: fα − fαeq Λf α |t + Λf α |t +δt fα − fαeq − + δt , fα (x + eα δ t , t + δ t ) − fα (x, t ) = − 2τ 2τ 2 (x,t ) (x+eα δ t ,t +δ t )
(1)
gα − gαeq gα − gαeq Λg α |t + Λg α |t +δt gα (x + eα δ t , t + δ t ) − gα (x, t ) = − − + δt , 2τ 2 τ 2 (x,t ) (x+eα δ t ,t +δ t )
(2)
and
Λf α = Λg α =
(eα − u) · [∇(ρ cs2 ) − ρ∇(ϕ − κ∇ 2 ρ)] cs2
(eα − u) · ∇(ρ cs2 ) cs2
gα = fα +
Γα (u) =
p cs2
fαeq
ρ
[Γα (u) − Γα (0)] +
Γα (u),
(3)
(eα − u) · (∇|∇ρ|2 − ∇ · ∇ρ∇ρ) cs2
[Γα (u)),
(4)
− ρ Γα (0),
(5)
,
(6)
where α is the discrete particle velocity direction, eα is the particle velocity, fαeq and gαeq are the equilibrium distribution functions, τ is the relaxation parameter, cs is the lattice speed of sound, ϕ is the chemical potential, and κ is the surface tension parameter. In this paper, we use a nineteen-velocity model in three dimensions (D3Q19). The discrete velocities of the model {eα } are
(0, 0, 0)c , eα = (±1, 0, 0)c , (0, ±1, 0)c , (0, 0, ±1)c , (±1, ±1, 0)c , (±1, 0, ±1)c , (0, ±1, ±1)c ,
α = 0, α = 1, . . . , 6, α = 7, . . . , 18,
(7)
C. Ming, L. Jing / Computers and Mathematics with Applications (
)
–
3
where c = δx /δt . For the athermal D3O19 model, the equilibrium distribution functions are
eα · u
fαeq = wα ρ 1 +
g α = wα eq
p cs2
cs2
+
+
ρ eα · u cs2
(eα · u)2
+
2cs4
−
ρ(eα · u)2 2cs4
u·u 2cs2
−
,
ρu · u
(8)
2cs2
,
(9)
where cs = 1/3 and wα are the corresponding integral weights: 1/3, wα = 1/18, 1/36,
α = 0, α = 1, . . . , 6, α = 7, . . . , 18.
(10)
The chemical potential ϕ is given by
ϕ(ρ) = 4β(ρ − ρ1 )(ρ − ρ2 )(ρ − 0.5(ρ1 + ρ2 )),
(11)
where ρ1 and ρ2 are the densities of the two fluids, respectively. The parameters β and κ depend on the surface tension (σ ) and the interface width (W ), as follows:
β= κ=
12σ
,
W (ρ2 − ρ1 )4
(12)
1.5W σ . (ρ2 − ρ1 )2
(13)
The density of fluid, ρ , the pressure, p, and the velocity, u, are calculated using
ρ=
18
fα ,
(14)
α=0
p=
18
gα cs2 +
α=0
ρu =
18
α=0
1 2
gα eα +
u · ∇(ρ cs2 )δ t , 1 2
(κ∇|∇ρ|2 − κ∇ · ∇ρ∇ρ)δ t .
(15)
(16)
The kinematic viscosity relates to the relaxation parameter by ν = τ cs2 δ t and is calculated by linear interpolation as
ν = ν1 +
ρ − ρ1 (ν2 − ν1 ), ρ2 − ρ1
(17)
where ν1 and ν2 are the kinematic viscosity of the two fluids, respectively. The corresponding macroscopic order parameter equation for Eq. (1) and macroscopic dynamical equations for Eq. (2) are
∂ρ + u · ∇ρ = ∇ · [τ (ρ∇ϕ − ∇ p)δ t ], ∂t 1 ∂p + ∇ · u = 0, ρ cs2 ∂ t ρ
∂u + (u · ∇)u = −∇ p + ∇ · Πν + ∇ · Πσ , ∂t
(18) (19) (20)
where
∇ · Πν = ρν(∇ u + ∇ uT ),
(21)
∇ · Πσ = κ(|∇ρ| I − ∇ρ∇ρ).
(22)
2
For low-Mach number flow, the time derivative of pressure in Eq. (19) is small and the incompressible condition is approximately satisfied.
4
C. Ming, L. Jing / Computers and Mathematics with Applications (
)
–
Fig. 1. Schematic sketch of a droplet impact on a moving plane wall with a thin liquid film.
Eqs. (1) and (2) are solved using the following three steps. Pre-streaming collision step:
eq δt ¯fα (x, t ) = fα (x, t ) − fα − fα + Λf α , 2τ (x,t ) 2 (x,t ) eq δt gα − gα . + g¯α (x, t ) = gα (x, t ) − Λ gα 2τ 2 (x,t ) (x,t )
(23)
(24)
Streaming step: fα (x + eα δ t , t + δ t ) = f¯α (x, t ),
(25)
gα (x + eα δ t , t + δ t ) = g¯α (x, t ).
(26)
Post-streaming collision step:
2τ + Λ , f α 2τ + 1 (x+eα δ t ,t +δ t ) 2τ + 1 (x+eα δ t ,t +δ t ) gα − gαeq 2τ g¯α (x + eα δ t , t + δ t ) = gα (x + eα δ t , t + δ t ) − + Λg α . 2τ + 1 (x+eα δ t ,t +δ t ) 2τ + 1 (x+eα δ t ,t +δ t ) f¯α (x + eα δ t , t + δ t ) = fα (x + eα δ t , t + δ t ) −
fα − fαeq
(27)
(28)
2.2. Flow configuration, boundary and initial conditions We consider that a liquid drop of diameter d0 , density ρ2 and dynamic viscosity µ2 normally impacts onto a wall with a thin liquid film of the same liquid, at an impact velocity u0 . The thickness of the film is h0 . The surrounding gas has density ρ1 and viscosity µ1 . The liquid–gas surface tension is σ . The drop is initially placed at a distance z0 = 0.5d0 + h0 + W from the wall; here, the interface width W = 5δx [24]. A schematic of a drop impingement onto a liquid film is shown in Fig. 1. We will limit ourselves to the flow stage of the crown spike formation. At that stage of drop splashing, gravity effects are typically not important [13,23]; therefore, the gravity is ignored in the present study. The computational domain is lx × ly × lz , where lx and ly are the lateral dimensions of the domain in the horizontal direction, and lz is the dimension in the vertical direction. The no-slip boundary condition is used at z = 0, while pressure boundary conditions are applied at z = lz . The periodic boundary conditions are used at the rest of the boundaries. The no-slip conditions in the LBM are realized with the bounce-back scheme in which all particles colliding with the moving wall not only reverse their momenta, but also gain the momentum imposed by the wall, gα¯ = gα − 6wα ρ2
eα · uw , c2
(29)
where gα¯ is the distribution function of the velocity eα¯ := −eα . A fixed pressure is imposed through the equilibria at z = lz . The initial velocity field is assigned by
−u0 k, u = uw i, 0,
(ρ = ρ2 , z > h0 ), (ρ = ρ2 , z ≤ h0 ), (ρ = ρ1 ),
(30)
and the initial condition of fα and gα can be set to its equilibrium value corresponding to the given velocity field: fα = fα(eq) (u) gα =
gα(eq) (u)
(t = 0),
(31)
(t = 0).
(32)
C. Ming, L. Jing / Computers and Mathematics with Applications (
)
–
5
Fig. 2. Comparison of time evolution of (a) the crown radius, and (b) the crown height at different mesh resolutions for Re = 2000, We = 300 and H = 0.1.
The dimensionless variables and parameters, position (R), well velocity (Uw ), time (T ), film thickness (H), density ratio (δ ), viscosity ratio (λ), Weber number (We) and Reynolds number (Re), are defined according to R=
r d0
,
Uw =
|uw | u0
,
T =
u0 t − W d0
,
H =
h0 d0
,
ρ2 ρ2 d0 u20 ρ2 d0 u0 µ2 , δ= , We = , Re = , µ1 ρ1 σ µ2 where r = xi + yj + zk is a spatial vector in the Cartesian coordinate system. λ=
(33)
2.3. Validation Numerical simulations of the 3D splash of a drop impact on a thin film liquid wall require massive computational resources. In this study, an LBM code with Message Passing Interface (MPI) is developed. In what follows, we will validate our LBM-MPI code by comparison with the existing results for a drop impact on a thin liquid film [5,27]. The method used in this paper has been previously verified and applied to 2D drop impact dynamics [22,24]. The influence of W on the results was studied by Lee and Lin [24]. The numerical result does not strongly depend on the interface thickness W in this particular case when W > 4. We find that the computational domain size of lx × ly × lz = 6d0 × 6d0 × 3d0 is sufficiently large for the present calculation (T ≤ 2). Therefore, we will use the computational domain size of lx × ly × lz = 6d0 × 6d0 × 3d0 unless otherwise stated. We conduct a grid convergence study using four different grid densities defined by the initial diameter of the drop, d0 = 120δx , 160δx , 200δx and 240δx , for a drop impact on a thin film (Uw = 0, H = 0.1, Re = 2000 and We = 300). In Fig. 2, we show a comparison of the time evolution of the crown radius Rc and the crown height Hc with different grid sizes. As the size increases, both the radius and height increase. When the size increases from d0 = 120δx to d0 = 160δx , the changes of the height with time are more obvious. When the size increases further from d0 = 160δx to d0 = 200δx , the values of Rc and Hc vary by less than 3% and 5%, respectively. It is found that the shapes obtained with d0 = 200δx and 240δx resolutions are slightly different but with the same flow features. However, it takes about 3 times the CPU time for the latter. Therefore, the size d0 = 200δx will be adopted for all subsequent calculations. We also check the LBM code performance using Np = 2i processors (i = 0, 1, . . . , 8). The execution time decreases with increasing number of processors. The speedup is linear up to Np = 128. Then the speedup performance slightly drops below the ideal linear curve when Np = 256. All simulations in this paper are carried out on the SGI Altix UV 1000 machine. A typical simulation on a 1200 × 1200 × 600 grid for a run of 104 time steps takes about 20 CPU hours using 256 processors. The present code performance in terms of million node updates per second (MNUPS) is about 1200 × 1200 × 600 × 104 /(20 × 60 × 60) = 120 MNUPS on 256 processors, or an equivalent performance of 0.486 MNUPS per processor. Here, it should be mentioned that if one is interested in the further evolution behavior of droplets beyond the present limitation (T ≤ 2), a relatively larger computational domain and a finer grid system will be necessary. To further validate our code, we simulate the case of drop splashing on a thin liquid film (Uw = 0, We = 300, Re = 2000 and H = 0.15). The flow is similar to, as close as possible, the experiment conducted by Lohse et al. (2008) [27]. In numerical calculations of such a flow, there are no destabilizing effects except for the numerical factors. Since the computational domain and the boundary conditions are symmetric, the calculation result is also symmetric. Fig. 3(a) shows the crown shapes at different times without any perturbation. It is found that as time increases the rim will remain axisymmetric. In real flows, once the drop impacts the surface, the instability takes place immediately following the moment of impact, which grows with time and eventually causes the crown spike formation [7]. Therefore an initial perturbation to flow is required
6
C. Ming, L. Jing / Computers and Mathematics with Applications (
)
–
Fig. 3. Crown shapes for a drop impact onto a stationary thin film (a) without and (b) with perturbation for Re = 2000, We = 300 and H = 0.15 at T = 0.3 and 0.6. The perturbation amplitude is ϵ = 0.02, the wave number n = 20.
to better reflect the dynamics of real drop impact. In our simulation, the perturbation is introduced into impacting velocity in the azimuthal direction by a cosine-wave perturbation for a short period of time. The perturbation is of the form u∗z = uz + ϵ uz exp
h0 − z d0
cos(nθ ).
(34)
The perturbation is specified in terms of its amplitude (ϵ) and wave number (n). We have checked the amplitude effects in the range of 0.01 ≤ ϵ ≤ 0.05. It is found that a higher ϵ does not significantly affect the overall characteristics of the crown spike, although it enhances the rim instability. The wave number determines the number of crown spikes, which can be estimated by [28]
n = 1.14 We .
(35)
Fig. 3(b) shows the crown shapes at different times for the drop impact with an initial perturbation of ϵ = 0.02 and n = 20. The perturbation is introduced into the fluid at T = 0.1 and withdrawn at T = 0.2. At the moment of introduction of the perturbation, only the equilibrium distributions of f and g are modified. We can see that the instability grows with time and induces the crown spike formation. A comparison of our numerical results and the experimental pictures of Lohse et al. (2008) [27] is shown in Fig. 4. Our simulation clearly captures the typical flow features of drop–film collision: a thin liquid sheet jets after the impact; the sheet grows into a crown and propagates radially from the drop; the rim grows at the edge of the crown becomes unstable, and followed by the spikes. The above perturbation parameters are adopted in all subsequent computation. 3. Results and discussion In this section we will study the normal impact of a liquid drop on a moving wall which is covered by a thin film of the same liquid. The ratio of the wall moving velocity to the drop impact velocity varies in the range of 0 ≤ Uw ≤ 1. The influence of the moving wall on the crown formation and evolution will be investigated systematically by varying Uw in the specified range. In order to isolate the effects of the film thickness, surface tension, density and viscosity the simulations are carried out for a set of fixed values: ϵ = 0.02, n = 20, H = 0.15, δ = 100, λ = 50, Re = 2000 and We = 300. 3.1. Flow patterns We firstly observe the flow patterns of drop splashing on the moving wall. Fig. 5 shows the crown formation and evolution of a drop impact on the moving wall at Uw = 0, 0.5 and 1. The azimuthal symmetry of the flow is the most notable feature in the case of Uw = 0, as shown in Fig. 5(a). When the wall moves at a constant velocity along the x direction, asymmetric splashing is observed in Fig. 5(b) and (c), and the splash is enhanced in the direction opposed to the motion and attenuated in the direction of motion. It is found that as Uw increases, the suppression effect becomes more obvious. Comparing Fig. 5(c) with Fig. 5(b), one can observe that when Uw = 1 the crown appears on both upward and downward sides, but the regular spikes do not occur on the downward side, and the crown height also falls. The present results suggest that when the wall velocity is sufficiently large, the splash can be suppressed completely on the downward side, which is similar to the experimental observation of Bird et al. (2009) for drop splashing on a moving dry surface [29]. 3.2. The influence of the moving wall on the ejecta formation We next study the effects of the moving wall on the ejecta which forms at the early stage of drop splashing. Based on the computational flow field data, we plot the interface (ρ = (ρ1 + ρ2 )/2) on the xz plane of symmetry for a given time
C. Ming, L. Jing / Computers and Mathematics with Applications (
)
–
7
Fig. 4. Evolution of the crown for a drop impact onto a thin liquid film at different times. (a) The experimental images of Lohse et al. (2008) [27] vs. (b) the LBM results for Re = 2000, We = 300 and H = 0.15.
Fig. 5. Formation and evolution of the crown for a drop impact onto a moving wall with a liquid film at different wall velocities. (a) Uw = 0, (b) Uw = 0.5, (c) Uw = 1.
T by use of the TECPLOT software, as shown in Fig. 6. Using the probe function of the software, we can obtain the x and z coordinates of all points on the interface curve and find out the local extreme value to determine the crown height (Hc ) by a point-by-point comparison. Fig. 6 shows the time evolution of the interface shapes near the neck of the drop impact at different Uw . For the case of Uw = 0, as shown in Fig. 6(a), once the drop impacts the wall, a layer of displaced liquid beneath it moves radially, and since it fails to accelerate the bulk of surrounding liquid, the inertial force squeezes the liquid layer radially and symmetrically from the neck so that the liquid layer is defected upwards and forms an ejecta. The ejecta sets in the middle of the neck between the oncoming drop and the liquid layer. Very shortly after initial impact, the ejecta experiences a very large acceleration in an approximately impulsive manner as can be seen from Fig. 6(a). When Uw = 0.5,
8
C. Ming, L. Jing / Computers and Mathematics with Applications (
)
–
Fig. 6. Interface shapes near the neck of a drop impacting onto a moving thin film at different times. (a) Uw = 0, (b) Uw = 0.5, (c) Uw = 1.
Fig. 7. Time evolution of the crown height on (a) the upward and (b) downward sides.
the evolutions of interface shapes near the neck are found to be basically the same as those in the case of Uw = 0. Although the ejecta is not uniform, the ejecta occurs on both sides (see Fig. 6(b)). However, when the wall velocity is increased to Uw = 1, as shown in Fig. 6(c), the ejecta only occurs at the upward side. The present results indicate that the moving wall influences the flux ejected into a splash jet on the upward and downward sides. The jetting velocity increases and decreases respectively on the upward and downward sides, resulting in a difference in crown height, as shown in Fig. 6(b) and (c). Fig. 7 plots the time evolution of the crown height for different Uw . We note that the height on the upward and downward sides obviously depends on the wall velocity. The larger the wall velocity, the higher the crown height on the upward side, while the lower the crown height on the downward side.
C. Ming, L. Jing / Computers and Mathematics with Applications (
)
–
9
Fig. 8. The velocity fields on the xz plane of symmetry for a drop impacting onto a moving thin film at different times. (a) Uw = 0 and (b) Uw = 1.
3.3. Velocity field around the ejecta Fig. 8 shows a comparison of the velocity fields on the xz plane of symmetry for the cases of the stationary wall and moving wall. When Uw = 0, at the time of impact and shortly thereafter, the inertial motion of the drop is transformed into the radial motion of liquid away from the axis; the radial movement of the liquid is indicated by the velocity vectors. When T = 0.2, the liquid in the film at radial distances greater than the position where the crown is formed has not started moving, as shown in Fig. 8(a). Initially, this discontinuity may aid in the formation of the crown. Subsequently, as liquid is fed into the crown, the velocity discontinuity persists. The velocity vectors in the crown point away from the wall in the radially outward direction as a result of its growth. As time increases, a characteristic feature of the flow appears near the tip of the crown; a vortex flow is formed near the tip of the crown due to the interaction between ejecta and surrounding fluid. We note that the velocity distribution around the drop is symmetric. We found that the maximum ejecta velocity can be several times higher in magnitude than the impact velocity, in agreement with the early observation of Toroddsen (1998), (2006) [30,31]. When Uw > 0, as shown in Fig. 8(b), at T = 0.2 the velocity discontinuity appears. As time increases, the velocity discontinuity gradually disappears on the downward side. This is because the moving wall causes the increase in relative velocity on the upward side and the decrease on the downward side. This rise and fall of the relative velocity explains the splash enhancement downward and attenuation upward. 3.4. The effect of the moving wall on the evolution of crown We now study the effect of the moving wall on the evolution of crown by examining the vorticity field. Fig. 9 shows the evolution of vorticity with time at Uw = 0 and 1 respectively, which allows us to quantitatively analyze the effect of moving wall on drop splashing. The vorticity distributions around the crown are directly correlated to the velocity field shown in Fig. 8. For the case of Uw = 0, Fig. 9(a) clearly illustrates that there is a pair of counter-rotating vortices around the crown due to the interaction of the ejecta with the surrounding fluid. The one above the ejecta is generated by the flow of fluid from the drop; another below the ejecta is generated by the flow of fluid from the initial layer. The strengths of vorticity are equal, which drive the liquid radially and symmetrically from the neck. The vortex expands radially rather rapidly in the initial stage of the impact during 0.2 ≤ T ≤ 0.4. The radial expansion slows down during T > 0.4. This vortex pair will assist in stripping the spike from the rim and lead to the secondary drop formation. For the cases of the moving wall, when Uw = 0.5 (not shown here), the patterns of vorticity contours’ evolution are found to be basically the same as those in the case of a normal impact, apart from the lack of the perfect rotational symmetry. A pair of counter-rotating vortices appears around the crown but their strengths are not uniform along the neck. When the wall velocity increases to Uw = 1, there is an obvious difference between the pair of counter-rotating vortices on the downward side, as shown in Fig. 9(b). As time increases, the viscous effect diffuses vorticity and significantly weakens the strength of the vorticity generated by the flow of the initial layer. When T > 0.2, the vortex on the downward side is quickly dissipated due to viscous. Thus, one vortex appears on the downward side. It means that the flow does not form the crown spikes on the downward side.
10
C. Ming, L. Jing / Computers and Mathematics with Applications (
)
–
Fig. 9. Evolution of the vorticity around the ejecta of a drop impacting onto a moving thin film at different times. (a) Uw = 0 and (b) Uw = 1.
4. Conclusions In this paper, a numerical study of a liquid drop impact on a moving wall with a pre-existing thin film of the same liquid is presented. A lattice Boltzmann model for incompressible two-phase flow at high density ratio is employed to investigate the dynamic behavior of drop splashing on the wall. In our simulation, an initial perturbation to flow is introduced to reflect the dynamics of drop evolution. The LBM-MPI code is validated by comparing the numerical results with the experimental observations for drop splashing. The qualitative features of the drop splashing on the moving wall at the early stage have been obtained by the present numerical simulation. We study the effects of the moving wall on the formation and evolution of ejecta and crown. It is found that an obvious effect is that the moving wall can act to enhance or suppress the splash. Further study will undoubtedly deepen our understanding about some important phenomena accompanying drop–wall interactions. Acknowledgments The authors gratefully acknowledge Professor Detlef Lohse for helpful discussions and his experimental photographs. They would also like to thank Dr. Z. Shang for his MPI-code development assistance. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17]
A.M. Worthington, A Study of Splashes, Green, London, 1908. H.E. Edgerton, J.R. Killian, Flash! Seeing the Unseen by Ultra-High-Speed Photography, Boston, Branford, 1954. G.E. Cossali, A. Coghe, M. Marengo, The impact of a single drop on a wetted solid surface, Exp. Fluids 22 (1997) 463–472. A.B. Wang, C.C. Chen, Splashing impact of a single drop onto very thin liquid films, Phys. Fluids 12 (2000) 2155–2158. C. Josserand, S. Zaleski, Droplet splashing on a thin liquid film, Phys. Fluids 16 (2003) 1650–1657. R. Rioboo, C. Bauthier, J. Conti, J. De Coninck, M. Voue, Experimental investigation of splash and crown formation during single drop impact on wetted surfaces, Exp. Fluids 35 (2003) 648–652. R.L. Vander-Wal, G.M. Berger, S.D. Mozes, Droplets splashing upon films of the same fluid of various depths, Exp. Fluids 40 (2006) 33–52. B. Ching, M.W. Golay, T.J. Johnson, Droplet impacts upon liquid surfaces, Science 226 (1984) 535–537. R. Krechetnikov, G.M. Homsy, Crown-forming instability phenomena in the drop splash problem, J. Colloid Interface Sci. 331 (2009) 555–559. D. Peregrine, The fascination of fluid mechanics, J. Fluid Mech. 106 (1981) 59–80. M. Rein, Phenomena of liquid drop impact on solid and liquid surfaces, Fluid Dyn. Res. 12 (1993) 61–93. A. Prosperetti, H. Oguz, The impact of drops on liquid surfaces and the underwater noise of rain, Annu. Rev. Fluid Mech. 25 (1993) 577–602. A.L. Yarin, Drop impact dynamics: splashing, spreading, receding, bouncing, Annu. Rev. Fluid Mech. 38 (2006) 159–192. J. Fukai, Z. Zhao, D. Poulikakos, C.M. Megaridis, O. Miyatake, Modeling of the deformation of a liquid droplet impinging upon a flat surface, Phys. Fluids A 5 (1993) 2588–2599. J. Fukai, Y. Shiiba, T. Yamamoto, O. Miyatake, D. Poulikakos, Wetting effects on the spreading of a liquid droplet colliding with a flat surface: experimant and modeling, Phys. Fluids 7 (1995) 236–247. M. Pasandideh-Fard, Y.M. Qiao, S. Chandra, J. Mostaghimi, Capillary effects during droplet impact on a solid surface, Phys. Fluids 8 (1996) 650–659. D.A. Weiss, A.L. Yarin, Single drop impact onto liquid films: neck distortion, jetting, tiny bubble entrainment, and crown formation, J. Fluid Mech. 385 (1999) 229–254.
C. Ming, L. Jing / Computers and Mathematics with Applications ( [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31]
)
–
11
M. Rieber, A. Frohn, A numerical study on the mechanism of splashing, Int. J. Heat Fluid Flow 20 (1999) 455–461. M. Bussmann, S. Chandra, J. Mostaghimi, Modeling the splash of a droplet impacting a solid surface, Phys. Fluids 12 (2000) 3121–3132. I.V. Roisman, C. Tropea, Impact of a drop onto a wetted wall: description of crown formation and propagation, J. Fluid Mech. 472 (2002) 373–397. S.D. Howison, J.R. Ockendon, J.M. Oliver, R. Purvis, F.T. Smith, Droplet impact on a thin fluid layer, J. Fluid Mech. 542 (2005) 1–23. S. Mukherjee, J. Abraham, Crown behavior in drop impact on wet walls, Phys. Fluids 19 (2007) 052103. G. Coppola, G. Rocco, L. de Luca, Insights on the impact of a plane drop on a thin liquid film, Phys. Fluids 23 (2011) 022105. T. Lee, C.L. Lin, A stable discretization of the lattice Boltzmann equation for simulation of incompressible two-phase flows at high density ratio, J. Comput. Phys. 206 (2005) 16–47. S. Chen, G. Doolen, Lattice Boltzmann method for fluid flows, Annu. Rev. Fluid Mech. 30 (1998) 329–364. L.-S. Luo, Theory of the lattice Boltzmann method: lattice Boltzmann models for nonideal gases, Phys. Rev. E 62 (2000) 4982–4996. W. van-Hoeve, T. Segers, H. Droes, D. Lohse, M. Versluis, A splash of red, in: The 61st Annual Meeting of the APS DFD, November 23–25, 2008 in San Antonio, Texas, 2008. N.Z. Mehdizadeh, S. Chandra, J. Mostaghimi, Formation of fingers around the edges of a drop hitting a metal plate with high velocity, J. Fluid Mech. 510 (2004) 353–373. J.C. Bird, S.S.H. Tsai, H.A. Stone, Inclined to splash: triggering and inhibiting a splash with tangential velocity, New J. Phys. 11 (2009) 063017. S.T. Thoroddsen, J. Sakakibara, Evolution of fingering pattern of an impacting drop, Phys. Fluids 10 (1998) 1359–1374. S.T. Thoroddsen, T.G. Etoh, K. Takehara, Crown breakup by marangoni instability, J. Fluid Mech. 557 (2006) 63–72.