Modelling of droplet impacts on dry and wet surfaces using depth-averaged form

An efficient time-adaptive multigrid algorithm is used to solve a range of normal and oblique droplet impacts on dry surfaces and liquid films using the Depth-Averaged Form (DAF) method of the governing unsteady Navier–Stokes equations. The dynamics of a moving three-phase contact line on dry surfaces is predicted by a precursor film model. The method is validated against a variety of experimental results for droplet impacts, looking at factors such as crown height and diameter, spreading diameter and splashing for a range of Weber, Reynolds and Froude numbers along with liquid film thicknesses and impact angles. It is found that, while being a computationally inexpensive methodology, the DAF method produces accurate predictions of the crown and spreading diameters as well as conditions for splash, however, underpredicts the crown height as the vertical inertia is not included in the model.

The Depth-Averaged Form (DAF) method used in this study is often used in modelling thin-film liquid flow to determine the concourse and velocity of rivers [37][38][39]. The DAF method can be thought of as a first-order accurate long-wave approximation with no Reynolds number limitation. Similar to VOF, this is based upon the N-S equations, and is being used after Veremieiev et al. used an efficient strategy to solve the governing DAF, paving the way for its applicability to droplet impacts [38]. Veremieiev has used this approach repeatedly, to investigate both electrified thin-film flow and droplets of bio-pesticide flowing over foliage [40,41]. The method of depth-averaging 2D N-S equations has also been referred to as the Integral Boundary Layer (IBL) approximation and has also been presented by Shkadov and Ruyer-Quil and Manneville [42][43][44].
In this study, there are 2D axisymmetric, 3D axisymmetric and 3D oblique simulations that are carried out for a variety of fluids, all of which are directly compared with experimental data from Vander Wal et al., Cossali et al. and Šikalo et al. [6,8,11,15]. This allows for a robust examination of the DAF method, as multiple parameters are considered and the accuracy of each determined.
Due to the computational demands of simulations, two main techniques are used to make the solver used more efficient; adaptive time-stepping and a multigrid solution. Adaptive time-stepping enables the algorithm used for the solution to be optimised in a controlled manner, subject to a specific error tolerance [45]. A multigrid solution, rather than reducing all errors on the same fine grid, which can be unacceptable in terms of computational time, will change to a coarse grid to reduce low frequency errors [46]. Different methods have been used to input both these techniques into various algorithms, increasing the efficiency of many computationally-demanding problems, of which droplet impact modelling is included [47][48][49].
The adaptive time-stepping scheme used in this algorithm is based on Heun's second-order predictor-corrector method [45,47], which automatically adjusts the time step subject to the error tolerance, allowing it to optimise. Multigrid methods use a smoother on a sequence of coarser computational grids to reduce high frequency error components on a particular grid [46,48,49]. They are often used to solve the discrete analogues within a wide variety of flow problems where advancements in efficiency and resources are being made [50]. This study combines both adaptive time-stepping and the multigrid method to reduce the computational resource required for the calculations as would otherwise be required. All impact simulations are modelled using C++ and run within the Linux-based Hamilton cluster environment.
The underlying theoretical concepts defining the droplet impact and the DAF method of the governing N-S equations are presented in Sect. 2, along with the conditions used for the various simulations carried out. Section 3 describes the numerical conditions used for the solution method, including the spatial discretisation of the impact domain and temporal discretisation used for the solution. In Sect. 4, a comparison between the theoretical predictions of the model and experimental data is undertaken. These discussions first establish the accuracy of the DAF method for modelling droplet impacts, before identifying the parameters that most affect behaviours and outcomes of impacts. Finally, conclusions drawn from the study are outlined in Sect. 5. Figure 1 shows a schematic of a droplet impacting upon a film in two stages. The liquid is assumed to be incompressible with constant density ρ, dynamic viscosity μ, and surface tension σ . For the droplet simulations the Cartesian X, Y, Z coordinate system is adopted, with the solution domain bounded vertically at the bottom of the liquid film where Z = 0, and from above by the free surface Z = H (X, Y, T ), where T is time and H is the height of the free surface above the solid impact plane. The domain is also bounded horizontally, with a length L in the X plane and a width W in the Y plane, where the origins of both the X and Y axes are located at the centre of each plane.

