# A Lagrangian meshfree model for solidification of liquid thin-films

Anand S Bharadwaj***Corresponding author<EMAIL_ADDRESS>Elisa
Thiel  Pratik<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

In this paper, a new method to model solidification of thin liquid films is
proposed.  This method is targeted at applications like aircraft icing and
tablet coating where the formation of liquid films from impinging droplets on
a surface form a critical part of the physics of the process. The proposed
model takes into account the (i) unsteadiness in temperature distribution,
(ii) heat transfer at the interface between the solid and the surface, (iii)
volumetric expansion/contraction and (iv) the liquid thin-film behaviour, each
of which are either partly or fully ignored in existing models. The liquid
thin-film, modeled using the Discrete Droplet Method (DDM), is represented as
a collection of discrete droplets that are tracked in a Lagrangian sense. The
height of the liquid film is estimated as a summation of Gaussian kernel
functions associated with each droplet. At each droplet location, a solid
height is also computed. The evolution of the solid height is governed by the
Stefan problem. The flow of the liquid thin-film is solved just as in the case
of DDM, while also taking into consideration the shape of the solidified
region lying beneath the droplet. The results presented in this work show the
reliability of the proposed model in simulating solidification of thin-films
and its applicability to complex problems such as ice-formation on aircraft
wings.  The model has been verified for canonical problems that have
analytical solutions. For the more complex problems of icing, the results of
the model are compared with data from literature, without considering a
background air flow. The comparison can be improved by coupling this model
with suitable air flow solvers, as shown in the final test case.

###### keywords:

Solidification , Thin-film flows , Stefan problem , Droplets , Smoothed
Particle Hydrodynamics , Meshless methods

[inst1]organization=Fraunhofer ITWM,addressline=Fraunhofer-Platz 1,
city=Kaiserslautern, postcode=67663, country=Germany

[inst2]organization=University of Luxembourg,addressline= 2 avenue de
l’universite, postcode=L-4365, state=Esch-sur-alzette, country=Luxembourg

[inst3]organization=Indian Institute of Technology Delhi,addressline= Hauz
Khas, postcode=110016, state=New Delhi, country=India

##  1 Introduction

Liquid thin-films are commonly found in lubrication [1], medicinal coatings
[2, 3], and in the automotive industry [4], among many other applications. In
several of these applications, the liquid thin-films are subject to freezing
due to exposure to low temperatures. For instance, aircraft icing [5] is a
very commonly observed case of freezing of water films, which forms layers of
ice on the aircraft surface. This is also observed quite commonly on wind-
turbines in cold weather [6]. Pharmaceutical thin-film freezing [7] is yet
another application, used in the production of pulmonary drugs, where a liquid
thin-film is solidified over a cryogenic surface. Spray-coating of tablets too
involves deposition and freezing of liquid drops over a substrate [8]. Thus,
the freezing of thin-films is an important phenomenon to study and numerical
models assist with better understanding of the physics of the process.

Freezing of thin-films has been an important area of focus for modelling and
simulation in aircraft icing research. Messinger [9] proposed one of the first
icing models. Subsequently, Myers [5, 10, 11, 12, 13] enhanced the model to
accommodate more physics involved in the phenomenon of icing. In the last two
decades, this model has been frequently used as a base for ice and water film
predictions[10, 14, 15, 16]. In this model, the evolution of thickness of the
ice and water film is tracked using a coupled pair of equations representing
mass and energy balance in the water and ice layers. These equations are
solved in an Eulerian framework over the entire surface.  The model assumes
that the surface has high thermal conductivity, such that water incident on
the surface immediately reaches the same temperature as the surface.
Additionally, the temperature distribution in the water and ice layers are
assumed to be quasi-steady and thus, linear.  Gori et al. [17] propose the use
of an analytical solution to the unsteady Stefan problem. In this model, the
temperature of the ice at the surface is assumed to be equal to the
temperature of the surface. Liu et al. [18] solve the unsteady heat equation
numerically using a finite-difference method in the ice layer, while making a
quasi-steady assumption in the water layer. Chauvin et al. [19] proposed a
three-layer icing model with an ice layer sandwiched between a running film
and a static film. They assume the temperature of the running layer to be
uniform, while accounting for unsteady variation of temperature in the ice
layer. Myers et al. [20] too proposed a three-layer model and solved the
unsteady heat equation in each of the layers by assuming a cubic polynomial
solution for the temperature in the individual layers. In the work of Malik et
al. [21], a mixture of ice and water is considered on a substrate and the
unsteady heat equation is solved in the mixture. In this method, the interface
between ice and water is captured using the enthalpy approach for the phase-
change. Impinging droplets on the thin liquid film have been handled using
both Eulerian [22, 23, 16] and Lagrangian approaches[24, 25, 26]. The two
approaches have their respective pros and cons. Since the Eulerian approach
models the droplets as a continuous field, the resolution of individual
droplets is not necessary and therefore, can be computationally efficient.
However, in the Lagrangian approach, the droplet physics can be captured more
accurately and phenomena like splash and rebound can be easily incorporated in
the model. Some of the prominent icing codes that handle freezing of water
films are LEWICE[25, 27, 28], ONERA[29, 30], GlennICE[31, 32, 33], FENSAP-
ICE[16, 34]. Makonnen [35] proposed an alternate model that predicts the rate
of ice accretion based on collision, sticking and accretion efficiencies.
These efficiencies are derived from empirical data fits. The model has been
used for icing predictions on wind turbines and wind farms [6, 36, 37]. In
recent times, machine learning-based models have gained popularity in
predictions of thin-film behaviour [38] and have been also been used in icing
predictions [39, 40]. High fidelity approaches have also been used to model
and solve freezing of impinging droplets over surfaces [41, 42, 43]. In these
appraoches, the complete Navier–Stokes equations are solved in the interior of
each droplet. Such approaches capture the physics very accurately but become
impractical when the number of droplets is large due to high computational
costs.

The model introduced in the present work, called Discrete Droplet Method with
Solidification (DDM-S), relaxes some of the assumptions  common to many of the
Messinger-family of models stated above.  While many of these models, with a
few exceptions [19, 21, 20], solve a quasi-steady thermal problem as a part of
the solidification process, the DDM-S model allows unsteady temperature
distributions in the solid and liquid films. This is accounted for by
evaluating the solid-liquid interface movement based on the solution of the
unsteady energy equation in the solid and liquid regions. The present model
also allows for the fluid in contact with the surface to have a different
temperature than that of the surface. In other words, it does not assume a
perfect thermal contact, as required by the extended Messinger’s model[5].
This reflects the physics more accurately since the liquid need not
necessarily attain the same temperature as the surface, especially on
impingement. The temperature of the liquid changes with time, owing to the
heat transfer occurring between the liquid and the surface. The DDM-S model
also conveniently incorporates volumetric change based on the density
difference between the solid and liquid phase. The change in height of the
solid layer and that of the liquid thin-film, is scaled using the density
ratios in order to account for the volumetric change. The DDM-S model does not
rely on empirical relations and attempts to capture the behaviour from
physical principles. The basis of this model is the discrete droplet method
(DDM) [44] that governs the evolution of liquid thin-films. The DDM, owing to
its ability to capture liquid thin-film behaviour accurately, forms a good
basis for the solidification model.  In the DDM, a thin-film is modelled as an
aggregation of droplets on a surface and the dynamics of the thin-film are
predicted based on the motion of these droplets. The model can capture the
dynamics of fully developed thin-films as well as developing thin-films.
Several of the advantages of using the DDM [44] are passed on to DDM-S. For
instance, the DDM efficiently captures partially wetted surfaces since the
descritization based on droplets is used only in regions where the thin-film
is present and not over the entire surface, thus, reducing computational
costs.

The layout of the paper is as follows. Section Methodology describes in detail
the different sub-parts of the proposed model, DDM-S. Some of the details that
are not directly relevant to the model are described in appendixes. Section
Results presents five test cases, which are arranged in an order of increasing
complexity.  The first three test cases validate specific parts of the model
by comparing the numerical solutions obtained from DDM-S with analytical
solutions. The last two test cases compare the numerical solution of the
entire DDM-S model with numerical and experimental data from literature.
Finally, we present the conclusions.

##  2 Methodology

In this section, we discuss in detail the proposed thin film solidification
model. The section begins with a brief overview of DDM which is the existing
liquid thin-film model that forms a base of the present work. We then proceed
to introduce the new solidification model (DDM-S) in the following
subsections.

###  2.1 Terminology

Typical droplet-based configurations for the DDM model and the DDM-S model are
shown in Figs. 1 and 2 respectively. In this work, the words – substrate and
surface, are used interchangeably and denote the underlying surface in
$\mathbb{R}^{3}$ over which liquid thin-films form. The liquid-thin film may
be referred to as liquid-film or thin-film. The solidified region is sometimes
also referred to as solid region or solid layer. The thickness of the film is
referred as its height.

###  2.2 Discrete Droplet Method (DDM)

(a)

(b)

Figure 1:  Discrete Droplet Model (DDM): (a) A surface in $\mathbb{R}^{3}$
with an aggregation of droplets (b) Liquid film height determined from droplet
locations on a surface.

In the discrete droplet method (DDM) [44], a thin liquid film forming over a
surface in $\mathbb{R}^{3}$ is discretized using droplets of prescribed
diameters as shown in Fig. 1(a). The surface is given by a triangulated mesh
with a set of consistently oriented surface normals. The DDM is a pseudo-2D
model, in which the droplets move along the surface and the film height
accounts for the third dimension. The aggregation of droplets represents the
liquid thin-film as shown in Fig. 1(b) and is used to estimate the film height
at a given location on the surface, as shown in the equation below.

$H_{l,i}=\sum_{j\in S_{i}}W_{ij}V_{j}$ (1)

