\section{Introduction}
\label{sec:Intro}

Since the days of Worthington \cite{worthington1876forms}, the problem of droplet impact has offered the fluid dynamics
research community exciting opportunities and challenges over the course of its history.
For the first time in a systematic manner, in his book entitled \textit{A study of splashes} \cite{worthington1908},
Worthington makes use of early photographic technology (alongside careful
sketchwork) to provide a comprehensive visual interpretation of splashing phenomena.
The framework has since captivated the interest of theoreticians and experimentalists alike, as it
incorporates one of the most invitingly simple geometrical configurations,
while at the same time giving rise to diverse and rich phenomena of immense scope.

A plethora of application areas benefit from understanding the outcomes of droplet impact events.
We emphasise in particular the role of droplet splashing (or absence thereof) in printing technologies (\cite{vanDam1,jung2012}),
combustion (\cite{Moreira1}), granular material interactions at all scales (\cite{thoroddsen2001,marston2012}),
electronics (\cite{kim2007}) and spray-cooling in nuclear reactors (\cite{sawan1975}).
The design of superhydrophobic coatings in relation to droplet impact dynamics (\cite{Tsai2,Deng1})
is yet another prime example of the widespread applicability
of this canonical problem.

Recent reviews provide an excellent insight into the state-of-the-art in the field within each decade (Rein \cite{rein1993} in the 1990's,
Yarin \cite{Yarin1} in the 2000's and more recently Josserand \& Thoroddsen \cite{josserand2016drop}).
The area has witnessed a very strong surge in the past decade, fueled in part by the development of
progressively more powerful imaging technologies, with both frame rates and resolutions capable of capturing
details beyond the scope of previous equipment (see \cite{thoroddsen2008} as well). Furthermore, the
improvement of numerical algorithms and usage of high performance computing has enabled computational studies
that complement and inform both experimental and analytical work. We focus particularly on the volume-of-fluid
package \textit{Gerris} \cite{popinet1, popinet2}, which is one of the most popular open-source tools due
to its strengths in dealing with interfacial flows on a range of very different scales. Comparisons with experiments,
as well as analytical work have been consistently robust, be it in cases of liquid-liquid impact \cite{Thoraval2,Agbaglah1} or
impacts of liquid on solid surfaces \cite{Visser1,philippi2016drop,wildeman2016spreading}.

In the case of normal (perpendicular) impact at low-to-moderate velocities (and depending on specific fluid properties),
an axisymmetric assumption can be used
in analytical and computational investigations. The reduction in dimensionality
is a significant advantage that has led to very efficient (axisymmetric) computations
and good agreement with experiments. Visser et al. \cite{Visser1} for example, while
innovating experimental technology enabling the time-resolved investigation of micron-sized drop impacts, have managed
to conduct successful comparisons with direct numerical simulations at impact speeds of up to $50$ m/s,
a regime which is commonplace in combustion, inkjet printing or aircraft-related applications.
In the respective scenario, the small drops spread onto the surface in what is known as pancaking motion, with
the axisymmetric approximation remaining valid in the absence of splashing events.

In cases where spreading and later retraction rather than splashing occurs, the vast majority of efforts
have been dedicated towards identifying quantities such as the maximal spreading radius
(\cite{stow1981,clanet2004maximal,fedorchenko2005effect,Roisman2009,Schroll1} and most recently \cite{wildeman2016spreading}),
as well as the resulting minimal film thickness, retraction dynamics and the role of the internal boundary layer - see \cite{bartolo2005retraction} and
in particular \cite{Eggers1} for a comprehensive investigation of the above.

At higher speeds however, there is still an ongoing debate as to how the splashing phenomena are first initiated,
and the splashing threshold in particular. Up until the groundbreaking experimental investigation of Xu et al. \cite{Xu1},
there have been numerous attempts to characterise the transition from spreading to splashing dynamics in the classical impact problem in terms of drop-related
parameters only (density, viscosity, impact velocity, surface tension). The Chicago group discovered, however,
that decreasing the ambient air pressure may completely suppress splashing. As such, a host of additional modelling, experimental
and numerical efforts have been initiated, with the work of Riboux et al. \cite{riboux2014experiments}
proposing a model deducing a threshold splashing velocity as a function of a generalised set of key parameters containing
the liquid density and viscosity, the drop radius, gas density and viscosity, the interfacial tension coefficient,
as well as the nanometric mean free path of the gas molecules.

Once the drop splashes, there is very little attention dedicated to the ensuing dynamics, with the sizes and velocities
of secondary drops being prohibitively small experimentally and computationally, although advances have taken place recently in terms of simplified models.
In particular, Riboux \& Gordillo \cite{riboux2015diameters} have proposed a one-dimensional approach to predicting sizes
and velocities of ejected droplets for $\mathcal{O}(1)$ mm sized impacting drops and low speeds, finding reasonable agreement with experiments.

As underlined by Josserand \& Thoroddsen \cite{josserand2016drop}, there are several exciting challenges
lying ahead, two of which are of great importance in the context of the present work. First of all,
gaining an improved understanding of splashing, particularly in difficult high speed conditions of industrial relevance, is moving more and more within reach,
and further detailed investigation using the available tools is needed. Secondly, oblique
impacts are rarely analysed (exceptions being \cite{Mundo1,Sikalo1, bird2009inclined}) due to the additional flow complexity.
Most often, qualitative rather than quantitative phenomena are explored in detail. The exceptions tend to focus on large
scale effects at the level of the entire drop, as opposed to details at the level of the splashing itself and the interesting local structures arising.
Both of these themes lie at the heart of the present work, which focuses on the modelling
and computation of oblique three-dimensional drop impact in aerodynamic conditions.

In aircraft-oriented research and design involving drop impact, the relevant scales are often dictated by the size of the parts
that are most affected by phenomena such as water impingement, retention and finally icing and its prevention. The wings or nacelles
are several metres long, while computing accurate air flows around them requires domains that span tens of metres in all dimensions.
This becomes highly prohibitive in terms of accurate resolution of the intricate and sensitive physical effects
pertaining to drop impact, which often happen at sub-micron scales in the order of tens to hundreds of microseconds.
As such, particle-trajectory calculations of various degrees of complexity have thus far proven to be the only
tractable solution in industrial setting.

There are several important limitations of current models, as pointed out by Gent et al. \cite{Gent1} in a relatively recent review:
\begin{itemize}
 \item droplets are assumed to be spherical and non-deformable as they approach the solid surface, hence topological transitions
 such as the emergence of secondary drops either before or after impact are not considered;
 \item phenomena related to multiple drops such as collisions are completely ignored;
 \item aerodynamic drag, gravity and buoyancy are assumed to be the sole forces affecting the drop trajectories;
 \item whereas the local velocity of the air flow is embedded into the ordinary differential equations governing
 the updates in drop trajectories, the liquid mass is assumed not to affect the surrounding air flow;
 \item once on the surface, empirical models translate the drop contribution towards liquid film formation
 and its movement further downstream along the surface of interest.
\end{itemize}
Many of these assumptions become inaccurate in the context of the large supercooled droplets (larger than several
tens of microns) found in the atmosphere. The difficulties outlined above have yet to be overcome,
and most modelling is performed at a highly coarse-grained level
(\cite{potapczuk1993ice,bragg1996aerodynamics,rutkowski2003numerical,wright2004,wright2005validation,wright2006,honsek2008,bilodeau2015}), with semi-empirical relations of varying complexity
being proposed in order to match with the rich but ultimately limited experimental data available
by NASA experiments conducted by Papadakis et al. \cite{papadakis2003,papadakisExp}.
The focus here is primarily on the final water retention values rather than the more fundamental problem
of the detailed impact process, making it ideal from an engineering standpoint but offering limited
insight into the underlying physics. In the past few years, the group at INTA/Madrid \cite{vargas2012mechanism,sor2015modeling} have looked
in more detail into the deformation of large-scale drops prior to impact, with results
that indicate regimes far more complex than captured by the typical assumptions mentioned above.
Several studies focusing on recent numerical advances in the high speed regime ($>50$ m/s impact velocity) have emerged,
particularly for impacts onto liquid, but also onto solid surfaces \cite{ming2014lattice, cheng2015numerical, guo2016investigation, Cherdantsev2017, Xie2017}.
These offer exciting opportunities to study short timescale phenomena beyond the reach of traditional particle methods,
however up to this point there have been few attempts to integrate the drop impingement process
into a framework that includes a more realistic model for the movement and effect of the air flow around the bodies of interest.

The present work bridges the relevant scales in the problem of drop impact onto aircraft surfaces
and proposes a suitable model for the air flow around the solid bodies of interest in which we then
accurately resolve the drop impingement process. While the drops are initialised as spherical
sufficiently far away from the body, we characterise their deformation prior to impact and
the spreading/splashing thereafter, depending on drop sizes and angles of impingement.
We focus on the asymmetric features of the drop spreading when droplets are very small (less than a few tens of microns),
phenomena which to our knowledge have yet to be observed. As the drop size increases, we quantify
the sizes and positions of the secondary drops emerging as a result of the impingement
and provide useful metrics for practitioners looking to improve water retention calculation methodology
and a deeper understanding of the physics involved in the impact process under challenging conditions.
All flow parameters have been carefully chosen to match with previous experimental studies or known flight-specific values,
while many of the quantified metrics are also compared to classical theoretical results where applicable.

The investigation is structured as follows. We introduce the proposed mathematical model in Section~\ref{sec:Model},
followed by a detailed description of the computational framework in Section~\ref{sec:Numerics}.
We then analyse our findings in Section~\ref{sec:Results}, focusing
on both pre-impact dynamics in subsection~\ref{subsec:preimpact} and post-impact dynamics in subsection~\ref{subsec:postimpact}.
These results are discussed and placed into context in Section~\ref{sec:Conclusions},
followed by concluding remarks.

\section{Mathematical Model}
\label{sec:Model}

In the present section we elaborate on how we adapt the classical problem of drop impact to the high speed flow conditions of interest
around aircraft surfaces.
First we discuss some useful assumptions allowing us to reduce geometrical complexity in the problem in subsection~\ref{subsec:Scales},
after which we expand on the mathematical model itself, outlining the relevant equations, initial and boundary conditions.

\subsection{Scale Transition}
\label{subsec:Scales}

The full model geometry discussed in previous paragraphs (aircraft wings/fuselage components) is far too complicated - and specific - from many points of view.
To begin with, our aim is to present a general methodology, applicable to a number of surfaces rather
than a specific specialised geometry. Secondly, the multi-scale modelling of both the background air flow
around the large scale body and the splashing dynamics within the much smaller drop impact regions
is beyond reach in terms of theoretical and current computational resources.
We thus employ several simplifications
to enable a closer inspection of a much more amenable problem, which still preserves the main physical
characteristics we wish to address.

Based on the disparity between the two scales in the problem (the impacting
droplet diameter and the solid body it impinges upon), we assume the curvature of the body
to have negligible effects.
To justify this approximation, the radius of curvature of the
leading edge of a typical NACA airfoil or nacelle lipskin, the most sensitive regions to water retention and icing,
is estimated to be of $R_b = \mathcal{O}(10^{-1})$ m for standard commercial aircrafts.
For a reasonably large droplet
of radius $R = 100\ \mu$m, we find $R/R_b \approx 10^{-3}$.
Thus, from the perspective of modelling the local droplet impact, the surface can be considered as approximately flat.
From a different viewpoint, we zoom in sufficiently close to the surface of the solid
body, such that in the respective region the droplet diameter is the representative lengthscale
and hence the details of the impact can be carefully studied.

\subsection{Governing Equations}
\label{subsec:Equations}

The framework of studying these fluids as incompressible in laminar flow conditions is a natural choice in the context of our problem,
as the primary target flight regimes of take-off and landing are characterised by relatively low velocities compared to
those reached at higher altitudes. Furthermore, most droplet impingement events are concentrated close to the leading
edge of the geometries of interest, where the flow has yet to enter the transition from laminar to turbulent state. Even in such circumstances,
a complex and likely empirical turbulence model would prevent the inspection of the detailed liquid dynamics, which is
the main goal of the present investigation.

The model fluids are assumed to be incompressible, immiscible and viscous.
Subscript $1$ is used to refer to the fluid inside the drop (taken to be water),
whereas subscript $2$ decorates quantities in the surrounding (air) flow.
Let $\rho_{1,2}$ and $\mu_{1,2}$ denote the constant densities and dynamic viscosities of the two fluids in the system.
The constant surface tension coefficient at the interface is given by $\sigma$.
Velocity vectors $\textbf{U}_{1,2} = (U_{1,2}, V_{1,2})$
and pressures $P_{1,2}$ are used in the formulation of the dimensional momentum and continuity equations
\begin{align}
 \rho_1 (\textbf{U}_{1t}+(\textbf{U}_1 \cdot \nabla)\textbf{U}_1) &= - \nabla P_1 + \mu_1 \Delta \textbf{U}_1, \label{eq:Droplets1} \\
 \rho_2 (\textbf{U}_{2t}+(\textbf{U}_2 \cdot \nabla)\textbf{U}_2) &= - \nabla P_2 + \mu_2 \Delta \textbf{U}_2, \label{eq:Droplets2} \\
  \nabla \cdot{{\textbf{U}}_{1,2}} &= 0. \label{eq:Droplets3}
 \end{align}
 Gravitational forces are assumed to be negligible.
 There are two lengthscales in the problem: the droplet diameter $D$, the natural choice for the reference lengthscale, and the size of the (finite)
 computational domain $L$. We scale lengths by $D$, velocities by a reference background velocity $U_{\infty}$
 and pressures by $\rho_1 U_{\infty}^2$. The emerging non-dimensional parameters are
\begin{equation}
\textrm{Re}= \rho_1 U_{\infty} D/ \mu_1,\quad \textrm{We} = \rho_1 U_{\infty}^2 D/\sigma,\quad K = \textrm{We} \sqrt{\textrm{Re}} = \sqrt{\rho_1^3 D^3 U_{\infty}^5/ (\sigma^2 \mu_1)}.
\end{equation}
The Reynolds number $\textrm{Re}$ and Weber number $\textrm{We}$ appear directly from the non-dimensionalisation procedure,
while the splashing parameter $K$ is introduced as an intrinsic element of a drop impact problem.
The expression, originally introduced by \cite{stow1981}, has been used to classify the possible outcomes of the impact.
This parameter has been controversial in the literature and
cannot independently account for the classification of the complicated impact process (see \cite{Xu1, mandre2012}),
however it serves as an indicator of the force of the splashing and permits comparisons with
previous investigations.

We also introduce density and viscosity ratios
\begin{equation}
r=\rho_1/\rho_2,\qquad m=\mu_1/\mu_2,
\end{equation}
and non-dimensionalise the governing equations~\eqref{eq:Droplets1}-\eqref{eq:Droplets3}, resulting in
\begin{align}
 \textbf{u}_{1t}+(\textbf{u}_1 \cdot \nabla)\textbf{u}_1 &= - \nabla p_1 + \textrm{Re}^{-1} \Delta \textbf{u}_1, \label{eq:Droplets4} \\
 \textbf{u}_{2t}+(\textbf{u}_2 \cdot \nabla)\textbf{u}_2 &= - r \nabla p_2 + rm^{-1}\textrm{Re}^{-1} \Delta \textbf{u}_2, \label{eq:Droplets5} \\
  \nabla \cdot{{\textbf{u}}_{1,2}} &= 0. \label{eq:dDroplets6}
 \end{align}
 The non-dimensional timescale is $D/U_{\infty}$.
 The typical fluid properties in the case of water and air at near freezing temperature (close to $0^{\circ}$ C)
 are given as follows. Water has density $\rho_1 = 999.8\ \textrm{kg/m}^3$ and dynamic viscosity $\mu_1 = 1.16 \times 10^{-3} \ \textrm{kg/ms}$, while
 the air density is $\rho_2 = 1.21\ \textrm{kg/m}^3$ and its dynamic viscosity $\mu_2 = 1.81 \times 10^{-5} \ \textrm{kg/ms}$.
 The constant surface tension coefficient is $\sigma = 7.2 \times 10^{-2}\ \textrm{N/m}$ and
 a representative value for the velocity of the background flow is $U_{\infty} = 78.44 \ \textrm{m/s}$.
 This value has been selected to coincide with classical experimental investigations \cite{papadakisExp},
 as well as subsequent numerical investigations in the aerospace engineering community (e.g. \cite{bilodeau2015}).
 We underline the large density ($r=826.28$) and viscosity ($m=64.09$) ratios, which pose significant numerical challenges - these
 are touched upon in Section~\ref{sec:Numerics}. Tables~\ref{tab:dropDeformationParameters} and~\ref{tab:dropImpactParameters}
 in the results discussion indicate the values of the key dimensionless groups in the problem and highlight
 the violent high speed impact regime investigated here.

To define the interfacial conditions governing the motion of the drop, we
assume a sharp interface $y=S(x,t)$; subsequently this is relaxed
in the context of the volume-of-fluid methodology employed in the direct numerical simulations.
The prescribed interfacial conditions are, in order, the kinematic condition, the continuity of normal
and tangential stresses, and continuity of velocity components:
\begin{align}
w_i &= S_t + u_i S_x + v_i S_y,\qquad i=1,2, \\
\left[\textbf{n}\cdot \boldsymbol{\mathcal{T}} \cdot \textbf{n}\right]^1_2 &= \textrm{We}^{-1} \kappa,\\
\left[ \textbf{t} \cdot \boldsymbol{\mathcal{T}} \cdot \textbf{n} \right ]^1_2 &= 0,\\
\left[\textbf{u}\right]^1_2 &= 0,
\end{align}
where $[(\cdot)]^1_2=(\cdot)_1-(\cdot)_2$ represents the jump across the interface,
$\textbf{n}$, and $\textbf{t}$ are the unit normal and tangent to the interface, respectively,
and $\kappa$ is the interfacial curvature. The stress tensor $\boldsymbol{\mathcal{T}}$ is given by
\begin{equation}
 \mathcal{T}_{ij}=-p\delta_{ij} + \mu \left( \dfrac{\partial u_i}{\partial x_j} + \dfrac{\partial u_j}{\partial x_i} \right),
\end{equation}
where the appropriate subscript is used in different fluid regions.
The initial and boundary conditions for the finite computational domain are described in the following subsection.

\subsection{Background Flow}
\label{subsec:Background}

One of the most important features of the model is the interaction
between the liquid drop and the air around it. In typical experimental conditions,
droplets are formed at the tip of an injection device and fall under gravity, with the height
of the device being varied in order to adjust the terminal velocity and hence fix the relevant
dimensionless parameters. In order to reach velocities beyond $\mathcal{O}(1)$ m/s it is necessary
to have some form of ejection mechanism that ensures not only reproducibility of the shapes, but
also a stability of the dynamics in early stages as the drop travels through the quiescent air flow
and may become immediately sheared and violently deformed and broken up. As such, most of the
investigations concerning velocities above $10$ m/s are restricted to very small drops (well below $100$ microns),
such that surface tension is strong enough to preserve the approximately spherical shape of the drop.

In flight conditions, leading edge droplet impact can be locally embedded in a stagnation-point
flow which develops into boundary layers on either side of the geometry.
As such, most droplets encounter a developing boundary layer structure with a strong shear component.
In an effort to reproduce the same type of air flow environment while preserving generality,
we proposed an oblique-stagnation point flow model for the air flow, with the liquid drop
being seeded sufficiently far away from the body on the dividing streamline of the flow. The reasons
behind this choice are twofold:
\begin{enumerate}
 \item Far away from the surface the drop should retain its shape and setting a uniform velocity field
 inside the drop with a zero (as in most desktop experiment setups) or purely horizontal (along the body) air motion would produce instantaneous
 breakup of the drop. The choice for stagnation-point background flow and the initial position of the drop ensures that the air
 flow undergoes only small changes until sufficiently close to the surface, which is when we expect the drop to start
 deforming in real life conditions.
 \item The stagnation-point flow has the same characteristics in the near vicinity of the point of zero
 velocities as on the aircraft surface, in that boundary layers are developing on either side of it and
 growing as we move further downstream. As such the liquid drop is subjected to the shear flow naturally occurring
 above the solid surface. This is best represented in the highly oblique impact cases,
 in which the air flow streamlines near the surface have strong deviations from their far-field orientation. Oncoming drops
 depart from their host streamlines close to the surface and their final impingement points
 are well within the boundary layer growth region. The choice in initial positioning of the drop thus retains generality, while
 ensuring suitable conditions for the early stages of the drop dynamics.
\end{enumerate}

There are however several points to be made prior to advancing to the mathematical description of the model.
First of all, the dynamics of drops in uniform flow has been extensively investigated and the deformation
characteristics for large enough drops are very rich (see \cite{Jalaal1} for a recent computational study).
Therefore, even in the case of tailoring the initial position of the drop to a region of uniform air flow,
the drop is anticipated to suffer significant deformations as it moves towards the solid body. The size
of the finite computational domain can then be used to alleviate (or enhance) this effect. Secondly, it should be
noted that there still remains a fundamental difference to the practical scenario in which a solid body is moving
through high liquid water content clouds (with stationary water drops of varying sizes) as opposed
to drops impinging onto a static solid surface, as in the present case. Here we are enhancing the inertial
contribution in the pre-impact drop dynamics and our choice in initial position of the drop does ultimately
affect the liquid volume impinging onto the surface. Previous experimental results have however been used
as guidance in order to best account for the complex flow dynamics, while retaining
a suitable well-controlled flow environment.

With the above properties in mind, we underline that the background air flow poses its own non-trivial
challenges. The history of the problem dates back to Hiemenz \cite{Hiemenz1}, who
was the first to present a solution for the two-dimensional normal stagnation point flow. Howarth \cite{howarth1951cxliv} then
extended the formulation to three dimensions.
The oblique case was first touched upon by Stuart \cite{Stuart1}, then later rediscovered independently by Tamada \cite{Tamada1} and Dorrepaal \cite{Dorrepaal1}.
There have been a number of corrections, extensions and generalisations on the main problem, such as \cite{wang2008similarity} and \cite{Tooke1}, including
extension to two-fluid systems (air flow impinging onto liquid films above a solid surface), as studied by \cite{tilley1998oblique} and \cite{Blyth1}.
As far as we know there is no general analytical solution to the three-dimensional stagnation-point flow problem at an arbitrary angle.
As such, we attempt to recreate this type of flow numerically using a combination of suitable boundary conditions that preserve
its main characteristics in the interior of the domain.

The three-dimensional computational box is selected to be of size $4L/D \times L/D \times 2L/D$ in $(x,y,z)-$directions,
where $L/D$ is taken to be of size $20$, i.e. $20$ drop diameters. The flat solid surface is taken to be in the $(x,z)-$plane,
with no-slip and impermeability conditions prescribed in this region,
such that at $\textbf{u}_{1,2} = 0$ at $y=0$.

In order to model the oncoming flow at an arbitrary angle of incidence
$\theta_i$
we prescribe inflow conditions given by
\begin{equation}
 u_2(x,L/D,t) = \cos(\theta_i), \qquad v_2(x,L/D,t) = -\sin(\theta_i) \textrm{ at } y= L/D.
\end{equation}

Laterally we impose typical free outflow conditions on all remaining four sides of the box.
The main reason for doing so relates to the movement of the secondary drops resulting as a consequence
of the splash which cause perturbations in the flow field, making it difficult to fix velocities at the boundaries.
The initial conditions are set to
\begin{equation}
 u_2(x,y,0) = \cos(\theta_i), \qquad v_2(x,y,0) = -\sin(\theta_i),
\end{equation}
prompting the need for the convergence of the background flow to a steady state prior
to the inclusion of the liquid droplets into the computational domain. For all
cases considered in the present investigation, an evolution of the flow spanning $100$ dimensionless time units
proved more than sufficient for this purpose, with a tolerance of $10^{-6}$
in the components of the velocity field selected to verify flow convergence to steady state.
We have confirmed this for all angles of incidence using a root mean square norm
of the velocity field, presented in Fig.~\ref{fig:fig2}b. Time $t=0$ is taken to be
the time at which the drop is seeded inside the domain and by this convention
the direct numerical simulations begin at $t=-100$.

\begin{figure}[!hbtp]
\centering
\includegraphics[width=1.0\textwidth]{1_BackgroundView.png}
\caption{Snapshots of the converged background velocity field obtained as a result of imposing uniform flow boundary conditions at an angle
of $\theta_i = 60^{\circ}$ on the upper boundary, hitting a flat solid surface at the bottom, with outflow conditions on all lateral boundaries.
The two cross-sections through a) the central $x-y$ plane and b) the $y-z$ plane illustrate the vertical velocity field (negative above, zero due to no-slip on the surface),
as well as streamlines of the flow. The water drop, shown in white, is initialised on the dividing streamline near the upper border of the geometry.}
\label{fig:fig1}
\end{figure}

