## Abstract

A computational study was conducted of axisymmetric droplet impingement on a flat surface at low droplet Reynolds numbers. The study was motivated by the problem of deposition of melted volcanic ash particles within aircraft gas turbine engines. The computations were performed using the combined level-set volume-of-fluid method for droplet Reynolds numbers between 0.05 and 10. The computational predictions were validated using existing experimental data. The computations indicate that contact radius increases over short time in proportion to the square root of time, in agreement with short-time analytical predictions. Typical assumptions made in development of approximate droplet impingement models were evaluated for low Reynolds number droplet impingement. The droplet shape was well approximated by a truncated spherical cap through most of the impingement process. The surface area over which the droplet spreads increases with increase in Reynolds number. The axial velocity component was found to be approximately independent of radial location over most of the droplet, and the radial velocity component was observed to vary log-normally in the axial coordinate and linearly in radius. The energy dissipation rate was distributed throughout the droplet for low Reynolds numbers cases, but became increasingly localized near the contact line as the Reynolds number increased past unity.

## 1 Introduction

Droplet impact on a flat surface is a ubiquitous fluid mechanics problem for many diverse applications, including areas such as ink-jet printing [15], thermal spray coating [68], plasma spraying [9,10], raindrop erosion [11,12], pesticide deposition [13,14], spray impingement cooling [15,16], lithography [17], additive manufacturing [18], fire suppression [19], icing from freezing rain [2024], precipitation effects on erosion and performance degradation of turbine blades [25,26], and blood splatter forensics [27,28]. Previous research on droplet impingement has been extensive, including a wide range of studies examining physical processes such as droplet spreading and oscillation [2934], rim formation [3537], splash formation [3840], air entrapment during droplet impact [4145], effect of surface contact angle and surface roughness on droplet spreading [4648], anisotropic effects caused by surface motion [49], effects of fluid viscoelasticity [50,51], and effects of supersonic flow [52]. Reviews of droplet impingement physics are given by Josserand and Thoroddsen [53] and Yarin [54], and Attané et al. [55] present a nice overview of energy-based theoretical methods for modeling of droplet impacts.

This paper was motivated by the problem of molten volcanic ash impingement in gas turbine engines [56]. Large volcanic eruptions eject into the atmosphere millions of tons of small ash particles measuring 1–100 μm in diameter. Particles on the smaller end of this scale, with diameters less than about 10 μm, can remain suspended for several weeks and can spread over distances of 1000 km or more from their source [57]. Volcanic ash is an amorphous glass composed primarily of silica (SiO2), along with a range of other minerals [5861]. Significant damage has occurred when volcanic ash particles are ingested into aircraft gas turbine engines at sufficiently high particle concentrations over sufficiently long exposure time durations [6264]. Particles passing through the engine combustion section are exposed to high temperatures, with turbine inlet temperature often in excess of 1000 °C. Exposure to high temperatures causes ash particles to soften, and eventually to fully melt into highly viscous liquid droplets [60]. Upon impact on a downstream surface, the droplets spread on the surface while exchanging heat with the underlying substrate, resulting in formation of a glassy deposit on the turbine nozzle blades and guide vanes that can plug film cooling holes and cause cracking of ceramic coating used for heat protection of hot section components [62].

Studies of ash particle deposition in scaled gas turbine engines were presented by Bonilla [65] and Shinozaki et al. [66]. Bonilla [65] used coal fly ash with particle diameters in the same range as volcanic ash particles. He studied in particular the effect of particle Stokes number on the film cooling system for nozzle guide vanes. Shinozaki et al. [66] used an optical fiber viewing system to examine ingestion of volcanic ash particles with diameters ranging from 25 to –100 μm into a small turbojet engine. Further experimental studies in which softened and partially melted particles at high temperatures were shot at high velocity toward a flat plate at different angles of attack were reported by Lawrence [67], Dean et al. [68], and Taltavull et al. [69]. Lawrence [67] found that for temperatures below the material softening point, the temperature had a small effect on the coefficient of restitution of the impacting particle. Taltavull et al. [69] and Shinozaki et al. [66] note that very fine particles (e.g., with diameters of 2–3 μm or less) have sufficiently low Stokes number that they are unlikely to impact on the engine surfaces, but instead are swept through the engine following the gas stream. On the other hand, relatively large particles (e.g., diameters greater than 100 μm) rarely enter combustion chambers and when they do so, they do not usually heat up sufficiently during passage through the combustion chamber to fully melt. It is thus the intermediate size particles (with nominal diameter around 10 μm) that are of most concern for deposition in gas turbine engines.

Thermal analyses recently presented by Seksinsky [70] reveal that (1) spreading of molten volcanic ash droplets on engine surfaces occurs on a much more rapid time scale than the heat exchange between the droplet and the surface and (2) the temperature profile across volcanic ash particles with diameter of 10 μm or less is nearly uniform as particles are advected through the gas turbine engine. Consequently, this analysis confirms that it is reasonable to model the molten volcanic ash deposition process as a short time-scale isothermal droplet impingement problem followed by a longer time-scale heat exchange and solidification problem. Giordano et al. [71] reported on a comprehensive model for estimating the viscosity of molten volcanic ash particles as a function of chemical composition. The resulting high values of fluid viscosity indicate that impacting droplet Reynolds numbers will be low even at the high particle velocities typical of gas turbine engines.

Droplet dynamics following impact onto a flat surface is governed by three independent dimensionless parameters – the Weber number $We=ρdV02/σ$, the Ohnesorge number $Oh=μ/ρσd$, and the equilibrium (static) contact angle $θe$. Here, $ρ$ and $μ$ are the droplet density and viscosity, $σ$ is the droplet surface tension in air, and $d$ and $V0$ are the diameter and speed of the droplet at the moment of impact. Other dimensionless parameters can also be formed from these variables. For instance, the droplet Reynolds number $Re=ρdV0/μ$ can be related to the above parameters by $Re=We/Oh$. The capillary number $Ca=μV0/σ$ can be written as $Ca=We/Re=OhWe$. For very small droplets, gravity can be assumed to be negligible (i.e., the Bond number $Bo=ρgd2/σ≪1$). The droplet dynamics upon impact was divided into four regimes by Schiaffino and Sonin [72] based on the values of We and Oh. An extended version of this plot is shown in Fig. 1, which also includes lines indicating Re = 1 and Ca = 1 conditions. In Region I, the flow is nearly inviscid and driven by the impact pressure. The droplet fluid is observed to spread rapidly into a thin sheet, and for conditions with Oh < 1 to undergo a period of oscillations in which the flow alternates from spreading radially outward to rebounding back inward. In region II, the flow is again nearly inviscid, but the spreading is driven more by capillary force than by the impact pressure difference. The droplet again experiences oscillations that are damped over time by viscous effects. Region III is characterized by highly viscous droplets whose outward spreading is driven by capillary force and resisted by viscous shear. Cases with Oh > 1 in regimes III and IV are overdamped and the droplet does not exhibit oscillations. Region IV is characterized by droplets that are highly viscous so that the droplet spreading is driven by the impact pressure and resisted by viscous shear. The droplet is observed to maintain approximately a spherical cap form. The low Reynolds number problems under investigation in this paper are characterized by region IV collisions, for which Schiaffino and Sonin [72] estimate the typical contact line spreading velocity $Uc$ as $Uc/V0=Re$ and the time scale $τc$ of spreading as $τcV0/d=1/Re$.

Fig. 1
Fig. 1
Close modal

Much of the existing literature for droplet impingement has been developed for relatively high droplet Reynolds number, of order 102–104, although some experimental studies at low Reynolds number are available. Experimental data for droplet spreading for Reynolds numbers as low as nine were reported by Rioboo et al. [73] using glycerin and silicon oil, with droplet diameter of the order of 1 mm. Rioboo et al. [73] proposed that spreading factor of the droplet immediately after impact is dominated by a kinematic phase, during which contact line motion is determined only by geometric considerations and fluid movement within the droplet in primarily downward. Gordillo et al. [74] reported experimental results for droplet spreading and impact force using silicon oils with a wide range of viscosities, resulting in Reynolds numbers ranging from 0.07 to 104. However, the authors noted the presence of significant viscoelastic effects in the material response for the high viscosity (low Re) cases, which may significantly impact the spreading data. Experiments using silicon oil to achieve very high droplet viscosities were also reported by Langley et al. [45]. This paper focused mostly on characterizing a thin air film that forms underneath the droplet following impact for cases with low Reynolds number. To our knowledge, there have been no reported numerical studies of droplet impact at Reynolds numbers less than about 10.

