Analytical consideration of liquid droplet impingement on solid surfaces

In industrial applications involving spray-cooling, combustion, and so on, prediction of the maximum spreading diameter of a droplet impinging on a solid surface permits a quantitative estimation of heat removal and energy consumption. However, although there are many experimental studies regarding droplet impingement behaviour, theoretical models have an applicability limit for predicting the maximum spreading diameter. In the present study, we have developed an analytical model for droplet impingement based on energy conservation that considers adhesion energy in both horizontal and vertical directions at the contact line. The theory is validated by our experiment and existing experimental data possessing a wide range of Weber numbers. We demonstrate that our model can predict βm (i.e., the maximum spreading diameter normalised in terms of initial droplet diameter) for various Newtonian liquids ranging from micro- to millimetre-sized droplets on different solid surfaces and can determine the transition between capillary and viscous regimes. Furthermore, theoretical relations for scaling laws observed by many researchers are derived.

Current knowledge of the detailed behaviour of droplet impingement derived from experimental studies has improved gradually with advancements in high-speed video technology [15][16][17][18][19][20][21][22][23][24] . Theoretical approaches, on the other hand, employ models that attempt to predict the maximum spreading diameter of droplet impingement based on energy balance, momentum balance, and empirical considerations 7,9,10,[25][26][27] . Existing models can be classified into two main categories: 1) β m as a function of Re and We (or a single parameter of Re or We) 2,17,19,22,[26][27][28] and 2) β m as a function of Re, We, and cos θ d , 9,10,25 , where θ d is the dynamic contact angle (advancing contact angle). It is only in the latter case that the interactions between the solid and liquid are considered. Due to the difficulty in predicting θ d theoretically, however, experimental data is often used instead. Although the studies found in the literature have generally reported scaling laws of β m ∝ We 1/2 in the capillary region and β m ∝ Re 1/5 in the viscous region 7 , most models have an applicability limit in the extent of the We (or Re) number for predicting the experimental data. Therefore, the quantitative prediction of β m for a wide range of We (or Re) numbers is indeed a challenging problem, especially since there are few models that can accurately predict solid surface properties. Thus, the important open question that persists is to determine the effect of different types of solids on β m as well as its theoretical prediction without the applicability limit.
In this work, we present a theoretical model derived using an energy balance approach to predict β m . Particularly, our model considers the adhesion energy at the contact line in the vertical direction in addition to the horizontal direction 29 . The derived equation can predict β m in a wide range of We (or Re) numbers for Newtonian liquid droplets on solid surfaces quantitatively without the use of arbitrary fitting parameters. We validate our model by comparing it to existing experimental data that employ micro-to millimetre-sized droplets 8,17,22,30,31 . In addition to these results, the transition point from the capillary regime to the viscous regime is theoretically determined.