Here, $H_{l,i}$ is the height of the liquid film at droplet $i$, $S_{i}$ is
the neighborhood of droplet $i$ determined by proximity, $W_{ij}$ is the
weight associated with droplet $j$ depending on its position with respect to
droplet $i$, and $V_{j}$ is the volume of the droplet $j$. As shown in Fig. 1,
in regions where there is more clustering of the droplets on the surface, the
reconstructed height function, $H_{l}$, is higher. The derivatives of height
and velocity are estimated using an SPH-like approach that uses Gaussian
kernel functions. The neighbourhoods $S_{i}$ required for the derivative
computations are obtained from search algorithms commonly used in meshfree
methods[45, 46].  The model solves the momentum equations at the droplet
locations and accounts for the shear forces at the interface of the thin-film
and the surface and the pressure gradients as a function of the gradients of
the film-height. For further details of the model, the work of Bharadwaj et al
[44] may be referred.

###  2.3 An overview of the new solidification model (DDM-S)

In the proposed solidification model, DDM-S, we consider a liquid thin-film
and a solidified region over a surface in $\mathbb{R}^{3}$. The basic
components of the model are presented below.

  * a.

The discretization of the liquid thin-film is still done using droplets, just
as in the case of the DDM. (See Sec. 2.2)

  * b.

The liquid height $H_{l}$ is computed as in the case of the DDM. (See Sec.
2.2)

  * c.

In addition, we compute a solid height, $H_{s}$, as an attribute/property of a
droplet. Consequently, a solidified region cannot be represented unless liquid
droplets are present. Droplets are not located directly on the surface.
Instead, they are at a distance $H_{s}=H_{s}(\vec{x},t)$ away from the
surface, as shown in Fig. 2. (See Sec. 2.4)

  * d.

For the evolution of the solid height $H_{s}$, an interface speed is evaluated
in a 1-D subdomain, spanning between the surface and the edge of the liquid
thin-film, by solving the Stefan problem. This interface speed governs the
actual phase change process. (See Sec. 2.5.1)

  * e.

The droplet volume is then updated based on the change in liquid height. This
step accounts for the volumetric change during solidification. (See Sec. 2.6)

  * f.

Unlike in the DDM, the liquid droplets move on the solidified region, along
the solid-liquid interface, as shown in Fig. 2. In DDM, droplets move along
the surface, with the velocity tangential to the surface. Here, in DDM-S,
droplets moving along the solid-liquid interface, with the velocity tangential
to this interface. Thus, the normals to the solid-liquid interface are
computed at each droplet location and at each time step. (See Sec. 2.7)

  * g.

As droplets move over the solid region, the solid height property, $H_{s}$, is
updated using its derivatives in an advection equation. (See Sec. 2.8)

Each of these aspects of the model are expanded in detail in the subsections
that follow.

###  2.4 Solid height and solid-liquid interface

The solid height at a droplet location, denoted as $H_{s}$, is stored as a
property of the droplet. The droplets are initialized with a solid height
field at the beginning of every simulation. The solid-liquid interface, thus,
is located at a distance $H_{s}$ away from the surface, along the surface
normal, as shown Fig. 2. Unlike the liquid film height ($H_{l}$), the solid
height ($H_{s}$) is not obtained by summing the Gaussian functions of the
droplets in the neighbourhood. Instead, it is treated as a field property
which is initialized and evolved at the droplet locations. The evolution of
$H_{s}$ occurs due to

  1. a.

freezing/melting, as governed by the Stefan problem (Secs. 2.5.1, 2.6).

  2. b.

movement of droplets over the stationary solidified region that is captured by
the advection equation that accounts for the velocity of the droplet and the
gradient of solid height at the droplet location (Sec. 2.8).

Note that since the solid region is stored only as a property at the droplet
locations, without droplets, a solid region cannot be represented by this
model. However, this does not imply that a liquid thin-film must always be
present to model a solidified region. For instance, an aggregation of droplets
with zero volume each would represent a liquid thin-film of zero height, and
can store a non-zero solid height ($H_{s}$).

Figure 2: Modelling of solid height and solid-liquid interface: $H_{l}$
indicates the liquid height and $H_{s}$ indicates the solid height.

###  2.5 Determination of interface speed: Stefan problem for phase-change

In applications relating to freezing of liquid thin-films, the temperature
gradients are typically much higher in the surface normal direction in
comparison to the tangential directions, thus, resulting in solidification
predominantly in the surface normal direction. We, thus, adopt a pseudo-1-D
approach that accounts for the surface normal gradients of temperature, to
solve the phase-change problem at each droplet location. This is done by
assuming a temperature profile at the droplet location, in a manner similar to
assuming a velocity profile in a liquid thin-film required for viscous force
calculations.
At a given droplet location, we consider a 1-D domain that lies between the
surface and the edge of the liquid thin-film, as shown in Fig 3. The solid-
liquid interface is located at a distance $H_{s}$ from the surface. The energy
equation is solved on either sides of this interface using a numerical
discretization of $N_{l}$ points on the liquid side and $N_{s}$ points on the
solid side, see Fig 3.

Figure 3: Evaluation of interface speed using a 1-D subdomain: The red points
are ghost points used in the evaluation of energy equation on the liquid side
($\Omega_{l}$) and the solid side ($\Omega_{s}$)

####  2.5.1 Stefan problem

The Stefan problem was first solved for a phase-change problem by J Stefan
[47]. We first derive the expression for the speed of the solid-liquid
interface from Rankine–Hugoniot (R-H) jump conditions for conservation laws
[48]. Consider a conservation law of the form

$\frac{\partial Q}{\partial t}+\nabla\cdot\vec{F}=0$ (2)

Here, $Q$ is the conserved quantity, and $\vec{F}$ is the flux. This equation
is applicable in the two phases that are separated by a phase-change
interface. The speed of this interface, $v_{I}$, is given by

$v_{I}=\frac{[F]\cdot\hat{n}}{[Q]}$ (3)

In the above equation, the pair of square brackets, $[\text{ }]$, denotes the
jump in the enclosed quantity across the solid-liquid interface and $\hat{n}$
is the unit normal vector of the interface.

In the present model, we use the R-H conditions to derive the solid-liquid
interface speed assuming the flow is incompressible in each phase.
Furthermore, the energy equation with only conduction is considered in each
phase. The R-H condition is applied to both the mass and energy conservation
equations. The continuity equation can be written as

$\frac{\partial\rho}{\partial t}+\nabla\cdot(\rho\vec{V})=0$ (4)

which is solved in both solid and liquid phases. Here, $\rho$ is the density
of the phase and $\vec{V}$ is the velocity of the phase. On applying the R-H
condition, Eq. 3, across the interface, the interface speed is estimated as

$v_{I}=\frac{[\rho\vec{V}]\cdot\hat{n}}{[\rho]}=\frac{(\rho_{L}\vec{V}_{L}-\rho_{S}\vec{V}_{S})}{\rho_{L}-\rho_{S}}\cdot\hat{n}$
(5)

Here the subscripts $L$ and $S$ denote the liquid phase and the solid phase
respectively. In the present model, the velocity of the solid phase is zero
since it is assumed to be stationary. Therefore,

$v_{I}=\frac{\rho_{L}\vec{V}_{L}}{\rho_{L}-\rho_{S}}\cdot\hat{n}\implies\vec{V_{L}}\cdot\hat{n}=v_{I}\frac{\rho_{L}-\rho_{S}}{\rho_{L}}$
(6)

The R-H condition is now considered for the energy equation, which can be
written as

$\frac{\partial\rho E_{t}}{\partial t}+\nabla\cdot(\rho
E_{t}\vec{V})=-\nabla\cdot\vec{q}$ (7)

which is also solved in both solid and liquid phases. Here, $E_{t}$ is the
internal energy and $\vec{q}$ is the heat transfer due to conduction in the
considered phase. Once again, applying the R-H condition, Eq. 3, at the solid-
liquid interface we arrive at

$v_{I}=\frac{[\rho E_{t}\vec{V}+\vec{q}]\cdot\hat{n}}{[\rho
E_{t}]}=\frac{(\rho_{L}E_{t,L}\vec{V}_{L}+\vec{q}_{L}-\rho_{S}E_{t,S}\vec{V}_{S}-\vec{q}_{S})}{\rho_{L}E_{t,L}-\rho_{S}E_{t,S}}\cdot\hat{n}$
(8)

Setting the solid phase velocity to zero and regrouping,

$v_{I}(\rho_{L}E_{t,L}-\rho_{S}E_{t,S})=\rho_{L}E_{t,L}(\vec{V}_{L}\cdot\hat{n})+(\vec{q}_{L}-\vec{q}_{S})\cdot\hat{n}$
(9)

Substituting for $\vec{V_{L}}\cdot\hat{n}$ from Eq. 6 and regrouping,

$v_{I}\left(\rho_{L}E_{t,L}-\rho_{S}E_{t,S}-\rho_{L}E_{t,L}\frac{\rho_{L}-\rho_{S}}{\rho_{L}}\right)=(\vec{q}_{L}-\vec{q}_{S})\cdot\hat{n}$
(10) $v_{I}\rho_{S}(E_{t,L}-E_{t,S})=(\vec{q}_{L}-\vec{q}_{S})\cdot\hat{n}$
(11)

The interface is assumed to be at the melting temperature. Therefore, the
difference in the energy $(E_{t,L}-E_{t,S})$ at the interface, is only due to
latent heat, $L$. Substituting for the heat fluxes using the Fourier’s law,

$v_{I}=\frac{(\vec{q}_{L}-\vec{q}_{S})}{\rho_{S}L}\cdot\hat{n}=\frac{1}{\rho_{S}L}\left(-\kappa_{L}\frac{\partial
T}{\partial n}\Big{|}_{L}+\kappa_{S}\frac{\partial T}{\partial
n}\Big{|}_{S}\right)$ (12)

Here $\kappa$ denotes the thermal conductivity. Using Eq. 12, the interface
speed is computed from the jump in the temperature gradients across the
interface. To evaluate the temperature gradients, we solve the energy
equation, Eq. 7, in the two phases, as a 1-D heat equation at each droplet
location. This is the case since the underlying flow is assumed to be
incompressible. The solution to the 1-D heat equation is explained in Sec.
2.5.3.