Focusing on the mid-$(x,y)$-plane (at $z=0$, see Fig.~\ref{fig:fig1}a), we find similar flow properties to the classical
case of oblique stagnation point-flow in two dimensions.
Using the typical definition for the stream function $\psi(x,y)$, where the horizontal velocity component $u=\psi_y$
and the vertical velocity component $v=-\psi_x$, sufficiently far away from the wall the flow takes the form
\begin{equation}
 \psi(x,y) = k x y + \frac{1}{2} \zeta y^2.
\end{equation}
This is effectively a superposition of irrotational stagnation-point flow of strength $k$ and a uniform shear flow parallel to
the solid surface (in the $x$-direction), where $k$ and $\zeta$ are scale constants (see \cite{Blyth1} for a recent exposition on this scenario).
$\psi = 0$ denotes the dividing streamline
onto which the liquid drop is initialised just below the upper boundary of the three-dimensional domain, with its center at $y = 19.25$ and $z=0$
and with $x$ varying as a function of the angle of incidence $\theta_i$ of the background flow.
Fig.~\ref{fig:fig1}a provides a visualisation of the converged flow field at the instance
of the initialisation of the drop for the case when $\theta_i = \pi/3$.

We underline that, despite the flow being essentially two-dimensional in the upper part of the domain (below the inflow boundary),
due to the presence of the solid surface and the lateral outflow condition, it develops its full three-dimensional structure
close to the impingement region, with a single stagnation point being present in the flow irrespective of the impingement angle.
This is best observed in Fig.~\ref{fig:fig2}a, but also in Fig.~\ref{fig:fig1}, where streamlines
are drawn on top of velocity fields (illustrated in colour) plotted in different two-dimensional cross-sectional planes to indicate the deflection in the air flow.