Theory
The energy conservation approach 9, 10 considers both kinetic and surface energies prior to droplet impingement as well as surface energy and viscous dissipation after impingement. We now proceed to derive the theoretical equation expressing β m as a function of θ d , Re, and We. Although some empirical and semi-empirical models exist in the literature 17,22,23,26 , those models lack a quantitative prediction of β m . Recently, the importance of the work done by the adhesion force at the contact line, not only in the horizontal direction, but also in the vertical direction 29 is revealed.
From an energy conservation standpoint, the contribution of the adhesion force in the vertical component must also be considered. Let E kine , E surf , E grav , E sprd , E vis , and E def be the kinetic energy, initial surface energy, gravitational potential of the droplet, adhesion energy, viscous dissipation, and deformation energy after the impingement, respectively. Then, the following energy conservation holds:  [m] in the liquid film that is a distance from the wall, respectively. In Eq. (5), the first and second terms correspond to the adhesion energy in the horizontal and vertical directions at the contact line, respectively, and θ is the simple averaged contact angle [°] of the static and dynamic contact angles. Of course, there exist a range of differences in droplet shapes ranging from a spherical cap to a flattened sphere. A quite low We number implies a spherical cap shape after impingement (gently depositing on the solid surface), whereas a large We number implies a flattened shape. Although it is well-known that the dynamic contact angle depends on We and Re, an exact determination of the contact angle is very difficult because of the scale differences that exist such as micro-and macro-contact angles 32 .
In the present study, we assume that the exact value of the contact angle in the spreading process exists in the range between the contact angle at a quite low We value (θ lowWe ) and that at a large We value (θ highWe ). Here, with respect to θ lowWe , the fluid motion is negligible small in the quite low We situation at constant temperature. In such the case, we assume that the value of θ st can be used as θ lowWe . Here, θ st is determined by measuring the static contact angle of the droplet. Then, we consider θ highWe to be θ d at the maximum spreading diameter. Finally, the simple averaged values of θ st and θ d are used to give θ = (θ st + θ d )/2. Moreover, in the deformation term of Eq. (7), the exact evaluation of the surface area is also very difficult. Therefore, the deformed surface (S def ) is defined as the harmonic average of the droplet surface of the spherical cap (S cap ) and of the disk (S disc ), given as  In most energy conservation approaches 9, 10 , the initial impinging velocity (u) is used to evaluate the viscous dissipation. However, since shear stress occurs in the liquid film that spreads along the solid surface, the radial liquid velocity along the solid surface (u R ) to evaluate the viscous dissipation term is more appropriate. When the liquid velocity reaches zero (i.e., kinetic energy is zero), the droplet diameter realizes its maximum spreading diameter. The treatment of the dissipation term is very difficult in this kind of analytical approach because the exact velocity distribution or profile is not known during the spreading process. However, it may be important to postulate a velocity profile for the evaluation of the viscous dissipation term. When the droplet impinges on the solid surface, the droplet shape initially becomes a bell shape owing to a recoil force from the solid surface, and then reaches the maximum spreading diameter caused by the surface tension 33 . This surface tension acts on the top of the bell-shaped droplet and pushes the liquid toward the solid surface, which then generates the radial liquid flow. This situation is like a wall jet along a solid surface [34][35][36] . Based on this assumption, the dissipation term can be evaluated using the velocity profile of wall jet.
The wall jet type velocity profile is non-linear and the peak of the velocity (i.e., maximum velocity) is located near the wall. From an experimental study of the wall jet 35 , it was found that the velocity peak is located at around one quarter of the effective height of the wall jet flow. However, in the case of the droplet impingement, the height is restricted by the droplet volume. If the liquid flow is confined by the wall and the liquid film surface (i.e., between parallel plates), however, the velocity profile becomes a parabolic shape where the peak of the velocity is located at half of the height. Thus, we postulate that the peak of the velocity profile in the case of the droplet impingement would be somewhere between these two situations of the wall jet flow and the parallel plate flow. Consequently, the effective height in Eq. (6) can be obtained by taking the harmonic average as follows: Therefore, in the present model, the maximum velocity is characterized at the effective height of Eq. (11). As a next step, to evaluate u R in Eq. (6) we need to calculate the initial radial velocity. However, the exact calculation is very difficult because the droplet shape is very complicated, as mentioned above. Thus, we estimate the initial radial-mean velocity u R 0 , which is defined by assuming a cylindrically-shaped droplet of diameter d 0 on the solid surface before spreading, as shown in Fig. 1. This assumption is used only for the analytical evaluation of the initial radial-mean velocity of the liquid. The droplet volume and velocity before impingement are denoted as V 0 and u, respectively. At the moment of impingement on the solid surface, the liquid flows out from the cylindrical surface with an initial velocity of u R 0 . Here, an equivalent height l of the cylindrical droplet is calculated as Because of mass conservation, Combining Eqs (12) and (13) yields the following relation: R 0 Because it is difficult to obtain a detailed velocity profile in the liquid, if the relationship between the maximum velocity u max and the mean velocity u mean satisfies u max ≈ u 2 R mean , then u R max ,0 as characterised by the effective height of Eq. (11) becomes Scientific RepoRts | 7: 2362 | DOI:10.1038/s41598-017-02450-4 As mentioned before, because the liquid velocity is zero when the droplet diameter reaches the maximum spreading diameter, the liquid velocity decreases from u R max ,0 to zero. For the sake of simplicity, the motion of the liquid film can also be characterised by a velocity that changes from u R max ,0 to zero. If the radial velocity in the radial direction of the spreading liquid film is known, the radial-mean velocity can be evaluated by where, r represents the spreading radius of liquid film, and u r ( ) R max represents the maximum velocity that changes from u R max ,0 at r = r 0 to zero at r = r max . Again, because it is difficult to determine the exact function of u r ( ) R max in the case of droplet impingement, we assume that u r ( ) R max linearly decays with respect to the spreading radius as (r m − r)u R max ,0 /(r m − r 0 ). Eventually, u R of Eq. (16) can be calculated as 3 u/8 in the present study. Then, in the viscous dissipation term (Eq. (6)), t m is the time when the kinetic energy E kine is completely converted into the adhesion energy E sprd , the viscous dissipation E vis and so on. Therefore, t m is given by r m /u. The viscous dissipation term can be calculated by substituting t m = r m /u, Eq. (11) and the relation of u R = 3 u/8 into Eq. (6) as follows: Finally, by substituting Eqs (2-5), (7) and (17) Figure 1. Simplified model to evaluate the radial velocity u R in the viscous dissipation term. In the schematic, u, V 0 , d 0 , u R 0 and l represent the initial impinging velocity, droplet volume, droplet diameter, initial radial velocity and initial radial height.
In Eq. (18), the first, second, and fifth terms represent the non-dimensional kinetic energy ⁎ E ( ) kine , the viscous dissipation ⁎ E ( ) vis , and the gravitational potential ⁎ E ( ) grav respectively, while the third and fourth terms combined represent the adhesion energy ⁎ E ( ) sprd . The sixth and seventh terms represent the initial surface energy * E ( ) surf and the surface energy of the deformed surface ⁎ E ( ) def , respectively. By using the definition of the Ohnesorge number (Oh = μ l /(ρ l d 0 σ lg ) 1/2 = We 1/2 Re −1 ), Eq. (18) can be solved for We 1/2 . From Eq. (18), we can derive two limiting solutions for the capillary and viscous regions. In the capillary region, the viscous dissipation is negligible. Thus, we obtain the following relation: From this relation, we obtain a scaling law of β m ∝ We 1/2 . In the viscous region, on the other hand, the kinetic energy and the viscous dissipation dominate, which leads to the following relation derived from Eq. (18):