Problem formulation
Stage (1) in Fig. 1a shows the droplet with an initial in-flight diameter D and initial velocity U 0 = (U t , 0, −U 0 ), where U 0 and U t are the normal and tangential velocity components, respectively. The droplet also has an angle of impact θ , where tan θ = U t /U 0 and is accelerating within a normal gravitational field g = (0, 0, −g 0 ), where  (1) shows the droplet before impact, with Stage (2) showing the droplet at initial contact with the liquid film where it is modelled as a spherical cap, b shows the formation of a crown after impact where 1 is the crown height, 2 is the lower external crown diameter, 3 is the residual shape of the impacting drop and 4 is the cross-section of the free rim g 0 = 9.81 ms −2 is the acceleration due to gravity. When θ = 0 a normal impact takes place and when θ > 0, the impact is oblique. The initial volume of the droplet is V 0 = (1/6)π D 3 .
In Stage (2), the droplet has impacted with a liquid film of thickness H * and the initial droplet impact site has an in-plane radius of R D . The total initial height of the film and droplet is H T and the initial droplet height upon impact is H D , where H T = H D + H * . Upon impact, the droplet is modelled as a spherical cap with volume Once the droplet is in contact with the film, they are considered as one volume of fluid. Immediately after impact, a crown of height H , which is a priori unknown, will form around the impact site and propagate away, capillary waves will also form and fingering or splashing may occur. Figure 1b shows formation of a crown after impact using the same Cartesian X, Y, Z coordinate system. Crown height, shown as 1, is the distance between the liquid film and the bottom of the free rim, shown as 4, which is unstable liquid at the top of the crown from which secondary droplets may form. Crown diameter, shown as 2, is taken to be the lower external diameter of the crown, and both the crown diameter and height are defined to be consistent with [15].

Governing equations
For comparison between the results of this study and experimental data, dimensionless variables are used both to define the simulation model and in the analysis. The variables that define the droplet impact model are transformed into dimensionless scales and are shown below: From Eq. (1), U = (U, V, W ) and p are the velocity and pressure within the fluid and (X, Y, Z , L , W, R, D, H, R D , H D , U, V, W, T, P) → (x, y, z, l, w, r, d, h, r d , h d , u, v, w, t, p) represent the change between dimensional and dimensionless variables, respectively.
The three dimensionless groups governing the droplet impact and its flow regime are defined below: The Reynolds number Re is the ratio of inertia to viscous forces within the fluid, the Weber number We is the ratio of inertia to surface tension forces and the Froude number Fr is a measure of the inertia with respect to gravitational forces upon the fluid flow. Comparatively, Fr is not as influential as Re and We in most practical situations as the effects of g upon small volumes of fluid is minimal. However, as velocity decreases and the film becomes thicker, the effect is no longer negligible; hence, Fr has been considered in this study.
Solving Eqs. (3) and (4) to the boundary conditions (5)-(8), although possible, is extremely challenging and would require a lot of computer resource and time in order to obtain a mesh-independent solution. To bypass this it is possible to use the long-wave approximation based on the assumption that the ratio E = H 0 /L 0 1, showing that the characteristic length scale in the vertical direction H 0 , is much smaller than the characteristic length scale in the horizontal direction L 0 . In addition to aiding the theoretical predictions of the splash, this assumption effectively reduces the dimensionality of the problem by one. Using this assumption, the resulting depth-averaged form of Eqs.
(3) and (4) retains the inertia terms necessary for modelling droplet impacts. The ratio of scales H 0 L 0 is true for most of the times of droplet impact (with the exception of the initial moments) and is used for the DAF derivation as presented in Appendices A and B and [38]. Please note, however, that for consistency with the previous studies and convenience of presenting results, the in-flight diameter is used as both the horizontal and vertical scale, such as H 0 = L 0 = D, in the main text: Re ∂v ∂t where the over-bar denotes the depth-averaged components of velocity as per Eq. (30). The dynamics of droplet impacts on dry surfaces are modelled by a precursor film model [51] that assumes the presence of a thin uniform liquid layer of thickness h * preceding and surrounding the droplet and covering the substrate. The application of the precursor film removes the triple point present at the contact line, removing the stress singularity [52]. Experimental evidence suggests that the physical thickness of the film is extremely small, around 1-100 nm as per [53]. However, the spatial resolution used must be at least of a comparable scale to the precursor film thickness, and using a precursor film that is too large leads to a systematic overprediction of the spreading rate. A disjoining pressure is then introduced across the film at the free-surface condition (7), to account for long-range inter-molecular interactions at the interfaces between the gas, the solid and the liquid droplet and the effect of the static equilibrium contact angle θ e , between the three phases: where n and m are constants representing the disjoining pressure attraction and repulsion, respectively, set to n = 6 and m = 3 for all trials in this study [54]. Note that for impacts on wet surfaces the disjoining pressure (13) is set to zero. Equations (9)-(12) are solved subject to symmetry boundary conditions at the edges of the computational domain, namely: 2