\begin{figure}[!hbtp]
\centering
\includegraphics[width=1.0\textwidth]{2_BackgroundValidation.png}
\caption{a) Magnitude of velocity vector on an $x-z$ plane immediately above the surface, at $y=0.01$. A single stagnation point is visible in the center
of the computational box, with streamlines aiding the visualisation of the flow as it increases in velocity towards the lateral boundaries. b) Root-mean-square
norm of the velocity vector for different angles of incidence $30^{\circ} \le \theta_i \le 90^{\circ}$ of the background air flow.}
\label{fig:fig2}
\end{figure}

Once the background air flow has reached its steady state, the initially spherical liquid drop
is prescribed to enter the computational domain at a desired location $(x_i,y_i = 19.25, z_i=0)$.
The drop then inherits the local velocity field of the background air flow which is
an approximately uniform flow directed towards the surface at an angle $\theta_i$,
and is advected towards the solid surface.
The droplet shape is subject to physical deformations up to the time of its impact.
Full hydrodynamic coupling determines its trajectory and shape, with no further assumptions being made beyond this point.

\section{Numerical Methodology}
\label{sec:Numerics}

The numerical computations in the present study have been carried out using the open-source package \textit{Gerris} \cite{popinet1,popinet2}
(\href{http://gfs.sourceforge.net/}{http://gfs.sourceforge.net/}), which has been used extensively with great success by
the multi-phase flow community over the last decade. The package is ideal for our purposes since it
accurately solves the incompressible Navier-Stokes equations (and a variety of additional multi-physics extensions) using
the finite volume method and a volume-of-fluid approach to account for fluid-fluid interfaces. The schemes
are second order accurate in both space and time, with strong adaptive mesh refinement capabilities ensuring
the computational cost remains relatively low even in challenging multi-scale contexts such as those
in the present problem.
In the following paragraph we elaborate on some of the specific
measures used to ensure a good numerical performance, and also underline the
overall features of our extensive computational effort.

The large density ratio (recall that for water-air flows $r=826.281$) between the fluids may cause convergence issues
for multi-grid Poisson solvers as the one used in
\textit{Gerris} \cite{tryggvason1}, causing slow convergence or leading to a breakdown of the numerical solution altogether.
A smoothing operator/filter has been proposed (\cite{popinet2,fuster2013}) in order to alleviate
this. Spatial filtering consists of averaging over the corners
of a computational cell (four in 2D and eight in 3D), which are in turn obtained
by averaging the centered values of the corner neighbours. Applying the filter effectively
smoothes the representation of the interface over a larger number of cells and can be applied
any number of times, although previous investigations on drop impact argue
that a single iteration of the filtering operator is sufficient \cite{thoravalThesis}.
As a result of this manipulation, the errors are maintained at a reasonable (and controllable) magnitude, while the convergence
properties of the solver are much improved.

The qualities of the package, in particular in terms of adaptive mesh refinement (AMR),
become evident in the study of the problem of drop impact at high velocities. The background
air flow requires a strong level of refinement close to the surface of the solid body
to account for the presence of the developing boundary layers around the stagnation point of the flow.
At the same time, capturing the evolution of the fluid-fluid interface demands an appropriate
resolution, enabling possible topological transitions. Splashing
entails the creation and subsequent tracking of a large number of secondary droplets, which may
or may not coalesce with other bodies of fluid. In addition, suitable choices
for refinement with respect to sharp changes in vorticity were also implemented.
We also note the more stringent treatment required during touchdown,
in which a reduced timestep and an extended local refinement region is necessary to avoid
numerical artefacts.
The computational gain when comparing to the case of a uniform mesh is remarkable.
The large domain would require $\mathcal{O}(10^{10-11})$ grid cells at the finest resolution,
however with the use of the adaptive mesh refinement this is decreased several
orders of magnitude down to $\mathcal{O}(10^{6})$ degrees of freedom, which becomes significantly more tractable.
Many of the results presented would have reached considerable runtimes (as well as challenging
memory and data storage requirements)
without the usage of adaptive mesh refinement, and possibly making many of the calculations presented here unrealisable.

We also employed the functionality to selectively eliminate droplets and bubbles whose dimensions are below
a threshold number of grid cells (chosen to be $16$), thus fixing the minimum lengthscale that computations can account for.
Note that this is already well within the sub-micron scale. This feature becomes useful
when secondary droplet break-off is violent and causes the fragmentation of the fluid
into droplets of a very small size which suffer from geometrical reconstruction errors as a result of
them spanning a small number of grid cells in each dimension. In practice the technique works by replacing
the connected volumes (under
a specified size) containing the drop fluid phase (water in our case) with the background fluid (air).
Furthermore, in our implementation complete droplet removal takes place if the droplets
are found within one spatial unit of lateral boundaries in order to limit numerical artefacts
when encountering the outflow region or sufficiently high above the surface of the solid
body ($y>10.0$) to avoid high speed satellite droplets reaching the inflow boundary
and causing numerical instabilities. In practice, the mentioned situation can be avoided
by prescribing a larger computational domain that demands increased computational costs.
The selective removal of droplets
ensures that a geometry of manageable size can still be used reliably.
The flow in the vicinity of the impact region is unaffected by this treatment restricted to the near-boundary regions,
hence no flow information is artificially lost.

Many of the problems of interest require the treatment of a triple contact point
between the solid surface and the liquid-gas interface. We note that the default mesh-dependent static contact
angle model with a selected value of $90^{\circ}$ is used here. The limitations of this basic method, as well
as proposed improvements have been recently discussed by
Afkhami et al. \cite{afkhami1}, who introduced a versatile dynamic contact angle model (implemented
in an extension of \textit{Gerris}).
In a related context, Pasandideh-Fard et al. \cite{Fard1} note that
the inertially dominated stages of the flow are unaffected
by changes in the contact angle, which had been altered with the use of surfactants in their investigation.
In general, the suitability of the static contact angle model in the inertia-dominated
spreading regime has been studied extensively (\cite{yokoi2009numerical,guo2016investigation}) and
the present choice is not restrictive. We have experimented numerically
with both grid sizes and different imposed static contact angle values in two dimensions, confirming
that in the early stages of the impact we are in a regime which is insensitive to the choice of
contact angle at the wall.

The runs in the present study have been performed at
multiple resolution levels, varying from $2^{10}$ to $2^{12}$ grid cells per spatial
dimension in each computational box.
As the interfacial shape is set to be resolved at this level, this would translate to up to approximately
$200$ cells per diameter for the initial spherical drop.
Before impingement we do not require such levels of refinement away from interfaces.
On the other hand, right before, during and after impingement, the entire liquid volume demands
a strong refinement level. Despite these stringent requirements, with the chosen settings
and drop sizes, each finest resolution cell spans from $0.097$ microns
for the smallest impinging drops studied to $1.15$ microns in the case of the largest drops of initial diameter of just over $230$ microns. These levels
have been selected to provide as much detail at the micron and sub-micron levels as possible.
Many features, such as for example the minimal film thickness arising as a result of the spreading
of a drop on the surface, have well-established sizes which are useful guidelines for
what scales need to be captured and are used for comparisons and validation.
From a more general standpoint, for the decision on a suitable level of refinement
and mesh adaptivity setting we have relied on three main criteria: a) mass conservation; b) changes
in defined metrics such as velocity and vorticity norms, secondary drop size distribution etc. and c) comparisons
to available analytical predictions and experimental data in the literature. Once all three criteria have been met,
the configuration in question was propagated towards full parameter studies over the variables of interest.
We emphasise that for the top two levels of refinement, volume conservation is accurate to within
$1\%$ across the entire set of tests in the present work, with only the most challenging of test cases (the largest initial drop diameter prescribed) causing
errors of the order of $3-4\%$ due to a combination of the difficult conditions (high Re/We) and the selective drop removal mechanism introduced above,
with smaller scale features being more frequent in this scenario.
A typical computation under these conditions requires in excess of $2\times 10^{3}$ CPU hours for the lowest resolution
tested and approximately $10^4$ CPU hours for the more refined cases.