Methods
In this work, we used two liquids -purified water (Wako Pure Chemical Industries, Ltd., Osaka, Japan) and pure ethanol (99.5% pure, Kenei Pharmaceutical Co. Ltd., Osaka, Japan) in order to understand droplet impingement for high-and low-surface energy liquids, whose densities (ρ), dynamic viscosities (μ), and surface tensions (σ lg ) are given as follows: ρ w = 998.   respectively. We released 1.1-µL droplets using a microsyringe from ten different heights, z = 1.5 to 700 mm. The droplet can be regarded as a free-falling object in all experiments. A high-speed video camera (HX-5, NAC image technology, Ltd., Japan) with a microscope (Leica Microsystems, Wetzlar, Germany), captured images of the droplet behaviour after striking the solid surfaces; the frame rate is 20,000 fps. We measured the impinging velocity u, droplet diameter d 0 , and the maximum spreading diameter d max based on the captured images. The values of We in our experiments ranged from 0.3 to 230. For high We number conditions, we used existing experimental data 8,17,22,30,31 to verify the present model. All experiments were performed three times and conducted within temperature and humidity ranges of 21.0-23.0 °C and 51.0-55.0%, respectively.  22)) and viscous (Eq. (23)) regimes obtained by the present model. Here, V 0 is the droplet volume and θ is the averaged contact angle of θ st and θ d . Red circles (⚬) represent the experimental data obtained in this study. White diamond (⬦) and inverted triangle (▿) markers represent water on an Al substrate and a superhydrophobic surface, respectively 17 . Black circle (⦁) and triangle (▴) markers represent existing experimental results for high-viscosity liquids (silicone oils) of μ = 20 mPa s and 300 mPa s, respectively 17 . Figure 2 shows images of the water droplet impinging on the SR substrate from varying heights of z = 10, 100 and 700 mm, where the time advances from the left image to the right image. As the time advances, the droplet deforms and the droplet diameter reaches its maximum spreading state. From these images, it can be seen that the droplet shape initially becomes bell shaped and then reaches a disk shape at its maximum spreading diameter. The top surface of the bell-shaped droplet pushes the liquid down into the droplet.  Figure 3-b1,b2 and b3 illustrate the relationships between β m and Re for the same combinations of liquids and solids shown in Fig. 3-a1,a2 and a3. In each figure, the solid red lines represent our model given by Eq. (18), while the solid black and blue lines correspond to existing models developed by Rosiman 27 , β m = 0.87Re 1/5 − 0.4Re 2/5 We −1/2 , and Pasandideh-Fard et al. 10 , β m = (We + 12) 1/2 /(3(1 − cosθ d ) + 4(WeRe −1/2 )) 1/2 , respectively. We used the value of θ d measured in our study for the model given by Pasandideh-Fard et al. 10 . The red circles denote experimental data that we collected in this study, while the white inverted triangles and diamond markers signify existing experimental data 17,22 . Our model provides a better fit for the experimental data in each figure as compared to the two existing models, especially at low We numbers, where the two existing models deviate significantly from the experimental data. This result indicates the importance of considering adhesion energy in the vertical direction, in addition to the horizontal, in the capillary region. In each figure, the dashed blue and black lines are the limiting solutions corresponding to the capillary and viscous regions, respectively. Existing experimental data in Fig. 3-b1,b2 and b3 for high viscous liquids (black circles and triangle markers) 17 are in good agreement with the limiting solutions in the viscous regime (Eq. (23)). Figure 4 depicts the transition between different droplet impingement conditions by taking a water droplet on SR as an example. Figure 4- Fig. 4-a2, which is the capillary region, ⁎ E kine is comparable to ⁎ E sprd . Between lines We [-] Non-dimensional energy  A and B, which is the capillary-to-viscous region, the effect of viscous dissipation gradually appears until line B, at which point ⁎ E vis exceeds ⁎ E sprd in the viscous region. While the scaling law of β m ∝ We 1/4 17 may potentially correspond to the intermediate region between lines A and B, we did not observe a scaling law of β m ∝ We 1/4 in our theoretical model. From Fig. 4-a2, the droplet condition of ⁎ E vis > ⁎ E sprd signals the onset of the viscous region. Thus, we can define the ratio of E r = ⁎ ⁎ E E / vis sprd (Eq. (A8) in Supplementary Information) in order to pinpoint the onset of the viscous region. The reason of the choice of ⁎ def is that the spreading process is directly affected by the wetting behaviour. When E r > 1, the droplet falls under the viscous region. Figure 4-b1 presents the relationship between E r and We. The intersection of E r = 1 and line B represents the transition point T. Figure 4-b2 displays the relationships between E r and impact numbers P = We/Re 4/5 and P = We/Re 2/5 7, 17, 37 . Two types of impact numbers are displayed in the figure. Although there is no consensus with respect to P, P = 1 is generally accepted as the boundary separating the capillary and viscous regions. However, our model indicates that both values of P cannot predict the transition point T. Strictly speaking, the point Q1, which is evaluated by P = We/Re 4/5 , is in the viscous region after the point T, whereas the point Q2, which is evaluated by P = We/Re 2/5 , is in the intermediate region between lines A and B. Here, the values of P = We/Re 4/5 and We/Re 2/5 are 0.2 and 6.5, respectively, when E r = 1. Figure 5-a1 and a2 compare our model with existing experimental data 19,22 . vapour layer is thick and prevents the liquid from touching the surface. However, this vapour layer has a non-uniform thickness, being thicker at the centre and thinner at the perimeter of the contact area 19,38 . Thus, the liquid may come in contact with the surface of the solid at this perimeter, and the interaction between the solid and liquid at this perimeter would affect the spreading behaviour of the droplet. Our model captures the trend set by the experimental data, although, strictly speaking, we should also account for the temperature dependency of the physical properties. This result indicates that our model has the potential to predict and understand droplet impingement behaviour, including thermal effects. In this case, we used the same value of θ used for the case of a water droplet on SR, owing to the lack of available contact angle data in literature 19 . In the case for the model given by Pasandideh-Fard et al. 9 , the value of θ d = 128.2° is used. We validated our model against existing experimental data for micro-sized droplet impingement on solid surfaces, as shown in Fig. 6. In this figure, both the white and black circles represent the existing experimental data from Visser et al. for micro-sized water droplets, d 0 = 48 μm and 50 μm, respectively, on a glass plate 8,30 . In the We <100 region, we used data reported in ref. 30, and in the We >100 region, we used data reported in ref. 8 because the accuracy of the data in the We >100 region in ref. 30 is reportedly poor. The solid red, blue, and black lines represent the theoretical results evaluated by our model, by Pasandideh-Fard et al. 10 , and by Rosiman 27 , respectively. From the image presented in ref. 30, we estimated the value of θ d to be approximately 90°. The static contact angle θ st is estimated as 79.2° using the ratio of θ st /θ d = 0.88 which is averaged value of the water droplets on SR and PC in our experiment. Thus, the value of θ is estimated as 84.6°. Although the result by Roisman shows good agreement with the experimental data in the high We number region, the results demonstrate that our model shows fairly good agreement with the experimental data from the low We number to the high We number region.
The white and black triangles represent the existing experimental data 31 for water droplets, d 0 = 213 μm and 618 μm, respectively, on a coal surface. The corresponding static contact angle is reported to be 57°. The dynamic contact angle is estimated as 64.8° using the ratio of θ st /θ d = 0.88. Thus, the value of θ is estimated as 60.9°, which we have adopted in our model (red and green dashed line). Our model was also successful in predicting the experimental data for droplets impinging on a coal surface.