.2 Initial conditions
The initial conditions of the droplet correspond to the moment the droplet impacts upon the liquid film or dry surface. Upon impact, the droplet is modelled as a spherical cap with in-plane radius r d , height above the film h d , and in-plane centre coordinates x d and y d : The relationship between r d and h d is determined via conservation of volume of the droplet before and after impact V 0 = V , which gives The droplet contact angle can now be determined from Eq. (15) for a spherical cap, as follows: The relationship between the radius of a spherical cap r d and the angle θ c at which it meets the surface (referred to as the contact angle for dry surfaces) is now determined by combining Eqs. (16) and (17) as follows: Note that Stage (1) in Fig. 1 cannot be modelled as the initial condition for DAF since the droplet free surface is multivalued. Therefore, the droplet is modelled starting at the earliest point, depicted as Stage (2) in Fig. 1, when θ c = 90 • and h d = r d = 2 −2/3 = 0.6300. Please note for good comparison between the DAF predictions and experimental data, time-shifts, denoted by t s , are applied that correspond to how long the droplet in the experimental data takes to reach d(t s ) = 2r d = 1.2600 and h(t s ) = h d = 0.6300. Following [1], the kinematic discontinuity at the contact line is taken as the initial condition for averaged velocities u andv. For axisymmetric impacts only one component,ū, is required. Oblique impacts require consideration of a tangential velocity, u t , and an impact angle, θ : As stated previously, the theoretical predictions of the DAF model are to be directly compared to the experimental data taken from [6,8,11,15]. The fluid properties for 85% glycerin, isopropanol, water and hexadecane are shown below in Table 1, taken from the respective experimental studies. For impacts on liquid films, the initial radial velocity u 0 = 1, while for impacts on dry surfaces u 0 is determined from the moving contact line velocity as per the Hoffman-Voinov-Tanner law [8], as follows: where c T ≈ 72 rad 3 is the experimental constant, initial contact angle θ c = 90 • , and the static equilibrium contact angle θ e is determined from Eq. (18) using experimental data for the radius achieved by the spreading droplet at its equilibrium; both θ e and u 0 are provided in Table 1. Hexadecane impacts are modelled for multiple values of We and Re; however, only one dry surface impact is modelled, which is reflected in the value of u 0 provided. Both θ e and u 0 are included for only dry surface impacts and reflect the value corresponding to the particular substrate used in each study, which was smooth glass for glycerin and isopropanol and aluminium for hexadecane, while no dry impacts were modelled for water. The separate conditions and related dimensionless values that were used to fully define the impacts for all simulations are shown in Table 2, again, taken from [6,8,11,15]. Oblique impacts were simulated using the same fluid and impact conditions as those modelled by [4], but with no directly comparable experimental data available.