In the next section we describe, in turn, our computational results for pre-impact deformation and post-impact dynamics,
concentrating on both fundamental phenomena and aspects related to the larger scale system itself.

\section{Results}
\label{sec:Results}

Once the drop is initialised in the steady background flow, it travels towards the surface
guided by an initially uniform (but $\theta_i$-dependent) velocity field, with streamlines deforming as the solid surface is approached.
Analysing the deviation from the initially spherical shape as a function of time is
one of the primary goals of the present work, since, as noted earlier, this effect
is often overlooked in standard water retention calculation models.

In order to provide a suitable
validation framework for the present results, we have tailored the parameters to coincide
with a subset of the data of the only experimental investigation of the pre-impact deformation and break-up phenomena we are aware of - see
Vargas et al. \cite{vargas2012mechanism} and Sor et al. \cite{sor2015modeling}.
Therein, an experimental setup consisting of a monosize droplet dispenser,
a rotating arm with a model wing fixed at its end,
as well as associated motor and camera equipment are used to capture the drop dynamics as
the solid body approaches the liquid droplets at velocities of up to $100$ m/s. As a result of the very violent
regime, the size variation for the drops is restricted to $D=300\ \mu$m and above. Very few pixels
per drop diameter are visible below this threshold and the resulting images can no longer
be comprehensively analysed. As a consequence,
in the results that follow we have selected three values within the respective range, as well as one smaller drop,
typical of the sizes found in the high liquid water content regions aircraft travel through.
The drop sizes, as well as all other associated dimensionless parameters are
summarised in Table~\ref{tab:dropDeformationParameters}, where we underline
that we have used a reference velocity of $U_\infty = 90$ m/s (the same as in the main series of experiments \cite{sor2015modeling})
and the physical properties of water and air at relevant near freezing temperatures.

\begin{center}
\captionof{table}{\label{tab:dropDeformationParameters} Relevant dimensionless parameters in the case of pre-impact deformation studies
in high speed conditions, matching in drop diameter to a subset of the studies performed by Sor et al.\cite{sor2015modeling}.}
\begin{tabular}{|c|c|c|c|c|c|c|}
\hline
 $D$ [m]& $\textrm{Re}= \rho_l U_{\infty} D/ \mu_l$ &  $\textrm{We} = \rho_l U_{\infty}^2 D/\sigma$ & $\textrm{Oh} = \sqrt{\textrm{We}}/\textrm{Re}$& $\textrm{Ca} = \mu_l U_\infty/\sigma$ & $\textrm{St} = \mu_g/(\rho_l D U_\infty)$  \\  \hline \hline
 $128 \times 10^{-6}$ & $8653.717$ & $10936.183$ & $0.012$ & $1.263$ & $1.803 \times 10^{-6}$ \\
 $362 \times 10^{-6}$ & $24473.794$ & $30928.893$ & $0.007$ & $1.263$ & $6.376 \times 10^{-7}$ \\
 $634 \times 10^{-6}$ & $42862.943$ & $54168.282$ & $0.005$ & $1.263$ & $3.640 \times 10^{-7}$ \\
 $1048 \times 10^{-6}$ & $70852.309$ & $89539.999$ & $0.004$ & $1.263$ & $2.202 \times 10^{-7}$ \\
\hline
\end{tabular}
\end{center}

Following impact itself, depending on the relevant parameters, the drop is anticipated to either spread due to its momentum
and subsequently recede under surface tension effects, or, in the cases of the larger drops, to splash and
break up into secondary droplets which move away from the surface but may later re-impinge.
Droplets found in the atmosphere typically lie within the interval of $20-250\ \mu$m in diameter
and as a consequence water catch studies reported in the literature \cite{papadakis2003,papadakisExp,wright2005validation,wright2006,honsek2008,bilodeau2015}
are found in this regime.
We consider four test cases ($D=20, 52, 111, 236\ \mu$m) for complete analysis of pre- and post-impact dynamics,
in order to facilitate comparisons with results in the field and provide further insight under
flight conditions of practical interest. The complete list of parameters is provided in
Table~\ref{tab:dropImpactParameters}, where the same water-air configuration is used,
however this time with reference velocity $U_{\infty} = 78.44$ m/s, in agreement with
datasets discussed in previously mentioned studies.

\begin{center}
\captionof{table}{\label{tab:dropImpactParameters}  Relevant dimensionless parameters in the case of long-time drop impact direct numerical simulations
in high speed conditions, matching in median volumetric diameter to a subset of the studies performed by
Papadakis et al. \cite{papadakis2003}.
The splashing parameter $K = \textrm{We} \sqrt{\textrm{Re}}$ varies between $ 6.283 \times 10^4$ and $2.547 \times 10^7$.
}
\begin{tabular}{|c|c|c|c|c|c|c|}
\hline
 $D$ [m]& $\textrm{Re}= \rho_l U_{\infty} D/ \mu_l$ &  $\textrm{We} = \rho_l U_{\infty}^2 D/\sigma$ & $\textrm{Oh} = \sqrt{\textrm{We}}/\textrm{Re}$& $\textrm{Ca} = \mu_l U_\infty/\sigma$ & $\textrm{St} = \mu_g/(\rho_l D U_\infty)$  \\  \hline \hline
 $20 \times 10^{-6}$ & $1352.143$ & $1708.779$ & $0.031$ & $1.263$ & $1.154 \times 10^{-5}$ \\
 $52 \times 10^{-6}$ & $3515.573$ & $4442.824$ & $0.019$ & $1.263$ & $4.438 \times 10^{-6}$ \\
 $128 \times 10^{-6}$ & $8653.717$ & $10936.183$ & $0.012$ & $1.263$ & $2.079 \times 10^{-6}$ \\
 $236 \times 10^{-6}$ & $15955.291$ & $20163.588$ & $0.009$ & $1.263$ & $9.779 \times 10^{-7}$ \\