####  2.5.2 Temperature profile at droplet locations

For each droplet, a temperature profile $T(n)$ needs to be reconstructed, as a
function of the normal coordinate $n$ which denotes the distance from the
surface in the normal direction. The temperature profile serves as an initial
condition for the heat equation described in the previous subsection. At the
solid-liquid interface, the temperature is equal to the melting temperature.
At the solid – surface interface and ambient-medium – liquid-film interface,
the temperature satisfies a convective heat transfer condition. Additionally,
the height averaged temperature distribution in the liquid thin-film should
equal the droplet temperature. These conditions are expressed mathematically
as

$T(n=H_{s})=T_{m}$ (13) $\frac{\partial T}{\partial
n}\Big{|}_{n=0}=\frac{h_{s}}{\kappa_{s}}(T(0)-T_{s})$ (14) $\frac{\partial
T}{\partial
n}\Big{|}_{n=H_{s}+H_{l}}=-\frac{h_{a}}{\kappa_{l}}(T(H_{s}+H_{l})-T_{a})$
(15) $\frac{1}{H_{l}}\int_{H_{s}}^{H_{l}}T(n)dn=T_{avg}$ (16)

In the above equations, $T_{m}$ is the melting/freezing temperature, $T_{s}$
is the surface temperature, $T_{avg}$ is the average droplet temperature and
$h_{s}$ and $h_{a}$ are the heat transfer coefficients at the surface-solid
interface and ambient medium-liquid interface respectively.  It is noted that
spatial and temporal variations in $T_{s}$ can be accommodated through Eq. 14.
In the liquid region, there are three conditions that are to be satisfied by
the temperature profile – temperature at the solid-liquid interface (Eq. 13),
heat transfer between the liquid-film and the ambient medium (Eq. 15) and the
average temperature of the droplet (Eq. 16). Thus, the temperature profile in
the liquid region is assumed to follow a quadratic curve of the form

$T(n)=a_{l}n^{2}+b_{l}n+c_{l}\text{ , }H_{s}<n\leq H_{l}$ (17)

The constants $a_{l}$, $b_{l}$ and $c_{l}$ can be evaluated from conditions of
Eqns. 13,15 and 16.
The temperature profile in the solid region is assumed to follow a linear
curve, owing to the presence of two conditions – the temperature at the solid-
liquid interface (Eq. 13) and the heat transfer between the solid region and
the underlying surface (Eq. 14).

$T(n)=a_{s}n+b_{s}\text{ , }0\leq n\leq H_{s}$ (18)

The constants $a_{s}$ and $b_{s}$ can be evaluated from conditions of Eqns. 13
and 14. The temperature profile is, thus, continuous at the solid-liquid
interface with its value equal to the melting temperature, but the normal
derivative of temperature could have a discontinuity. It is this discontinuity
in the derivative that causes the movement of the solid-liquid interface, as
seen in Eq. 12. The derivation of the constants involved in the temperature
profile $(a_{l},b_{l},c_{l})$ and $(a_{s},b_{s})$ are elaborated in A.

It is reiterated that the temperature profile obtained as described here,
serves as an initial condition. As the simulation progresses, the change in
the temperature profile is given by the heat equation described in Sec. 2.5.3.

####  2.5.3 Numerical solution of the energy equation in a 1-D sub-domain with
an evolving interface

At each droplet location, the temperature gradients on the liquid side and the
solid side are needed, as shown in Eq. 12. The heat equation is, thus, solved
at each droplet location in a 1-D subdomain discretized by $N$ points, see Fig
3.

$\frac{\partial T}{\partial t}=\frac{\kappa}{\rho
c}\frac{\partial^{2}T}{\partial n^{2}}=\alpha\frac{\partial^{2}T}{\partial
n^{2}}$ (19)

Here, $c$ denotes specific heat and $\alpha$ denotes the thermal diffusivity.
The initial condition at each point is provided by the analytical form
described in Sec.2.5.2. The physical extent of the 1-D sub-domain begins at
the surface and ends at the edge of the liquid film. There are two ghost
points on either side, as shown by the red points in Fig 3, to facilitate the
imposition of boundary conditions. The solid-liquid interface is at the
melting temperature and is used as a boundary condition for the heat equation
on both the solid and liquid phases. The temperatures at the two ghost points
($g1$ and $g2$) are determined based on the convective heat transfer boundary
condition at the surface and edge of the liquid film, as shown below:

$Q_{s}=\kappa\frac{\partial T}{\partial n}\Big{|}_{s}=h_{s}(T_{g1}-T_{s})$
(20) $Q_{a}=\kappa\frac{\partial T}{\partial
n}\Big{|}_{a}=-h_{a}(T_{g2}-T_{a})$ (21)

Here, $h_{s}$ and $h_{a}$ are the convective heat transfer coefficients at the
surface and at the interface with the ambient fluid respectively. $T_{s}$ and
$T_{a}$ are the surface temperature and the temperature of the ambient fluid
respectively. $Q_{s}$ and $Q_{a}$ denote the heat transfer rate due to the
convective boundary condition at the respective boundary.
Eq. 19 is solved using a finite difference method.

$T_{k}^{m+1}(1+2\lambda)-T_{k+1}^{m+1}\lambda-T_{k-1}^{m+1}\lambda=T_{k}^{m}$
(22)

where $\lambda=\frac{\alpha\Delta t}{\Delta n^{2}}$, and $T_{k}^{m}$ denotes
the temperature $T$ at point $k$ of the 1-D sub-domain at time level $m$..
Near the boundaries and in the vicinity of the solid-liquid interface, the
above equation is modified to impose the boundary conditions. For instance, at
the first and last points of the 1-D sub-domain, the discretization is shown
below.

$T_{1}^{m+1}(1+2\lambda)-T_{2}^{m+1}\lambda=T_{1}^{m}+T_{g1}\lambda$ (23)
$T_{N}^{m+1}(1+2\lambda)-T_{N-1}^{m+1}\lambda=T_{N}^{m}+T_{g2}\lambda$ (24)

For the point that lies on the solid-liquid interface, the melting temperature
is forced.

$T_{N_{s}}=T(n=H_{s})=T_{m}$ (25)

These equations lead to a system of linear equations with a tridiagonal matrix
of the form

$\setcounter{MaxMatrixCols}{11}\begin{bmatrix}1+2\lambda&-\lambda&0&.&.&.&.&.&.&.&.\\\
-\lambda&1+2\lambda&-\lambda&0&.&.&.&.&.&.&.\\\ .&&&&&&&&&&\\\ .&&&&&&&&&&\\\
.&&&&&&&&&&\\\ .&&.&.&0&1&0&.&.&&\\\ .&&&&&&&&&&\\\ .&&&&&&&&&&\\\
.&&&&&&&&&&\\\ .&.&.&.&.&.&.&.&-\lambda&1+2\lambda&-\lambda\\\
.&.&.&.&.&.&.&.&0&-\lambda&1+2\lambda\\\
\end{bmatrix}\begin{bmatrix}T_{1}^{m+1}\\\ T_{2}^{m+1}\\\ .\\\ .\\\ .\\\
T_{N_{s}}\\\ .\\\ .\\\ .\\\ T_{N-1}^{m+1}\\\ T_{N}^{m+1}\\\
\end{bmatrix}=\begin{bmatrix}T_{1}^{m}+T_{g1}\lambda\\\ T_{2}^{m}\\\ .\\\ .\\\
.\\\ T_{m}\\\ .\\\ .\\\ .\\\ T_{N-1}^{m}\\\ T_{N}^{m}+T_{g2}\lambda\\\
\end{bmatrix}$ (26)

The above system of equations is solved using the LSQR algorithm. The number
of points in the 1-D sub-domain may be chosen such that it sufficiently
resolves the temperature gradients.  The use of too few points may adversely
affect the solution. In the test cases presented in this paper, the total
number of points in the 1-D subdomain is 20 with 10 points discretising the
solid region and 10 points discretising the liquid region. For these test
cases, using a larger number of points does not change the solution.

By solving for the temperature field from the heat equation, as described
here, this model accommodates for unsteadiness in the temperature distribution
and relaxes assumptions imposed by Myers’ model [5] about the thermal problem
being quasi-steady and the temperature profile being linear.

####  2.5.4 Temperature interpolation due to change in droplet neighborhood

Figure 4: A schematic representing a scenario of temperature projection when a
droplet moves or neighbourhood changes between time-steps. The primed
quantities are before and non-primed are after the droplet movement/change in
neighbourhood has occurred.

In Eq. 22, $T^{m}$ denotes the temperature of the droplet at its location at
the beginning of time level $m+1$. However, this temperature is not directly
available since the previously stored temperature profile corresponds to the
droplet’s previous location and neighborhood, before the Lagrangian motion.
Therefore, the temperature profile at the end of time level $m$ needs to be
projected to the current droplet locations, to accommodate its new location
and new neighborhood before solving Eq. 22. We denote the known temperature
profile before the droplet movement as $T^{k^{\prime}}$ and the projected
temperature profile as $T^{k}$ (corresponds to $T^{m}$ of Eq. 22).
This process is illustrated in Fig. 4. We first consider the solid region. Let
the normal coordinate before projection be denoted as $n^{\prime}$ and after
the projection be denoted as $n$ . The projected temperature profile of the
1-D sub-domain is given as

$T^{k}(n)=T^{k^{\prime}}(n^{\prime})+\frac{\partial
T}{\partial\tau}(n^{\prime})\Big{|}_{k^{\prime}}\Delta\tau$ (27)

where $\frac{\partial T}{\partial\tau}$ is the variation of temperature in the
tangential direction and $\Delta\tau$ is the distance moved by the droplet on
the tangent plane. We assume a linear variation of $\frac{\partial
T}{\partial\tau}$ w.r.t $n^{\prime}$. The gradient at the wall $\frac{\partial
T}{\partial\tau}\Big{|}_{n^{\prime}=0}=q_{1}$ and at the solid-liquid
interface $\frac{\partial
T}{\partial\tau}\Big{|}_{n^{\prime}=H^{\prime}_{s}}=0$, serve as boundary
conditions for the linear profile. This gives