Spatial discretisation
Equations (9)- (12) are solved subject to the attendant boundary conditions, shown in Eq. (14), on a rectangular computational domain (x, y) = (0, l) × (0, w). This domain consists of regular rectangular cells divided into n x and n y cells, respectively, with sides of length Δx = l/n x , and width Δy = w/n y . Each cell is subdivided using a staggered arrangement of unknownsū,v, p and h, as shown in Fig. 2. Similarly, for axisymmetric flow, the computational domain r = [0, l] is divided into n x cells with a staggered arrangement for unknownsū, p and h. For all axisymmetric impacts, the values of l and w were set to l = w = 10 with the centre of the impact at (x d , y d ) = (0, 0), for oblique impacts, the domain dimensions were also set to l = w = 10; however, the centre of impact was set at (x d , y d ) = (−2.5, 0) to accommodate for increased fluid movement in the x-direction. The unknown variables required to predict the droplet impact; thickness h, pressure p and the velocity components (ū,v), are located at cell centres (i, j), and cell faces (i +1/2, j) and (i, j +1/2), respectively. This staggered mesh arrangement avoids the checkerboard instability which results in a periodic pattern of high and low Pseudo-densities [55]. The instability occurs if central differencing is applied to first-order pressure term derivatives and to the terms in the continuity equation when pressure and velocity components are collocated. Further details concerning the spatial discretisation scheme are available in [38].

Temporal discretisation and multigrid solver
An automatic adaptive time-stepping procedure is utilised to decrease computational resource required, as per [38]. This uses estimates of the local truncation error obtained from the difference between an explicit predictor stage and the current solution stage to optimise the time step size. The local truncation error is subject to a tolerance, ε, which then uses the β-method to advance the solution in time, with β = 0.75. This value of β gives a mix of first and second-order differentiation, giving a balance between efficiency and reliability. Fully explicit second-order equations are used to predict values for the unknowns for both 3D and axisymmetric impacts. The equations are solved using a combination of the Full Approximation Storage (FAS) method and a customised multigrid strategy [47]. The discretised equations are solved using a fixed number of FAS V-cycles on intermediate grid levels with up to 10 V-cycles on the first grid level so that residuals are reduced below a specified tolerance. The benefits of using at least one V-cycle at the intermediate level are clear as the relative residuals are reduced by almost two orders of magnitude. Each V-cycle runs through a pre-smoothing, coarse-grid correction and post-smoothing stage. This time-stepping scheme is both efficient and accurate, ensuring small time steps are avoided when the solution varies slowly and enabling the error to be controlled throughout the solution process.