\hline
\end{tabular}
\end{center}

For the smallest $20\ \mu$m drops we consider an extensive parameter study in terms of impingement angles $10^{\circ} \le \theta_i \le 90^{\circ}$
in increments of $10^{\circ}$. This enables a detailed analysis of the effects related to the competition between inertial and capillary regimes,
while noting the influence of the background flow on the drop dynamics. For the more challenging larger droplets we focus on two
specific cases, namely $\theta_i = 60^{\circ}$ and $\theta_i = 90^{\circ}$, guiding us towards results in both
symmetric and asymmetric impact, described in full detail in subsection~\ref{subsec:postimpact}.

\subsection{Pre-impact dynamics}
\label{subsec:preimpact}

In the present subsection we describe qualitative and quantitative features related
to the motion of droplets prior to them impacting the solid surface.
Intuitively we expect the most deformation
and possible break-up to happen close to the solid surface as the air flow slows down and the droplet
encounters developing boundary layers. We note however that, particularly for large drops,
a rich dynamics characterised by so-called bag break-up and rupture can be observed
even in the case of simple uniform flow and in the absence of any streamline deflection \cite{Jalaal1}.
These strongly time-dependent morphological changes underline the importance of one of
the parameters in the presented model, namely the initial position of the drop relative to the solid surface.
If prescribed too far away from the surface, the initial spherical drop may become completely fragmented
by the time it reaches the surface, while seeding it too close to the surface may not allow sufficient
time for its natural dynamics to occur before impingement. As such, the comparison to the experimental results
from INTA/NASA \cite{vargas2012mechanism, sor2015modeling} serves as an important validation step.
The authors focused on describing
and modelling the change in shape, as well as the consequences thereof in terms of predicting
the drag coefficient of the evolving shape. They found that for moderate-sized droplets (with diameters in the hundreds
of microns) the approximation of the shape as an oblate spheroid proves to be reasonably accurate,
quantifying this deformation as $a(t)/R$, where $a(t)$ denotes the evolving
major semi-axis of the spheroid, normalised by the initial drop radius. This value was reported
to increase smoothly from $1.0$ as the drop is sufficiently far away from the surface to
values in the range of $1.3$ for $D=362\ \mu$m, to $2.0$ for $D \approx 1$ mm, increasing monotonically as a function
of the size of the initial drop. As they approach the surface,
the larger drops suffer considerable deformations in which the symmetric framework postulated before
is no longer applicable. Finally, when close to within $10$ mm of the solid surface, the drops violently
break up into a cloud of secondary droplets which can only be described qualitatively in the experiments.

Example evolutions of the droplet shapes are shown in panels a) and b) of Fig.~\ref{fig:fig3},
in which we analyse the deformation of a relatively small drop ($D=362\ \mu$m),
as well as a large drop ($D=1048\ \mu$m) alongside their experimental counterparts.
In the former case, we find that the proposed mild deformation into an oblate spheroidal
shape is recovered and good qualitative agreement with the experiments is found. The same
applies for the latter larger drop case, in which the flattening of the shape is much more pronounced
and asymmetric features arise in the latter stages. Note how the center of gravity
of the shape shifts towards
the lower part of the drop in the third subimage, only to develop secondary structures around the
edges which ultimately rupture from the main shape and break off into smaller droplets prior to impact.
It should be noted that there is a difference in timescales when comparing the experimental
and computational results; in the experimental data the deformation takes place
over a distance of several hundred drop diameters, whereas in all computational results this evolution
takes place within the prescribed distance of roughly $20$ initial drop diameters. The flow field
and its extensional nature is effectively scaled down to the size of the computational box.

\begin{figure}[!ht]
\centering
\includegraphics[width=1.0\textwidth]{3R1_DeformationSnapshots.png}
\caption{Pre-impact drop deformation visualisation for spherical drops of diameter a) $D=362\ \mu$m and b) $D=1048\ \mu$m.
Inside each panel the left images are experimental results by Sor et al \cite{sor2015modeling}, while the right images are the
corresponding DNS results. The images are reproduced with permission by Instituto Nacional de T\'{e}cnica Aeroespacial.
c) Quantification of the drop deformation in terms of the drop semiaxis $a$ normalized by the initial radius $R=D/2$, with the corresponding parameters
described in Table~\ref{tab:dropDeformationParameters}. The timestep
at which the drop first touches the solid surface is also highlighted with an open circle.}
\label{fig:fig3}
\end{figure}

From a quantitative perspective, for comparison purposes we use the same semiaxis deformation metric $a(t)/R$
in Fig.~\ref{fig:fig3}c to uncover an excellent agreement with the experimental data. We mark the time of impact
with an open circle and note that the obtained values are within $10\%$ of their experimental counterparts,
while the evolution of this measurement in time also shows the same features. Notably, for the larger drop
we plot the full extent of the liquid volume (accounting for the shedding of secondary droplets).
If these are to be excluded, at a distance of half a diameter above the solid surface, the deformation
is found to be $1.36$, $1.72$ and $1.85$ for the $362\ \mu$m, $634\ \mu$m and $1048\ \mu$m drops, respectively,
with approximately $1.3$, $1.7$ and $1.94$ being the equivalent values in the experiment. The inclusion
of secondary drops becomes visible around $t \approx 15$ in both cases
and causes an increase in this metric to just below $2.0$ and $3.0$ for the two largest droplets,
indicating the complexity of the flow in the respective regimes as the drops approach the surface.

The computational framework developed here can be used to access information on the flow field
and drop shape at distances very close to the surface that are beyond the frame-restricted capabilities
of current powerful video technology.
Consequently, we consider the case
of a smaller drop of initial diameter $D=128\ \mu$m and find very small deviations
from the imposed shape during its entire evolution. A small initial flattening of the shape into an oblate spheroid
suffers corrections prior to impact and ultimately impinges almost undeformed.

\begin{figure}[!h]
\centering
\includegraphics[width=1.0\textwidth]{4_ImpactSnapshots.png}
\caption{Initially spherical droplets of diameter a) $D=128\ \mu$m, b) $D=362\ \mu$m, c) $D=634\ \mu$m, and d) $D=1048\ \mu$m, at the moment of impact
onto a flat solid surface, having been deformed by the background stagnation point flow. The smallest drop retains its shape,
while the edges of the largest drop break up into a large number of secondary droplets even before impact. This compares favourably
to e) previous experimental investigations of drop deformation prior to impacting a moving solid body ($D_{\textrm{exp}}=1048\ \mu$m as well) by Vargas et al \cite{vargas2012mechanism}.
The last image is reproduced with permission by Instituto Nacional de T\'{e}cnica Aeroespacial.}
\label{fig:fig4}
\end{figure}

For completeness, all four cases are illustrated in Fig.~\ref{fig:fig4} at the last computed timestep before touchdown, with the last
$D=1048\ \mu$m case being placed side-by-side with its experimental counterpart.
For the smallest drop, deformation is hardly visible
(as confirmed by Fig.~\ref{fig:fig3}a), with an approximately spherical liquid volume impinging onto the surface. A strong
flattening of this shape with the beginning of breakup features becoming visible around the edges, takes place for slightly larger drops
and this ultimately leads to progressively smaller liquid fragments/drops being shed from the sides. In the largest drop volume case,
the cloud of droplets behind the main liquid volume becomes visible and resembles the experimental result.

In what follows we focus on the impingement process itself and in particular on the spreading or splashing
characteristics of the flow, as well as the associated secondary drop formation and dynamics.

\subsection{Post-impact dynamics}
\label{subsec:postimpact}

Once the drop approaches the region very close to the wall, the gradually thinner
air film below is forced to move away laterally and the pressure underneath the droplet continues
to grow until impact takes place. We note the presence of either a single or multiple
air bubbles entrained under the surface. In the classical context with a quiescent air flow
and small to moderate impact velocities, the size and evolution of the air bubble is well studied,
and its effect on the splashing process itself has been shown to be negligible \cite{riboux2014experiments}.
In the present case however, there are two fundamental differences from the traditional impact problem
due to the very high impact velocity, as well as the strong pre-impact drop deformation,
particularly in the oblique impingement cases. In the first instance and on the basis of Fig.~\ref{fig:fig5},
we will provide a qualitative assessment of the results. The studied parameter space consists of
the two cases outlined in Table~\ref{tab:dropImpactParameters}
of normal impact and oblique impact at $60^{\circ}$, and four different drop diameters, producing eight
studies in total. The time sequence of top views of the liquid drop on the left hand side of Fig.~\ref{fig:fig5}, illustrates the drop
shape at
three key times
in its evolution, namely:
\begin{enumerate}[i.]
 \item when the center of mass of the drop lies one initial diameter above the surface (left column);
 \item when the drop reaches its maximum spread on the surface and before retraction under capillary forces takes place (middle column);
 \item ten time units later, which serves as an indication of how the longer timescale of the impact develops
 into either retraction for the smaller drops or violent rupture and splashing for the larger drops (right column).
\end{enumerate}
For the $90^{\circ}$ impact case we concentrate on the second of the above time instances, namely when the smallest
drop reaches its maximum diameter - results are shown in the right hand side column of Fig.~\ref{fig:fig5}.
In each image a reference
lengthscale of $20\ \mu$m is added as a visual aid to the extent of the drop atomisation (or lack thereof).

\begin{figure}[!ht]
\centering
\includegraphics[width=0.9\textwidth]{5R1_TopView.png}
\caption{Splashing dynamics for drops of sizes $D = 20,\ 52,\ 111$ and $236\ \mu$m (each row represents a different drop
size, with the complete list of parameters defined in Table~\ref{tab:dropImpactParameters}) at an angle
of incidence of $\theta_i = 60^{\circ}$. The left column illustrates the drop shapes as their center of mass is at $y=D$ above the surface,
the images in the second column are plotted at the dimensionless timestep at which the
drop illustrated on the top row (smallest drop, with initial diameter $D=20\ \mu$m, and relatively regular spreading behaviour)
reaches its maximum spread $t_s^{\textrm{max}}$, while the third column shows
the drop shapes ten time units later, once either retraction or more pronounced splashing has occurred. The rightmost column is used
to visualise the splashing for the $\theta_i = 90^{\circ}$ impact case at $t_s^{\textrm{max}}$.}
\label{fig:fig5}
\end{figure}