$\frac{\partial
T}{\partial\tau}(n^{\prime})=q_{1}\Big{(}1-\frac{n^{\prime}}{H^{\prime}_{s}}\Big{)}$
(28)

Substituting this in Eq. 27,

$T^{k}(n)=T^{k^{\prime}}(n^{\prime})+q_{1}\Big{(}1-\frac{n^{\prime}}{H^{\prime}_{s}}\Big{)}\Delta\tau$
(29)

Secondly, we consider the liquid region. We proceed along similar lines as in
the solid region,

$T^{k}(n)=T^{k^{\prime}}(n^{\prime})+\frac{\partial
T}{\partial\tau}(n^{\prime})\Big{|}_{k^{\prime}}\Delta\tau$ (30)

In this case, however, the normal coordinate is measured starting from the
solid-liquid interface. The boundary conditions for the linear profile of
$\frac{\partial T}{\partial\tau}$ is $\frac{\partial
T}{\partial\tau}\Big{|}_{n^{\prime}=0}=0$ and $\frac{\partial
T}{\partial\tau}\Big{|}_{n^{\prime}=H^{\prime}_{l}}=\gamma$. This gives,

$\frac{\partial T}{\partial\tau}(n^{\prime})=\frac{\gamma
n^{\prime}}{H^{\prime}_{l}}$ (31)

Substituting this in Eq. 30, we get

$T^{k}(n)=T^{k^{\prime}}(n^{\prime})+\frac{\gamma
n^{\prime}}{H^{\prime}_{l}}\Delta\tau$ (32)

Usually, the tangential derivative, $\gamma$ at $n^{\prime}=H^{\prime}_{l}$ is
unknown. The normal derivative, on the other hand, can be estimated from the
convective heat transfer occurring between the liquid-film and the
surroundings. Assuming this to be equal before and after projection,

$\frac{\partial T}{\partial n}\Big{|}_{n=H_{l}}=\frac{\partial T}{\partial
n}\Big{|}_{n^{\prime}=H^{\prime}_{l}}=q_{2}$ (33)

Differentiating Eq. 32 w.r.t $n$, we get

$\frac{\partial}{\partial n}T^{k}(n)=\frac{\partial}{\partial
n^{\prime}}T^{k}(n)\frac{\partial n^{\prime}}{\partial
n}=\Big{[}\frac{\partial}{\partial
n^{\prime}}T^{k}(n^{\prime})+\frac{\gamma}{H^{\prime}_{l}}\Delta\tau\Big{]}\frac{H^{\prime}_{l}}{H_{l}}$
(34)

In the above equation, $\frac{\partial n^{\prime}}{\partial n}$ denotes the
scaling of the normal coordinate which is given by the ratio of the liquid
heights $\frac{H^{\prime}_{l}}{H_{l}}$. Now, setting $n=H_{l}$ and
substituting from Eq. 33, we get

$q_{2}=\Big{[}q_{2}+\frac{\gamma}{H^{\prime}_{l}}\Delta\tau\Big{]}\frac{H^{\prime}_{l}}{H_{l}}$
(35)
$\gamma=q_{2}\Big{(}1-\frac{H^{\prime}_{l}}{H_{l}}\Big{)}\frac{H^{\prime}_{l}}{\Delta\tau}$
(36)

Substituting for $\gamma$ in Eq. 32, we get the below equation that is used to
project the temperature in the liquid region of the 1-D sub-domain.

$T^{k}(n)=T^{k^{\prime}}(n^{\prime})+q_{2}n^{\prime}\Big{(}1-\frac{H^{\prime}_{l}}{H_{l}}\Big{)}$
(37)

###  2.6 Updating heights, droplet repositioning and updating droplet volume

The movement of the interface, given by the Stefan problem (see Sec. 2.5.1),
corresponds to a change of the heights of the solid and liquid parts of the
thin film, as shown below.

$v_{I}=\frac{\partial H_{s}}{\partial t}=-\beta\frac{\partial H_{l}}{\partial
t}$ (38)

In this equation $\beta$, which is the ratio of density of the solid to that
of the liquid, accounts for volumetric expansion/contraction.
When the solid height changes based on the Stefan condition, the droplet needs
to be re-positioned to the updated solid-liquid interface by moving it by the
appropriate distance.

$\Delta\vec{x}=\Delta H_{s}\vec{n}_{s}$ (39)

Here, $\vec{n}_{s}$ is the unit normal to the surface. The volume of the
droplet also needs to be updated depending on whether $H_{s}$ increases or
decreases. For instance, if $H_{s}$ increases, the height of the liquid-film
$H_{l}$ will reduce owing to mass conservation. This reduction in height is
accounted for by a reduction in droplet volume. To derive an expression for
the updated volume, a solid fraction is defined at the droplet location as

$f=\frac{H_{s}}{H_{s}+H_{l}}$ (40)

The volume of the droplet relates to the solid fraction as

$V=(1-f)V_{0}$ (41)

where $V_{0}$ corresponds to the volume of the droplet when no solid is
present (i.e. $f=0$). Therefore,

$\Delta V=-V_{0}\Delta f$ (42)

The change in solid fraction can be estimated by using the old and new heights
in Eq. 40.

###  2.7 Updating solid normal

For the movement of the liquid thin-film over the solidified region, the
surface normals of the conventional DDM need to be revisited. The normal to
the solidified region, associated with a given droplet is referred to as the
solid normal. The distinction between the solid and surface normal is
illustrated in Fig.5.  While the surface normal is used for the solution of
the heat equation in the 1-D subdomain, the solid normal is used for the
movement of droplets over the solidified region.

Figure 5:  Solid normal vs Surface normal: While DDM requires the surface
normal for droplet movement, the DDM-S requires the solid normal.

The solid normal is calculated based on the position of droplets in the
neighbourhood of the given droplet. This procedure is elaborated in [49]. For
a given droplet $i$, the quantity $J$ shown below, needs to be minimized. $J$
denotes sum of the projections of the solid normal onto the solidified region.
The solidified region is represented by the droplets in the neighbourhood of
droplet $i$. When this sum, $J$, is minimized, the resulting vector
$\vec{n}_{i}$ is the required solid normal.

$J=\sum_{j\in S_{i}}W_{ij}(\vec{n}_{i}\cdot\vec{r}_{ij})^{2}$ (43)

Here, $\vec{n}_{i}$ is the unit solid normal at the droplet $i$ and
$\vec{r}_{ij}$ is the vector joining droplets $i$ and $j$. This can be recast
as

$J=\vec{n}^{T}C\vec{n}$ (44)

where the matrix $C$ is a $3\times 3$ matrix is given by

$\begin{bmatrix}\sum_{j\in S_{i}}W_{ij}\Delta x^{2}&\sum_{j\in
S_{i}}W_{ij}\Delta x\Delta y&\sum_{j\in S_{i}}W_{ij}\Delta x\Delta z\\\
\sum_{j\in S_{i}}W_{ij}\Delta xW_{ij}\Delta y&\sum_{j\in S_{i}}W_{ij}\Delta
y^{2}&\sum_{j\in S_{i}}W_{ij}\Delta y\Delta z\\\ \sum_{j\in S_{i}}W_{ij}\Delta
xW_{ij}\Delta z&\sum_{j\in S_{i}}W_{ij}\Delta z\Delta y&\sum_{j\in
S_{i}}W_{ij}\Delta z^{2}\\\ \end{bmatrix}$ (45)

such that $\vec{r}_{ij}=\Delta x\text{ }\hat{i}+\Delta y\text{ }\hat{j}+\Delta
z\text{ }\hat{k}$. A Gaussian function is used to evaluate the weight,
$W_{ij}$, between droplets $i$ and $j$. This minimization leads to an Eigen
value problem.

${C}\vec{n}=\lambda\vec{n}$ (46)

The normalized Eigen vector corresponding to the smallest Eigen value of the
matrix $C$, is the required solid normal.

###  2.8 Advection of solid height

The solid height, as mentioned earlier, is treated as a property of the liquid
droplet. Since the solidified region is stationary and the droplets move, the
solid height cannot be advected along with the droplet with no change.
Instead, the solid height needs to be revised as

$\frac{\partial H_{s}}{\partial t}-\vec{V}_{\text{drop}}\cdot\nabla H_{s}=0$
(47)

Note the negative sign preceding the convective term. In other words, the
solid region moves at a velocity $-\vec{V}_{\text{drop}}$ w.r.t the droplets.
The derivative of the solid height is given by an SPH-based formulation as

$\nabla H_{s}=-\sum_{j\in S_{i}}\nabla W_{ij}H_{j}\Delta A_{j}$ (48)

From the consistency condition $\sum_{j\in S_{i}}\nabla W_{ij}\Delta A_{j}=0$,
the above equation may be rewritten as

$\nabla H_{s}=-\sum_{j\in S_{i}}\nabla W_{ij}(H_{j}+H_{i})\Delta
A_{j}=-2\sum_{j\in S_{i}}\nabla W_{ij}\frac{(H_{j}+H_{i})}{2}\Delta A$ (49)

$(H_{i}+H_{j})/2$ can be interpreted as an average, $H_{ij}$, at the midpoint
of the line joining droplets $i$ and $j$.

$\nabla H_{s}=-2\sum_{j\in S_{i}}\nabla W_{ij}H_{ij}\Delta A$ (50)

Upwinding can now be performed on the term $H_{ij}$, to render the marching of
the equation stable.

$\displaystyle H_{ij}=H_{i}\text{, if
}(\vec{x}_{i}-\vec{x}_{j})\cdot\vec{V}_{\text{drop}}>0$ (51)
$\displaystyle=H_{j}\text{, if
}(\vec{x}_{i}-\vec{x}_{j})\cdot\vec{V}_{\text{drop}}<0$ (52)