Axisymmetric impacts
Axisymmetric impacts on both dry surfaces and liquid films were modelled. Simulations for 85% glycerin and isopropanol impacting on dry surfaces were compared with [8] and [11], respectively. Water impacting upon a film of depth h * = 0.67 was also modelled and the numerical results compared with [15]. For full data defining the impacts, refer to Table 2.
As stated in Sect. 3.2; for each prediction, the time-stepping scheme requires a specific error tolerance, ε, to be defined to adjust the time step according to this value, optimising the time step in a controlled manner. To determine the optimum value of ε, the computational time required to solve different sizes of mesh grid was compared for a variety of values, using a water droplet impact with We = 297, u 0 = 1 and h * = 0.67. Figure  To ensure the accuracy of results, grid independency studies were carried out to ensure the computational domain was divided into enough cells such that the results were independent of both n x and n y . Two studies were undertaken for both dry and wet surfaces, with the assumption that dry surfaces would require a finer grid due to the length scale of the precursor films compared with liquid films. Both studies were 2D axisymmetric models, and as all 3D models used n x = n y , the results are valid across both 2D and 3D simulations. Figure 4 shows the results of the grid independency study for liquid films, conducted for a water droplet impact with We = 297, u 0 = 1, h * = 0.67 and a time-shift of t s = 0.35, using both crown height and diameter to determine at what number of cells the result was independent. From the studies, for height and diameter, between n x = 2048 and n x = 4096, there is a 4.5% height change and 0.2% diameter change in the final values; therefore, n x = 2048 is considered to be sufficient to ensure accurate results within these error ranges.
In addition to the grid independency study, the figure also compares both the evolution of crown diameter (a) and height (b), as defined in Fig. 1b, over time against experimental data from [15], highlighting the difference in accuracy of the DAF predicting horizontal and vertical action in the droplet impact. In (a) the numerical prediction of crown diameter evolution is almost identical to the experimental data from [15], which is found to have the relation d = C · (t − t 0 ) n where for this case; C = 2.06, n = 0.45 and t 0 = 0. In (b) DAF predicts the crown to form and breakdown much faster than experimental data, as well as underestimating the maximum crown height. This could be due to two related reasons; foremost of which is that vertical inertia is not taken into account, so there is no momentum to carry the fluid to the correct height. This is due to all second order, O(E 2 ), or higher, terms in the DAF equations being factored out, meaning the vertical DAF equation only consists of gravity, surface tension and pressure forces as seen in Eq. (12). As horizontal inertia is considered within the DAF equations, the diameter shows strong agreement. A further reason for incorrect predictions of the crown height is the initial droplet and velocity profile used in the model. From Fig. 2 in [11], it can be seen that when the droplet impacts and reaches h = 0.6300, the fluid is not a perfect spherical cap and there is already spreading of the droplet. Though [11] shows this for a dry surface impact, for a liquid film impact, this will cause the fluid outside and below the impact site to have been  disturbed before the point at which the DAF prediction begins. This forces the fluid down further, generating more vertical momentum than the spherical cap predicts and therefore taking longer to form and dissipate the crown. Figure 5 shows the results of the dry surface grid independency study, carried out for a glycerin impact with We = 51, θ e = 63.8 • , according to the relation in Eq. (17), u 0 = 0.02 and a precursor film thickness of h * = 0.001. Though two values of precursor film thickness are compared in the discussion below, the fact that h * = 0.001 is a smaller scale than h * = 0.01, allows the results of the grid independency study to be applied to both precursor films. From the figure, it can be seen that between n x = 8192 and n x = 16384 there is only 1.1% change in the final spreading diameter; therefore, n x = 8192 was used for all dry surface impact models. Figure 6 shows the diameter of the same glycerin impact over time where We = 51, θ e = 63.8 • and u 0 = 0.02, but compares two values of precursor film, h * = 0.01 and h * = 0.001, with each other and experimental data from [8]. A time-shift of t s = 0.47 is required to align the initial condition of the DAF impact where h d = r d = 0.6300 with the corresponding instant of the experimental data. A precursor film is required for dry impacts due to the nature of DAF, which requires a depth of fluid at every point in the computational domain to be able to be solved. From the figure it can be seen that decreasing the thickness of h * by a factor of 10 results in a more accurate prediction, with d = 1.499 at t = 100 for h * = 0.001 and d = 1.467 for h * = 0.01. This is compared with d = 1.560 from the experimental data and shows an increase in accuracy of 2.1% for the lower value of h * . As the thicker precursor film is closer in scale to the droplet thickness h, it will have a greater effect upon the spreading rate and diameter, increasing the error from the experimental data. The reason for the underprediction of both film thicknesses is that until total equilibrium the droplet may not be a perfect spherical cap and may have a slight rim, as seen in Fig. 1a in [8]. This would lead to an overprediction of the static equilibrium contact angle; however, once the droplet reached complete equilibrium, the prediction should match the experimental data for both precursor films. The agreement of the prediction for h * = 0.001 is strong, with the diameter at t = 100 within 4.5% of the experimental data. Figure 7 shows the diameter of an isopropanol impact over time with We= 93, θ e = 6.5 • and u 0 = 0.40. As with the glycerin impact, two values of precursor film are compared with the experimental data, with a time-shift of t s = 0.21. The results of this simulation show less agreement with experimental data than glycerin. As isopropanol has a much lower viscosity and lower static equilibrium contact angle it takes a lot longer to reach equilibrium. At t = 10, h * = 0.001 is 4.2% above the experimental data and h * = 0.01 is 7.1% above, giving an increase in accuracy of the prediction by 3.0%. This discrepancy can again be explained by the physical droplet not being a perfect spherical cap at both the initial impact and at equilibrium. As the initial modelled droplet starts at a lower height than the actual impact, there is less vertical acceleration, reducing the spreading rate. The experimental data from [11] also did not reach times past t = 8, so though a line of best fit was used to predict the droplet diameter at equilibrium, it is possible the spreading rate could have changed beyond this time, and would otherwise align more with the DAF predictions. Decreasing h * results in more computational resource being required; however, for this 2D axisymmetric case, this is not a limiting factor. In contrast, for 3D cases, the simulation takes much longer and a compromise between accuracy and time has to be considered.

