Nomenclature
h,	film thickness;	n,	direction vector;
D,	drop diameter;	a,	1.5 times grid spacing;
U0,	impact velocity;	fi,	difference of position of crown between data set i and case of infinite Re;
ρ,	density;
σ,	surface tension coefficient;	||fi||,	2-norm of function fi;
μ,	dynamic viscosity;	Dev,	function indicating global dissipation;
V,	velocity;	Vr,	radial velocity;
p,	pressure;	Vz,	vertical velocity;
Fst,	surface tension;	Φ,	dissipation function;
φ,	level-set (LS) function;	v,	kinematic viscosity;
(r, z),	radial and vertical coordinates;	V,	volume;
|d|,	distance from interface;	ev,	dissipation rate;
α,	liquid volume fraction;	E0,	initial kinetic energy;
Ω,	fluid volume;	H,	film thickness, H = h/D;
∂Ω,	position of interface;	Re,	impact Reynolds number, Re = ρlDU0/μl
κ	curvature of surface
δs	numerical thickness of surface
We	impact Weber number, We = ρlDU02/σ;	Rc	crown radius scaled by D;
K	combined non-dimensional number, K = We · Oh-0.4 = (Re · We2)0.4;	λ	initial loss factor;
		Re*	crown Reynolds number;
T	non-dimensional time, T = tU0/D;	Rec	critical impact Reynolds number with respect to threshold c.
T0	time shifting;

Subsripts

i, material i (l for liquid, w for water, g for gas).

Supersripts

T, tensor transposition.

1 Introduction

Drop impact is a pervasive problem in natural phenomena and industrial production, such as raindrop erosion, fuel atomization, ink-jet printing, corrosion of turbine blades, and aircraft icing. An in-depth study of the mechanism helps us to gain a better understanding of the relevant physical processes, with important scientific and applied value^[1].

There are several factors where a drop impacts on a thin planar film of the same Newtonian liquid^[1-2]. Through dimensional analysis, we get three major independent non-dimensional numbers, with the effect of air and gravity neglected: the non-dimensional thickness of the liquid film H=h/D, the impact Reynolds number and the impact Weber number The impact velocity is U₀, the diameter of the drop is D, the thickness of the liquid film is h, the density of the liquid is ρ_l, the viscosity of the liquid is μ_l, and the surface tension is σ. The effect of air is negligible, since the density ratio of liquid to air is much larger than one. In addition, gravity is neglected because of the large Froude number and the small Bond number over an extremely short time^[1]. Other mechanisms are not considered here, such as surface property of the wall^[3] (possibly important when the film is extremely thin), heat effect when there exists temperature gradient^[4], non-Newtonian effect^[5].

In a high-speed impact, a crown-like sheet is a typical phenomenon^[1], as shown in Fig. 1. The end rim of the initial jetting becomes unstable and develops into fingers of liquid. Then, because of a Rayleigh-Plateau-like instability (with a Rayleigh-Taylor type involved)^[6], the top rims of the fingers break up into secondary droplets. This whole process is called splashing (strictly speaking, crown splashing, considering other types of splashing^[7]). The case with no secondary droplet but only jetting in a moderate impact is a crown or a spreading, whereas the case with no jetting when the impact speed is low enough is called deposition. The thresholds between these cases are believed to be controlled by ^[7-10]. We are not concerned here with the secondary droplets or the involving mechanism, but rather with the spreading motion of the crown.

The motion of the crown is a priority concern of researchers. The most important result was given by Yarin and Weiss^[8], who treated the crown as a kinematic discontinuity in a shallow water approach. In their theory, the governing equation of the quasi-one-dimension flow is of a Burgers type, which leads to a square-root law of the motion of the crown (all variables are non-dimensionalized with D and U₀):

(1)

in which R_c is the radius of the crown, H is the thickness of the film, T is the time, and T₀ is the shift time to be determined empirically. This formula generally over-predicts the experimental data by about 15%, which is believed to be the result of neglecting the effects of viscosity^[1].

There have been many researches on the issue of viscosity. Roisman and Tropea^[11] used the solution of Stokes's first problem as an estimation, and drew a conclusion that viscosity has no significant effect on the flow. Trujillo and Lee^[12] extended the work of Yarin and Weiss's theory to viscous flow and found that the effect of viscosity on the movement of the crown can be ignored. Josserand and Zaleski^[13] carried out numerical simulations using the volume of fluid (VOF) method and found that the position of the jetting depends approximately on the square root of time for different Re numbers, which indicates that viscous effects do not change the basic form of (1). Some scholars^{[9, 14]} found experimentally that viscosity has little effect on the motion of the crown under common experimental conditions. Though, many studies view the liquid viscosity as having little effect on this problem, all are limited to low or moderate viscosity (or high Re).