The evolution of solid height using the method described above makes it
susceptible to dispersive and dissipative numerical errors. While smooth solid
profiles can be satisfactorily captured, for profiles that vary sharply these
errors need to be tackled by using sufficiently small temporal and spatial
resolutions or higher-order discretizations.

###  2.9 Surface point and position correction

For each droplet moving on the solid, a corresponding point on the surface is
stored. This is the surface point that lies closest to the droplet such that
the line joining the point and the droplet is normal to the surface. When a
droplet moves over the solidified region, the solid height at the new
location, $H_{s}$, may not necessarily match with the height obtained by the
advection step, $H_{s}^{\text{adv}}$, described in Sec. 2.8.  To ensure that
the actual height , $H_{s}$, matches the height from the advection equation,
$H_{s}^{\text{adv}}$, a correction is introduced in the position as shown
below, where $\hat{n}_{s}$ is the normal to the surface.

$\Delta\vec{x}_{c}=\hat{n}_{s}(H_{s}^{\text{adv}}-H_{s})$ (53)

###  2.10 Impinging droplets

The present model incorporates handling of impinging droplets, which is an
important aspect in applications like aircraft or wind-turbine icing. When a
liquid droplet moving in the bulk-flow meets the surface, it becomes a part of
the thin-film model (DDM). At the location of this droplet, the solidified
height is initialized by performing an interpolation from neighbouring
droplets and an initial temperature profile is prescribed.
The solid height at the new droplet location is interpolated as

$H_{s,i}=\sum_{j\in S_{i},j\neq i}W_{ij}H_{s,j}$ (54)

For the temperature profile at the new droplet location, the following
boundary conditions are used to determine a quadratic temperature profile.

$T(n=H_{s})=T_{m}$ (55) $\frac{\partial T}{\partial
n}\Big{|}_{n=H_{s}+H_{l}}=-h_{a}[T(n=H_{s}+H_{l})-T_{a}]$ (56)
$\frac{1}{H_{l}}\int_{H_{s}}^{H_{l}}T(n)dn=T_{avg}$ (57)

Here, $T_{m}$ is the melting temperature, $T_{a}$ is the ambient temperature
and $h_{a}$ is the convective heat transfer coefficient at the edge of the
liquid film.

###  2.11 Summary of algorithm

The entire DDM-S model is implemented in an in-house solver,
MESHFREE‡‡‡www.meshfree.eu. A summary of steps of the algorithm is given
below, with links to the appropriate sections above.

Algorithm 1 Discrete Droplet Method with Solidification

1:Initialize solid height $H_{s}$ at the droplet locations.

2:Move the droplets from the surface to the solid-liquid interface

3:while Time-stepping loop do

4: Lagrangian movement of droplets [50]

5: Recomputing surface points - Sec. 2.9

6: Advection of solid height - Sec. 2.8

7: Position correction - Sec. 2.9

8: Impinging droplet treatment - Sec. 2.10

9: Stefan problem - Evaluation of interface speed - Sec. 2.5.1

10: Update heights, droplet volume and position - Sec. 2.6

11: Update solid normals - Sec. 2.7

12: Compute new  droplet velocity using DDM [44]

13:end while

##  3 Results

In this section, we present the results of the proposed model in stages of
increasing complexity. The test cases are tabulated in Table. 1. While test
cases 1-3 verify individual parts of the model, the last two test cases
involving a circular conductor and an aircraft wing incorporate all aspects of
the model. In all the test cases, the ratio of densities of the solid and
liquid, $\beta$, is set to unity implying no volumetric change. A validation
of the underlying DDM model (without solification), can be found in our
earlier work [44].

Table 1: List of test cases

Test case | Remarks
---|---
1: Thin-film of constant velocity sliding over a solid layer of varying height | Verifies advection of solid height
2: Sliding thin-film of constant height over a solid layer of varying height | Verifies DDM solution in the presence of a solid region
3: Solidification of a stationary thin film | Verifies the Stefan condition-based phase-change model
4: Solidification on a circular conductor with droplet impingement | Validates the complete model
5: Solidification on an aircraft wing with droplet impingement | Validates the complete model

###  3.1 Test Case 1: Thin-film of constant velocity sliding over a solid
layer of varying height

In this test case, we consider a stationary, wavy solidified region over which
the droplets of the thin-film slide at a constant velocity of $10$m/s
magnitude. The purpose of this test case is to verify the movement of droplets
over the solid and the advection of solid-height as the movement occurs. A
flat surface, as shown in Fig. 6, is considered and a solidified region is
assumed, whose height is given by

$H_{s}=0.1+0.05\text{sin}\left(\frac{4\pi x}{L}\right)$ (58)

This equation represents the analytical solid height i.e. the wavy solid
region. $L$ is taken to be 5 units.

(a) A sketch of the solid profile given by Eq. 58, over which droplets slide
at a constant velocity.

(b) The numerical domain of computation: a flat plate initialized with a wavy
solidified region.

Figure 6: Test Case 1: Domain of computation.

Fig 7 shows a droplet of the thin-film being tracked as the film slides. The
other droplets are hidden to enhance visualization. The coloured contour of
the droplet represents the solid height $H_{s}$ stored at its location. In
Fig. 7(a) and (b), the droplet is located close to a crest. The solid height
$H_{s}$ is correspondingly high, as seen by the red colour of the droplet. As
the droplet continues moving, it reaches a trough and takes on lower values of
$H_{s}$, shown by the blue colour in Fig. 7(c) and (d). The trend repeats in
Fig. 7(e) and (f), as the droplet goes on to climb the next crest. Therefore,
it is seen from the frames in the figure that as the thin-film slides, the
liquid droplet moves in accordance to the topology of the solidified region
and also, updates the value of the solid height as per Eq. 47.  Since Eq. 47
is an advection equation, it is likely to introduce numerical errors. The
error is quantified and a convergence with respect to droplet diameter is
reported. We consider three cases with different diameters such that
$d_{1}$=3.375e-3, $d_{2}$ = 2$d_{1}$ and $d_{3}$ = 2$d_{2}$. In these three
cases, the height of the liquid film remains the same, which implies that
smaller diameters would correspond to a higher number of droplets. It is,
therefore, evident that cases with smaller diameter of droplets will result in
smaller numerical errors in comparison to cases with larger droplet diameters.
Fig. 8(a) shows the position of droplets moving over the surface at $t=0.25$s,
for the case with finest diameter of droplets ($d_{1}$). It is noted that some
droplets are located below the analytical surface, which is attributed to the
dissipation errors due to the numerical solution of the advection of solid-
height. Fig. 8(b) shows the convergence of error with droplet diameter. An
error of 2.5%, as defined in Eq. 59, was obtained for the case with finest
diameter.

$\text{error}=\frac{1}{A_{s}}\int|H^{adv}_{s}-H_{s}|dA$ (59)

Here, $H^{adv}_{s}$ denotes the solid height obtained from the advection
equation, $H_{s}$ denotes the actual solid height in Eq. 58 and $A_{s}$
denotes the area of the surface.

(a) t = 0s

(b) t = 0.03s

(c) t = 0.12s

(d) t = 0.15s

(e) t = 0.24s

(f) t = 0.27s

Figure 7: Test Case 1: A droplet is tracked as it slides in the liquid thin-
film over the wavy solidified region. The colour of the droplet represents its
stored solid height. Through frames (a)-(f), the droplet is seen following the
topology of the solidified region. As it moves, the solid height stored is
higher in the crests (red colour) and lower in the troughs (blue colour). This
is achieved by treating solid height as an advected property as per Eq. 47.
Only one droplet is shown for easier visualization.

(a)

(b)

Figure 8:  Test Case 1: Error quantification: (a) Droplets positions at
$t=0.25$s for $d_{1}$=3.375e-3m (b) Error convergence w.r.t droplet diameter;
an error of 2.5% is obtained for the case with finest diameter ($d_{1}$).

###  3.2 Test Case 2: Sliding thin-film of constant height over a solid layer
of varying height

In this test case, the DDM model that solves the liquid thin-film flow is
verified in the presence of a solidified region of linearly varying height. As
seen in Fig. 9, a linearly varying solid height is considered.

$H_{s}(x)=0.05+x\tan\left(\frac{\pi}{36}\right)$ (60)

This solid profile effectively creates a ramp of 5 degrees on which the liquid
thin-film slides. The below analytical expression for velocity of a thin-film
of constant height, can be used for comparing the velocity obtained
numerically.

$V_{\text{drop}}=g\sin\theta\frac{\rho
H_{l}^{2}}{\eta}\left(1-\exp\left[-\frac{\eta t}{\rho
H_{l}^{2}}\right]\right)$ (61)

where $V_{\text{drop}}$ is the magnitude of the velocity of the droplet in the
direction along the solid ramp,  $g=\|\vec{g}\|$ (=$10^{3}m/s^{2}$) is the
acceleration due to gravity, $\theta$ is the angle of inclination of the plane
set to $5^{o}$, $\rho$ (=$1000$ kg/m3) is the density of the liquid film,
$\eta$ (=$0.001$ Ns/m2) is the dynamic viscosity of the liquid film, and
$H_{l}$ is the height of the liquid film. The droplets on the surface are
distributed uniformly over the surface to achieve a film height $H_{l}\approx
3.6e-4$m.

(a) A sketch of the solid profile given by Eq. 60, over which droplets slide
as governed by the DDM.

(b) The numerical domain of computation: a flat plate initialized with an
inclined solidified region.

Figure 9: Test Case 2: Domain of computation.

(a)

(b)

Figure 10:  Test Case 2: Comparison of velocity from the solution of the DDM
model with Solidification (DDM-S) with the analytical solution; (a) Comparison
of simulations with three different diameters,$d_{1}=6.75e-3$,$d_{2}=2d_{1}$,
$d_{3}=2d_{2}$. (b) Error convergence plot.