3D normal impacts
3D normal impacts on both dry surfaces and liquid films were modelled and compared with the photographic compendium from [6], recorded at the frame rate of 2,000 frames/s; the most viscous fluid hexadecane was considered. The same adaptive time-stepping scheme was used, as per Fig. 3, with ε = 0.1. A grid size of n x ×n y = 2048×2048 was used for 3D impacts. This size of mesh is not grid independent for dry surface impacts, as shown in Fig. 5;  [11], for We = 93, θ e = 6.5 • , u 0 = 0.40 and applying a time-shift of t s = 0.21 however, due to the time constraints, performing a simulation with n x ×n y = 4096×4096 or n x ×n y = 8192×8192 was infeasible. Any 3D results for dry surface impacts therefore have a potential ±5.2% error range. Each 3D DAF prediction corresponds to the same time as the frame shown in the adjacent experimental figure and is viewed from a high angle. For full data defining the impacts, refer to Table 2.
The dry surface hexadecane impacts were solved in 3D to compare with [6] and provided qualitative confirmation of the analysis of the glycerin and isopropanol impacts. Figure 8 shows a 3D hexadecane droplet impacting on a dry surface for We = 106, θ e = 11.2 • and u 0 = 0.31, with no time-shift. θ e was again calculated using the relation in Eq. (18), having assumed Frame 12 in Fig. 8a shows the equilibrium state and measuring d = 3 in this frame. No time-shift was required for this impact as the initial droplet shown in Frame 1 mirrors the prescribed spherical cap initial condition. A precursor film of h * = 0.01 was used to reduce computational time needed for the calculation. The figure shows strong agreement with the experimental data, showing similar spreading rates and diameters throughout the impact, reaffirming the results found in Sect. 4.1. The depth of the droplet across the diameter is also accurate in the numerical prediction. In Frame 3, the residual shape of the initial impact is visible, which is reflected in the DAF prediction, in Frame 12, an outer rim is formed which is thicker than the internal diameter. The accuracy of the droplet height prediction in this figure is due to minimal vertical inertia in a dry surface impact compared to a liquid film impact. As there is little vertical movement after impact, the consideration of gravitational and pressure forces are sufficient to describe the droplet height. Figures 9 and 10 show 3D hexadecane impacts on thin films, where h * = 0.1, for We = 269 and We = 566, respectively, with u 0 = 1. Figure 9 has a time-shift of t s = 0.51 and Fig. 10 a time-shift of t s = 0.17. For the same thickness of film, faster velocities require smaller time-shifts due to the time required to reach the condition of h d = r d = 0.6300. As with water impacts, the crown diameter is predicted accurately in both cases. The thickness of the crown, however, as predicted by DAF is larger than the experimental data, no rim is shown in the predictions and the crown height is again underestimated. The inaccuracies of the predictions of both crown thickness and height are both due to the lack of vertical inertia in the model. The kinematic discontinuity of the crown is not as prominent as in the physical situation, so more fluid forms the crown due to surface tension and viscosity, described by Re and We, which also explains the lack of a rim. However, a rim would be unable to be predicted by DAF Fig. 9 Hexadecane droplet impinging upon a thin liquid film: comparison of DAF prediction with the complementary experimental data, Fig. 4 from [6], for We = 269, u 0 = 1, h * = 0.1 and applying a time-shift of t s = 0.51 for times t = 1.12, 2.20, 3.29 and 4.37: a free-surface profiles, b cross-section profiles through the centre of the impact anyway due to the multivalued free surface. In both figures the DAF prediction also shows capillary waves being produced, which only begin to form in later frames in the physical situations. In the experimental data, the capillary waves are only produced as the crown falls, and again due to the lack of vertical inertia in the DAF predictions, the crown begins to fall for the DAF predictions quicker than the physical situation, leading to capillary waves being produced earlier. Figure 10 also shows splash occurring after impact, between Frames 7 and 13, with the DAF method diverging within this timeframe. Though Fig. 10 shows some slight instabilities in the crown, Fig. 11 shows the point at which the method diverges, corresponding to approximately Frame 10, for the same impact but with a grid size of n x × n y = 512 × 512. This figure highlights the creation of inertia-driven instability within the crown and how this causes points on the crown to become too thin for the mesh grid used. The coarser mesh better shows the instabilities due to the resolution of the figures produced, the instabilities are still present for a grid of n x × n y = 2048 × 2048; however, this resolution makes it more difficult to observe the exact gridpoints the instabilities cause the DAF method to diverge.  Without knowing the exact time at which splash occurred in the physical impact, to further investigate whether splash is predicted correctly by DAF, different cases were explored. Tables 3 and 4 compare the occurrence of splash in experimental data with divergence of DAF. From the tables, it can be seen that for all but one case, the method diverges when splash occurs. For all these cases, the divergence occurred within the correct timeframe when compared to [6]. The only case in which the method diverging and presence of splash in the experimental data do not correspond is where the splash was minimal. In this case, it may be that the precursor film of h * = 0.01 prevents the splashing due to surface tension or by being of a similar scale, and the number of cells in the computational domain would need to be increased to detect the splashing. By observation, not all splashes in these cases occur in the same manner as in Fig. 11, as some happen instantaneously with no crown formed. For all cases, it is assumed that the reasons the DAF method diverges and therefore the prediction of splash is based upon the fluid free surface become multivalued or too thin for the mesh to solve.