Recently, Gao and Li^[15] carried out detailed experiments and observed a visible influence of viscosity of the liquid on the spreading of the crown. They introduced a viscosity-related empirical parameter λ (named the initial loss factor) to reduce the difference between the experimental data and Yarin and Weiss's theory. The formula they deduced reads as follows:

(2)

with λ =0.26Re^0.05We^-0.07H^-0.34. The loss factor is around 0.5 in most situations, which suggests a loss of kinetic energy of about 75% in the initial stage of the impact (with the surface energy neglected).

The initial loss is related to the phenomena at earlier time during the formation of the crown. An initial jetting (called an ejecta sheet) at the neck of the drop is discovered at much earlier time during the impact than that of the crown by Weiss and Yarin^[16] and Thoroddsen^[17] based on boundary-integral simulations and experimental observations, respectively. The physics of ejecta sheet is not totally clear as yet, despite many studies on this issue^{[12, 17-21]}. Specifically, Deegan et al.^[7] have stated that the common splashing crown is a combination of an ejecta sheet and a so-called Peregrine sheet, which can be treated as a one-dimensional discontinuity in Yarin and Weiss's concept. The way of the combination is still under discussion because of the complex evolution of the ejecta sheet. Notably, Zhang et al.^[22] have pointed out that the ejecta sheet and Peregrine sheet form a single continuous structure for low Re number based on experimental observations.

It should be noted that many of the aforementioned studies focus on the motion of crown and discuss the effect of viscosity on it, regardless of the inner mechanism of the flow. Answers to questions such as ``why does or does not the liquid viscosity influence the crown's evolution" and ``where and how does the impact energy dissipate" remain unclear. Thus, the present study aims to bridge the gap between theories and experiments on this issue and discuss the influence of viscosity. In particular, the mechanism of viscous dissipation in the drop impact process is considered.

2 Model, method and validation 2.1 Coupled level-set volume of fluid (CLSVOF) method

To tackle two-phase flows with surface evolution, the LS method and the VOF method are two of the most common ways among other numerical methods^[2], such as the smoothed particle hydrodynamics (SPH) method^[23], the lattice Boltzmann method (LBM)^[24], and the boundary integral method (BIM)^[16] for potential flow. The CLSVOF method is the coupled version of these two separate methods.

In the LS method, a continuous phase function is set in order to characterize the phase distributions. The interface is defined clearly by a smooth LS function, which allows the surface tension to be calculated more accurately. However, solving the LS function does not guarantee its conservation if a re-initialization procedure is inevitable, which often results in loss of mass. The VOF method is conservative. One uses the volume fraction of the space occupied by one phase of the whole, but with a lower accuracy in calculation of the interface. This drawback leads to an interface reconstruction procedure, which can cause errors in calculating the surface tension. Sussman and Puckett^[25] successfully achieved the coupling of the two methods. They used the continuous LS function to calculate the normal of the interface in the reconstruction procedure, which ensures the accuracy of the interface. They then re-initialize the LS function by setting it as a zero iso-surface at the interface to reduce the non-conservation.

We use ``one fluid approach"^[26], considering a single fluid with variable viscosity and density, which yields the incompressible Navier-Stokes equation as follows:

(3)

(4)

(5)

where F_st is the surface tension term.

For interface tracking, we define the VOF function as the liquid volume faction, and the LS function as follows:

(6)

where |d| is the distance from the interface, Ω_i is the volume occupied by phase i, and ∂Ω is the position of the interface.

The global density and viscosity can then be defined as follows: , respectively. Substituting these into (3), we get

(7)

The LS function satisfies

(8)

The LS function requires re-initialization if the interface undergoes severe stretching and tearing. We implement the piecewise linear interface construction (PLIC) method here at every time step.

The surface tension is estimated using the continuum surface force (CSF) model^[27]:

(9)

(10)

(11)

To reduce the strengthening effect of the high ratio of density and viscosity on the numerical instability at the surface, we let

(12)

where a is equal to 1.5 times the grid spacing.

Those equations are solved by the finite volume method in two-dimensional axisymmetric geometry with the commercial software ANSYS FLUENT. The pressure and velocity are decoupled by the semi-implicit method for pressure-linked equations (SIMPLE) method. The domain is 20 mm×15 mm with a set of meshes shown in Fig. 2. The meshes are adaptive within four levels, depending on the gradient of fluid density. A set of meshes for boundary layer is set near the solid wall. The pressure term is discretized using pressure staggering option (PRESTO) scheme. The convection term and LS function are discretized using a second order upwind scheme, and the diffusion term is discretized by least-squares-cell-based interpolation. A first order implicit scheme is used for the time integration. The minimum time step is 2.5× 10^-8 s, which allows the convergence condition to be satisfied.