As per Eq. 61, the liquid film starts from rest and accelerates under the
influence of gravitational forces until it reaches a terminal velocity due to
the balancing viscous forces.
We present a convergence study by considering three cases of varying droplet
diameters such that the height of the liquid film remains the same. The case
with the smallest diameter corresponds to $d_{1}=6.75e-3$m and the other two
cases correspond to $d_{2}=2d_{1}$ and $d_{3}=2d_{2}$. In order to ensure that
these three cases have the same liquid film-height, the number of droplets are
scaled accordingly i.e. the case with small diameter of droplets would have a
higher number of droplets to maintain the same height of liquid film.  Fig.
10(a) shows the evolution of velocity at a fixed point on the ramp. It
compares the velocity obtained from the three simulations of DDM with
solidification (DDM-S) with that of the analytical solution.  The DDM-S
solutions converge as the diameter of the droplets is refined. This is also
quantified by the plot of error in Fig. 10(b), in which the error is defined
as

$\text{error}=\frac{1}{TV_{t}}\int|V_{\text{DDMS}}-V_{\text{ana}}|dt$ (62)

where $V_{\text{DDMS}}$ is velocity obtained from the DDM-S model,
$V_{\text{ana}}$ is the analytical solution, $T$ is the total simulation time
and $V_{t}$ is the analytical terminal velocity of the film. For the case with
the smallest diameter ($d_{1}$), the error is within 0.34% of the analytical
solution.

###  3.3 Test Case 3: Solidification of a stationary thin film

In this test case, a stationary film resting on a flat surface undergoes
directional solidification starting from the surface, as shown in Fig. 11. The
temperature in the liquid film is uniformly at the melting temperature. The
temperature in the solidified region varies from a forcing temperature at the
surface to the melting temperature at the solid-liquid interface. An
analytical solution for the interface location and the temperature profile is
shown below [51].

$H_{s}=v_{I}t$ (63)
$T(n,t)=T_{m}+\frac{L}{C_{v}}\left[1-\exp\left(\frac{v_{I}^{2}}{\alpha}t-\frac{v_{I}}{\alpha}n\right)\right]$
(64)

Here, $L$ is the latent heat (=$10^{3}$J/kg), $C_{v}$ is the specific heat
(=$10^{3}$J/kg), $v_{I}$ is the interface speed, $\alpha$ is the thermal
diffusivity (=0.5m2/s), $T_{m}$ (=0oC) is the melting temperature and $n$ is
the normal coordinate. In the numerical simulation, a Dirichlet boundary
condition is imposed at the surface that can be obtained by substituting $n=0$
in Eq. 64. Imposition of this boundary condition, pushes the solid-liquid
interface away from the surface causing a growth in the solid height
($H_{s}$). Fig. 12(a) compares $H_{s}$ of the DDM-S model with the analytical
solution and Fig. 12(b) compares the temperature profiles of the analytical
and numerical solution at $t=0.07$s.  These results correspond to 20 points
discretizing the 1-D sub-domain. Fig. 13 shows the convergence of error with
respect to the number of points ($N$) used to discretize the 1-D sub-domain.
Beyond $N$=20, the solution does not change significantly.

(a) A sketch of the domain: Stationary droplets are nearly uniformly
distributed over a flat plate and the solidification occurs in the direction
shown.

(b) The numerical domain of computation: a flat plate with droplets
initialized with $H_{s}=0$.

Figure 11: Test Case 3: Domain of computation.

(a) Comparison of solid-liquid interface location

(b) Comparison of temperature profiles

Figure 12: Test Case 3: Comparison of numerical solution of DDM-S and
analytical solution Figure 13:  Test Case 3: Convergence plot. The error drops
with the number of points discretizing the 1-D subdomain (N).

###  3.4 Test Case 4: Solidification on a circular conductor with droplet
impingement

Having verified the different components of the proposed thin film
solidification model in the previous test cases, we now validate the DDM-S
model as a whole. In this test case, a circular, cylindrical conductor is
considered on which water droplets impinge to form ice. At the beginning of
the simulation, there is no water film or ice present on the surface. The
icing conditions are tabulated in Table. 2.

Table 2: Test Case 4: Test conditions $V_{\infty}$ | 10$m/s$
---|---
$T_{a}$, $T_{s}$, $T_{avg}$ | $-5^{o}$
Liquid Water Content (LWC) | $1.8\text{g/m}^{3}$
Mean Volume Diameter | 26 $\mu$m
Icing time | 10 minutes
Cylinder diameter | 34.9 mm

Fig. 14 shows the geometry of the circular cylinder. A fictitious plane
upstream of the geometry is initialized with randomised positions of droplet
sources. The rate of droplet generation depends on the volume flux of water
through this plane. The volume flux of water droplets through the plane of
area $A$, is

$Q_{w}=Q_{a}\frac{\text{LWC}}{\rho_{w}}=AV_{\infty}\frac{\text{LWC}}{\rho_{w}}$
(65)

The subscripts $w$ and $a$ denote water droplets and air respectively, $\rho$
denotes density, LWC stands for Liquid Water Content and $V_{\infty}$ denotes
the freestream air velocity. The diameter of the droplets is related to the
number of droplets introduced per second, $N_{d}$, at the fictitious plane, as

$N_{d}\frac{\pi}{6}d^{3}=Q_{w}$ (66)

The above equations, Eqns. 65 and 66, are used in test cases which involve
impinging droplets with a given LWC, as done in the next test case too.  In
the numerical model, the heat transfer coefficients at the interfaces of the
thin-film with the ambient and the surface, $h_{s}$ and $h_{a}$, are chosen to
be $100$ W/m2K.

Fig. 15 shows the development of the ice profiles. From Fig. 15(a)-(d), an
increase in $H_{s}$ is noted by the darkening red coloured contours of
$H_{s}$. The increase in $H_{s}$ can also be noticed as a widening gap between
the droplets and the cylinder surface. The numerical solution of the model
DDM-S, is obtained by interpolating the solid (ice) heights onto the mid-
section of the cylindrical conductor.  Fig. 16 compares this interpolated ice
shape, with that reported in literature[52], where a glaze ice model was used
to predict the icing. We define an error as in Eq. 67. We note an error of
around 11 % with respect to the numerical data reported in literature.

$\text{error}=\frac{\int_{S}|H_{s,DDM-S}-H_{s,lit}|}{\int_{S}H_{s,lit}}$ (67)

Figure 14: Test Case 4: Geometry of circular conductor.

(a) t/T = 0.25

(b) t/T = 0.5

(c) t/T = 0.75

(d) t/T = 1

Figure 15:  Test Case 4: Formation of ice on a circular conductor. The
coloured contours of ice-height show an increasing trend with simulation time.
Figure 16: Test Case 4: Solidified region of ice on a circular cylinder;
results compared with literature [52].

###  3.5 Test Case 5: Solidification on an aircraft wing with droplet
impingement

In this test case, a GLC-305 wing is considered with droplet impingement. The
reference test case is taken from the experimental data of the Icing Research
Tunnel (IRT) of NASA Glenn Research Centre [53]. The icing condition chosen
for simulation corresponds to Ice-3 of the experimental data. Fig. 17 shows
the wing geometry. The details of the wing are tabulated in Table. 3a.

Figure 17: Test Case 4: GLC-305 wing geometry. Table 3: Test Case 3: Geometry
and icing conditions

Airfoil Section | GLC-305
---|---
Angle of attack | $6^{o}$
Wing twist | $4^{o}$
Root chord |  $25.2$ in
Tip chord |  $10.08$ in
(a) GLC-305 wing details

$V_{\infty}$ | 89.7 m/s
---|---
$T_{\text{a}}$, $T_{\text{s}}$, $T_{\text{avg}}$ | $-11.27^{o}$
Liquid Water Content (LWC) | $0.51\text{g/m}^{3}$
Mean Volume Diameter | 14.5 $\mu$m
Icing time |  $5$ minutes
(b) Icing conditions

The conditions of icing are tabulated in Table. 3b. The numerical values of
the parameters of the solidification model are derived from these icing
conditions, as done in the previous test case, and with the help of inputs
from the NASA Glenn Research Centre.  In the numerical model, as in the
previous test case, the heat transfer coefficients at the interfaces of the
thin-film with the ambient and the surface, $h_{s}$ and $h_{a}$, are chosen to
be $100$ W/m2K.
Fig. 18 shows the formation and evolution of the liquid thin-film (water) over
the lower surface of the wing and its solidification to form ice. The coloured
contours denote the ice height. It is noted that as the simulation time
increases, the water thin-film develops over the wing surface. The ice-height
at the leading-edge is higher than the other parts of the wing and it
increases with time.

(a) t = 2s

(b) t = 15s

(c) t = 25s

(d) t = 30s

Figure 18: Test Case 5: Formation of water film on the lower surface of the
wing. The coloured contours of ice-height show an increasing trend with
simulation time and larger ice accretion at the leading edge.

Fig. 19 shows the ice shape based on interpolated ice heights at the root
section of the wing. The predictions of the proposed DDM-S model matches well
with the experiment near the leading edge of the airfoil. However, further
downstream of the leading edge, the model over-estimates ice-accretion in
comparison to the experiment.  A possible source of this mismatch is that the
aerodynamics plays a prominent role in droplet incidence on the wing, which
has not been considered in the present model. The evolution of the thin-film,
therefore, does not account for the shear and inertial forces exerted by the
bulk-flow on the film, which is a limitation of the model in its present form.
It is also noted that the liquid water content (LWC) used in the experiment
may have been under-reported, as confirmed by NASA Glenn Research Centre.
However, this effect is difficult to quantify.

The role of aerodynamics will be actively captured in future extensions of the
model which will involve the coupling of the present model with a bulk-flow
model. For the purposes of the present work, its effect has been indirectly
considered by using the collection efficiency $\zeta$ which signifies the
amount of droplet incidence occurring at a given section on the surface [54].
We consider a distribution of collection efficiency $\zeta$ over the airfoil
using the freestream velocity vector ($\vec{V}_{\infty}$) and the outward
surface normal ($\hat{n}$).