3D oblique impacts
3D oblique impacts on liquid films were modelled. As no experimental data were available to compare the impacts with, attempts to compare the model with the analytical solution from [4] were made. For oblique impacts, a finer grid was required of n x × n y = 4096 × 4096 to obtain a solution. This finer grid increased the time taken for the simulations to run; therefore, only small time scales could be predicted. The results in [4] are shown for t = 5, 10 and 15, but due to the limitation of computational resource, no comparison could be made for t = 15. The predicted crown shape for an oblique water impact with We = 297, h * = 0.67, impact angle θ = 45 • and u 0 = u t = 1 for times t = 0.1, 0.3, 1 and 3 is shown in Fig. 12a-d. The figure shows the droplet moving laterally, with most action between t = 0.3 and t = 1, with the diameter of the impact doubling and crown height increasing. There could be the same issue with predicting crown height as with earlier simulations, in that the lack of consideration of vertical inertia predicts most movement at an earlier time than the physical reality. This could also mean that the crown is under-predicted and the actual impact would show more dramatic dynamics. At each time, instabilities within the crown can be seen, which could be a mechanism of splashing, as in Sect. 4.2. However, these instabilities smooth out and the solution does not diverge, possibly due to the larger grid size used. Figure 12e shows a comparison of the analytical solution presented of the shape of the external base of the crown for t = 5 and t = 10 with the predictions in Fig. 7b from [4]. The figure shows [4] to predict greater crown diameter than that of DAF at these times. As shown in Fig. 3 DAF predicts the crown diameter for normal impacts well, while it is expected that the analytical solution [4] overpredicts the crown diameter since it is obtained ignoring gravity, viscosity and surface tension, i.e. in the limits of We → ∞, Re → ∞ and Fr → ∞. Data from [15] also demonstrate that the crown diameter increases with impact velocity as well as We, Re and Fr.