The smallest drop size (initial diameter $D=20\ \mu$m) impingement is characterised by
inconsequential pre-impact deformation with the approximately spherical shape retained
up to very near the time of impact, followed by a strong spreading motion in a highly inertial regime,
finally followed by retraction due to surface tension. Intriguing corner-type features emerge particularly
for the oblique impact cases due to the directionality of the impact, which will be discussed in detail
in subsequent paragraphs.
Referring to the oblique impact scenario, the asymmetry becomes more visible for medium sized drops (at the order of $100\ \mu$m in initial diameter)
prior to impact and particularly after impact as fluid volumes have sufficient momentum to
overcome surface tension and push outside the typical nearly circular contour, instead spreading laterally outward
towards the front of the drop. The dynamics is however still dominated by one large fluid volume
from which small secondary drops are ejected as the drop increases in size. We underline that the imposed angle of incidence
has a clear influence on the angle and extent of lateral spread of the liquid mass.
The largest drops ($D=236\ \mu$m) experience violent splashing, with visible liquid threads forming in the forward and laterally outward directions
as the fluid mass disintegrates into hundreds of droplets.
Similar features are observed in the normal impact case in terms of fragmentation, with
traditional spreading motion transitioning to azimuthal instabilities, followed
by a rupture of the liquid rim into small drops, but with a main fluid mass still intact near the impact site.
Ultimately a violent fragmentation breaks up the liquid volume into thin filaments near the surface,
and numerous secondary drops are advected away from the impact region under the influence
of the background flow.

Conducting a systematic analysis of the drop's morphology during the early and intermediate
stages of the impact is most accessible for the smallest drops (below several tens of microns in initial diameter,
top row of Fig.~\ref{fig:fig5}, when no splashing occurs),
where early and very recent analytical results are available for comparison when
$\theta_i = 90^{\circ}$. Following this baseline, the generalisation to the predominantly three-dimensional
effects of the asymmetric impact are best constructed. Even in the normal impact case however, the presence of
the non-quiescent air flow at high speeds is anticipated to produce some modifications in the standard metrics
surrounding the characterisation of the impingement process, which will be emphasised in the following paragraphs.

In order to aid future comparisons, in Fig.~\ref{fig:fig6} we define several quantities of interest, namely
the time-dependent drop diameter
in the $x-$direction $D_x(t)$ (the direction of impact for the non-normal incidence cases),
the drop diameter in the transverse $z-$direction $D_z(t)$, as well as the height of the drop near
its center of mass $h_f$. The first two metrics are best observed from the top view ($x-z$ plane) presented in the
top part of panel a), while a cut through the $x-y$ plane provides information on the minimum thickness of the film.
The entrapment of a small air bubble due to impact cushioning, results in a small variation
in the drop's curvature just above this feature, which is why for the relevant local minimum
we select a point where this local adjustment is negligible.
Several notable studies (see Introduction) have addressed the topic of the maximum spread $D_m$
of the drop in normal impact conditions, with the recent investigation of Wildeman et al. \cite{wildeman2016spreading}
chosen as reference here. Plugging our parameters into their main result, we find $D_m \approx 4.112$,
which compares very well with the computed value for $\theta_i = 90^{\circ}$ in Fig.~\ref{fig:fig6}c.
Symmetry in this case is preserved and we find $D_x^{\textrm{max}}  = D_z^{\textrm{max}} \approx 3.9$, which alongside the good agreement
also indicates that the surrounding flow has a limited influence on the maximum spread.

As the angle of incidence is decreased down to $\theta_i = 30^{\circ}$, the flow speed in conjunction with
the increasingly pronounced directionality of the impact
enables the liquid mass to advance towards the front side (in the $x-$direction) of the impacting drop,
pushing more strongly towards the front edge and increasingly distorting it in this direction. Fig.~\ref{fig:fig6}c
indicates this monotonic increase in $D_x^{\textrm{max}}$ and decrease in $D_z^{\textrm{max}}$ as
$\theta_i$ is reduced, with the final aspect ratio being measured at almost a factor of two.
We point out that in this regime the drop is also subjected to a stronger air flow as it lies
further away from the dividing streamline and the background flow velocity has an increased magnitude.
For illustrative purposes, in Fig.~\ref{fig:fig6}b we expand on how the maximum diameter values
are obtained in the asymmetric cases, with the two diameters $D_x(t)$ and $D_z(t)$ being shown throughout
their evolution for an angle of incidence $\theta_i = 60^{\circ}$. The dynamics in the $x-$direction
is chosen as reference, as this is the dominant motion due to our choice in impact directionality.
The value of $D_z^{\textrm{max}}$ is then defined as the value of $D_z(t_x^{\textrm{max}})$,
where $t_x^{\textrm{max}}$ it the timestep at which $D_x$ reaches its maximum, despite
it not necessarily being the highest absolute value in the $z-$direction. The figure
shows negligible deformation up to the time of impact $t \approx 20.0$, followed by a sharp increase
in diameter in both directions but more strongly in $x$, with the rim finally retracting under
the effect of surface tension from all directions. The reference values (see vertical dashed line)
derived from similar studies of each incidence angle are then used to construct Fig.~\ref{fig:fig6}c.

\begin{figure}[!hbtp]
\centering
\includegraphics[width=1.0\textwidth]{6R1_SpreadingMetrics.png}
\caption{Spreading dynamics of microdroplets with an initial diameter of $D = 20\ \mu$m at angles of incidence ranging from $30^{\circ}$ to $90^{\circ}$.
a) Top view schematic of the spreading diameter in the $x-$direction (the impingement direction) and the $z-$direction, as well
as the minimal film thickness $h_f$ (in side view below), for an angle of incidence of $\theta_i = 60^{\circ}$.
b) Evolution in time of the spreading and later retracting liquid drop for the $60^{\circ}$ impingement
angle case. c) Summary of the maximum spread in both $x$ and $z$ for a collection of angles of incidence, indicating the transition from
symmetric spreading to a strongly asymmetric final shape in the direction of impact.}
\label{fig:fig6}
\end{figure}

Another key morphological metric we consider is the minimum film height, as extracted near
the drop center, sufficiently far away from the entrapped bubble.
In the normal impact case and in the strong inertial regime described here,
Eggers et al.\cite{Eggers1} estimate this thin film height to reach a minimum $h_f/R \approx \textrm{Re}^{-2/5}$,
which would give $h_f \approx 0.028$ in our case. This is the height at which the thinning film reaches the liquid boundary
layers within the drop itself and ceases its decrease. The result obtained in our
investigation is $h_f \approx 0.033$ and we found no evidence of significant
variation as a result of modifying the angle of incidence. The very slight
overestimation is perhaps counterintuitive given that the fast air flow pushing from above
would be expected to enhance the thinning effect. We note that even at these small
lengthscales the mesh is sufficiently fine with several gridpoints spanning
the thin film region; changes in the resolution did not result in meaningful
changes of this value.

One of the most salient features of the drop impact in the modelled high speed regime is the emergence
of a corner-type feature near the advancing front of the spreading liquid mass; this feature
becomes highly prominent, particularly as the
angle of incidence $\theta_i$ is $60^{\circ}$ or lower.
Above the respective angle, normal impact is characterised by approximately axisymmetric behaviour,
while in slightly oblique impacts ($\theta_i \approx 70^{\circ}-80^{\circ}$) the footprint can
 be described as elliptical, although a slight symmetry-breaking tilt to the front becomes observable
 on the lower side of this range.
A comprehensive analysis
of the corner-type property has been performed for angles varying in the range $10^{\circ} \le \theta_i \le 90^{\circ}$
and small drop size (initial diameter $D=20\ \mu$m), and the results are given in Fig.~\ref{fig:fig7}.
Herein we track the evolving drop shape from above, and concentrate on the moment where its spreading
diameter reaches its maximum value for each of the particular cases.
Two angular metrics are then defined as illustrated on the right hand side of Fig.~\ref{fig:fig7}:
$\varphi_t$ (angle from the tip of the advancing front in the impingement $x-$direction to the top part of the drop, the maximum in the $z-$direction)
is a more global measure of the deformation, whereas $\varphi_n$ is a local measure
of the corner angle near the tip of the advancing front, defined by a triangle whose base is fixed to be a quarter
of the initial radius $R/4$, as shown in the figure.
We emphasise that while the discussed feature is called a corner (or of corner-type) throughout this subsection,
the shape would be more accurately described as an apparent corner, since locally near the tip of the advancing
front surface tension always induces a smoothing of the shape.

\begin{figure}[!h]
\centering
\includegraphics[width=1.0\textwidth]{7R1_AngleSummarySketch.png}
\caption{a) Characterisation of the geometric feature arising at the leading (front) side of the drop due to the oblique impact.
Angle $\varphi_t$ is measured from the most advanced point of the drop in the direction of impact to the maximum in the spread
in the perpendicular direction of the same plane, while $\varphi_n$ represents the more local feature arising at $0.25R_0$ behind the
front, where $R_0$ represents the initial drop radius. Both angles are defined in panel b), while the three insets present top views
of the drop shape at the moment of maximum spread in the $x-$direction, the timestep at which all the angles in the figure are calculated.}
\label{fig:fig7}
\end{figure}