Conclusion
In this article, we have presented experimental and theoretical considerations of droplet impingement on solid surfaces. Our model accurately predicts the impinging behaviour of several kinds of liquids on various solid surfaces. According to the equations that we derived based on our theoretical considerations, β m observes a scaling law of β m ∝ We 1/2 (∝Re) in the capillary region and β m ∝ We 1/10 (∝Re 1/5 ) in the viscous region. In addition, the contribution of each energy component to the variation of β m indicates that impact numbers, such as P = We/ Re 4/5 and P = We/Re 2/5 , cannot predict the transition point between the capillary and viscous regimes. Instead, our theory proposed using E r instead of P to identify the boundary between the capillary and viscous regions. The present work, however, mainly considers a large-density ratio situation where the effect of the surrounding gas is negligible. To confirm the validity of our model, the low-density ratio situation must be considered because some vortex motions of the surrounding fluid may affect the impinging behaviour. In addition, we did not address non-Newtonian fluids in this study. To do so, the viscous term would have to be reconsidered to reflect the different expressions for the shear stress as compared to that of a Newtonian fluid. Nevertheless, the theoretical model that we presented in this paper has the potential of becoming a powerful tool to analyse droplet impingement behaviour. In particular, for inkjet droplets in precision engineering applications, such as soldering of electronics and microarrays for semiconductor components, our model can guide the development and precise fabrication of nano-and microstructures used in high-performance systems and devices.