Numerous simplifying assumptions are typically made in development of theoretical approximations for droplet spreading at high Reynolds numbers, such as the assumption of inviscid straining flow within the droplet or the assumption that energy dissipation occurs within a thin boundary layer under the droplet. The primary objective of this study is to examine the extent to which these simplifying assumptions remain valid for problems occurring at low Reynolds numbers. In cases where the simplifying assumptions break down, we attempt to develop alternative assumptions for low Reynolds number based on the computational data. The paper reports on a computational study of axisymmetric droplet impingement onto a flat surface over a range of Reynolds numbers from 0.05 to 10. The computational method used in the study is based on the coupled level-set—volume-of-fluid (CLSVOF) method, which is briefly described in Sec. 2. Computational results are summarized in Sec. 3, including validation of the computational results using the available low Reynolds number experimental data. Comparison of results at different Reynolds numbers and comparison of common assumptions for droplet velocity profile and shape with computational predictions are provided. Conclusions are given in Sec. 4.

## 2 Computational Method

Computations were conducted of droplet impingement onto a flat surface using the finite volume based formulation in ansysfluent v19.2 and employing the CLSVOF method to represent the droplet fluid and the external fluid phases [24,7577]. The flow field was computed in polar cylindrical coordinates ($r,θ,z$), with axisymmetric velocity profile given by
$u=u(r,z,t)er+w(r,z,t)ez$
(1)
The governing equations for the fluid flow are based on the assumption of an incompressible Newtonian viscous fluid with viscosity $μ$ and density $ρ$ that are functions of the level-set function $ϕ$, given by
$∇⋅u=0$
(2)
$ρ(ϕ)[∂u∂t+(u⋅∇)u]=−∇p+∇⋅[μ(ϕ)(∇u+∇uT)]−Fsf$
(3)
The magnitude of the level-set function $ϕ$ is the distance from the interface, and the sign of the level-set function is positive in the liquid phase and negative in the gas phase. The level-set function is evolved in time using the standard equation
$∂ϕ∂t+∇⋅(uϕ)=0$
(4)

The level-set function is regularly reinitialized using the geometrical interface-front construction method to preserve the property $|∇ϕ|=1$ [78].