The progressively more stretched shape of the drop, as well as the evolution near the tip of the advancing front capturing
the corner-type feature itself are both embedded in the above quantities, which are presented at
the bottom of Fig.~\ref{fig:fig7}, with examples of the underlying drop shapes depicted in the row above for $\theta_i = 40^{\circ},\ 60^{\circ}$
and $80^{\circ}$.
In the intermediate case small distortions of the liquid rim are already visible, while
at $40^{\circ}$ a pronounced outgrowth near the advancing front
selected by the direction of impact is clearly identifiable.
Due to the preserved axisymmetry, at $\theta_i = 90^{\circ}$ we compare the numerical
results with simple predictions.
We naturally expect $\varphi_t \approx 45^{\circ}$ and
based on the maximal spreading radius described in Fig.~\ref{fig:fig6}c, the definition of the angle $\varphi_n$, as well as using basic trigonometry,
we estimate $\varphi_n \approx 75^{\circ}$. We recover $\varphi_t = 44.93^{\circ}$ and $\varphi_t = 74.76^{\circ}$ by analysing the data,
which is well aligned with the anticipated axisymmetric evolution.
Both angle measurements are expected to decrease in oblique impact scenarios, with
the elongation of the liquid shape gradually reducing their values as $\theta_i$ decreases.
This is indeed
the case, with a smooth monotonic variation in $\varphi_t$ finalising at approximately $22^{\circ}$ for the $\theta_i = 30^{\circ}$ impingement case.
The local angle $\varphi_n$ naturally begins at a much higher value, but again, as the impingement angle $\theta_i$ decreases,
the deformation of the spreading drop is more is enhanced and this results in a steady drop from $\varphi_n \approx 75^{\circ}$
for the normal impact case down to $\varphi_n \approx 25^{\circ}$ at $\theta_i = 30^{\circ}$. The slope characterising
this decrease becomes markedly larger in absolute value below $\theta_i = 60^{\circ}$.
Note however that due to the shape irregularity this definition becomes less practical for angles below $\theta_i = 50^{\circ}$,
as high variation is induced depending on the choice of distance defining $\varphi_n$.
In particular, the top left hand side inset in Fig.~\ref{fig:fig7} reveals the formation of a small
finger-like extension that becomes smoothed out under surface tension. The strongly varying curvature
of this shape near the tip makes it difficult to design a universally useful local metric to describe the corner,
which is why the choice of a distance of $R_0/4$ behind the advancing tip should in some sense be interpreted with
an embedded variation, as would any other choice.

Finally, a physical interpretation of the corner formation process is proposed as
a combination of two different mechanisms, one related to the background flow, the other
to the liquid movement itself. In the latter case, during the spreading motion
liquid pushes into the front of the rim in a preferential direction given by the progressively
more oblique impact as defined by $\theta_i$. There is sufficient inertia to drive more liquid mass
towards the advancing front, however this is insufficient to overcome surface tension in the lateral direction
due to the relatively small size of the drops. As such, the liquid that accumulates in the rim in the lateral regions is also steered towards
the front of the drop, where at the meeting point the conditions for a localisation of the interfacial
shape into a corner are met before surface tension relaxes this feature. At the same time,
high speed air is pushing from above in the same direction as the primary impacting motion, further guiding liquid into this
front region. This becomes far more evident for the lower impingement angles in which the drop spreads
on the surface in a region several diameters away from the flow stagnation point and as such
the local shear forces become gradually stronger, supporting the complete manifestation
of the observed corner at the advancing front of the spreading drop.

\begin{figure}[!ht]
\centering
\includegraphics[width=1.0\textwidth]{8_SideViews10.png}
\caption{Side view ($x-y$ plane) of the late time impingement dynamics of a $D = 20\ \mu$m drop at an angle of $\theta_i = 10^{\circ}$,
resulting in the break-up of a liquid volume
in the leading region. The panels on the right
illustrate a magnification of this area.
The interface is coloured in the magnitude of the horizontal velocity at the respective points, while the adaptive grid underlying
each timestep is also shown. For reference, the smallest grid cell measures approximately $0.39\ \mu$m in dimensional terms.}
\label{fig:fig8}
\end{figure}

For impingement angles lower than $\theta_i \approx 30^{\circ}$, we find entirely different phenomena
captured in Fig.~\ref{fig:fig8}. The near-glancing
incidence of the drop, coupled with a more uni-directional background air flow,
both contribute to a pinch-off near the advancing front of the liquid volume, as opposed
to the creation and subsequent relaxation of a corner-type feature. For most of
its development, the region near the front of the drop advances on top of a very thin liquid film,
with the bulk of the liquid mass eventually catching up with significant horizontal velocity (used
to colour the liquid interface in each subfigure in Fig.~\ref{fig:fig8}). The advancing front of the strongly
elongated shape retains sufficient inertia to eventually detach from the surface and subsequently
break off into several small liquid drops that progress at very high speed towards the edge of
the finite computational domain, with the majority of the fluid quantity retained on the solid surface.
One of the additional causes underlying the observed fillamentation
is that the drop lands further away from the global stagnation point in the background flow, implying that locally
the flow is predominantly moving in the direction of the spreading in the front, promoting the lateral movement as opposed to
pushing down onto the liquid.

Having discussed at length the rich features already appearing in the impact of the smallest drops,
we proceed to a quantitative study of the splashing dynamics throughout
the entire duration of the direct numerical simulations for the full range of drop sizes considered; these results are
summarised in Fig.~\ref{fig:fig9}. Recall that the drops are initialised at $t=0$,
impacting the surface at $t \approx 20.0$ and engaging in either spreading motion or generation of secondary droplets
being tracked over $80$ subsequent dimensionless time units. Detailed results are presented for
the asymmetric case with an angle of incidence $\theta_i = 60^{\circ}$ and four different
droplet sizes ranging from $20\ \mu$m to $236\ \mu$m.

Focusing on the left panel of Fig.~\ref{fig:fig9}, we notice the effect of initial drop size
on the formation and break-up of secondary drops. The smallest
drop follows the described spreading and retraction motion detailed in subsection~\ref{subsec:postimpact}
and no secondary drops are ejected from the surface. Instead, a steady spherical cap
solution is observed during the final stages. At intermediate sizes, but still below
an estimated $100\ \mu$m threshold, small drops detach from the edges of the rim
in all directions however with a stronger preference towards the direction of impact,
with approximately $20-30$ secondary droplets being swept away by the air flow and
advected towards the lateral outflow boundaries; a few re-impinge
onto the surface and remain at such respective locations.
For the largest droplets of initial diameter $D=111\ \mu$m and $D=236\ \mu$m,
a very violent splashing motion ensues (also visible in Fig.~\ref{fig:fig5}),
with several hundreds very small drops moving away from the impingement region.
As expected, due to the larger drop size the effect of surface tension is weaker and
break-up into progressively smaller liquid fragments is more pronounced.
In such cases, even the secondary droplets are subject to subsequent break-ups, with
a minimum size being again restricted by surface tension.
We have made extensive verifications of the selected grid refinement and minimum cell size
in order to restrict numerical artefacts at this level. In fact, the final resolution
for these studies was partly selected in light of a convergence to a minimal secondary
drop size captured by our simulations. We note, however, that the selective droplet removal
procedure performed primarily to avoid instabilities around the boundaries, does affect
the secondary droplet count, with some of the smallest droplets ejected as part of the prompt splash away
from the surface at very high speed being removed from the finite computational domain
within several time units after the impact.

The size and position of each individual drop is tracked after impact and hence
statistical information on the secondary drops is compiled and studied dynamically.
Two particular points in time have been selected for visualisation purposes on the right hand side
of Fig.~\ref{fig:fig9}, which represent the early post-impact stage when the maximum
number of secondary drops is found in the domain ($t_1$), followed by the point in time halfway through
the evolution at $t=50.0$ when the main body of fluid no longer ejects secondary drops in the impact region and
the secondary drops are airborne ($t_2$). Both
of these apply to the case of the largest $D=236\ \mu$m initial diameter drop,
selected due to the impact and splash producing the richest secondary drop dataset.

\begin{figure}[!h]
\centering
\includegraphics[width=1.0\textwidth]{9R1_DropAnalysis60.png}
\caption{Secondary droplet ejection characteristics as a result of spherical drops of initial diameter $D = 20,\ 52,\ 111$ and $236\ \mu$m impinging
onto a solid surface at an angle of incidence $\theta_i = 60^{\circ}$. Evolution of the number of droplets in time for each case (left), with the two panels
on the right indicating the secondary drop size distribution (normalised by the initial droplet size) for the $236\ \mu$m drop at the two different times, $t_1$ and $t_2$, highlighted in
the left panel.}
\label{fig:fig9}
\end{figure}

\begin{figure}[!h]
\centering
\includegraphics[width=1.0\textwidth]{R1_SecondaryDropEvolution.png}
\caption{Secondary droplet size distribution evolution characteristics for the largest drop
in the test batch summarised in Table~\ref{tab:dropImpactParameters} and indicated as $D_4$ in Fig.~\ref{fig:fig9}, impinging
onto a solid surface at an angle of incidence $\theta_i = 60^{\circ}$.
a) Left panel: three-dimensional view, with timesteps $t_1$ and $t_2$ detailed in Fig.~\ref{fig:fig9} marked with vertical bars.
b) Right panel: top view of the same dataset in the form of a contour plot.}
\label{fig:figR1DistributionEvolution}
\end{figure}

Once the impact has taken place and sufficiently many secondary drops have formed,
the volumes of these liquid fragments follows an approximately log-normal distribution (see detail at $t_1$),
centered around a mean of $10^{-4}$ relative to the volume of the initial drop. If assumed to be spherical
(which is seldom the case) this translates into droplets with a radius of $1/20$ relative to the radius of the initial drop.
As further fragmentation takes place due to the interaction between the fast movement of the drop
and the surrounding air boundary layer flow, a second local maximum becomes prominent,
with tens of drops with volumes of the order of $10^{-6}$ relative to the initial volume being present
far away from the impact area. Under the action of surface tension, these drops are often
more regular (spherical) in shape if still airborne, with some of them re-impinging
far away from the impact region and becoming spherical caps as in the case of the previously
studied $D=20\ \mu$m initial diameter drops. In fact, the distribution highlighted
at $t_2$ in Fig.~\ref{fig:fig9} also includes the minimum droplet volume captured within
this computation, which is found to be of just less than $10^{-7}$ or of a radius of $1/200$ relative
to the initial drop - just above $1\ \mu$m in dimensional terms (recall the initial drop diameter is $D=236\ \mu$m in this
numerical experiment).