The physical parameters we used in the simulation are the same as those in the experiments of Levin and Hobbs^[28], namely an impact velocity of 4.8 m/s, a drop diameter of 2.9 mm, a film thickness of 0.5 mm. The density of the liquid phase is 1 000 kg/m³ and the viscosity is 0.001 kg/(m · s). The gaseous phase is set to be air, the density is 1.225 kg/m³ and the viscosity is 1.79 × 10^-5 kg/(m· s). The surface tension is 0.0723 N/m. We also use different phase labels to distinguish the drop and the film.

2.2 Validation

We first validate our numerical method based on the experiments of Levin and Hobbs^[28] (marked as `LH exp.' in Fig. 3) and Cossali et al.<sup>[9]</sup> (marked as `C exp.' in Fig. 4), as shown in Fig. 3 and Fig. 4, respectively. We compare our numerical results (marked as `N.D.' in both figures) with the experimental data, with the formula of Yarin and Weiss<sup>[8]</sup> (`YW') and Gao and Li^[15] (`GL'), and with the theory of Trujillo and Lee<sup>[12]</sup> (`TL').

In both cases, our numerical results are in overall good agreement with the experimental data. In contrast, the theoretical predictions differ from the experimental data more or less, from which we can conclude that our numerical results are credible.

3 Results

Six cases are used to discuss the influence of the liquid viscosity on the evolution of the crown. Only the viscosity of the liquid is changed; all other physical parameters are held fixed. The viscosities normalized by the viscosity of water are set as 0 (inviscid), 0.02, 1, 20, 50, 100, with Re numbers infinite, 6.96× 10⁵, 1.39× 10⁴, 696, 278, and 139, respectively. The We number is We=924, the film thickness is H≈ 0.17, and the non-dimensional factor K is of order 10³.

The evolution of the crown is shown in Fig. 5. Two different colors are used to distinguish between the drop material and the film. The wall of the crown is mainly composed of the material from the film, whereas the drop material is concentrated mostly on the inner wall. There is no apparent mixing between them, as was observed by Thoroddsen^[17] and Josserand and Zaleski^[13]. As can be seen from the evolution, after the ejecta-sheet phase, the drop material beneath the static liquid turns to expand horizontally, together with the agitated film discharge into the downstream zone. Thus, an obvious crown-like sheet takes shape at the crossing part. After the momentum of the drop has been transferred totally, the crown spreads by its inertia. A higher liquid viscosity suppresses the splashing; the height of the crown is lower, the wall is inclined more outwardly, and the jetting is slower. Rioboo et al.^[10] pointed out that a crown appears for K>400 and that splashing arises for K>2 100, given 0.08 < H < 0.14, which is consistent with our numerical results. Moreover, the ejecta sheet becomes the leading edge of the Peregrine sheet for small Re number, which is the same as the situation asserted by Deegan et al.^[7] and Zhang et al.^[22].

Note that the droplets displaced in the simulation are in fact spherical rings. Moreover, the assumption of axisymmetric flow is not valid after a certain time if splashing occurs, which was pointed out by Zhang et al.^[29] from their experiment. We do not concentrate on small-scale phenomena concerning splashing droplets, but rather on the motion of the crown. This is believed to be unaffected by the droplets emitted from the top rim, in which a two-dimensional geometry is a reasonable assumption.

To study the influence of viscosity on the motion of the crown, we take its intersection with a plane at a height of 20% above the film (the capillary waves are within 20% of the height of the film^[8]) and use the liquid-gas interface as the outer edge of the crown. In this way, we are able to capture the position of the crown automatically, as shown in Fig. 1. The curves of the motion of the crown are plotted in Fig. 6, which shows the position of the crown as a function of time for different viscosities. The spreading is insensitive to the variation in viscosity when the viscosity is low, but a distinct reduction of the spreading speed of the crown is observed at high viscosity. For T < 0.4, the viscosity has no visible influence on the position of the sheet, whereas the curves for different viscosities separate from each other at later time.

4 Discussion

From an overall perspective, the distinction among the curves in Fig. 6 can be a measurement of the effect of viscosity. A function f_i is defined by calculating the difference between each data set and the inviscid case. Rescaling its 2-norm with the case of water, we have a function D_ev(Re)=||f_i||/||f_w|| representing the effect of viscosity. A fitted curve D_ev(Re)=1 217× Re^-0.55 is gotten with R²=0.99, as shown in Fig. 7. D_ev increases significantly as Re decreases, especially for the cases of low Re (e.g., Re < 1 000), although no obvious turning point could be found in this curve.

To estimate the effect of viscosity the crown takes shape, Trujillo and Lee^[12] proposed a model based on the shallow water approximation. Their theory shows that is a metric of the ratio of inertia to the viscosity upstream of the crown; specifically, Re^*≈ 8 in the case of the highest viscosity in our simulation (Re=139, H=0.17). There is only a slight difference between their result and the formula of Yarin and Weiss, except for the case of an extremely thin film (H~0.01). This means that the effect of viscosity on the motion of the crown is not obvious unless the film is extremely thin. In their theory, the initial velocity at upstream of the crown is considered uniform and equals in magnitude to the impact velocity. However, this assumption is less tenable at high viscosity.

