## 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 [1–5], thermal spray coating [6–8], 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 [20–24], 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 [29–34], rim formation [35–37], splash formation [38–40], air entrapment during droplet impact [41–45], effect of surface contact angle and surface roughness on droplet spreading [46–48], 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 (SiO_{2}), along with a range of other minerals [58–61]. 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 [62–64]. 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=\rho dV02/\sigma $, the Ohnesorge number $Oh=\mu /\rho \sigma d$, and the equilibrium (static) contact angle $\theta e$. Here, $\rho $ and $\mu $ are the droplet density and viscosity, $\sigma $ 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=\rho dV0/\mu $ can be related to the above parameters by $Re=We/Oh$. The capillary number $Ca=\mu V0/\sigma $ 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=\rho gd2/\sigma \u226a1$). 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 $\tau c$ of spreading as $\tau cV0/d=1/Re$.

Much of the existing literature for droplet impingement has been developed for relatively high droplet Reynolds number, of order 10^{2}–10^{4}, 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 10^{4}. 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

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

**n**and mean curvature $\kappa $ are defined by

*a*as [80]

**n**at interface cells adjacent to the wall is set as

where $nW$ and $tW$ are unit vectors normal and tangent to the wall, respectively, and $\theta 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 *L _{e}* 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*=

*5*

*d*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.

*d*, the impact velocity $V0$, and the liquid density $\rho \u2113$ as

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\xd710\u22126$. 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 *L _{e}*

*=*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]).

Property | Symbol | Units | Range of values | Nominal value |
---|---|---|---|---|

Diameter^{a} | $d0$ | μm | 3–30 | 10 |

Impact velocity^{b} | $V0$ | m/s | 130–450 | 200 |

Viscosity^{c} | μ | Pa·s | 10–10^{5} | 25 |

Density^{d} | ρ | kg/m^{3} | 2350–2450 | 2400 |

Surface tension^{e} | σ | N/m | 0.35–0.37 | 0.36 |

Equilibrium contact angle^{f} | $\theta e$ | deg | 18–71 | 55 |

Property | Symbol | Units | Range of values | Nominal value |
---|---|---|---|---|

Diameter^{a} | $d0$ | μm | 3–30 | 10 |

Impact velocity^{b} | $V0$ | m/s | 130–450 | 200 |

Viscosity^{c} | μ | Pa·s | 10–10^{5} | 25 |

Density^{d} | ρ | kg/m^{3} | 2350–2450 | 2400 |

Surface tension^{e} | σ | N/m | 0.35–0.37 | 0.36 |

Equilibrium contact angle^{f} | $\theta e$ | deg | 18–71 | 55 |

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.

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

Viscosity varies from 10–10^{5} 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.

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

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

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

### 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 ($\varphi =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 3b–3e 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.

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(\varphi )$ 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 [41–43,45,87–90]. Experimentally, studies indicate that this air film beneath impacting droplets is sometimes very thin, only a few tens of nanometers thick [41].

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(\varphi )(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.

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\u2264Re\u226410$, 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.

*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

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.

*d*as $aGL=AGL/d2$ and $rc=Rc/d$, the potential energy can be written in dimensionless form as

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.

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.

where $\Pi h$ is the integral over the droplet of the local energy dissipation rate $\pi 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.

### 3.3 Velocity Profile Within the Droplet.

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\u2264Re\u22641000$.

*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*

where $\Phi (z)=1\u22120.5\u2009erfc[lnz\u2212\mu LN/2\u2009\sigma 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 $\Delta CLN=0.008$, $\Delta \mu LN=0.004$, and $\Delta \sigma 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$, $\mu LN$, and $\sigma LN$ are plotted in Fig. 16 as functions of Reynolds number. The data indicate that $\mu LN$ and $\sigma LN$ are approximately constant while $CLN$ decreases significantly in time. We do not see any clear trend in values of $\mu LN$ and $\sigma LN$ with Reynolds number, but rather the values are scattered around mean values of approximately −1.1 and 0.85, respectively.

### 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.

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.

## 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\u2264Re\u226410$). 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\u2264Re\u22647$, 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)

*a*=_{GL}air–liquid interfacial surface area (dimensionless)

*A*=_{GL}air–liquid interfacial surface area (m

^{2})- Bo =
bond number (dimensionless)

*c*=coefficient in velocity profile (1/m-s)

- $CLN$ =
log-normal coefficient (dimensionless)

*d*=droplet diameter (m)

*e*=_{c}total kinetic energy (dimensionless)

*e*=_{p}total potential energy (dimensionless)

*E*=_{c}total kinetic energy (J)

*E*=_{p}total potential energy (J)

*f*=log-normal distribution function (dimensionless)

**F**=_{sf}surface tension force (N)

*g*=gravity (m/s

^{2})*h*=height of top of droplet (dimensionless)

*H*=smoothed Heaviside function (dimensionless)

*L*=_{e}smallest mesh size (dimensionless)

**n**=unit normal (dimensionless)

- Oh =
Ohnesorge number (dimensionless)

*p*=pressure field (Pa)

*r*=radial coordinate (m)

*r*=_{c}contact radius (dimensionless)

*R*=_{c}contact radius (m)

- Re =
Reynolds number (dimensionless)

*s*=straining rate (dimensionless)

*t*=time (s)

*u*=radial velocity component (dimensionless)

**u**=velocity vector (dimensionless)

*U*=_{c}spreading velocity (m/s)

*V*=_{0}droplet impact speed (m/s)

*w*=axial velocity component (dimensionless)

- We =
Weber number (dimensionless)

*z*=axial coordinate (m)

*z*=_{m}center of mass (dimensionless)

### Greek Symbols

*δ*=Dirac delta function (dimensionless)

*θ*=polar coordinate (radians)

*θ*=_{e}static contact angle (radians)

*κ*=mean curvature (m

^{−1})*μ*=fluid viscosity (Pa·s)

- $\mu g$ =
gas phase viscosity (Pa·s)

- $\mu \u2113$ =
liquid phase viscosity (Pa·s)

- $\mu LN$ =
log-normal coefficient (dimensionless)

*π*=_{h}dimensionless dissipation rate (dimensionless)

*ρ*=fluid density (kg/m

^{3})- $\rho g$ =
liquid phase density (kg/m

^{3}) - $\rho \u2113$ =
liquid phase density (kg/m

^{3}) *σ*=liquid–gas surface tension (N/m)

- $\sigma LN$ =
log-normal coefficient (dimensionless)

- $\sigma SG$ =
solid–gas surface tension (N/m)

- $\sigma SL$ =
liquid–solid surface tension (N/m)

*τ*=_{c}spreading time (

*μ*s)*ϕ*=level-set variable (m)

- Ф =
cumulative normal distribution function (dimensionless)