$\zeta=-\frac{\vec{V}_{\infty}\cdot\hat{n}}{|\vec{V}_{\infty}|}$ (68)

The distribution is shown in Fig 20(a). The abscissa denotes the ratio of arc
length from the leading edge (s) to the chord (c), with positive values
signifying the upper surface and negative values signifying the lower surface.
The collection efficiency is high close to the leading edge (s/c=0) and decays
with distance from the leading edge. This decay is quicker over the upper
surface than the lower surface due to the positive angle of attack. The mass
of droplets impinging at a location is given by

$\dot{m}_{imp}=\text{LWC }V_{\infty}\zeta A$ (69)

Here, $A$ denotes the area of impingement. At given locations that discretise
the wing, we reproduce droplet impingement by placing droplet sources that
generate the volume of water proportional to the collection efficiency. Thus,
the change in volume of the droplets is estimated using the below relation

$\rho_{w}\frac{d\text{Vol}}{dt}=\dot{m}_{imp}$ (70)

where, $\rho_{w}$ and Vol denote the density of water and volume of a droplet
respectively. Having indirectly modelled the effect of aerodynamics through
the collection efficiency $\zeta$, the DDM-S solver proceeds to estimate the
ice height. As seen in Fig. 20(b), the ice heights show a relatively closer
match with the experimental data in comparison to the previous case (Fig. 19)
where $\zeta$ was not considered in the droplet impingement model. This
result, therefore, emphasizes the importance of coupling the aerodynamics of a
bulk-flow with the DDM-S model.

Figure 19: Test Case 5: Comparison of ice profiles at the root section
obtained from the model with experimental data [53].

(a)

(b)

Figure 20:  Test Case 5: (a) Distribution of collection efficiency $\eta$ over
the wing cross-section (b) Comparison of ice profiles at the root section
obtained from the model with experimental data [53].

##  4 Conclusions

The Discrete Droplet Method with Soldification (DDM-S) is based on a
previously proposed liquid thin-film model (DDM) that incorporates
solidification phase-change. The numerical model stores solid height as a
property at each droplet location. The evolution of the solid film thickness
depends on the movement of the solid-liquid interface, as dictated by the
Stefan problem. The movement of the droplet occurs over the solidified region
and the velocity is governed by the momentum equation as in DDM.

DDM-S accounts for unsteadiness in temperature distributions in the liquid and
solid regions of the thin-film. It also allows the liquid temperature to
evolve as per the heat transfer occurring with the surface without the
assumption of perfect thermal contact. It also conveniently incorporates
volumetric change in its framework. One of the unique aspects of the model is
the treatment of the solid-height as a property at a droplet location. This
treatment has some advantages and disadvantages. The advantage is that in this
framework, the solid region can also have its own velocity and this velocity
can be accounted for in the advection velocity of the solid-height advection
equation, Eq. 47. This is done without a new set of numerical points, but
simply using existing set of droplets as the numerical points. However, the
disdvantage with this approach is that solid-height cannot be modelled without
droplets and its evolution as per an advection equation makes it susceptible
to numerical errors such as dissipation and dispersion, that can be countered
only by using higher-order discretizations.

Some of the future directions of research would include coupling the DDM-S
with a bulk-flow Navier-Stokes model to make predictions in the presence of
aerodynamic loads. This would be very useful in aircraft wing and wind turbine
applications. Another interesting direction is to allow for the movement of
the solid phase that could be governed by a different set of governing
equations. This could help in problems like simulation of ice-sheet movement.

## Acknowledgements

Anand S Bharadwaj would like to acknowledge Fraunhofer-Gesellschaft and SERB-
NPDF for providing financial support during the course of the research. The
authors would like to acknowledge Dennis Eck, Jr. of NASA Glenn Research
Centre for providing valuable information about the test case of icing on
GLC-305 swept wing. Pratik Suchde would like to acknowledge funding from the
Institute of Advanced Studies, University of Luxembourg, under the AUDACITY
programme.

##  Appendix A Temperature profile at droplet location

Let us consider the temperature profile in the liquid thin-film of the form

$T(n)=a_{l}n^{2}+b_{l}n+c_{l}$ (71)

The constants $a_{l}$, $b_{l}$ and $c_{l}$ are determined by the conditions

$T(0)=T_{m}$ (72) $\frac{\partial T}{\partial
n}\Big{|}_{n=H_{l}}=-\frac{h_{a}}{\kappa}(T(H_{l})-T_{a})=\alpha$ (73)
$\frac{1}{H_{l}}\int_{0}^{H_{l}}T(n)dn=T_{avg}$ (74)

Here, $n$ is the normal coordinate such that $n=0$ at the solid-liquid
interface. Substituting Eq. 72 in Eq. 71, we get $c_{l}=T_{m}$.
Differentiating Eq. 71 and setting $n=H_{l}$,

$2a_{l}H_{l}+b_{l}=\alpha$ (75)

Eq. 74 results in

$a_{l}\frac{H_{l}^{2}}{3}+b_{l}\frac{H_{l}}{2}=T_{avg}-T_{m}=\beta$ (76)

Solving Eqns. 75 and 76 for $a_{l}$ and $b_{l}$,

$a_{l}=\frac{3}{2H_{l}^{2}}\left(\frac{\alpha H_{l}}{2}-\beta\right)$ (77)
$b_{l}=\left(\frac{3\beta}{H_{l}}-\frac{\alpha}{2}\right)$ (78)

However, $\alpha$ is in terms of $T(H_{l})$ as in Eq. 73, which is still
unknown. Thus, we substitute the above expressions for $a_{l}$ and $b_{l}$
back in Eq. 71 and solve for $T(H_{l})$.

$T(H_{l})=\frac{H_{l}\frac{h_{a}}{4k}T_{a}+\frac{3}{2}(T_{avg}-T_{m})+T_{m}}{1+H_{l}\frac{h_{a}}{4k}}$
(79)

In the above equation, all terms on the RHS are known. Thus, the order of
solving for the unknowns is as follows - $T(H_{l})$, $\alpha$, $\beta$,
$a_{l}$ and $b_{l}$.
A similar approach is adopted for the temperature profile in the solid region
as well. We begin with assuming a temperature profile that varies linearly.

$T(n)=a_{s}n+b_{s}$ (80)

The constants $a_{s}$ and $b_{s}$ are determined by the conditions

$T(H_{s})=T_{m}$ (81) $\frac{\partial T}{\partial
n}\Big{|}_{n=0}=\frac{h_{s}}{\kappa}(T(0)-T_{s})$ (82)

Substituting Eq. 81 in Eq. 80, we get

$T_{m}=a_{s}H_{s}+b_{s}$ (83)

Differentiating Eq. 80 and setting $n=0$,

$a_{s}=\frac{h_{s}}{\kappa}(T(0)-T_{s})$ (84)

Substituting for $a$ in Eq. 83,

$b_{s}=T_{m}-H_{s}\frac{h_{s}}{\kappa}(T(0)-T_{s})$ (85)

Since $T(0)$ is still unknown, we substitute $a_{s}$ and $b_{s}$ in Eq. 80 and
solve for $T(0)$.

$T(0)=\frac{T_{m}+\frac{h_{s}}{\kappa}T_{s}H_{s}}{1+\frac{h_{s}}{\kappa}H_{s}}$
(86)

We, therefore, solve for the unknowns in the following order - $T(0)$, $a_{s}$
and $b_{s}$.