Figure 8 shows the distributions of the radial velocity near the crown for different viscosities. The radial velocity increases linearly with the distance from the impact point, which is an example of a stagnation flow. After reaching its peak, the velocity diminishes from the inner wall of the crown. The peak velocity is the input upstream velocity of the crown. At T=1, the input velocity is about 75% of the impact velocity in magnitude in those cases in which the viscosity is lower than that of water. This velocity is significantly reduced for higher viscosities, e.g., 45% of the impact velocity in the case of 100 times the viscosity of water. Given that the upstream velocity varies with viscosity, different initial conditions should be proposed in shallow water flow.

We next consider the loss from an energy perspective, which has the potential to explain the difference in initial conditions caused by the viscous effect. The dissipation function of the liquid is

(13)

which represents the rate at which work is converted into heat by the liquid viscosity. The distribution of this function is shown in Fig. 9. It can be seen that there is a region of large dissipation at the neck and beneath the impact point during the early phase of the impact, whereas the viscous dissipation near the wall grows gradually into the integral part some time later. The value and the region of major dissipation grow with the viscosity. However, the dissipation function is so strongly dependent on the structure of the flow that an exact theoretical calculation is difficult.

A rough estimate is given below of the magnitude in the early stage (U₀t < h, namely, T < H ≪1) according to Josserand and Zaleski^[13]. The penetration depth is U₀t, the position of ejecta sheet is and the ejecta velocity is assuming the viscous length to be the magnitude of Thus, we have

(14)

with T=U_0t/D. The major dissipation region is the annulated column at the neck, with a volume of . The ratio of the dissipation to the initial kinetic energy is then with T < H≪ 1 and Re≫ 1.

(15)

If we take c as the threshold of e_ν below which the effect of viscosity can be neglected in this phase, solve e_ν(T=H, Re_c)=c, then we can get a critical Reynolds number

(16)

In fact Re_c=16H^1/2 c^-1 is good enough since c ≪ 1 and H≪ 1. With Re_c, we can estimate the dissipation in another way during T < H, i.e., the dissipation ratio is smaller than c when Re is larger than Re_c. We calculate some typical values with H=0.17 and We=924 in consistent with our simulation as follows: Re_0.01=670 (Oh_0.01=0.045); Re_0.1=69 (Oh_0.1=0.437); and Re_0.02=301 (Oh_0.02=0.1). These values indicate that the dissipation caused by viscous force in the extremely early stage is rather small.

Compared with the initial loss of 25%-55% at T=1 in our simulation, a non-ignorable part of dissipation probably arises during H < T < 1. This is because when the dissipation is gradually dominated by the viscous force near the wall rather than by the shear layer along the material interface between the drop and the film. However, because of the complex evolution of the morphology, few theories are available for that phase. Such theories would explain the appreciably different motion for different Re numbers. In a similar way, one could also estimate the dissipation in the later period of time after the crown takes shape and spreads outwards, referring to Pasandideh-Fard et al.^[30], Trujillo and Lee^[12], Roisman^[31] and Roisman et al.^[32]. Note that the transformation into surface energy is rather complicated with regard to the droplets breakup, and is not considered here.

Thus, a three-segment model is believed to be a better approach to the motion of the crown-like jetting with regard to the viscous effect from an energy perspective. The dissipation at the neck is negligible when the penetration is smaller than the film (i.e., T < H), which can be modeled by referring to Howison et al.^[19], Josserand and Zaleski^[13]. The loss due to the viscous force near the wall, together with the shear around the neck, requires further insights during the formation of the crown for H < T < 1. As a discontinuity, the crown spreads outward when the shallow-water approach is a proper approximation. This requires a reasonable initial condition, which can be determined by the loss in the previous period.

5 Summaries and conclusions

The effect of different viscosities on the motion of the crown is studied using the CLSVOF method. Various theories and models involving the influence of viscosity are reviewed. The main conclusions of this study are as follows:

(ⅰ) A three-segment model from an energy perspective might be a better approach to describing crown evolution induced by drop impact, compared with the existing models.

(ⅱ) The dissipation in the extremely early stage (T < H) is slight, and further insights are needed during the formation of the crown (H < T < 1).

(ⅲ) The influence of viscosity on the motion of the crown-like jetting is not obvious in the early stage (around T < 0.4), whereas a visible inhibition of the motion of the crown is seen later on (T>1).

(ⅳ) The initial condition in the models based on the shallow-water approach when time is sufficiently large should be proposed differently according to the viscosity of the liquid, which is dependent on the dissipation in the previous period.