The distributed surface tension force $Fsf$ is evaluated using an extension of the Brackbill et al. [79] equation as
$Fsf=2H(ϕ)δ(ϕ)σκ n$
(5)
where the interface unit normal n and mean curvature $κ$ are defined by
$n=∇ϕ|∇ϕ|ϕ=0,κ=∇⋅n$
(6)
A phase-smoothed Heaviside function $H(ϕ)$ and Dirac delta $δ(ϕ)$ can be defined for an interface with thickness a as [80]
$H(ϕ)={0forϕa (liquid phase)12[1+ϕa+1πsin(πϕa)]for−a<ϕ
(7)
$δ(ϕ)={0for|ϕ|>a1+cos(πϕ/a)2afor−a<ϕ
(8)
The fluid density and viscosity were specified in terms of the values ($ρℓ,μℓ$) in the liquid and the values ($ρg,μg$) in the gas as
$ρ(ϕ)=[1−H(ϕ)]ρg+H(ϕ)ρℓ$
(9a)
$μ(ϕ)=[1−H(ϕ)]μg+H(ϕ)μℓ$
(9b)
The contact line coincides with the intersection of the zero level-set surface $ϕ=0$ with the wall ($z=0$). A dynamic boundary condition is used in which the interface normal n at interface cells adjacent to the wall is set as
$n=nW cos θE+tW sin θE$
(10)

where $nW$ and $tW$ are unit vectors normal and tangent to the wall, respectively, and $θE$ is the equilibrium contact angle. For all computations reported in the paper, we set θE = 18 deg.

A block-structured mesh containing four blocks was generated using the ansys meshing utility, as shown in Fig. 2(a). The mesh is uniform in blocks 1 and 3, where the mesh size is set equal to a value Le in block 1 and 0.62 in block 3. In blocks 2 and 4, the mesh gradually varied from the small grid size in block 1 to the large grid in block 3. We set R = H =5d as a good compromise between computational cost and minimizing effect on simulation residuals. Tests with different grid sizes showed no significant effect on the computational results. A fully implicit coupled pressure and velocity solver was used to solve Eqs. (2) and (3) on the block-structured grid. The momentum Eq. (3) and the level-set Eq. (4) were solved using a second-order upwind scheme. The pressure interpolation was performed using the PRESTO! (PREssure STaggering Option) scheme, and gradients were computed using a cell-based least-squares method.

Fig. 2
Fig. 2
Close modal
The domain initial and boundary conditions are shown in Fig. 2(b). The droplet is entrained in a downward axisymmetric stagnation-point flow. A uniform downward velocity was set at the velocity inlet, and the no-slip condition was prescribed at the wall boundary. The density and viscosity of the gas phase were set to 1000 times smaller than the liquid phase. Runs were also performed with significantly larger viscosity differences between the gas and liquid phases, but they were found to eventually become unstable. The offset distance between the droplet surface and the wall at the initial time was set equal to $d/10$. This value was chosen to allow the droplet to deform as it encounters the thin slip-film when approaching the flat surface, while not un-necessarily wasting computation time as the droplet approaches the wall. The transient droplet solution was initialized using the downward axial velocity from the steady-state solution for the gas phase. All variables in the computations were nondimensionalized using the initial droplet diameter d, the impact velocity $V0$, and the liquid density $ρℓ$ as
$u′=u/V0,x′=x/d,p′=p/ρℓV02,t′=tV0/d$
(11)

Primes on dimensionless variables are dropped in the remainder of the paper for convenience.

## 3 Computational Results

Computations of droplet impingement were performed for parameter values typical of volcanic ash droplet impingement in gas turbine engines, as listed in Table 1. The nominal values in Table 1 lead to dimensionless parameter values of Re = 0.19, Ca = 14,000, We = 2667, Oh = 271, and Bo = $6.5×10−6$. The dimensionless parameter values used for the four computational cases examined in the paper are listed in Table 2. To maximize computational efficiency, a variable time stepping method was used with maximum Courant number of 0.75. The grid size near the contact line was selected as Le= 0.006. Lack of grid sensitivity for the impinging droplet simulations was verified by repeating one simulation with four values of the grid cell length near the contact line, ranging from $Le=0.004$ to 0.01 (see also related discussion by Afkhami et al. [81]).

Table 1

Typical values of input parameters for volcanic ash impingement in gas turbine engines

PropertySymbolUnitsRange of valuesNominal value
Diametera$d0$μm3–3010
Impact velocityb$V0$m/s130–450200
ViscositycμPa·s10–10525
Densitydρkg/m32350–24502400
Surface tensioneσN/m0.35–0.370.36
Equilibrium contact anglef$θe$deg18–7155
PropertySymbolUnitsRange of valuesNominal value
Diametera$d0$μm3–3010
Impact velocityb$V0$m/s130–450200
ViscositycμPa·s10–10525
Densitydρkg/m32350–24502400
Surface tensioneσN/m0.35–0.370.36
Equilibrium contact anglef$θe$deg18–7155
a

Based on a distal particle size distribution analysis, Dacre et al. [57] found that 90% of ash particles (by mass) ejected from Eyjafjallajökull volcano were in size range of 3–30 μm.

b

Velocity range within the turbine section of a gas turbine engine.2

c

Viscosity varies from 10–105 over the range of temperatures typical of GTEs, as recorded above, based on measurements of Giordano et al. [71]. Values are strongly dependent on temperature, but the dependency is well known and a good correlation exists. Viscosity nominal value based on measurements by Song et al. [60] at 1200 °C.

d

Based on data from Vogel et al. [82], Wilson et al. [83], and Shipley and Sarna-Wojcicki [84] for amorphous volcanic glass particles.

e

Based on data from Song et al. [60] and Li et al. [85]. Surface tension has negligible dependence on temperature.

f

Song et al. [60] found that the contact angle has a strong dependence on temperature. Nominal value used based on measurements at 1200 °C.

Table 2

Listing of parameter values in different computational cases examined

CaseReynoldsCapillaryWeberOhnesorge
10.05$5.3×104$26661032
20.19$1.4×104$2666271
312633266651
41025826665
CaseReynoldsCapillaryWeberOhnesorge
10.05$5.3×104$26661032
20.19$1.4×104$2666271
312633266651
41025826665

### 3.1 Results for a Typical Low Reynolds Number Case.

Typical results for volcanic ash transported in a gas turbine engine (case 2, with Re = 0.19) are plotted in Fig. 3 showing a time series of the velocity magnitude both inside and external to the droplet during collision. The droplet surface ($ϕ=0$) is indicated using a solid black line. The droplet first makes contact with the surface in frame 3a, at which time the velocity within the droplet is nearly unchanged from the initial state, save for a decrease in velocity magnitude within a small region surrounding the collision point. The external fluid exhibits a high-velocity jet away from the droplet along the impact surface, as the squeeze film is emptied of fluid by the impinging droplet. Frames 3b3e show the droplet flattening as the contact line moves outward along the impact surface. The velocity magnitude rapidly decreases both within the droplet and external to the droplet. A jet of external fluid is emitted from the vicinity of the contact line in frames 3b and 3c due to the corner flow in this region, but the magnitude decreases rapidly.

Fig. 3
Fig. 3
Close modal

The pressure field internal to the droplet is shown in Fig. 4(a) at t =0.36. A region of high pressure is observed along the impact surface, just inward of the contact line. This pressure field is very similar to that obtained in the droplet impact simulations of Philippi et al. [86] and Eggers et al. [35]. Figure 4(b) shows contours of the step function $H(ϕ)$ indicating the fluid phase at a time of 0.36, where the gray region of the plot indicates the region inside the droplet. Streamlines of the velocity field are also shown. The presence of a recirculating vortical flow in the external fluid streamlines just outside of the impacting droplet is similar to the flow structures present in the numerical computations of droplet impingement by Guo et al. [38]. The phase plot in Fig. 4(b) exhibits a thin strip along the bottom surface of the droplet in which the phase transitions back to the external fluid. This strip, which was a common feature observed in computations that we performed for low Reynolds number droplet impingement, is consistent with observations of air film formation beneath impinging droplets made by numerous researchers [4143,45,8790]. Experimentally, studies indicate that this air film beneath impacting droplets is sometimes very thin, only a few tens of nanometers thick [41].

Fig. 4
Fig. 4
Close modal

Theories for droplet impingement on surfaces typically utilize a balance between the rate of change of potential and kinetic energy in the droplet and the rate of energy dissipation [1,29,55]. Consequently, it is of interest to examine the contours of some of these energy terms within the droplet. A time series showing the contours of the kinetic energy $12H(ϕ)(u2+w2)$ for case 2 is presented in Fig. 5 for the same six times as used in Fig. 3. The kinetic energy contour levels are reset at each time frame since the value of kinetic energy decreases significantly during the time series. The value of the kinetic energy gradually increases with height and with radius within the droplet at any time during the spreading process.

Fig. 5
Fig. 5
Close modal
A time series of the contours of the rate of hydrodynamic energy dissipation $πh$ (nondimensionalized by $ρℓd2V03$) within the droplet is plotted in Fig. 6, which is defined by
$πh(r,z,t)=2H(ϕ)Re[(∂u∂r)2+(ur)2+(∂w∂z)2+12(∂u∂z+∂w∂r)2]$
(12)
Fig. 6
Fig. 6
Close modal

At the droplet contact time in frame 6a, the energy dissipation rate is dominated by a region surrounding the contact point, and it is small in the upper part of the droplet. As the droplet spreads on the impacting surface, the magnitude of the energy dissipation rate decreases but the region with significant values of energy dissipation rate spreads out over a wide region within the droplet. We note that most theories for droplet impingement at high Reynolds number assume that all energy dissipation occurs within a thin boundary layer lying between the droplet and the impact surface [1,29]. By contrast, our direct computation results for energy dissipation rate at low Reynolds number exhibit significant energy dissipation rate values within the body of the droplet and relatively small dissipation rate along the impact surface (except very close to the contact line).

### 3.2 Reynolds Number Effect.

The effect of Reynolds number on the droplet impingement computations was examined by repeating the above computations with four different Reynolds number values over the interval $0.05≤Re≤10$, with the same value of We (see Table 2). The computed contact radius $rc$ and the droplet center of mass height $zm$ are plotted in Fig. 7 as functions of time for Re values of 0.05, 0.19, 1, and 10. The maximum spread factor is observed to increase as the Reynolds number increases.

Fig. 7
Fig. 7
Close modal
A comparison of our computational results for time variation of the contact radius with experimental data and analytical predictions from other investigators is shown in Fig. 8, with both a long-time variation plot (a) and a short-time variation plot (b). Our computational data are plotted as solid curves, with increasing contact radius as the Reynolds number increases (as indicated by the arrow from A to D). The computational results are compared to experimental data for low Reynolds number droplet impingement obtained by Gordillo et al. [74] for Reynold's numbers of 0.11, 0.7, and 7, by Rioboo et al. [73] for a Reynolds number of 9, and by Langley et al. [45] for Reynolds numbers of 0.1 and 1. A short-time analytical result from Philippi et al. [86] based on approximate solution of the inviscid flow equations is also plotted in Fig. 8, which can be written in dimensionless form as
$rc=3t/2$
(13)
Fig. 8
Fig. 8
Close modal

Both the computational results and the experimental data exhibit the trend that the contact radius decreases with decrease in Reynolds number. The Philippi et al. [86] inviscid flow solution yields a higher contact radius than any of the finite Reynolds number computational or experimental results.

It might first be noted that there are significant differences between different experimental data in Fig. 8. For instance, the Langley et al. [45] data for Re = 0.4 (circles) are close to the Gordillo et al. [74] data for Re = 0.7 (squares), but significantly above their data for Re = 0.11 (deltas). Similarly, the Langley et al. [45] data for Re = 4 are close to the Rioboo et al. [73] for Re = 9 (diamonds) but above the Gordillo et al. [74] data for Re = 7 (x′s). Our computational predictions agree well with the experimental data for short times (Fig. 8(b)), but exhibit significant differences during longer times (Fig. 8(b)). For instance, the computational prediction for Re = 10 (line D) is slightly above the Gordillo et al. [74] data for Re = 7 (x′s) but well below the Rioboo et al. [73] for Re = 9 (diamonds). The computational prediction for Re = 1 (line C) is close to the Gordillo et al. [74] data with 0.7 (squares) up to a time of about t =3, after which the experimental data suddenly flattens out. A similar flattening out of the contact radius is observed in the Gordillo et al. [74] data for Re = 0.11 for t >0.2. Gordillo et al. [74] in fact discussed this flattening behavior observed for low Reynolds number cases and attributed it to presence of viscoelastic effects present in their silicon oils at high viscosities. These viscoelastic effects, which are not present in the computations, explain the differences between the Gordillo et al. [74] data and predictions of lines A and B after finite time. The data in Fig. 8(a) were replotted as a log–log plot in Fig. 9. It is noted that all of the experimental and computational data exhibit a region in with contact radius varies in proportion to $t1/2$ at short times, in agreement with the Philippi et al. [86] solution.

Fig. 9
Fig. 9
Close modal
We have also examined the total kinetic and potential energy of the droplet, $Ec$ and $Ep$, which can be written in dimensionless form as $ec=Ec/(ρℓd3V02)$ and $ep=Ep/(ρℓd3V02)$. The potential energy consists of both gravitational potential energy and interfacial potential energy. The potential energy is given by the sum of the air–liquid interfacial energy $σAGL$ and the liquid–solid interfacial energy $πRc2σSL$, minus the interfacial energy $πRc2σSG$ that was present before the surface was covered by the liquid. Here, $AGL$ is the area of the liquid–gas interface, $Rc$ is the contact line radius on the planar solid substrate, and $σ$, $σSL$, and $σSG$ denote the liquid–gas, solid–liquid and solid–gas surface tensions, respectively. The resulting interfacial potential energy is given by
$Ep=σAGL+πRc2(σSL−σSG)$
(14)
Young's equation defines the equilibrium contact angle $θe$ by
$σ cos θe=σSG−σSL$
(15)
so that Eq. (13) reduces to
$Ep=σ(AGL−πRc2 cos θe)$
(16)
Nondimensionalizing the interfacial area and the contact radius using the droplet radius d as $aGL=AGL/d2$ and $rc=Rc/d$, the potential energy can be written in dimensionless form as
$ep=1We(aGL−πrc2 cos θe)$
(17)
The dimensionless kinetic energy within the droplet is given by
$ec=12∫VH(ϕ)(u2+w2) dv$
(18)

Plots of the time variation of $ep$ and $ec$ are shown in Fig. 10 for the same Reynolds number values examined in Fig. 7. A log–log plot is used for the kinetic energy since it reduces rapidly following droplet impact. The potential energy flattens out at a value between 50 and 60% of the value for the initial spherical droplet, whereas the kinetic energy decreases by nearly eight orders of magnitude during the computation. The potential energy of the droplet is observed to decrease as the Reynolds increases. This trend is further examined in Fig. 11, where the product of potential energy and Weber number in the long-time steady-state configuration of the droplet is plotted versus Reynolds number both for our computational results and for values calculated from the droplet impingement videos given in the Gordillo et al. [74] supplementary information. The multiplication by We in Fig. 10 was used to normalize the data, as suggested by Eq. (16), since our computations and the Gordillo et al. [74] experiments were conducted at different Weber number values. Both our computations and the Gordillo et al.'s [74] experiments indicate a decrease in potential energy as the Reynolds number increases.

Fig. 10
Fig. 10
Close modal
Fig. 11
Fig. 11
Close modal

Figure 12 shows the droplet shape and contours of the local dissipation rate plotted at two different times for three different Reynolds numbers. The dissipation rate is largest in the low Reynolds number cases, and it occurs over a larger fraction of the droplet volume. This trend is expected because the dissipation defined in Eq. (11) scales with the inverse Reynolds number. The velocity gradients in Eq. (11) define the shape of the contours. For the Re =10 case, the dissipation rate is dominated by velocity gradients near the contact line, whereas dissipation rate is more spread out in the low Reynolds number cases.

Fig. 12
Fig. 12
Close modal
We also examined the energy budget assumption for the droplet of the form
$Πh=−ddt(ec+ep)$
(19)

where $Πh$ is the integral over the droplet of the local energy dissipation rate $πh$ defined in Eq. (12). A log–log plot showing time variation of both the left- and right-hand sides of Eq. (19) is given in Fig. 13. Both sides of this equation are found to be close to each other for short times, less than about $t=1$, but a difference develops for longer times. The absolute value of this difference is quite small, however, since both of these terms reduce to less than about 5% of their initial values by t =1. The difference in these values can readily be accounted for by factors such as energy dissipation by the moving contact line or work done on the external gas phase. For times greater than t =10, both terms are very small and seem to be dominated by numerical truncation error.

Fig. 13
Fig. 13
Close modal

### 3.3 Velocity Profile Within the Droplet.

An important feature of theoretical models for impacting droplets is the assumption of a velocity profile within the droplet. The droplet flow field is often assumed to have the form of an inviscid stagnation-point flow [1,29,34,48] such that
$u=sr,w=−2sz$
(20)
where the straining rate $s(t)$ is determined by an additional boundary condition, such as the kinematic condition $w(0,h,t)=h˙$ at the top of the droplet. This form of velocity profile is appropriate at sufficiently long times after impact for cases with high Reynolds numbers. Since the radial velocity in Eq. (20) does not satisfy the no-slip condition, it is assumed to have a thin boundary layer underneath it in which most of the viscous dissipation occurs. Alternatively, Madejski [90] developed a theory using a cylinder model for the lamellar film that forms at long time as the droplet spreads on the surface, in which the radial velocity was assumed to have the form
$u=czr,w=−cz2$
(21)

where $c=c(t)$ is a coefficient with dimensions of 1/m⋅s. This velocity profile is more appropriate for droplet impact cases with low Reynolds number at long time, as the radial velocity satisfies the no-slip condition on the substrate surface.

A similarity theory for droplet impact for short times after impact was developed by Philippi et al. [86] based on solution of the potential flow equations, which gives droplet velocity and pressure fields as a function of time after impact for cases with high droplet Reynolds numbers. Philippi et al. [86] also provide a solution for the viscous boundary layer flow lying under the self-similar inviscid flow within the droplet, where the boundary layer velocity field varies as a function of $zRe/t$, in agreement with the classical Stokes impulsive motion solution. Similar results are noted by Wildeman et al. [91] based on computational solution of droplet impacts with Reynolds number in the interval $100≤Re≤1000$.

In this subsection, we explore the direct computation results for droplet impingement to determine a reasonable approximate form for the velocity profile within the droplet that is appropriate for low Reynolds number collisions. One common feature of the velocity profiles in Eqs. (20) and (21) is that the axial velocity component w is approximately independent of radius r, at least near the central and middle parts of the droplet. This feature is examined in Fig. 14, in which the axial velocity is plotted as a function of height at three different radial locations (r =0, 0.15, and 0.3) and at different times during the droplet impact. It is observed that at each of these times, the velocity profile at each radial location is very similar, indicating that the axial velocity can reasonably be assumed to be independent of radius, or $w=w(z,t)$, for low Reynolds number droplet impingement problems. In this case, the radial velocity can be obtained from the axisymmetric continuity equation as
$u(r,z,t)=−r2∂w∂z$
(22)
Fig. 14
Fig. 14
Close modal
We seek a functional fit for the radial and axial velocity profiles, which satisfies the no-slip and no-penetration boundary condition on the substrate surface, which from Eq. (21) requires that $w=∂w/∂z=0$ on $z=0$. Satisfaction of these boundary condition can be ensured by approximating the gradient of the axial velocity as proportional to a log-normal distribution $f(z)$ such that
$∂w∂z=−CLNf(z)=−CLNσLN z2πexp(−(lnz−μLN)22σLN2)$
(23)
and from Eq. (22)
$u=CLNr2f(z)$
(24)
Here $CLN$, $μLN$, and $σLN$ are undetermined fitting coefficients. The axial velocity profile $w(z,t)$ can be obtained by integrating Eq. (23) to obtain
$w(z,t)=−CLN[Φ(z)−Φ(0)]$
(25)

where $Φ(z)=1−0.5 erfc[lnz−μLN/2 σLN]$ is the cumulative normal distribution function and erfc(⋅) is the complementary error function.

A plot demonstrating these functional fits for the case with Re = 0.19 is given in Fig. 15. The radial velocity profile at r =0.15 was fit using Eq. (24), which is shown in Fig. 15(a). The best-fit coefficient values were obtained by searching for values yielding minimum least-square error with increment sizes $ΔCLN=0.008$, $ΔμLN=0.004$, and $ΔσLN=0.002$. Corresponding plots for the axial velocity profile were obtained using Eq. (25) with the same coefficient values as used for the curves in Fig. 14(a). The resulting predictions for axial velocity profiles are shown in Fig. 15(b) to be in excellent agreement with axial velocity data taken along the symmetry axis ($r=0)$. This procedure was repeated for the velocity profiles obtained for computations with all four Reynolds number values listed in Table 2, and the best-fit values for the coefficients $CLN$, $μLN$, and $σLN$ are plotted in Fig. 16 as functions of Reynolds number. The data indicate that $μLN$ and $σLN$ are approximately constant while $CLN$ decreases significantly in time. We do not see any clear trend in values of $μLN$ and $σLN$ with Reynolds number, but rather the values are scattered around mean values of approximately −1.1 and 0.85, respectively.

Fig. 15
Fig. 15
Close modal
Fig. 16
Fig. 16
Close modal

### 3.4 Droplet Shape During Impact.

It is also necessary in many analytical theories for droplet impact to assume a family of shapes for the droplet. Common choices are a truncated sphere or a cylinder (6), although other options such as the rimmed cylinder or cylinder with spherical rims are sometimes employed. A comparison of the computed droplet shape over a series of times during droplet impingement with the popular cylinder and truncated sphere shapes is given in Fig. 17 for the Re = 0.19 case. The cylinders and truncated spheres shown for comparison have the same volume and center of mass as the computed droplet. It is apparent from these results that the truncated sphere is an excellent approximation of the droplet shape during the early stages of droplet impingement, and even for long times remains a reasonable approximation to the computed profile.

Fig. 17
Fig. 17
Close modal

In Fig. 18, a time series showing the droplet shapes (in gray shading) and velocity vectors during droplet impingement are compared for cases with two Reynolds numbers: Re = 0.19 (on the left) and Re = 10 (on the right). All other parameters except the fluid viscosity are the same for the two computations. The lower Reynolds number problem clearly shows less spreading and less droplet flattening at times t =0.614 and 1.752 after the initial impact has occurred. The higher Reynolds number case not only spreads more than the lower Reynolds number case, but it also spreads more rapidly along the substrate surface.

Fig. 18
Fig. 18
Close modal

## 4 Conclusions

Small volcanic ash particles (with diameter from 3–30 μm) can be suspended in the atmosphere for periods of several weeks following a volcanic eruption. Jet aircraft flying through volcanic ash clouds ingest suspended particles into their gas turbine engines. These small particles can heat up and melt in the combustion portion of the engine, and then deposit on the downstream sections of the engine. Such deposits are of particular concern when they block film cooling holes on the turbine blades, leading the blades to heat up and the material to fail. Scaling studies of molten volcanic ash particle impingement in gas turbine engines indicate that this collision occurs at high Weber numbers (We = O(1000)) but at low droplet Reynolds numbers (Re = O(0.1 − 1)). A survey of the literature on droplet impingement on flat surfaces reveals that nearly all existing work was done for high droplet Reynolds numbers, and does not adequately capture the parameter range involved in the volcanic ash droplet impingement problem. This paper presents a computational study of droplet impingement at high Weber number (We = 2776) and at relatively low Reynolds number ($0.05≤Re≤10$). This Reynolds number range was selected because it extends into the range of interest for molten volcanic ash impingement while also extending up into the bottom part of the range covered by previous moderate Reynolds number studies. The computations were conducted using an axisymmetric finite volume approach with the CLSVOF method used for interface tracking. Validation of the computations was performed by comparing predictions for contact radius variable with time to experimental data from Gordillo et al. [74] for $0.11≤Re≤7$, experimental data Rioboo et al. [73] for Re = 9, and the experimental data of Langley et al. [45] for Re = 0.1 and 1. The comparison of the computational predictions with the Gordillo et al. data was good for short times, but the presence of viscoelastic effects caused a spreading of the droplet to terminate prematurely (as discussed by Gordillo et al. [74]). The results were also compared over short times to the inviscid-flow analytical solution of Philippe et al. [86], which seems to form a high Reynolds number upper bound for the data. At all Reynolds number values examined, the contact radius is observed to increase in proportion to the square root of time for short times after impact, in agreement with the Philippe et al. [86] prediction. In general, the results show that the droplet spreads less rapidly, and reaches a smaller value of the maximum spread radius, as the Reynolds number decreases.

The computational results were used to examine the time and spatial variation of the kinetic energy and dissipation rate, as well as the time variation of the potential energy of the droplet. The dissipation rate decreases with increase in Reynolds number. At low Reynolds number, the dissipation rate occurs throughout the middle part of the droplet, but as the Reynolds number approaches a moderate value of about ten, the dissipation rate becomes increasingly dominated by the region around the contact line. Spatial patterns of dissipation rate differ significantly from assumptions made in high Reynolds number approximate models (e.g., Kim and Chun [29]), which typically assume that all dissipation occurs in a thin boundary layer underneath the essentially inviscid droplet flow. Droplet kinetic energy increases with Reynolds number, corresponding with the observed decrease in dissipation rate. The droplet potential energy was found to decrease with increase in Reynolds number in the computed data. This trend was also observed in the experimental data of Gordillo et al. [74].

The computed droplet shape was shown to be well approximated by a truncated sphere throughout the low Reynolds number impingement process. The axial velocity was found to be nearly independent of radius within the central part of the droplet, from which it follows that the radial velocity increases linearly with radius. The radial velocity profile was shown to fit closely to a log-normal form, from which the axial velocity profile was shown to be proportional to the cumulative probability distribution function, which can be written in terms of the complementary error function. The velocity profiles were shown to approach the form proposed by Madejski [90] at long times, in which the radial velocity varies linearly in height z and the axial velocity varies quadratically in z.

## Acknowledgment

The authors gratefully acknowledge support and encouragement from James Uhlman and John Grant at New England Research and Development LLC, from John Lekki and James Zakrajsek at NASA Glenn Research Center, from Allan Aubert and Michael Smith at the Naval Air Systems Command, and from Ryan Lundgren at Pratt & Whitney.

## Funding Data

• NASA EPSCoR (Grant No. NNX15AK55A; Funder ID: 10.13039/100000104).

• Naval Air Systems Command, U.S. Department of Defense (Contract N6833517C0495; Funder ID: 10.13039/100010464).

## Nomenclature

• a =

interface thickness in volume-of-fluid method (dimensionless)

•
• aGL =

air–liquid interfacial surface area (dimensionless)

•
• AGL =

air–liquid interfacial surface area (m2)

•
• Bo =

bond number (dimensionless)

•
• c =

coefficient in velocity profile (1/m-s)

•
• $CLN$ =

log-normal coefficient (dimensionless)

•
• d =

droplet diameter (m)

•
• ec =

total kinetic energy (dimensionless)

•
• ep =

total potential energy (dimensionless)

•
• Ec =

total kinetic energy (J)

•
• Ep =

total potential energy (J)

•
• f =

log-normal distribution function (dimensionless)

•
• Fsf =

surface tension force (N)

•
• g =

gravity (m/s2)

•
• h =

height of top of droplet (dimensionless)

•
• H =

smoothed Heaviside function (dimensionless)

•
• Le =

smallest mesh size (dimensionless)

•
• n =

unit normal (dimensionless)

•
• Oh =

Ohnesorge number (dimensionless)

•
• p =

pressure field (Pa)

•
• r =

•
• rc =

•
• Rc =

•
• Re =

Reynolds number (dimensionless)

•
• s =

straining rate (dimensionless)

•
• t =

time (s)

•
• u =

•
• u =

velocity vector (dimensionless)

•
• Uc =

•
• V0 =

droplet impact speed (m/s)

•
• w =

axial velocity component (dimensionless)

•
• We =

Weber number (dimensionless)

•
• z =

axial coordinate (m)

•
• zm =

center of mass (dimensionless)

### Greek Symbols

Greek Symbols

• δ =

Dirac delta function (dimensionless)

•
• θ =

•
• θe =

•
• κ =

mean curvature (m−1)

•
• μ =

fluid viscosity (Pa·s)

•
• $μg$ =

gas phase viscosity (Pa·s)

•
• $μℓ$ =

liquid phase viscosity (Pa·s)

•
• $μLN$ =

log-normal coefficient (dimensionless)

•
• πh =

dimensionless dissipation rate (dimensionless)

•
• ρ =

fluid density (kg/m3)

•
• $ρg$ =

liquid phase density (kg/m3)

•
• $ρℓ$ =

liquid phase density (kg/m3)

•
• σ =

liquid–gas surface tension (N/m)

•
• $σLN$ =

log-normal coefficient (dimensionless)

•
• $σSG$ =

solid–gas surface tension (N/m)

•
• $σSL$ =

liquid–solid surface tension (N/m)

•
• τc =

•
• ϕ =

level-set variable (m)

•
• Ф =

cumulative normal distribution function (dimensionless)

## References

1.
Bechtel
,
S. E.
,
Bogy
,
D. B.
, and
Talke
,
F. E.
,
1981
, “
Impact of a Liquid Drop Against a Flat Surface
,”
IBM J. Res. Dev.
,
25
(
6
), pp.
963
971
. 10.1147/rd.256.0963
2.
Son
,
Y.
, and
Kim
,
C.
,
2009
, “
Spreading of Inkjet Droplet of Non-Newtonian Fluid on Solid Surface With Controlled Contact Angle at Low Weber and Reynolds Numbers
,”
J. Non-Newtonian Fluid Mech.
,
162
(
1–3
), pp.
78
87
.10.1016/j.jnnfm.2009.05.009
3.
Zable
,
J. L.
,
1977
, “
Splatter During Ink Jet Printing
,”
IBM J. Res. Dev.
,
21
(
4
), pp.
315
320
.10.1147/rd.214.0315
4.
Lim
,
T.
,
Han
,
S.
,
Chung
,
J.
,
Chung
,
J. T.
,
Ko
,
S.
, and
Grigoropoulos
,
C. P.
,
2009
, “
Experimental Study on Spreading and Evaporation of Inkjet Printed Pico-Liter Droplet on a Heated Substrate
,”
Int. J. Heat Mass Transfer
,
52
(
1–2
), pp.
431
441
.10.1016/j.ijheatmasstransfer.2008.05.028
5.
van Dam
,
D. B.
, and
Le Clerc
,
C.
,
2004
, “
Experimental Study of the Impact of an Ink-Jet Printed Droplet on a Solid Substrate
,”
Phys. Fluids
,
16
(
9
), pp.
3403
3414
.10.1063/1.1773551
6.
Alavi
,
S.
,
Passandideh-Fard
,
M.
, and
Mostaghimi
,
J.
,
2012
, “
Simulation of Semi-Molten Particle Impacts Including Heat Transfer and Phase Change
,”
J. Therm. Spray Technol.
,
21
(
6
), pp.
1278
1293
.10.1007/s11666-012-9804-8
7.
Attinger
,
D.
,
Zhao
,
Z.
, and
Poulikakos
,
D.
,
2000
, “
An Experimental Study of Molten Microdroplet Surface Deposition and Solidification: Transient Behavior and Wetting Angle Dynamics
,”
ASME J. Heat Transfer
,
122
(
3
), pp.
544
556
.10.1115/1.1287587
8.
Pasandideh-Fard
,
M.
,
Chandra
,
S.
, and
Mostaghimi
,
J.
,
2002
, “
A Three-Dimensional Model of Droplet Impact and Solidification
,”
Int. J. Heat Mass Transfer
,
45
(
11
), pp.
2229
2242
.10.1016/S0017-9310(01)00336-2
9.
Bertagnolli
,
M.
,
Marchese
,
M.
,
Jacucci
,
G.
,
Doltsinis
,
I.
,
St
., and
Noelting
,
S.
,
1995
, “
Modelling the Impact of Particles on a Rigid Substrate Under Plasma Spraying Conditions
,”
J. Therm. Spray Technol.
,
4
(
1
), pp.
41
49
.10.1007/BF02648527
10.
McDonald
,
A.
,
Lamontagne
,
M.
,
Moreau
,
C.
, and
Chandra
,
S.
,
2006
, “
Impact of Plasma-Sprayed Metal Particles on Hot and Cold Glass Surfaces
,”
Thin Solid Films
,
514
(
1–2
), pp.
212
222
.10.1016/j.tsf.2006.03.010
11.
Abuku
,
M.
,
Janssen
,
H.
,
Poesen
,
J.
, and
Roels
,
S.
,
2009
, “
Impact, Absorption and Evaporation of Raindrops on Building Facades
,”
Build. Environ.
,
44
(
1
), pp.
113
124
.10.1016/j.buildenv.2008.02.001
12.
Zhao
,
R.
,
Zhang
,
Q.
,
Tjugito
,
H.
, and
Cheng
,
X.
,
2015
, “
Granular Impact Cratering by Liquid Drops: Understanding Raindrop Imprints Through an Analogy to Asteroid Strikes
,”
Proc. Natl. Acad. Sci. U. S. A.
,
112
(
2
), pp.
342
347
.10.1073/pnas.1419271112
13.
Bergeron
,
V.
,
Bonn
,
D.
,
Martin
,
J. Y.
, and
Vovelle
,
L.
,
2000
, “
Controlling Droplet Deposition With Polymer Additives
,”
Nature
,
405
(
6788
), pp.
772
775
.10.1038/35015525
14.
Wirth
,
W.
,
Storp
,
S.
, and
Jacobsen
,
W.
,
1991
, “
Mechanisms Controlling Leaf Retention of Agricultural Spray Solutions
,”
Pest. Sci.
,
33
(
4
), pp.
411
420
.10.1002/ps.2780330403
15.
Kim
,
J.
,
2007
, “
Spray Cooling Heat Transfer: The State of the Art
,”
Int. J. Heat Fluid Flow
,
28
(
4
), pp.
753
767
.10.1016/j.ijheatfluidflow.2006.09.003
16.
Li
,
S. C.
,
Libby
,
P. A.
, and
Williams
,
F. A.
,
1995
, “
Spray Impingement on a Hot Surface in Reacting Stagnation Flows
,”
AIAA J.
,
33
(
6
), pp.
1046
1055
.10.2514/3.12526
17.
Banine
,
V. Y.
,
Koshelev
,
K. N.
, and
Swinkels
,
G. H. P. M.
,
2011
, “
Physical Processes in EUV Sources for Microlithography
,”
J. Phys. D: Appl. Phys.
,
44
(
25
), p.
253001
.10.1088/0022-3727/44/25/253001
18.
Suli
,
L.
,
Zhengying
,
W.
,
Jun
,
D.
,
Pei
,
W.
, and
Bingheng
,
L.
,
2017
, “
A Numerical Analysis on the Metal Droplets Impacting and Spreading Out on the Substrate
,”
Rare Met. Mater. Eng.
,
46
(
4
), pp.
893
898
.10.1016/S1875-5372(17)30118-2
19.
Manzello
,
S. L.
, and
Yang
,
J. C.
,
2002
, “
On the Collision Dynamics of a Water Droplet Containing an Additive on a Heated Solid Surface
,”
Proc. R. Soc. A: Math., Phys. Eng. Sci.
,
458
(
2026
), pp.
2417
2444
.10.1098/rspa.2002.0980
20.
Fu
,
S. P.
,
Sahu
,
R. P.
,
Diaz
,
E.
,
Robles
,
J. R.
,
Chen
,
C.
,
Rui
,
X.
,
Klie
,
R. F.
,
Yarin
,
A. L.
, and
,
J. T.
,
2016
, “
Dynamic Study of Liquid Drop Impact on Supercooled Cerium Dioxide: Anti-Icing Behavior
,”
Langmuir
,
32
(
24
), pp.
6148
6162
.10.1021/acs.langmuir.6b00847
21.
Ju
,
J.
,
Yang
,
Z.
,
Yi
,
X.
, and
Jin
,
Z.
,
2019
, “
Experimental Investigation of the Impact and Freezing Processes of a Hot Water Droplet on an Ice Surface
,”
Phys. Fluids
,
31
(
5
), p.
057107
.10.1063/1.5094691
22.
Mishchenko
,
L.
,
Hatton
,
B.
,
,
V.
,
Taylor
,
J. A.
,
Krupenkin
,
T.
, and
Aizenberg
,
J.
,
2010
, “
Design of Ice-Free Nanostructured Surfaces Based on Repulsion of Impacting Water Droplets
,”
ACS Nano
,
4
(
12
), pp.
7699
7707
.10.1021/nn102557p
23.
Antonini
,
C.
,
Innocenti
,
M.
,
Horn
,
T.
,
Marengo
,
M.
, and
Amirfazli
,
A.
,
2011
, “
Understanding the Effect of Superhydrophobic Coatings on Energy Reduction in Anti-Icing Systems
,”
Cold Regions Sci. Technol.
,
67
(
1–2
), pp.
58
67
.10.1016/j.coldregions.2011.02.006
24.
Blake
,
J.
,
Thompson
,
D.
,
Raps
,
D.
, and
Strobl
,
T.
,
2015
, “
Simulating the Freezing of Supercooled Water Droplets Impacting a Cooled Substrate
,”
AIAA J.
,
53
(
7
), pp.
1725
1739
.10.2514/1.J053391
25.
Zhou
,
Q.
,
Li
,
N.
,
Chen
,
X.
,
Xu
,
T.
,
Hui
,
S.
, and
Zhang
,
D.
,
2008
, “
Liquid Drop Impact on Solid Surface With Application to Water Drop Erosion on Turbine Blades: Part II—Axisymmetric Solution and Erosion Analysis
,”
Int. J. Mech. Sci.
,
50
(
10–11
), pp.
1543
1558
.10.1016/j.ijmecsci.2008.08.002
26.
Corrigan
,
R.
, and
DeMiglio
,
R.
,
1985
, “
Effects of Precipitation on Wind Turbine Performance
,” NASA, Washington, DC, Report No.
NASA-TM-86986
.https://ntrs.nasa.gov/citations/19850019074
27.
Hulse-Smith
,
L.
,
,
N.
, and
Chandra
,
S.
,
2005
, “
Deducing Drop Size and Impact Velocity From Circular Bloodstains
,”
J. Forensic Sci.
,
50
(
1
), pp.
1
10
.10.1520/JFS2003224
28.
Laan
,
N.
,
De Bruin
,
K. G.
,
Bartolo
,
D.
,
Josserand
,
C.
, and
Bonn
,
D.
,
2014
, “
Maximum Diameter of Impacting Liquid Droplets
,”
Phys. Rev. Appl.
,
2
(
4
), pp.
1
7
.10.1103/PhysRevApplied.2.044018
29.
Kim
,
H.
, and
Chun
,
J.
,
2001
, “
The Recoiling of Liquid Droplets Upon Collision With Solid Surfaces
,”
Phys. Fluids
,
13
(
3
), pp.
643
659
.10.1063/1.1344183
30.
Lee
,
J. B.
,
Derome
,
D.
,
Guyer
,
R.
, and
Carmeliet
,
J.
,
2016
, “
Modeling the Maximum Spreading of Liquid Droplets Impacting Wetting and Nonwetting Surfaces
,”
Langmuir
,
32
(
5
), pp.
1299
1308
.10.1021/acs.langmuir.5b04557
31.
Clanet
,
C.
,
Béguin
,
C.
,
Richard
,
D.
, and
Quéré
,
D.
,
2004
, “
Maximal Deformation of an Impacting Drop
,”
J. Fluid Mech.
,
517
, pp.
199
208
.10.1017/S0022112004000904
32.
Gao
,
X.
, and
Li
,
R.
,
2014
, “
Spread and Recoiling of Liquid Droplets Impacting Solid Surfaces
,”
AIChE J.
,
60
(
7
), pp.
2683
2691
.10.1002/aic.14440
33.
Roisman
,
I. V.
,
Rioboo
,
R.
, and
Tropea
,
C.
,
2002
, “
Normal Impact of a Liquid Drop on a Dry Surface: Model for Spreading and Receding
,”
Proc. R. Soc. A: Math., Phys. Eng. Sci.
,
458
(
2022
), pp.
1411
1430
.10.1098/rspa.2001.0923
34.
Bathel
,
B. F.
,
Stephen
,
N.
,
Johnson
,
L.
,
Ratner
,
A.
, and
Huisenga
,
M.
,
2007
, “
Prediction of Postcontact Parameters of Fluid Droplet Impact on a Smooth Surface
,”
AIAA J.
,
45
(
7
), pp.
1725
1733
.10.2514/1.24553
35.
Eggers
,
J.
,
Fontelos
,
M. A.
,
Josserand
,
C.
, and
Zaleski
,
S.
,
2010
, “
Drop Dynamics After Impact on a Solid Wall: Theory and Simulations
,”
Phys. Fluids
,
22
(
6
), p.
062101
.10.1063/1.3432498
36.
Roisman
,
I. V.
,
Berberović
,
E.
, and
Tropea
,
C.
,
2009
, “
Inertia Dominated Drop Collisions. I. On the Universal Flow in the Lamella
,”
Phys. Fluids
,
21
(
5
), p.
052103
.10.1063/1.3129282
37.
Visser
,
C. W.
,
Frommhold
,
P. E.
,
Wildeman
,
S.
,
Mettin
,
R.
,
Lohse
,
D.
, and
Sun
,
C.
,
2015
, “
Dynamics of High-Speed Micro-Drop Impact: Numerical Simulations and Experiments at Frame-to-Frame Times Below 100 Ns
,”
Soft Matter
,
11
(
9
), pp.
1708
1722
.10.1039/C4SM02474E
38.
Guo
,
Y.
,
Lian
,
Y.
, and
Sussman
,
M.
,
2016
, “
Investigation of Drop Impact on Dry and Wet Surfaces With Consideration of Surrounding Air
,”
Phys. Fluids
,
28
(
7
), p.
073303
10.1063/1.4958694
39.
Mani
,
M.
,
Mandre
,
S.
, and
Brenner
,
M. P.
,
2010
, “
Events Before Droplet Splashing on a Solid Surface
,”
J. Fluid Mech.
,
647
, pp.
163
185
.10.1017/S0022112009993594
40.
Riboux
,
G.
, and
Gordillo
,
J. M.
,
2014
, “
Experiments of Drops Impacting a Smooth Solid Surface: A Model of the Critical Impact Speed for Drop Splashing
,”
Phys. Rev. Lett.
,
113
(
2
), pp.
1
5
.10.1103/PhysRevLett.113.024507
41.
Kolinski
,
J. M.
,
Rubinstein
,
S. M.
,
Mandre
,
S.
,
Brenner
,
M. P.
,
Weitz
,
D. A.
, and
,
L.
,
2012
, “
Skating on a Film of Air: Drops Impacting on a Surface
,”
Phys. Rev. Lett.
,
108
(
7
), pp.
1
5
.10.1103/PhysRevLett.108.074503
42.
Liu
,
Y.
,
Tan
,
P.
, and
Xu
,
L.
,
2013
, “
Compressible Air Entrapment in High-Speed Drop Impacts on Solid Surfaces
,”
J. Fluid Mech.
,
716
, p.
R9
.10.1017/jfm.2012.583
43.
Wang
,
Y.
, and
Bourouiba
,
L.
,
2017
, “
Drop Impact on Small Surfaces: Thickness and Velocity Profiles of the Expanding Sheet in the Air
,”
J. Fluid Mech.
,
814
, pp.
510
534
.10.1017/jfm.2017.18
44.
Xiong
,
W.
, and
Cheng
,
P.
,
2018
, “
Numerical Investigation of Air Entrapment in a Molten Droplet Impacting and Solidifying on a Cold Smooth Substrate by 3D Lattice Boltzmann Method
,”
Int. J. Heat Mass Transfer
,
124
, pp.
1262
1274
.10.1016/j.ijheatmasstransfer.2018.04.056
45.
Langley
,
K.
,
Li
,
E. Q.
, and
Thoroddsen
,
S. T.
,
2017
, “
Impact of Ultra-Viscous Drops: Air-Film Gliding and Extreme Wetting
,”
J. Fluid Mech.
,
813
, pp.
647
666
.10.1017/jfm.2016.840
46.
Tang
,
C.
,
Qin
,
M.
,
Weng
,
X.
,
Zhang
,
X.
,
Zhang
,
P.
,
Li
,
J.
, and
Huang
,
Z.
,
2017
, “
Dynamics of Droplet Impact on Solid Surface With Different Roughness
,”
Int. J. Multiphase Flow
,
96
, pp.
56
69
.10.1016/j.ijmultiphaseflow.2017.07.002
47.
Kim
,
S. J.
,
Kim
,
J.
,
Moon
,
M. W.
,
Lee
,
K. R.
, and
Kim
,
H. Y.
,
2013
, “
Experimental Study of Drop Spreading on Textured Superhydrophilic Surfaces
,”
Phys. Fluids
,
25
(
9
), p.
092110
.10.1063/1.4821985
48.
Šikalo
,
Š.
,
Wilhelm
,
H. D.
,
Roisman
,
I. V.
,
Jakirlić
,
S.
, and
Tropea
,
C.
,
2005
, “
Dynamic Contact Angle of Spreading Droplets: Experiments and Simulations
,”
Phys. Fluids
,
17
(
6
), p.
062103
.10.1063/1.1928828
49.
,
H.
, and
Amirfazli
,
A.
,
2017
, “
Asymmetric Spreading of a Drop Upon Impact Onto a Surface
,”
Langmuir
,
33
(
23
), pp.
5957
5964
.10.1021/acs.langmuir.7b00704
50.
Izbassarov
,
D.
, and
,
M.
,
2016
, “
Effects of Viscoelasticity on Drop Impact and Spreading on a Solid Surface
,”
Phys. Rev. Fluids
,
1
(
2
), pp.
1
18
.10.1103/PhysRevFluids.1.023302
51.
Venkatesan
,
J.
, and
Ganesan
,
S.
,
2019
, “
Computational Modeling of Impinging Viscoelastic Droplets
,”
J. Non-Newtonian Fluid Mech.
,
263
(
1
), pp.
42
60
.10.1016/j.jnnfm.2018.11.001
52.
Forney
,
L. J.
,
1990
, “
Droplet Impaction on a Supersonic Wedge—Consideration of Similitude
,”
AIAA J.
,
28
(
4
), pp.
650
654
.10.2514/3.10442
53.
Josserand
,
C.
, and
Thoroddsen
,
S. T.
,
2016
, “
Drop Impact on a Solid Surface
,”
Annu. Rev. Fluid Mech.
,
48
(
1
), pp.
365
391
.10.1146/annurev-fluid-122414-034401
54.
Yarin
,
A. L.
,
2006
, “
Drop Impact Dynamics: Splashing, Spreading, Receding, Bouncing
,”
Annu. Rev. Fluid Mech.
,
38
(
1
), pp.
159
192
.10.1146/annurev.fluid.38.050304.092144
55.
Attané
,
P.
,
Girard
,
F.
, and
Morin
,
V.
,
2007
, “
An Energy Balance Approach of the Dynamics of Drop Impact on a Solid Surface
,”
Phys. Fluids
,
19
(
1
), p.
012101
.10.1063/1.2408495
56.
Rose
,
W. I.
,
1987
, “
Interaction of Aircraft and Explosive Eruption Clouds—A Volcanologist's Perspective
,”
AIAA J.
,
25
(
1
), pp.
52
58
.10.2514/3.9579
57.
Dacre
,
H. F.
,
Grant
,
A. L. M.
, and
Johnson
,
B. T.
,
2013
, “
Aircraft Observations and Model Simulations of Concentration and Particle Size Distribution in the Eyjafjallajökull Volcanic Ash Cloud
,”
Atmos. Chem. Phys.
,
13
(
3
), pp.
1277
1291
.10.5194/acp-13-1277-2013
58.
Kueppers
,
U.
,
Cimarelli
,
C.
,
Hess
,
K. U.
,
,
J.
,
,
F. B.
, and
Dingwell
,
D. B.
,
2014
, “
The Thermal Stability of Eyjafjallajökull Ash Versus Turbine Ingestion Test Sands
,”
J. Appl. Volcanol.
,
3
(
1
), pp.
1
11
.10.1186/2191-5040-3-4
59.
Song
,
W.
,
Hess
,
K.-U.
,
Damby
,
D. E.
,
,
F. B.
,
Lavallée
,
Y.
,
Cimarelli
,
C.
, and
Dingwell
,
D. B.
,
2014
, “
Fusion Characteristics of Volcanic Ash Relevant to Aviation Hazards
,”
Geophys. Res. Lett.
,
41
(
7
), pp.
2326
2333
.10.1002/2013GL059182
60.
Song
,
W.
,
Lavallee
,
Y.
,
Hess
,
K. U.
,
Kueppers
,
U.
,
Cimarelli
,
C.
, and
Dingwell
,
D. B.
,
2016
, “
Volcanic Ash Melting Under Conditions Relevant to Ash Turbine Interactions
,”
Nat. Commun.
,
7
(
1
), p.
10795
.10.1038/ncomms10795
61.
,
F. B.
,
Vasseur
,
J.
,
Aulock
,
F. W.
,
von
,
Hess
,
K.-U.
,
Scheu
,
B.
,
Lavallée
,
Y.
, and
Dingwell
,
D. B.
,
2014
, “
Nonisothermal Viscous Sintering of Volcanic Ash
,”
J. Geophys. Res.: Solid Earth
,
119
(
12
), pp.
8792
–87
50
.10.1002/2014JB011453
62.
Chen
,
W. R.
, and
Zhao
,
L. R.
,
2015
, “
Review—Volcanic Ash and Its Influence on Aircraft Engine Components
,”
Procedia Eng.
,
99
, pp.
795
803
.10.1016/j.proeng.2014.12.604
63.
Dunn
,
M. G.
,
2012
, “
Operation of Gas Turbine Engines in an Environment Contaminated With Volcanic Ash
,”
ASME J. Turbomach.
,
134
(
5
), pp.
1
18
.10.1115/1.4006236
64.
Guffanti
,
M.
,
,
T.
, and
Budding
,
K.
,
2010
, “
Encounters of Aircraft With Volcanic Ash Clouds: A Compilation of Known Incidents, 1953–2009
,”
U.S. Geol. Survey Data Ser.
,
545
(1.0), pp.
1
12
.https://pubs.usgs.gov/ds/545/DS545.pdf
65.
Bonilla
,
C.
,
2012
, “
The Effect of Particle Size and Film Cooling on Nozzle Guide Vane Deposition
,”
M.S. thesis
,
The Ohio State University
,
Columbus, OH
.10.1115/1.4007057
66.
Shinozaki
,
M.
,
Roberts
,
K. A.
,
Van De Goor
,
B.
, and
William Clyne
,
T.
,
2013
, “
Deposition of Ingested Volcanic Ash on Surfaces in the Turbine of a Small Jet Engine
,”
,
15
(
10
67.
Lawrence
,
M. J.
,
2013
, “
An Experimental Investigation of High Temperature Particle Rebound and Deposition Characteristics Applicable to Gas Turbine Fouling
,”
M.S. thesis
,
The Ohio State University
,
Columbus, OH
68.
Dean
,
J.
,
Taltavull
,
C.
, and
Clyne
,
T. W.
,
2016
, “
Influence of the Composition and Viscosity of Volcanic Ashes on Their Adhesion Within Gas Turbine Aeroengines
,”
Acta Mater.
,
109
, pp.
8
16
.10.1016/j.actamat.2016.02.011
69.
Taltavull
,
C.
,
Dean
,
J.
, and
Clyne
,
T. W.
,
2016
, “
Adhesion of Volcanic Ash Particles Under Controlled Conditions and Implications for Their Deposition in Gas Turbines
,”
,
18
(
5
), pp.
803
813
70.
Seksinsky
,
D.
,
2020
, “
Modeling Volcanic Ash Particle Impingement in Gas Turbine Engines
,” M.S. thesis,
University of Vermont
,
Burlington, VT
.
71.
Giordano
,
D.
,
Russell
,
J. K.
, and
Dingwell
,
D. B.
,
2008
, “
Viscosity of Magmatic Liquids: A Model
,”
Earth Planet. Sci. Lett.
,
271
(
1–4
), pp.
123
134
.10.1016/j.epsl.2008.03.038
72.
Schiaffino
,
S.
, and
Sonin
,
A. A.
,
1997
, “
Molten Droplet Deposition and Solidification at Low Weber Numbers
,”
Phys. Fluids
,
9
(
11
), pp.
3172
3187
.10.1063/1.869434
73.
Rioboo
,
R.
,
Marengo
,
M.
, and
Tropea
,
C.
,
2002
, “
Time Evolution of Liquid Drop Impact Onto Solid, Dry Surfaces
,”
Exp. Fluids
,
33
(
1
), pp.
112
124
.10.1007/s00348-002-0431-x
74.
Gordillo
,
L.
,
Sun
,
T. P.
, and
Cheng
,
X.
,
2018
, “
Dynamics of Drop Impact on Solid Surfaces: Evolution of Impact Force and Self-Similar Spreading
,”
J. Fluid Mech.
,
840
, pp.
190
214
.10.1017/jfm.2017.901
75.
Griebel
,
M.
, and
Klitz
,
M.
,
2017
, “
CLSVOF as a Fast and Mass-Conserving Extension of the Level-Set Method for the Simulation of Two-Phase Flow Problems
,”
Numer. Heat Transfer, Part B: Fundam.
,
71
(
1
), pp.
1
36
.10.1080/10407790.2016.1244400
76.
Sussman
,
M.
, and
Puckett
,
E. G.
,
2000
, “
A Coupled Level Set and Volume-of-Fluid Method for Computing 3D and Axisymmetric Incompressible Two-Phase Flows
,”
J. Comput. Phys.
,
162
(
2
), pp.
301
337
.10.1006/jcph.2000.6537
77.
Sun
,
D. L.
, and
Tao
,
W. Q.
,
2010
, “
A Coupled Volume-of-Fluid and Level Set (VOSET) Method for Computing Incompressible Two-Phase Flows
,”
Int. J. Heat Mass Transfer
,
53
(
4
), pp.
645
655
.10.1016/j.ijheatmasstransfer.2009.10.030
78.
Sethian
,
J. A.
,
1999
,
Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry, Fluid Mechanics, Computer Vision, and Materials Science
,
Cambridge University Press
,
Cambridge, UK
.
79.
Brackbill
,
J. U.
,
Kothe
,
D. B.
, and
Zemach
,
C.
,
1992
, “
A Continuum Method for Modeling Surface Tension
,”
J. Comput. Phys.
,
100
(
2
), pp.
335
354
.10.1016/0021-9991(92)90240-Y
80.
Croce
,
R.
,
Griebel
,
M.
, and
Schweitzer
,
M. A.
,
2009
, “
Numerical Simulation of Bubble and Droplet Deformation by a Level Set Approach With Surface Tension in Three Dimensions
,”
Int. J. Numer. Methods Fluids
,
62
(
9
), p. 993.10.1002/fld.2051
81.
Afkhami
,
S.
,
Zaleski
,
S.
, and
Bussmann
,
M.
,
2009
, “
A Mesh-Dependent Model for Applying Dynamic Contact Angles to VOF Simulations
,”
J. Comput. Phys.
,
228
(
15
), pp.
5370
5389
.10.1016/j.jcp.2009.04.027
82.
Vogel
,
A.
,
Diplas
,
S.
,
Durant
,
A. J.
,
Azar
,
A. S.
,
Sunding
,
M. F.
,
Rose
,
W. I.
,
Sytchkova
,
A.
,
,
C.
,
Krüger
,
K.
, and
Stohl
,
A.
,
2017
, “
Reference Data Set of Volcanic Ash Physicochemical and Optical Properties
,”
J. Geophys. Res.: Atmos.
,
122
(
17
), pp.
9485
9514
.10.1002/2016JD026328
83.
Wilson
,
T. M.
,
Stewart
,
C.
,
Sword-Daniels
,
V.
,
Leonard
,
G. S.
,
Johnston
,
D. M.
,
Cole
,
J. W.
,
Wardman
,
J.
,
Wilson
,
G.
, and
Barnard
,
S. T.
,
2012
, “
Volcanic Ash Impacts on Critical Infrastructure
,”
Phys. Chem. Earth
,
45–46
, pp.
5
23
.10.1016/j.pce.2011.06.006
84.
Shipley
,
S.
, and
Sarna-Wojcicki
,
A.
,
1983
, “
Distribution, Thickness, and Mass of Late Pleistocene and Holocene Tephra From Major Volcanoes in the Northwestern United States: A Preliminary Assessment of Hazards From Volcanic Ejecta to Nuclear Reactors in the Pacific Northwest
,” United States Geological Survey Miscellaneous Field Studies Map MF-1435,
Report
.https://pubs.usgs.gov/mf/1983/1435/report.pdf
85.
Li
,
Z.
,
Giese
,
R. F.
,
Wu
,
W.
,
Sheridan
,
M. F.
, and
van Oss
,
C. J.
,
1997
, “
The Surface Thermodynamic Properties of Some Volcanic Ash Colloids
,”
J. Dispersion Sci. Technol.
,
18
(
3
), pp.
223
241
.10.1080/01932699708943732
86.
Philippi
,
J.
,
Lagrée
,
P. Y.
, and
Antkowiak
,
A.
,
2016
, “
Drop Impact on a Solid Surface: Short-Time Self-Similarity
,”
J. Fluid Mech.
,
795
, pp.
96
135
.10.1017/jfm.2016.142
87.
Hicks
,
P. D.
, and
Purvis
,
R.
,
2010
, “
Air Cushioning and Bubble Entrapment in Three-Dimensional Droplet Impacts
,”
J. Fluid Mech.
,
649
, pp.
135
163
.10.1017/S0022112009994009
88.
Thoroddsen
,
S. T.
,
Etoh
,
T. G.
, and
Takehara
,
K.
,
2003
, “
Air Entrapment Under an Impacting Drop
,”
J. Fluid Mech.
,
478
, pp.
125
134
.10.1017/S0022112002003427
89.
Thoroddsen
,
S. T.
,
Etoh
,
T. G.
,
Takehara
,
K.
,
Ootsuka
,
N.
, and
Hatsuki
,
Y.
,
2005
, “
The Air Bubble Entrapped Under a Drop Impacting on a Solid Surface
,”
J. Fluid Mech.
,
545
(
1
), pp.
203
212
.10.1017/S0022112005006919
90.
,
J.
,
1976
, “
Solidification of Molten Metal Droplets Impinging on a Cold Surface
,”
Int. J. Heat Mass Transfer
,
19
(
9
), pp.
1009
1013
.10.1016/0017-9310(76)90183-6
91.
Wildeman
,
S.
,
Visser
,
C. W.
,
Sun
,
C.
, and
Lohse
,
D.
,
2016
, “
On the Spreading of Impacting Drops
,”
J. Fluid Mech.
,
805
, pp.
636
655
.10.1017/jfm.2016.584