## References

  * [1] Jianbin Luo, Shizhu Wen, and Ping Huang.  Thin film lubrication. part i. study on the transition between ehl and thin film lubrication using a relative optical interference intensity technique.  Wear, 194(1-2):107–115, 1996.
  * [2] KA Overhoff, KP Johnston, J Tam, J Engstrom, and RO Williams III.  Use of thin film freezing to enable drug delivery: a review.  Journal of Drug Delivery Science and Technology, 19(2):89–98, 2009\.
  * [3] Ruifen Jiang and Janusz Pawliszyn.  Thin-film microextraction offers another geometry for solid-phase microextraction.  TrAC Trends in Analytical Chemistry, 39:245–253, 2012.
  * [4] Ehsan Rabiei, Uwe Haberlandt, Monika Sester, and Daniel Fitzner.  Rainfall estimation using moving cars as rain gauges–laboratory experiments.  Hydrology and Earth System Sciences, 17(11):4701–4712, 2013.
  * [5] Tim G Myers.  Extension to the messinger model for aircraft icing.  AIAA journal, 39(2):211–218, 2001.
  * [6] Jennie Molinder, Heiner Körnich, Esbjörn Olsson, Hans Bergström, and Anna Sjöblom.  Probabilistic forecasting of wind power production losses in cold climates: A case study.  Wind Energy Science, 3(2):667–680, 2018.
  * [7] Sawittree Sahakijpijarn, Chaeho Moon, John J Koleng, Dale J Christensen, and Robert O Williams Iii.  Development of remdesivir as a dry powder for inhalation by thin film freezing.  Pharmaceutics, 12(11):1002, 2020.
  * [8] Ki-Soo Seo, Rajiv Bajracharya, Sang Hoon Lee, and Hyo-Kyung Han.  Pharmaceutical application of tablet film coating.  Pharmaceutics, 12(9):853, 2020.
  * [9] Bernard L Messinger.  Equilibrium temperature of an unheated icing surface as a function of air speed.  Journal of the aeronautical sciences, 20(1):29–42, 1953.
  * [10] Tim G Myers and Jean PF Charpin.  A mathematical model for atmospheric ice accretion and water flow on a cold surface.  International Journal of Heat and Mass Transfer, 47(25):5483–5500, 2004.
  * [11] TG Myers and DW Hammond.  Ice and water film growth from incoming supercooled droplets.  International Journal of Heat and Mass Transfer, 42(12):2233–2242, 1999.
  * [12] Tim G Myers, Jean PF Charpin, and S Jon Chapman.  The flow and solidification of a thin fluid film on an arbitrary three-dimensional surface.  Physics of Fluids, 14(8):2788–2803, 2002.
  * [13] TG Myers, JPF Charpin, and CP Thompson.  Slowly accreting ice due to supercooled water impacting on a cold surface.  Physics of fluids, 14(1):240–256, 2002.
  * [14] Giulio Gori, Pietro M Congedo, Olivier Le Maître, Tommaso Bellosta, and Alberto Guardone.  Modeling in-flight ice accretion under uncertain conditions.  Journal of Aircraft, 59(3):799–813, 2022.
  * [15] Ningli Chen, Xian Yi, Qiang Wang, and Yu Liu.  Numerical study on wind-driven thin water film runback on an airfoil.  AIAA Journal, pages 1–9, 2023.
  * [16] Héloïse Beaugendre, François Morency, and Wagdi G Habashi.  Fensap-ice’s three-dimensional in-flight ice accretion module: Ice3d.  Journal of aircraft, 40(2):239–247, 2003.
  * [17] Giulio Gori, Gianluca Parma, Marta Zocca, and Alberto Guardone.  Local solution to the unsteady stefan problem for in-flight ice accretion modeling.  Journal of Aircraft, 55(1):251–262, 2018.
  * [18] Tong Liu, Kun Qu, Jinsheng Cai, and Shucheng Pan.  A three-dimensional aircraft ice accretion model based on the numerical solution of the unsteady stefan problem.  Aerospace Science and Technology, 93:105328, 2019.
  * [19] Rémi Chauvin, Lokman Bennani, Pierre Trontin, and Philippe Villedieu.  An implicit time marching galerkin method for the simulation of icing phenomena with a triple layer model.  Finite Elements in Analysis and Design, 150:20–33, 2018.
  * [20] Tim G Myers, Sarah L Mitchell, G Muchatibaya, and MY Myers.  A cubic heat balance integral method for one-dimensional melting of a finite thickness layer.  International Journal of Heat and Mass Transfer, 50(25-26):5305–5317, 2007.
  * [21] Yasir A Malik, Lokman Bennani, Stephan Bansmer, Pierre Trontin, and Philippe Villedieu.  Experimental and numerical investigation of accretion inception and heat transfer physics in ice crystal icing.  International Journal of Heat and Mass Transfer, 214:124364, 2023\.
  * [22] Ellen Norde, Jacco M Hospers, Edwin van der Weide, and Harry W Hoeijmakers.  Splashing model for impact of supercooled large droplets on a thin liquid film.  In 52nd Aerospace Sciences Meeting, page 0738, 2014.
  * [23] Ellen Norde, Edwin Theodorus Antonius van der Weide, and Hendrik Willem Marie Hoeijmakers.  Eulerian method for ice crystal icing.  AIAA journal, 56(1):222–234, 2018.
  * [24] Chao Wang, Shinan Chang, Mengyao Leng, Hongwei Wu, and Bo Yang.  A two-dimensional splashing model for investigating impingement characteristics of supercooled large droplets.  International Journal of Multiphase Flow, 80:131–149, 2016.
  * [25] Gary A Ruff and Brian M Berkowitz.  Users manual for the nasa lewis ice accretion prediction code (lewice).  Technical report, 1990.
  * [26] Philippe Villedieu, Pierre Trontin, Didier Guffond, and David Bobo.  Sld lagrangian modeling and capability assessment in the frame of onera 3d icing suite.  In 4th AIAA Atmospheric and Space Environments Conference, page 3132, 2012.
  * [27] William B Wright, Peter Struk, Tadas Bartkus, and Gene Addy.  Recent advances in the lewice icing model.  2015\.
  * [28] Edward T Rocco, Yiqiang Han, Richard Kreeger, and Jose Palacios.  Super-cooled large droplet experimental reproduction, ice shape modeling, and scaling method assessment.  AIAA journal, 59(4):1277–1295, 2021.
  * [29] Thierry Hedde and Didier Guffond.  Onera three-dimensional icing model.  AIAA journal, 33(6):1038–1045, 1995.
  * [30] Pierre Trontin, Ghislain Blanchard, Alexandros Kontogiannis, and Philippe Villedieu.  Description and assessment of the new onera 2d icing suite igloo2d.  In 9th AIAA Atmospheric and Space Environments Conference, page 3417, 2017.
  * [31] William B Wright.  User manual for the nasa glenn ice accretion code lewice. version 2.2. 2.  Technical report, 2002.
  * [32] William Wright, Phil Jorgenson, and Joseph Veres.  Mixed phase modeling in glennice with application to engine icing.  In AIAA Atmospheric and Space Environments Conference, page 7674, 2010.
  * [33] William Wright, Mark Potapczuk, and Laurie Levinson.  Comparison of lewice and glennice in the sld regime.  In 46th AIAA Aerospace Sciences Meeting and Exhibit, page 439, 2008\.
  * [34] Yves Bourgault, Héloïse Beaugendre, and Wagdi G Habashi.  Development of a shallow-water icing model in fensap-ice.  Journal of Aircraft, 37(4):640–646, 2000.
  * [35] Lasse Makkonen.  Models for the growth of rime, glaze, icicles and wet snow on structures.  Philosophical Transactions of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences, 358(1776):2913–2939, 2000.
  * [36] Lukas Strauss, Stefano Serafin, and Manfred Dorninger.  Skill and potential economic value of forecasts of ice accretion on wind turbines.  Journal of Applied Meteorology and Climatology, 59(11):1845–1864, 2020.
  * [37] Neil Davis, Andrea N Hahmann, Niels-Erik Clausen, and Mark Žagar.  Forecast of icing events at a wind farm in sweden.  Journal of Applied Meteorology and Climatology, 53(2):262–281, 2014\.
  * [38] Cristina P Martin-Linares, Yorgos M Psarellis, Georgios Karapetsas, Eleni D Koronaki, and Ioannis G Kevrekidis.  Physics-agnostic and physics-infused machine learning for thin films flows: modeling, and predictions from small data.  arXiv preprint arXiv:2301.12508, 2023.
  * [39] Jennie Molinder, Sebastian Scher, Erik Nilsson, Heiner Körnich, Hans Bergström, and Anna Sjöblom.  Probabilistic forecasting of wind turbine icing related production losses using quantile regression forests.  Energies, 14(1):158, 2020.
  * [40] Markus Kreutz, Abderrahim Ait Alla, Michael Lütjen, Jan-Hendrik Ohlendorf, Michael Freitag, Klaus-Dieter Thoben, Florian Zimnol, and Andreas Greulich.  Ice prediction for wind turbine rotor blades with time series data and a deep learning approach.  Cold Regions Science and Technology, 206:103741, 2023.
  * [41] Koji Fukudome, Yusuke Muto, Ken Yamamoto, Hiroya Mamori, and Makoto Yamamoto.  Numerical simulation of the solidification phenomena of single molten droplets impinging on non-isothermal flat plate using explicit moving particle simulation method.  International Journal of Heat and Mass Transfer, 180:121810, 2021\.
  * [42] Huimin Liu, Enrique J Lavernia, and Roger H Rangel.  Numerical simulation of substrate impact and freezing of droplets in plasma spray processes.  Journal of Physics D: Applied Physics, 26(11):1900, 1993.
  * [43] Truong V Vu, Khoa V Dao, and Binh D Pham.  Numerical simulation of the freezing process of a water drop attached to a cold plate.  Journal of Mechanical Science and Technology, 32:2119–2126, 2018\.
  * [44] Anand S Bharadwaj, Joerg Kuhnert, Stéphane PA Bordas, and Pratik Suchde.  A discrete droplet method for modelling thin film flows.  Applied Mathematical Modelling, 112:486–504, 2022.
  * [45] JM Domínguez, AJC Crespo, M Gómez-Gesteira, and JC2869628 Marongiu.  Neighbour lists in smoothed particle hydrodynamics.  International Journal for Numerical Methods in Fluids, 67(12):2026–2042, 2011.
  * [46] Christian Drumm, Sudarshan Tiwari, Jörg Kuhnert, and Hans-Jörg Bart.  Finite pointset method for simulation of the liquid–liquid flow field in an extractor.  Computers & Chemical Engineering, 32(12):2946–2957, 2008.
  * [47] Joseph Stefan.  Über einige probleme der theorie der wärmeleitung.  Sitzungber., Wien, Akad. Mat. Natur, 98:473–484, 1889.
  * [48] Timothy G Myers, Matthew G Hennessy, and Marc Calvo-Schwarzwälder.  The stefan problem with variable thermophysical properties and phase change temperature.  International Journal of Heat and Mass Transfer, 149:118975, 2020\.
  * [49] Felix R Saucedo-Zendejo, Edgar O Reséndiz-Flores, and Jörg Kuhnert.  Three-dimensional flow prediction in mould filling processes using a gfdm.  Computational Particle Mechanics, 6(3):411–425, 2019.
  * [50] Pratik Suchde and Jörg Kuhnert.  Point cloud movement for fully lagrangian meshfree methods.  Journal of Computational and Applied Mathematics, 340:89–100, 2018\.
  * [51] Christian Stolze, Tobias Janoschka, Ulrich S Schubert, Frank A Müller, and Stefan Flauder.  Directional solidification with constant ice front velocity in the ice-templating process.  Advanced Engineering Materials, 18(1):111–120, 2016.
  * [52] Zhou Chao and Yin Jiaqi.  Glaze icing process of moveable overhead conductor and its aerodynamic characteristics with numerical method.  International Journal of Heat and Mass Transfer, 176:121436, 2021\.
  * [53] Michael Papadakis, Hsiung Wei Yeong, Mario Vargas, and Mark Potapczuk.  Aerodynamic performance of a swept wing with ice accretions.  In 41st Aerospace Sciences Meeting and Exhibit, page 731, 2003.
  * [54] Yihua Cao, Chao Ma, Qiang Zhang, and John Sheridan.  Numerical simulation of ice accretions on an aircraft wing.  Aerospace Science and Technology, 23(1):296–304, 2012.