The full evolution of the secondary drop size distribution in this case
is illustrated with two different visualisation techniques in Fig.~\ref{fig:figR1DistributionEvolution}.
We are concentrating on the case of angle of incidence $\theta_i = 60^{\circ}$, however
we have found the qualitative behaviour described below to be consistent with variation in $\theta_i$.
The observed log-normal distribution evolves (with the number of drops increasing, but the structure
being retained) over roughly $15$ time units. At this stage ($t \approx 50$), most drops in the system follow one of two
general tendencies:
\begin{enumerate}
 \item the larger drops in the distribution are the ones which are detached from the main liquid mass
 but lie on the solid surface (after early dynamics or later re-impingement) as approximately spherical caps.
 After they reach this configuration they will only increase in size as a result
 of coalescence with neighbouring spherical caps or incoming smaller secondary drops that re-impinge onto the surface.
 This region, the right hand side local maximum centered around $10^{-3}$ in normalised volume) remains relatively steady in both number of drops and extent of variation.
 \item the fragmentation process around the left hand side second local maximum in the distribution (with mean of approximately $10^{-7} - 10^{-6}$)
 is rich and spans roughly
 two orders of magnitude in normalised volume. These are primarily airborne drops that
 continue to travel, break up or coalesce as a result of the interaction with the background air flow.
\end{enumerate}
Rather remarkably, we find that the separation between the two types of drops (the two local maxima in the distribution)
 is relatively well preserved when changing both initial drop size and angle of impingement.

 The same type of data analysis as presented above may be useful in furthering our understanding
 of the dynamics very close to the splashing threshold, with possible comparisons
 with the work of Riboux \& Gordillo \cite{riboux2015diameters} becoming an interesting
 future line of investigation. In its current form, the high speed regime and presence
 of the background flow in the model make such direct comparisons difficult, however
 early calculations indicate that the first drops ejected as part of the splashing
 process lie encouragingly close to the previously mentioned predictions.

\begin{figure}[!h]
\centering
\includegraphics[width=1.0\textwidth]{R1_DropSurfaceArea.png}
\caption{Normalised liquid-gas surface area as a function of time
for four different drop sizes (summarised in Table~\ref{tab:dropImpactParameters}) and
two impingement angles: a) $\theta_i = 60^{\circ}$ and b) $\theta_i = 90^{\circ}$.
The insets in each subplot concentrate on the evolution of the surface area during the five
time units just before the moment of touchdown, marked by a vertical black dashed line in the full scale images.}
\label{fig:figR1SurfaceArea}
\end{figure}

The present work was in part motivated by the question of how extending our understanding of liquid droplets impinging onto
solid surfaces in high speed conditions acts as one of the early building blocks in the broader context of icing
prevention and aircraft safety and design.
As a result of the deformation and splashing dynamics forming the subject of the present work,
a useful metric to discuss is the evolution of the liquid-gas surface area,
which may in the future be considered in view of coupling to thermodynamic effects.
We illustrate our findings for the four cases introduced in Table~\ref{tab:dropImpactParameters}
and two angles of incidence, $\theta_i = 60^{\circ}$ and $\theta_i = 90^{\circ}$ in Fig.~\ref{fig:figR1SurfaceArea}.
To aid the discussion around the quantitative information, we introduce two simple approximations
of what could be anticipated in light of previously elucidated dynamics and focus on the normal impact scenario for clarity.

For the smallest drops studied ($D_1$ here), we find that spreading and retracting behaviour is still characteristic
and as such we can use a straightforward analogy in terms of a flat cylinder of appropriate radius and height to estimate
the corresponding surface area. In particular, considering the time of maximum spread $t_s^{\textrm{max}}$, we find
a value of the resulting $D_x^\textrm{max} \approx D_z^{\textrm{max}} \approx 3.9$ (see Fig.~\ref{fig:fig6}c).
Assuming (and having verified that) volume conservation holds, we find an approximate height for the cylinder,
which in this case is roughly $1/24$. As we do not account for the bottom of the cylinder (the side adhering
to the solid surface), the adjusted formula for surface area (accounting for the top surface and the side)
gives a normalised result $S_1^{\textrm{max}} \approx 3.83$, in
excellent agreement with the corresponding maximum in Fig.~\ref{fig:figR1SurfaceArea}b.
Now looking to the other end of the spectrum i.e. the largest drops considered in which case strong splashing is observed, the drop size distributions
in Fig.~\ref{fig:fig9} indicate that most of the drops have volumes with means of
either $10^{-3}$ or $10^{-6}$ relative to the initial drop volume. Accounting again for conservation of volume,
we would anticipate the normalised surface area to vary between $10$ and $100$ should we naively assume all drops to be airborne
and perfectly spherical. Instead we of course have a combination of these drops, many of which are also
on the surface (with the area in contact with the surface not counted), effects which would
all bring our estimate closer to the former value of $S_4^{\textrm{max}}  \approx 10$. For both impingement angles
we notice that this estimate does not deviate significantly from the obtained value.

Comparing the surface area evolution for both scenarios, we naturally find an increase
with drop diameter, however this increase appears to be much more ordered in the normal impact case in
terms of both distance between maxima and profile shapes. It should also be noted that, as the insets present,
the pre-impact deformation and hence surface area generated is more pronounced for the case with oblique incidence.
Once touchdown occurs, the initial strong increase in surface area is dominated by the spreading mass of liquid
in both cases, with fragmenting secondary drops only adding negligible features to the evolution for large initial diameter cases.
Once a maximum is reached, capillary forces cause smooth retraction for the smaller drops. This
is in contrast to the larger drop cases, in which, despite the apparent smoothness of the curves, a closer inspection reveals prominent bumps indicating
individual fragmentation or coalesence events for small drops. Ultimately capillary retraction as well as
drops exiting the finite computational domain through lateral boundaries lead to a decrease in overall liquid-gas surface area.

While much of our attention is dedicated to modelling and quantifying
aspects related to single drop impingement dynamics at high speeds, we emphasise that the produced datasets
encoding sizes and locations in space and time of all secondary drops becomes
central in determining further re-impingement events and elucidating retention properties in general.
From a broader perspective, throughout the present section we have described not only fundamental flow features,
but also intricate dynamics and a level of physical detail which is invariably omitted once highly simplified coarse-graining
procedures are considered instead.
In a practical context, this information may provide an accurate physical foundation for engineering models on larger scales
and is the subject of ongoing work in our group.

\section{Conclusions}
\label{sec:Conclusions}

In the present work, three-dimensional drop impact at high velocities has been investigated numerically
through the use of high accuracy direct numerical simulations
in order to advance our understanding of drop deformation and splashing dynamics in an aerodynamic context.
A model that encompasses the transition from larger lengthscales of typical engineering applications (of the order
of an airfoil chord or nacelle diameter) down to the local impact region, whilst accounting
for the surrounding (non-quiescent) air flow has been proposed. Using an oblique-stagnation
point flow model for the background flow provides a suitable framework for the impinging
drop to naturally interact with structures such as the growing boundary layers it would encounter
on aircraft surfaces. Pre-impact deformation and break-up is compared to available experiments
in the relevant parameter range of $\mathcal{O}(100)$ m/s impact velocities and $D=\mathcal{O}(100-1000)\ \mu$m diameter drops.
The impingement of the smallest droplets is described by a regular spreading motion,
while beyond a certain size (estimated to be of approximately $50\ \mu$m), violent splashing is observed.
Topological changes such as droplet break-off or coalescence, as well as potential
subsequent re-impingement of smaller fragments have been taken into account.

Comparison with existing analytical results is possible for the tractable spreading dynamics of the smallest
drops in the tested range (under $50\ \mu$m in diameter). Variation in the angle of impingement (not previously performed in this regime)
reveals intriguing corner-type features at the advancing impact front which warrant further investigation.
These emerge as a result of the asymmetric impact pushing more fluid mass in a preferred direction,
however with surface tension preventing break-up above a $30^{\circ}$ angle of incidence, and the action
of the strong background flow induced shear. Below this angle,
fluid filaments have been shown to form and detach from the surface near the advancing droplet front,
breaking into large droplets
and being carried away by the surrounding air flow.
For larger drops, splashing and secondary drop ejection is captured numerically in detail
and examined in order to advance understanding of the water retention process.
The average number of secondary drops resulting from the impact increases with the initial droplet diameter,
with their sizes well-represented by a log-normal distribution, with a secondary local
maximum emerging as a result of further break-up at the final stages of the simulated dynamics.

These findings provide detailed insight into the highly complex fluid dynamical
processes occurring during aircraft flight through high liquid water content regions.
The present approach and simulations are a significant advance of
standard particle-based methods which are common industrial practice
and which rely heavily on semi-empirical arguments. The modelled background air flow
provides a reliable local description of the flow in which a full interaction between
the air and the liquid is permitted and the deformation of the drops is captured in detail before
impact, while the emergence and movement of secondary drops in an active flow region
is also treated realistically. Significant efforts have been made in order to ensure
a highly accurate resolution of the flow. Nevertheless, we emphasise that these results
originate from very intensive and resource heavy computational efforts (both in terms
of runtime and data storage and processing requirements). Further advances in this area
of research in terms of both algorithms and raw computing power will facilitate more understanding
of these violent regimes, with access to smaller grid sizes, larger domains (the issue of locality)
as well as sensitivity to contact angle dynamics, being only a small subset
of the possible future directions within this methodology.

In conclusion, this research provides a renewed perspective on the modelling of water
catch on aircraft surfaces, with possible ramifications towards other areas involving
high speed drop impact, such as inkjet printing, combustion and agricultural sprays. The presented results have shown very favourable agreement
with recent experimental and analytical results, where possible, in an incredibly
challenging regime, whilst new phenomena and detailed quantification of practical
information beyond the capabilities of present video technology and analytical treatments
has also been provided.
We believe that the proposed numerical framework is a valuable tool not only
from the fundamental perspective in the study of drop impact, but also in an engineering context as
a means of using scale transition to include detailed physical and fluid-related processes
in water retention estimation and associated phenomena such as icing.

\section*{Acknowledgement}

The authors acknowledge the support of Innovate UK through the SANTANA (System Advances in Nacelle Technology AerodyNAmics, project reference 113001) program.
The authors were also partly supported by EPSRC grants EP/K041134/1 and EP/L020564/1.
We are grateful to Richard Newman and Hui Yao from Bombardier Aerospace Ltd. for fruitful practical discussions.
Finally, we would like to thank the anonymous referees for their insight and thoughtful comments.

\section*{References}