Conclusions
This study used the DAF method to predict a range of droplet impacts, for different fluids, film thicknesses and impact angles, comparing the predictions to experimental data to validate the accuracy of DAF in droplet impact modelling. For dry surface impacts, there is strong agreement between the theoretical prediction and experimental data. Wet surface impacts contrasted the accuracy of the DAF in predicting horizontal and vertical movement in droplet impacts. While predictions of droplet and crown diameters were accurate, crowns were predicted as developing and breaking down too quickly, with overall height underestimated. This is due to the lack of consideration of vertical inertia in the DAF model used, caused by ignoring all terms of order O(E 2 ) or higher in the DAF equations. The vertical DAF equation therefore only consists of gravity and pressure terms.
In all but one case, DAF also predicted splash within the correct time range, by the fact the method diverges. In the exceptional case, there was minimal splashing and therefore to make DAF a greater predictor of splash, a finer mesh is required. The splash shown in Fig. 11 occurred due to instabilities within the crown; however, different types of splash occurred for different cases, all of which are assumed to cause instabilities where the fluid becomes too thin for the mesh to solve, or the free surface became multivalued.
To progress the work, further adjustments need to be made to the model. Foremost, the long-wave assumption used in the DAF predictions will have to be altered so that vertical inertia is incorporated into the model and crown height can be modelled as accurately as the diameter. To better predict splash, finer meshes will be required, though this increases the computational resource required for each simulation. Therefore, further optimisation of the model is required. In its current state, this model can therefore accurately predict spreading and splashing, which are married in terms of applicability. The DAF prediction could also be used to investigate cases in which spray density is high enough that the assumption that behaviour of a single drop represents that of a spray no longer applies. Further work could therefore investigate how multiple droplets interact when impacting both near to or on top of each other. The model could be particularly useful in fields where thin-film deposition is of importance, such as spray painting, surface cleaning and ink-jet printing.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Appendix A. DAF derivation
As discussed in Sect. 2.1, solving Eqs. (3) and (4) to the boundary conditions (5)-(8), although possible, is computationally challenging and would require a lot of computer resource and time in order to obtain a mesh-independent solution. To bypass this it is possible to use the long-wave approximation based on the assumption that the ratio E = H 0 /L 0 1, showing that the characteristic length scale in the vertical direction H 0 is much smaller than the characteristic length scale in the horizontal direction L 0 . In addition to aiding the theoretical predictions of the splash, this assumption effectively reduces the dimensionality of the problem by one. Formulating (3) and (4) in terms of scales H 0 and L 0 leads to: The boundary conditions (7) and (8) where ξ = z/ h. Using relations (30) and (34) then leads to the following analytical expressions for these terms: The final form of the DAF equations, using these friction and dispersion terms are found in Eqs. (9)-(12).

Appendix B. Axisymmetric DAF equations
The axisymmetric version of the DAF equations leads to two extra inertial pressure terms in the governing DAF equations, as seen below: E Re ∂ū ∂t +ū ∂ū ∂r