# Predicting Energy Budgets in Droplet Dynamics: A Recurrent Neural Network
Approach

Diego A. de Aguiar<EMAIL_ADDRESS>ORCID: 0000-0001-5159-867X
Departamento de Matemática e Computação, Faculdade de Ciências e Tecnologia,
Universidade Estadual Paulista “Júlio de Mesquita Filho”, Presidente Prudente,
Brazil  Hugo L. França<EMAIL_ADDRESS>ORCID: 0000-0002-5361-7704
Van der Waals-Zeeman Institute, Institute of Physics, University of Amsterdam,
Amsterdam, The Netherlands  Cassio M<EMAIL_ADDRESS>ORCID:
0000-0002-5361-7704 Departamento de Matemática e Computação, Faculdade de
Ciências e Tecnologia, Universidade Estadual Paulista “Júlio de Mesquita
Filho”, Presidente Prudente, Brazil

(October 2023)

Neural networks in fluid mechanics offer an efficient approach for exploring
complex flows, including multiphase and free surface flows. The recurrent
neural network, particularly the Long Short-Term Memory (LSTM) model, proves
attractive for learning mappings from transient inputs to dynamic outputs.
This study applies LSTM to predict transient and static outputs for fluid
flows under surface tension effects. Specifically, we explore two distinct
droplet dynamic scenarios: droplets with diverse initial shapes impacting with
solid surfaces, as well as the coalescence of two droplets following
collision. Using only dimensionless numbers and geometric time series data
from numerical simulations, LSTM predicts the energy budget. The marker-and-
cell front-tracking methodology combined with a marker-and-cell finite-
difference strategy is adopted for simulating the droplet dynamics. In
addition, as the code must handle various drop shapes, we have also employed a
method designed for extracting particles from interfaces in digital images.
Using a recurrent neural network (RNN) architecture fed with time series data
derived from geometrical parameters, as for example droplet diameter
variation, our study shows the accuracy of our approach in predicting energy
budgets, as for instance the kinetic, dissipation, and surface energy trends,
across a range of Reynolds and Weber numbers in droplet dynamic problems.
Finally, a two-phase sequential neural network using only geometric data,
which is readily available in experimental settings, is employed to predict
the energies and then use them to estimate static parameters, such as the
Reynolds and Weber numbers. While our methodology has been primarily validated
with simulation data, its adaptability to experimental datasets is a promising
avenue for future exploration. Finally, we hope that our strategy can be
useful for diverse applications, spanning from inkjet printing to combustion
engines, where the prediction of energy budgets or dissipation energies is
crucial. This analysis provides insights for developing machine learning
modeling using numerical and/or experimental data from droplet dynamics.

keywords: Droplets $|$ Numerical solution $|$ Prediction $|$ Surface Tension
$|$ LSTM $|$ Energy budget

## I. Introduction

Droplet dynamics, a phenomenon of remarkable complexity and practical
significance, extends its influence across various domains of science and
engineering, ranging from raindrop splashes on leaves to the precise
deposition of microdroplets in advanced manufacturing processes. Achieving a
fundamental understanding of the dynamics governing droplet impact represents
a significant scientific challenge, needing the collaborative engagement of
experimental, analytical and numerical methodologies. This concerted effort is
essential to unravel the complexities inherent in the complex phenomenon of
droplet impact dynamics.

Advancements in these methodologies have unveiled a wealth of intricate
physical mechanisms governing liquid-drop impacts on rigid surfaces. These
insights not only contribute to the fundamental understanding within the realm
of fluid dynamics but also create innovations across a diverse array of
applications. The implications of droplet impact dynamics reverberate through
fields such as aerospace, energy, materials science, biotechnology, forensic
analysis, and studies related to food and cosmetics, among others.

In recent decades, the exploration of the droplet dynamics has yielded a rich
body of literature, including excellent review articles that provide
comprehensive insights into the multifaceted aspects of droplet dynamics.
Notable contributions in this context have been made, as well detailed in [1,
2, 3, 4]. These articles clarify some important aspects as for instance: the
dimensionless numbers and their influence on the maximum spreading of the
droplet, different regimes considering pre- and post-impact and their
consequences in terms of energy mechanisms, etc.

Concerning droplet dynamics and its applications, the complexities of the
inkjet printing process stand out as a pertinent challenge. This process
involves various fluid dynamics challenges that significantly impact the
outcome of the printing process. A comprehensive overview of these challenges,
particularly focusing on fundamental aspects, has been extensively discussed
in the review paper by [5]. A fundamental question of inkjet printing
challenge lies in the collision dynamics of droplets during their trajectory
from ejection to eventual impact on the substrate. The understanding of this
collision process is important, influencing the spatial distribution, size,
and morphology of the printed droplets. As droplets traverse through the air,
interactions can be observed, governed by fluid dynamics principles, as for
instance coalescence, fragmentation, bouncing and other complex behaviors.
After this, a different dynamic happens when the droplets impact on the
substrate. Factors such as droplet velocity, angle of impact, and the
properties of the drop and substrate all contribute to the final result.
Moreover, the subsequent spreading of the droplets upon contact with the
substrate adds another layer of complexity to the overall process. The
dynamics of spreading govern how the ink interacts with the substrate surface,
influencing factors as for instance the shape of the droplet. In essence, the
challenges within inkjet printing underscore the profound interplay of fluid
dynamics phenomena, encompassing droplet collisions, substrate interactions,
and spreading dynamics. In this work, we have focused on two stages described
in the review paper [5]: the impact of a droplet with different shapes in a
substrate and the coalescence of droplets during the collision process.

Over the past decades, the simulation of droplet dynamics has been adopted for
advancing our understanding of complex fluid phenomena [6]. These numerical
investigations have opened new horizons, allowing for innovative approaches in
the study of droplet behaviors and interactions. In these computational
frameworks, the governing equations are rigorously solved on fixed or adaptive
meshes, employing classical discretization schemes such as finite-element
methods [7], finite volume methods [8], finite difference methods [9], among
others. A critical computational consideration in simulating the impact and
collisions of droplets is the treatment of the moving interface. This involves
the solution of the adequate boundary conditions, coupled with the accurate
representation and dynamic updating of the interface throughout the transient
problem. Traditional mesh-based frameworks, commonly employed in the
simulation of free surface or multiphase fluid flows, include the well-
established front-tracking and front-capturing methodologies [10]. Beyond
mesh-based approaches, mesh-free schemes have gained widespread application in
simulating droplet dynamics. Smoothed Particle Hydrodynamics (SPH), for
instance, offers an efficient and effective alternative of solving fluid
dynamics problems without relying on a fixed grid allowing the solution of
complex droplet interactions and deformations (see [11]). It is crucial to
highlight the significance of both analytical papers within the context of
impacting droplets [12] and experimental studies focusing on droplet dynamics
[13]. These results play an important role in not only validating but also
enriching the fidelity of computational simulations. Analytical works provide
essential theoretical foundations, offering insightful equations and
principles that guide numerical model development, as for instance the
spreading of the droplet in a substrate [14]. Simultaneously, experimental
investigations have been adopted to establish different scaling laws in the
spreading of droplets [15, 16, 17, 18, 4, 19] and in droplet collisions [20,
21, 22, 23, 24]. The synergy between these analytical, experimental, and
computational methods enhanced the understanding of droplet interactions [25,
26, 27].

More recently, the employment of Machine Learning (ML) strategies in the
context of droplet dynamics received more attention due to the advances of new
computation technologies, the increasing data availability and the
implementation facilities on open-source codes (e.g. TensorFlow, Pytorch,
MXNet, etc). The ML-based methods open a promising avenue to reduce both time
and costs associated with traditional numerical simulations. This approach not
only introduces efficiency gains but also serves as a robust alternative for
predicting outcomes, reflecting the growing synergy between data-driven
machine learning and traditional computational fluid dynamics. For example, in
the context of droplet dynamics, different machine learning models have been
adopted in order to predict physical or geometrical properties observed in
complex flows, such as predicting the maximum spreading of a droplet by using
the dimensionless numbers and properties related to each regime [28, 29, 30,
31, 32, 33]. The idea of applying related features to predict fluid properties
and/or characteristics was also explored by predicting drop morphologies
adopting image data [34], classifying splashing and nonsplashing drops using
images from when half of the drop impacted the surface [35], using an
explainable artificial intelligence (XAI) to interpret how said classifier
works and uncover a possible correlation between the spreading of a drop and
the normalized impact force on a solid surface [36], predicting interface
curvature through neighboring particles from 2D front-tracking representation
[37], predicting size, velocity and shape characteristics of droplets from
droplet-on-demand liquid metal jetting (DoD-LMJ) [38], etc.

The connection between droplet dynamics and the energy budget has been
extensively documented. When a droplet approaches a surface or interacts with
another droplet, it initiates a process of energy transfer involving kinetic,
surface, and dissipative components. Examining the energy budget in droplet
dynamics not only enhances our understanding of fundamental fluid dynamics
principles but also facilitates the development of simple analytical models
based on energy balance. For instance, expressions derived from energy
conservation can be employed to construct equations predicting the maximum
diameter of a droplet spreading on a rigid surface. Seminal works by [15],
[39], [40] and [19] exemplify the application of such energy balance
principles.

Exploring the energy budget in droplet dynamics is of practical significance,
offering valuable insights for optimizing processes such as inkjet printing.
Achieving precise control over droplet impact dynamics is crucial for
obtaining desired outcomes in applications like these. Recent literature has
presented several papers delving into the energy budget in complex droplet
dynamics. For instance, in the context of impinging droplets on rigid
surfaces, [41] utilized a Volume-of-Fluid code to investigate the effects of
low and high velocities on surface energy and viscous dissipation.
Simultaneously, [42] conducted an energy-budget analysis using the diffusive-
interface method, with a focus on a two-dimensional case. Additional insights
into the connection between splash formation and energy dissipation were
provided by [43], who employed experimental results involving a ferrofluid
droplet with various shapes. In droplet collision studies, [44] applied a
front-tracking numerical scheme to study the energy budget using unequal-size
droplets while [45] combined an energy-based strategy and binary droplet
collisions through experimental setups to reanalyze the bouncing regime
boundary on the regime map proposed by [22]. More recently, [46] explored the
energy balance during coalescence and breakup after collision, utilizing an
experimental process for droplet atomization. Additionally, [47] investigated
electro-mechanical energy conversion in droplets.

Therefore, based on the literature review described in the previous statement,
it is evident that numerous studies have focused on various aspects, such as
the maximum spreading diameter, different scaling laws, and boundaries on the
regime maps of droplet collisions, using arguments of energy conservation.
However, these results are typically presented in the context of static
analyses, as there are few works [48, 41] that deal with the transient
variation of properties (such as the spreading of a droplet impacting a
surface or the diameter and height of a coalesced drop after a drop-drop
collision) and their relation to the transient energy budget. Notably, recent
investigations using learning algorithms are also limited to static data
analysis. In this work, we propose a machine learning model based on a
recurrent neural network to predict the temporal variation of the energy
budget. We have adopted the Long Short-Term Memory (LSTM) [49] with a
geometric input data set ${z(t_{1}),\dots,z(t_{n})}$, which learns a mapping
to produce temporal predictions in terms of the energy budget. The term
“geometric” refers to the time variation of the droplet diameter during
spreading on a surface or the diameter and height of the coalesced droplet
after the collision of two droplets. Expanding the potential of our framework
utilizing LSTM, incorporating fluid parameters such as the Reynolds and Weber
numbers, as well as flow parameters like impact parameter on the collision,
enhances the accuracy of predictions for the transient behavior of kinetic,
surface energies, and viscous dissipation. Although our tests have exclusively
employed numerical results, the versatility of this non-intrusive methodology
in predicting energy budgets during intricate droplet dynamics can be further
extended to experimental investigations. This extension is feasible since it
relies solely on temporal series data, specifically variations in droplet
diameter, which can be easily captured by experimental devices.

## II. Methods

### A. Mathematical model and overview of the numerical method

The equations that govern the flow behavior of isothermal and incompressible
fluids are the continuity (1) and momentum (2) equations:

$\displaystyle\nabla\cdot\mathbf{u}=0,$ (1)
$\displaystyle{\rho\left(\frac{\partial\mathbf{u}}{\partial
t}+\nabla\cdot(\mathbf{u}\mathbf{u})\right)=-\nabla p+\nabla\cdot\bm{\tau}},$
(2)

where $\mathbf{u}$ and $p$ are the velocity and pressure fields and $\rho$ is
the fluid’s density. The stress tensor $\bm{\tau}$ assumes a form of a
Newtonian liquid as

$\bm{\tau}=2\eta\bm{\dot{\gamma}},$ (3)

where $\eta$ is the viscosity of the fluid and
$\bm{\dot{\gamma}}=\frac{1}{2}(\nabla\mathbf{u}+(\nabla\mathbf{u})^{T})$ is
the strain-rate tensor with $\nabla\mathbf{u}$ being the velocity gradient.

A dimensionless version of the previous equations can be obtained by scaling
the variables as follows

$\bar{\mathbf{x}}=\frac{\mathbf{x}}{D},\hskip
10.0pt\bar{t}=\frac{tU}{D},\hskip
10.0pt\mathbf{\bar{u}}=\frac{\mathbf{u}}{U},\hskip 10.0pt\bar{p}=\frac{p}{\rho
U^{2}},\hskip 10.0pt\bar{\boldsymbol{\tau}}=\frac{D\boldsymbol{\tau}}{\eta
U},$ (4)

using $U$ and $D$ as characteristic velocity and length, respectively. To
simplify the nondimensionalization of the governing equations for isothermal
and incompressible flows, we redefined them by using the previous
dimensionless terms without its bars, as follows

$\displaystyle\nabla\cdot\mathbf{u}=0,$ (5)
$\displaystyle\frac{\partial\mathbf{u}}{\partial
t}+\nabla\cdot(\mathbf{u}\mathbf{u})=-\nabla
p+\frac{1}{Re}\nabla^{2}\mathbf{u},$ (6)

where $Re=\frac{\rho UD}{\eta}$ is the dimensionless Reynolds number.

In our tests, there are two main boundary conditions to be considered: solid
walls condition, to define behavior across domain boundaries, and the free
surface condition, for the contact between fluid and the external environment
(air). In summary, a no-slip condition is applied for solid wall contours,
while the free surface condition is defined as

$\displaystyle\left\\{\begin{array}[]{ll}\mathbf{n}\cdot(\mathbf{T}\cdot\mathbf{n})=\gamma\kappa,\\\
\mathbf{m}\cdot(\mathbf{T}\cdot\mathbf{n})=0,\end{array}\right.$ (9)

where $\mathbf{n}$ and $\mathbf{m}$ are, respectively, the normal and
tangential vectors to the free surface, $\gamma$ is the surface tension
coefficient, $\kappa$ is the interface curvature and
$\mathbf{T}=\boldsymbol{\tau}-p\mathbf{I}$ is the total stress tensor. The
nondimensionalized forms of (9) can be written as

$\displaystyle\left\\{\begin{array}[]{ll}\mathbf{n}\cdot(\bar{\mathbf{T}}\cdot\mathbf{n})=\frac{1}{We}\bar{\kappa},\\\
\mathbf{m}\cdot(\bar{\mathbf{T}}\cdot\mathbf{n})=0,\end{array}\right.$ (12)

where $\bar{\kappa}=D\kappa$,
$\mathbf{\bar{T}}=2\boldsymbol{\bar{\overset{.}{\gamma}}}-Re\cdot\bar{p}\mathbf{I}$
and $We=\frac{\rho U^{2}D}{\gamma}$ is the Weber number. Furthermore, an extra
dimensionless number is adopted for the collision regime, the impact parameter
$B=\frac{2b}{D_{1}+D_{2}}$, where $b$ is the center-to-center relative
distance between the drops and $D_{1}$, $D_{2}$ are the diameters for each
drop, as seen in [45].

To solve the governing equations numerically, we have employed the established
projection method. In summary, the method uses the Helmholtz-Hodge
decomposition to decouple velocity and pressure from equations (5) and (6). In
the Helmholtz-Hodge decomposition, as presented in [50], any vector field
defined in a region $\Omega$ with smooth boundary $\partial\Omega$ can be
uniquely decomposed as

$\tilde{\mathbf{u}}=\mathbf{u}+\nabla\psi,$ (13)

where $\mathbf{u}$ is a vector field such that $\nabla\cdot\mathbf{u}=0$ and
$\psi$ is a scalar field also defined in $\Omega$.

Initially, the projection method solves the momentum equation (6) for an
intermediate velocity field $\tilde{\mathbf{u}}$. For the purpose of
decoupling velocity and pressure, the unknown $p$ is approximated by some
known scalar field $\tilde{p}$, resulting in the following equation

${\frac{\partial\tilde{\mathbf{u}}}{\partial
t}+\nabla\cdot(\tilde{\mathbf{u}}\tilde{\mathbf{u}})=-\nabla\tilde{p}+\frac{1}{Re}\nabla^{2}\tilde{\mathbf{u}}},$
(14)

with appropriated boundary conditions. After solving equation (14) and
obtaining the velocity field $\tilde{\mathbf{u}}$, we still need to make sure
the continuity equation will be satisfied. By using the Helmholtz-Hodge
decomposition and equations (13) and (5), we can write the following Poisson
equation for a scalar field $\psi$

$\nabla^{2}\psi=\nabla\cdot\tilde{\textbf{u}}.$ (15)

Finally, by solving equations (14) and (15) we obtain $\tilde{\mathbf{u}}$ and
$\psi$ so that the decomposition (13) can be used to obtain a final velocity
field $\mathbf{u}$ that satisfies both the momentum and continuity equations.
The discretization process is achieved using the finite-difference method for
free surface problems as well discussed in [51]. It is important to mention
that we have employed an axisymmetric formulation to solve the droplet
impacting problem, while droplet collision simulations were performed using a
two-dimensional code.

### B. Energy budget

Since there will be no gravitational influence on the generated data, the
energy budget is given by the kinetic energy, the surface energy and the
viscous dissipation. The kinetic energy $E_{k}$ is defined as

$E_{k}=\int_{V}\frac{1}{2}\rho\left\|\mathbf{u}\right\|^{2}dV,$ (16)

calculated by the integral over the fluid volume V. Numerically, it is
calculated as the sum of the kinetic energy inside all cells fully inside a
droplet. The surface energy $E_{s}$ is computed as

$E_{s}=\int_{S}\gamma dS.$ (17)

calculated by the integral over the fluid surface $S$. Lastly, based on the
viscous dissipation, we can define

$E_{d}=\int_{t}\int_{V}2\hskip 2.84544pt\eta\hskip
2.84544pt\left\|\overset{.}{\bm{\gamma}}\right\|^{2}dV\hskip 2.84544ptdt,$
(18)

where $\int_{V}2\hskip 2.84544pt\eta\hskip
2.84544pt\left\|\overset{.}{\bm{\gamma}}\right\|^{2}dV$ is called the viscous
dissipation rate.

The law of conservation of energy defines that the total sum of all the
energies ($E_{t}$) in the system must remain constant through time so that

$E_{t}(t)=E_{k}(t)+E_{s}(t)+E_{d}(t)=E_{t}(0).$ (19)

In order to construct a nondimesional investigation, we have adopted the
following nondimensional energy quantities:

$\bar{E}_{k}=ReWe\int_{\bar{V}}\frac{1}{2}\left\|\bar{u}\right\|^{2}d\bar{V},$
$\bar{E}_{s}=Re\int_{\bar{S}}1\hskip 2.84544ptd\bar{S},$
$\bar{E}_{d}=We\int_{\bar{t}}\int_{\bar{V}}2\left\|\bar{\overset{.}{\bm{\gamma}}}\right\|^{2}d\bar{V}\hskip
2.84544ptd\bar{t}.$

The focus of this work is to use geometric time series data extracted from
numerical simulations to generate time series predictions for each type of
energy. From now on, we have omitted the bar over the nondimensional
variables.

### C. Problem descriptions and data preparation

To apply the neural network solution, we first need to generate adequate
datasets. In order to do so, numerical simulations are used to generate two
different flow regimes: drop impact on a solid surface, to generate diverse
spreading lengths with different initial drop shapes, and droplet collision,
to test the capacity of the model to predict features that oscillate through
time. Both problem descriptions are represented in Fig. 1 (left) while the
schematic to obtain the dataset is also described in the same figure (right).

Fig. 1: Dataset is extracted from numerical simulations of two droplet dynamic
problems.

In order to generate the dataset, each case is executed with varying fluid and
flow parameters, and/or initial drop shapes. In the droplet impact scenario,
we vary the droplet shape by extracting interfaces from experimental droplet
images [43]. Conversely, for the droplet collision problem, we adopt the
impact parameter B to induce variations. In both cases, modifications are also
conducted changing the fluid parameters based on Reynolds and Weber numbers.
In summary, samples within the dataset can be defined as individual data
points, each characterized by a specific combination of features encompassing
fluid and flow parameters as well as initial drop characteristics, as it will
be detailed in the next sections.

####  C.1 Spreading dataset

For the first dataset, we were inspired by the drop shapes from experiments
presented by [43]. To test the efficiency of the machine learning model, the
shapes were chosen in such a way that the train, validation and test subsets
would have sufficiently different shapes, as shown in Fig. 2.

Fig. 2: Different initial drop shapes and the train-validation-test splits.

To extract particles from these shapes and simulate each case, an image
interface extraction method was implemented, as described in Appendix A. The
initial diameter and height of each shape, as presented in Table 1, are
selected in such a way that each droplet has approximately the same volume.
Since we used an axisymmetric simulation for the spreading case, as shown in
Section A, this is necessary to make each droplet have approximately the same
mass.

$Shape$ |  $D(0)$ (in mm) |  $H(0)$ (in mm)
---|---|---
1 |  3.92 |  3.98
2 |  3.6771 |  4.64
3 |  2.62 |  4.64
4 |  4.64 |  3.6829
5 |  4.6086 |  3.1343
6 |  4.8143 |  2.5714
7 |  3.8229 |  4.24
8 |  4.5057 |  3.3857
9 |  3.5571 |  5.0086
10 |  4.6343 |  3.6229
Table 1: The initial spreading length and height for each drop shape.

The parameters chosen to configure the numerical simulations are presented in
(Table 2). For each set of parameters and one droplet shape, a numerical
simulation is run to generate one sample. Some combinations have fewer samples
due to potential numerical errors in the solver, resulting in the failure to
generate them. The Reynolds and Weber numbers were arbitrarily defined by
varying the impact velocity. In total, 132 samples were extracted.

Spreading dataset
---
$Samples$ |  Re |  We
10 |  1.5063 |  240.625
10 |  4.5188 |  2165.625
10 |  5.875 |  212.5
10 |  10.2438 |  184.375
10 |  17.6250 |  1912.5
10 |  20.4375 |  118.75
10 |  24.8063 |  90.625
9 |  30.6313 |  53.125
10 |  30.7313 |  1659.375
7 |  35 |  25
10 |  61.3125 |  1068.75
9 |  74.4188 |  815.625
9 |  91.8938 |  478.125
8 |  105 |  225
Table 2: Number of samples for each group of parameters used in the spreading
simulations.

To visualize the spreading behavior while varying Re and We, we have presented
in Fig. 3 the droplet morphology and the temporal evolution of the diameter
for shapes 1 and 10. As depicted in this figure, increasing the Reynolds
number results in a greater spread for both shapes. This behavior is
quantified by the subplots in Fig. 3 depicting the time variation of the
droplet diameter. Figs. 3 and 3 illustrate snapshots comparing two shapes for
two different sets of dimensionless numbers, while Figs. 3 and 3 compare the
evolution of the spreading diameter through time for both shapes. It is
evident from Figs. 3 and 3 that at lower Reynolds numbers, the droplets spread
more slowly compared to the second test involving higher $Re$.

Fig. 3: Effect of drop shape (shapes 1 and 10 from Fig 2) on impact when
varying Re and We: morphology for the first (a) and second (b) set of
parameters and diameter time variation for both shapes from the first (c) and
second set (d).

####  C.2 Collision dataset

The second dataset consists of two drops of 1mm in diameter colliding
horizontally with one another. The pair of drops merge after collision and
generates horizontal and vertical oscillation over time, helping to create a
time series that is harder to predict and that can be used to test the machine
learning model’s efficiency.

To define the numerical simulation parameters for each sample, cases similar
to those presented in the experimental coalescence map from [45] were selected
by varying the impact velocity. For each set of Reynolds and Weber numbers, 20
simulations are run by varying the impact parameter B in the intervals shown
in Table 3. Once again, certain samples were not successfully run due to
solver errors. In total, 72 samples were generated.

Collision dataset
---
$Samples$ |  Re |  We |  B (interval)
20 |  74.451 |  22.0 |  $[0.0,0.4]$
20 |  85.4788 |  29.0 |  $[0.0,0.4]$
9 |  97.8478 |  38.0 |  $[0.0,0.2316]$
5 |  104.0863 |  43.0 |  $[0.0211,0.1684]$
7 |  109.9715 |  48.0 |  $[0.0421,0.1684]$
6 |  115.5573 |  53.0 |  $[0.0421,0.2105]$
5 |  120.8853 |  58.0 |  $[0.0,0.2316]$
Table 3: Number of samples for each group of parameters used in the collision
simulations.

To visualize the difference in horizontal and vertical spreading after impact
resulting from changes in the physical parameters and the impact parameter,
snapshots from some simulated samples are shown in Fig. 4. Figs. 4, 4, and 4
illustrate snapshots of two coalescing droplets on three different sets of
dimensionless numbers, with pairs 4 and 4 varying only the Reynolds and Weber
numbers, and 4 and 4 varying only the impact parameter. For Figs. 4, 4, and 4,
we compare the time series of diameter and height oscillation for each
example.

Fig. 4: Droplet collision dynamics with different impact parameters and
varying Re and We: morphology for the first (a), second (b) and third (c) set
of parameters and time variation in diameter $D$ and height $H$ for the first
(d), second (e) and third set (f).

### D. Machine learning model

Since the transient features in the extracted data have a certain degree of
time-correlation, it would be beneficial to use a neural network architecture
that can capture it effectively. For this purpose, a Recurrent Neural Network
[52], which can be used to find correlations on data sequences, is a good
option.

The drawback of a simple RNN architecture is the likely introduction of the
vanishing gradient problem [53], considering that the data has a large number
of time steps. To avoid this, the LSTM [54], a type of RNN architecture
employed to efficiently learn long-term dependencies in data, is chosen to be
implemented.

The LSTM architecture introduces two new concepts to increase its capacity to
deal with more time steps: a cell state to capture long-term information and
three gates that control the intensity of the data passing through them. The
structure of an LSTM cell and all of its additional concepts are represented
in Fig. 5.

Fig. 5: LSTM Cell structure.

In contrast to the hidden state, which mostly carries information from the
previous time step, the cell state is capable of passing data through
multiples cells, thanks to its contents suffering only minor reductions and/or
increments between each time step.

The gate mechanisms are constructed by a sigmoid neural network layer
responsible for transforming a certain input into outputs with values in the
range $[0,1]$. These values are called keys, since they are used to tune the
data passing through the gates, with 0 erasing the data, 1 passing it
integrally and $]0,1[$ passing it partially.

The first gate, also known as the forget gate, is responsible for controlling
how much each information from the previous cell state $C_{t-1}$ will pass on
to $C_{t}$. Its keys ($f_{t}$) are defined as

$\displaystyle f_{t}=\sigma(W_{f}\hskip
2.84544pt\cdot[\bar{y}_{t-1},x_{t}]+b_{f}),$ (20)

where $\sigma$ is the sigmoid activation function, $W_{f}$ is a weight matrix,
$[\bar{y}_{t-1},x_{t}]$ is the concatenation of the previous cell’s output
with the current cell’s input and $b_{f}$ is a bias.

The second gate, known as the input gate, controls the intensity with which
the information obtained from the current cell will be added to the new cell
state. The new information candidate $\tilde{C}_{t}$ and the keys $i_{t}$ that
regulate its intensity are defined as follows

$\displaystyle i_{t}=\sigma(W_{i}\hskip
2.84544pt\cdot[\bar{y}_{t-1},x_{t}]+b_{i})$ (21)
$\displaystyle\tilde{C}_{t}=tanh(W_{C}\hskip
2.84544pt\cdot[\bar{y}_{t-1},x_{t}]+b_{C}),$ (22)

with that, all the mechanisms necessary to define the new cell state $C_{t}$
were described. While the forget gate controls how much of the old cell state
data will continue to pass, the input gate controls how much of the new
information will be added to the cell state. To update it, the expression is
defined as

$\displaystyle C_{t}=f_{t}*C_{t-1}+i_{t}*\tilde{C}_{t}.$ (23)

The third and final gate, also called the output gate, is responsible for
defining how much of the cell state information will be passed to the output
$\bar{y}_{t}$ of the cell and the hidden state. Before using the information
from the cell state, the hyperbolic tangent function is applied to reduce the
range of its values. The keys ($o_{t}$) and the final output of the cell
($\bar{y}_{t}$) are given by

$\displaystyle o_{t}=\sigma(W_{o}[\bar{y}_{t-1},x_{t}]+b_{o}),$ (24)
$\displaystyle\bar{y}_{t}=o_{t}*tanh(C_{t}).$ (25)

Now that the main features of the LSTM architecture and the motivation for
using it have been highlighted, the next step is to organize the overall
machine learning architecture that will be used to predict certain features.

The most interesting features to predict, whether in a numerical or in a
experimental sense, are divided in two types depending on the output type,
being static predictions, where each predicted feature has only one value, and
transient predictions, where each predicted feature has its values changing
through time. Our work focuses on the transient predictions, but an
exploration of static predictions is also provided in Section C.2.

For the input, there can be both static and transient features. Since the
input dimensions need to be the same for every feature, all of the static data
are repeated so that each have multiple time steps with the same value. For
example, for the input of the dimensionless number Re, the repetition would
generate $Re(0)=Re(1)=\cdots=Re(t_{max})$.

For the hidden layers, we have a number of LSTM layers, followed by a Dense
layer to generate values within the expected range for predictions and a
Reshape layer to control the output shape. The schematics for the input/output
relations of one of the explored neural network structures for the droplet
collision problem is shown in Fig. 6.

Fig. 6: Machine Learning architecture for transient energy predictions using
transient geometric data for the droplet collision problem.

To better visualize the explored predictions, all test cases are shown in Tab.
4, characterized by output type (static or transient), dataset used, input
features and output features.

Test cases
---
Type |  Dataset |  Input |  Output
Transient |  Spreading |  $t,D(t)$, Re, We |  $E_{k}(t),E_{s}(t),E_{d}(t)$
Transient |  Collision |  $D(t),H(t)$, Re, We, B |  $E_{k}(t),E_{s}(t),E_{d}(t)$
Transient |  Collision |  $D(t),H(t)$ |  $E_{k}(t),E_{s}(t),E_{d}(t)$
Static |  Collision |  $E_{k}(t),E_{s}(t),E_{d}(t)$ |  Re, We
Table 4: Description of the explored test cases.

## III. Results

A supervised LSTM network was employed to train models capable of predicting
the energy budget based on geometric time series data extracted from numerical
simulations. Initially, to test the predictive performance and adaptability of
the architecture, two models were trained by using each of the previously
defined datasets.

### A. Performance of the model formulations

Firstly, as the hyperparameter tuning is a crucial stage for applying machine
learning modeling, in the current work we have employed the open source
hyperparameter optimization framework Optuna [55]. After a detailed
investigation considering the datasets obtained by the simulations of droplets
with varying initial shapes impacting on solid surfaces and collision of
droplets, the hyperparameters were chosen as follows

  * •

Learning rate of 0.0001;

  * •

Maximum of 5000 epochs;

  * •

Batch sizes of 8 for spreading and 32 for collision;

  * •

3 LSTM hidden layers with 32 neurons each;

  * •

1 Dense output layer with $n\times m$ neurons, where $n$ represents the number
of time steps and $m$ represents the number of features;

  * •

1 Reshape layer (transient output only);

  * •

Adam as an optimization algorithm;

  * •

Mean square error applied as the loss function.

Additionally, an early stopping is applied based on training loss (collision)
or validation loss (spreading).

To measure the predictive capability of each model, we used the Coefficient of
determination - $R^{2}$ (26), the Root Mean Square Error - $RMSE$ (27) and the
Normalized Root Mean Square Error - $NRMSE$ (28) defined, respectively, as

$\displaystyle
R^{2}=1-\frac{\sum_{i=1}^{N}(y_{i}-\tilde{y}_{i})^{2}}{\sum_{i=1}^{N}(y_{i}-\bar{y})^{2}},$
(26) $\displaystyle
RMSE=\sqrt{\frac{1}{N}\sum_{i=1}^{N}(y_{i}-\tilde{y}_{i})^{2}},$ (27)
$\displaystyle NRMSE=\frac{RMSE}{y_{max}-y_{min}},$ (28)

where $N$ is the number of samples, $y_{max}$ and $y_{min}$ are, respectively,
the maximum and the minimum possible value for the feature, $y_{i}$ and
$\tilde{y}_{i}$ are, respectively, the expected and the predicted value, and
$\bar{y}$ is the average of all expected values.

The $R^{2}$ is employed to assess how well the variability of the dependent
variables (energies) can be explained by the variability of the independent
variables (geometric data and dimensionless numbers). A score approaching 0
indicates a poor explanatory capability, while a score nearing 1 signifies a
strong explanatory capability for the dependent variable. $RMSE$ and $NRMSE$
are both used to evaluate the quality of the predictions. While RMSE generates
a value dependent on the scale of the predicted variables, NRMSE is used to
provide a relative measure of error regardless of the data scale. In both
cases, a value close to zero indicates a well-fitted model.

While metrics (26) and (27) are used to individually measure the efficiency of
each model, (28) is used to compare models and check the adaptability of the
chosen architecture.

For the spreading dataset, 10 different shapes are split into 70% training (7
shapes, 95 samples), 20% validation (2 shapes, 27 samples) and 10% test (1
shape, 10 samples). To optimize the training process, validation loss is used
as an early stopping control parameter.

The $R^{2}$ score for each individual sample ranges from 0.98315 to 0.9994,
while the $RMSE$ ranges from 0.01122 to 0.05893. Overall, both metrics fall
within the expected optimal range and the majority of the predictions fall
within a $\pm 10\%$ error range. Additionally, we observe that these metrics
remain rather consistent across all samples, providing further evidence that
the model is well-fitted.

For the collision dataset, 72 samples are split into 70% training (50 samples)
and 30% test (22 samples). To optimize the training process, training loss is
used as an early stopping control parameter.

The $R^{2}$ score for each of the 22 test samples ranges from 0.95291 to
0.99997, with the $RMSE$ ranging from 0.00115 to 0.04618. Given the increased
variability in energy behavior, it is understandable how the range of the
metrics is larger than in the spreading case. Another contributing factor is
that there are less samples in the training dataset with higher Reynolds
values, resulting in lower accuracy for predictions within that range.

Despite this difference, both metrics still fall near the optimal ranges, even
presenting higher accuracy with fewer training samples than the spreading
dataset, as will be shown next. This also shows that predicting the behavior
of a new drop shape (spreading case) is more complex than predicting the
behavior of new flow parameters (collision case).

Previously, the metrics were analyzed for each sample individually. To assess
the overall prediction accuracy of the models, Table 5 compares the error
metrics for all training and test samples from both datasets.

Metric | Training samples | Test samples
---|---|---
Spreading | Collision | Spreading | Collision
$R^{2}$ |  0.9924 |  0.99996 |  0.99429 |  0.99559
$RMSE$ |  0.03677 |  0.00124 |  0.02957 |  0.0137
$NRMSE$ |  0.03674 |  0.0019 |  0.02961 |  0.02096
Table 5: Comparison of the evaluation metrics for each dataset.

Since the training datasets are employed in training the model, it is expected
that their metrics should also be optimal. Nevertheless, comparing these
metrics with those of the test datasets provides insights into how well the
model has acquired adaptability to new sets of data that were not used during
training.

The most important metric for comparison between both models is the $NRMSE$,
as it offers a normalized measure of error, disregarding the scale of the
data. With these insights, we can infer that the chosen architecture is not
only reasonably accurate but also highly versatile for different regimes.

### B. Prediction for the energy budget

To better visualize the model’s efficiency for the spreading case, we first
analyse the prediction behavior of two different samples shown in Fig. 7. Each
energy plot is normalized by the total energy in the first time step
($E_{t}(0)$).

In Figs. 7(a) and 7(b), we can observe pronounced fluctuations caused by the
sliding technique (as depicted in Fig. 1) applied to the dataset. Given that
each sample is divided into multiple sub-samples, the predictions can
oscillate between each one of them, as observed in the surface energy $E_{s}$
shown in Fig. 7(a) and the viscous dissipation $E_{d}$ from Fig. 7(b). While
it may appear to increase errors, some tests have shown that the technique, in
fact, enhances model accuracy by generating more samples and using fewer time
steps, which reduces training complexity.

Since most samples from the training and test datasets exhibit a lower range
of surface energy ($E_{s}$) values, as illustrated in Fig. 7(a), their
prediction errors end up being the most visually accentuated.

The predominance of samples within this lower range also influences the
training process. Samples with lower Weber numbers, which have higher surface
tension effects, exhibit higher surface energies compared to the majority. As
depicted in Fig. 7(b), this higher surface energy undermines accuracy, as the
model is primarily tuned for lower values. On the other hand, a higher
Reynolds number can impact the prediction of kinetic energy $E_{k}$, affecting
the capability of accurately predicting the curve in the few steps leading up
to when the drop stops moving, as somewhat depicted in Fig. 7(b).

(a) 1st sample - Re = 30.73, We = 1659.38.

(b) 4th sample - Re = 24.81, We = 90.62.

Fig. 7: Comparison between the expected and predicted energy time series from
two samples of the spreading test dataset.

In general, the main difference in energy behavior observed in the spreading
dataset lies in the steepness of its curves. Since the initial drop shape is
one of the key factors influencing this behavior, the analysis indicates that
the model is capable of generating accurate predictions for samples with a new
drop shape (see Fig. 2).

For the collision dataset, we illustrate in Fig. 8 the prediction behavior of
two different test cases. Again, the energy plot is normalized using the total
energy in the first time step.

The collision adds another layer of complexity to the problem, as the energies
now exhibit oscillations over time. In Fig. 8(a) we observe that the model can
produce accurate predictions even when there is an anomaly in the curves (as
seen at time $t\approx 6$). Furthermore, although fluctuations in viscous
dissipation $E_{d}$ are still detectable, they occur much less frequently.
Decreasing Re and We, as illustrated in Fig. 8(b), the coalesced droplet
experiences an oscillatory regime that can also be correctly captured by the
prediction. Overall, according to these figures, we can confirm the excellent
agreement between the expected and predicted results.

(a) 1st sample - Re = 115.56, We = 53.0, B = 0.17.

(b) 21st sample - Re = 74.45, We = 22.0, B = 0.27.

Fig. 8: Comparison between the expected and predicted energy time series from
three samples of the collision test dataset.

### C. Two-phase prediction with sequential neural networks

To further investigate the significance of our proposed methodology, we adopt
a two-phase prediction strategy tailored specifically for the droplet
collision problem. Initially, we consider only the drop diameter $D(t)$ and
height $H(t)$ as inputs to our model. Through this initial phase, our aim is
to accurately estimate the energy budget associated with the collision
process. Subsequently, leveraging the insights gained from the energy budget
prediction, we extend our model to predict static values, such as the
nondimensional parameters Re and We.

This two-phase prediction approach holds particular relevance in experimental
settings where obtaining precise fluid parameters is challenging. Often,
despite having access to video data from experiments, the lack of
comprehensive knowledge regarding the fluid properties impedes accurate
prediction of the energy budget. In such scenarios, our method offers a robust
solution by including geometric data to not only predict the energy budget but
also estimate essential dimensionless parameters crucial for understanding the
underlying physics of droplet dynamics. Using this innovative two-phase
prediction process, we bridge the gap between experimental observations and
theoretical insights, thereby enhancing our ability to analyze and interpret
droplet collision phenomena with greater precision and efficiency.

In summary, this process involves two separate neural networks. The output of
the first network, which predicts the energy budget, serves as input for
predicting the dimensionless numbers in the second network, as illustrated in
Fig. 9.

Fig. 9: Two-phase schematics for the energies and dimensionless numbers
predictions.

####  C.1 Energy budget prediction

For the first step, we employ an LSTM network to predict all three energies
using only the diameter $D(t)$ and height $H(t)$ of the droplet. Fig. 10
illustrates an example of how the energy prediction behavior is affected when
dimensionless numbers are omitted as inputs.

(a) With dimensionless numbers.

(b) Without dimensionless numbers.

Fig. 10: Comparison between energy predictions with and without the
dimensionless numbers as inputs - Re = 115.56, We = 53.0, B = 0.17.

Although the model experiences a decrease in accuracy when the dimensionless
numbers are removed, the overall behavior of the energies is still correctly
captured. Moreover, given the extremely limited availability of collectable
data from just a video, this increase in error remains tolerable. To assess
the differences in metrics across training and test samples for both
configurations, let us examine Table 6.

Metric | Training samples | Test samples
---|---|---
With | Without | With | Without
$R^{2}$ |  0.99996 |  0.99963 |  0.99604 |  0.9776
$RMSE$ |  0.00124 |  0.00401 |  0.02732 |  0.03087
$NRMSE$ |  0.0019 |  0.00615 |  0.02736 |  0.04725
Table 6: Comparison of the evaluation metrics for models with and without
dimensionless numbers Re, We and B as inputs.

In general, if the dimensionless numbers are known, it is advisable to choose
the model that has been trained with them included. However, upon comparing
the metrics, it is evident that utilizing only geometric data for the model
remains a viable option with decent accuracy.

####  C.2 Dimensionless numbers prediction

In a parallel test, we attempted to predict the dimensionless numbers directly
using only geometric data. However, this approach resulted in poor prediction
accuracy for the collision dataset, as shown in Appendix B.

Therefore, rather than directly using $D(t)$ and $H(t)$ to predict the
dimensionless numbers Re and We, the second step involves utilizing the
previously trained neural network to predict the energies and then training a
new neural network to predict the dimensionless numbers using these energies.

Given that each sample is divided into 100-time-step slices, the selected
prediction for a sample is determined by the median of the predicted values
for each of its 10 slices. This is done to mitigate potential noise caused by
high errors resulting from using approximated energies as inputs. The
prediction results for the collision test dataset are depicted in Fig. 11.

Fig. 11: Dimensionless numbers predictions for each of the 22 samples from the
droplet collision test dataset.

Despite the introduction of some noise through the use of approximated
energies, the overall predicted values still appear to align with the
expected. To quantitatively analyze the results, the previous metrics are
calculated both for training and test data in Table 7.

Metric |  Training samples |  Test samples
---|---|---
$R^{2}$ |  0.99366 |  0.99823
$RMSE$ |  2.50602 |  1.39455
$NRMSE$ |  0.02534 |  0.0141
Table 7: Comparison of the evaluation metrics for the prediction of the
dimensionless numbers by using predicted energies.

With this, we have demonstrated that, in certain regimes, it may be feasible
to predict energies and dimensionless numbers using only geometric data
extracted from experiment videos. As a follow-up to this study, new
architectures could be developed and tested with additional datasets to
enhance the neural network generalization, thereby expanding its applicability
to a wider range of regimes.

## IV. Conclusions

Initially, we proposed two LSTM models to predict the energy balance in two
moving free surface problems: the impact of non-spherical droplet in a rigid
surface and the collision of two equal-sized droplets. Since the first dataset
required generating data from different and unusual shapes, a method capable
of extracting particles from interfaces present in digital images was also
developed (as shown in Appendix A). Following this test, a two-phase
sequential prediction approach was applied in an attempt to predict important
properties using exclusively geometric data. This method is specially useful
for experimental observations, since it only uses data which can be easily
extracted from experiment images or videos. It provides an easy-to-use tool
for the challenging task of predicting fluid parameters such as the energies
and the dimensionless numbers.

In the context of one-phase prediction, both models shared the same overall
architecture, with the only variation being in the data used for training. The
first model was trained using data from 95 numerical simulations involving
seven different shapes, extracted from images of real experiments, impacting a
solid surface and spreading over time. The second model utilized data from 50
numerical simulations of two colliding drops that spread and retracted after
coalescence, causing oscillations in the energy budget. For the droplet impact
problem, the inputs to the neural network consist of simulation time, the
temporal series of the droplet diameter and nondimensional numbers. However,
for the droplet collision scenario which is more complex than the spreading
case, the model is initialized with the temporal series of diameter and height
of the coalesced droplet, in combination with nondimensional numbers and the
impact parameter. In both models, the output corresponds to the energy budget.

To evaluate the performance of each model, we conducted a comprehensive metric
analysis commonly employed in machine learning modeling, which includes the
Coefficient of Determination ($R^{2}$), the Root Mean Square Error ($RMSE$),
and the Normalized Root Mean Square Error ($NRMSE$). The $R^{2}$, $RMSE$ and
$NRMSE$ metrics for both models are, respectively, 0.99429 and 0.99559,
0.02957 and 0.0137, and 0.02961 and 0.02096. After analyzing the metrics, we
can observe that both models achieved performances sufficiently close to the
ideal values for each metric, with $R^{2}s$ close to one and $RMSEs$
approaching zero. The analyzed models are accurate in a similar way, even when
applied to considerably different problems, which can be evaluated by the
similarity of their $NRMSEs$. The spreading model showed promising results by
successfully handling a completely new shape. Similarly, the collision
dataset’s model performed well with samples featuring sets of parameters not
used during training. The accuracy of the models is also evaluated through
plots depicting the kinetic, surface, and viscous dissipation, confirming
excellent agreement between the simulation data and the predicted values. This
suggests that the machine learning architectures have the potential to be
highly adaptable and accurate across a range of fluid regimes.

According to the results for a two-phase prediction, the architecture remains
similar, but now follows a two-step scheme. Firstly, an LSTM network predicts
the energy budget using only geometric data. Then, the predicted energies
serve as inputs in a second LSTM network to predict the static dimensionless
numbers. The metrics for the first phase (energies) and the second phase
(static prediction for Reynolds and Weber numbers) are, respectively, 0.9776
and 0.99823, 0.03087 and 1.39455, and 0.04725 and 0.0141. Despite a slight
increase in error when compared to the one-phase prediction strategy, the
metrics for two-phase predictions still demonstrate high accuracy. After the
performance investigation concerning the metrics, we have also included the
results for predicting the Reynolds and Weber numbers. Results show excellent
agreement for low values of Re and We, whereas for larger values, the neural
network appears to have more difficulty in establishing accurate predictions.
This difficulty can be attributed to the fact that, for the training data,
there are fewer instances with higher Re and We values, as seen in Table 3.
Finally, this approach proves useful in the sense that directly attempting to
predict the dimensionless numbers yielded higher errors, as demonstrated in
Appendix B.

While the metrics indicate high versatility, there are still some questions
which must be explored. One notable task is the challenge of predicting
regimes beyond coalescence in the collision model, given that our current
model predominantly relies on coalescence data. For instance, capturing the
distinct behaviors associated with the bouncing regime poses a significant
challenge, as transitioning between regimes presents a complex task for
machine learning models. Future research could be focused on finding new
appropriate fluid parameters specific to the bouncing phenomena (or any other
regime, for that matter) and exploring interpretation and generalization of
the model.

We hope that this research will shed light on the advantages of leveraging the
performance offered by LSTM architectures for predicting different types of
transient data in complex flows of complex fluids, as for instance dynamics of
viscoelastic droplets [56, 57].

## Acknowledgements

The authors would also like to thank the financial support given by the São
Paulo Research Foundation (FAPESP) grants #2013/07375-0 and #2021/14953-6, the
National Council for Scientific and Technological Development (CNPq), grants
#305383/2019-1 and CAPES. The authors also acknowledge the Numerical
Simulation Laboratory at FCT/UNESP for their support with cluster resources.
C.M. Oishi would like to thank the insightful discussions on machine learning
modeling held with the research team during his visit at the AI Institute in
Dynamic Systems, University of Washington.

##  Appendix A Drop shape extraction from an image

One of the limitations observed in most flow solvers is that defining
interfaces with unusual formats is a difficult process. In the experimental
environment, creating spherical or elliptical drops that are strictly
symmetrical is a practically impossible task. However, when running numerical
simulations to reproduce such experiments, it is often chosen to use these
formats because they present parametric equations that easily define them.

In order to be able to define drops with shapes similar to those observed in
experiments, thus allowing the testing of network architectures without
introducing noise due to shape approximations, a method capable of extracting
interfaces present in digital images was also developed in this work. To
extract particles from an image, the steps shown on Fig. 12 are performed:

Fig. 12: Steps to extract particles from an image.

Step 1:

open an image that contains the interface. Preferably, the interface should be
outlined with another color to help the contour detection process;

Step 2:

change the image to grayscale or directly to black and white (recommended only
for cases that feature easy-to-extract interfaces with little to no noise);

Step 3:

for images converted to grayscale, apply a threshold. Based on a threshold
value, image pixels are converted to white or black. If a pixel has a value
lower than the threshold (L), it is changed to black (0). Otherwise, it is
changed to white (255). For the example in Fig. 12, a threshold of $L=110$ was
applied;

Step 4:

extract the contour(s) of the image that represents the interface(s). In the
example previously shown, there is only one continuous interface, so only one
contour was detected;

Step 5:

approximate the contour shape by collecting a set of discretized points
$(i,j)$. Since the points collected by the function are actually the indices
of the pixels from the image’s contour, they need to be adjusted so that they
fall in the ranges $[x_{min},x_{max}]$ and $[y_{min},y_{max}]$:

$\displaystyle x=\frac{j(x_{min}-x_{max})}{m}+x_{max}$ $\displaystyle
y=\frac{i(y_{min}-y_{max})}{n}+y_{max}$

With $n\times m$ being the dimensions of the input image and $(x,y)$ being the
new coordinates for the points extracted from the interface. As the list of
pixels returned is ordered, the method facilitates the creation of particles
for representation by front-tracking.

Step 6:

If the function returns an insufficient number of points or generates points
with large spacing (as seen in the previous example), it is possible to apply
linear interpolation to create new intermediate points.

To do this, the Euclidean distance between each pair of points is calculated
and, if the distance is greater than a chosen maximum tolerance, new
intermediate points are generated.

As the method involves extracting points through the pixels of an image’s
contours, the main restriction observed is that the extraction is limited by
the image’s resolution. The more pixels an image has, the more accurate the
extracted interface will be.

Furthermore, as the points are extracted from the lower left corner of each
pixel, aliasing behavior is generated in the extracted points, as can be seen
in Fig. 13(a).

This disturbance in the initialization of points can generate atypical
behavior in the earlier steps of the simulation. To help solve the aliasing
problem, a technique was applied directly to the solver.

The Trapezoidal Sub-grid Undulations Removal (TSUR), developed by [58], is a
mass conserving smoothing method used to reduce aliasing and abrupt
undulations in curves. The result of applying the method is shown in Fig.
13(b).

(a) Aliasing generated when extracting points from pixels.

(b) The new points after applying TSUR to reduce aliasing.

Fig. 13: Comparison before and after applying the TSUR technique.

In addition to the usefulness of TSUR in removing the initial aliasing of the
interface, the technique can also reduce numerical noise that may appear
during the simulation, generating more physically accurate approximations.

##  Appendix B Static predictions for the dimensionless numbers using
geometric data

To assess the feasibility of predicting dimensionless numbers using only
geometric data, we trained an LSTM network using the collision training
dataset. Fig. 14 illustrates the predicted values for each instance of the
collision test dataset while a detailed metric analysis is presented in Table
8.

Fig. 14: Direct dimensionless numbers predictions for each of the 22 samples
from the collision test dataset.

According to the results in Figs. 11 and 14 and analysing the difference in
metrics from Table 8, it is evident that the performance of the one-phase
neural network, operating solely on geometrical data as inputs, is inferior to
that of our two-stage prediction strategy presented in Section C.2. This
comparison highlights the crucial role that is undertaken by the two-phase
sequential neural network approach in enhancing the accuracy and robustness of
dimensionless number predictions.

Metric |  Two-phase |  Direct
---|---|---
$R^{2}$ |  0.99823 |  0.95478
$RMSE$ |  1.39455 |  7.04277
$NRMSE$ |  0.0141 |  0.07122
Table 8: Comparison of the evaluation metrics for the two-phase and direct
approaches.

## References

  * [1] Alexander L. Yarin  “DROP IMPACT DYNAMICS: Splashing, Spreading, Receding, Bouncing…”  In _Annual Review of Fluid Mechanics_ 38 Annual Reviews, 2006, pp. 159–192  DOI: 10.1146/annurev.fluid.38.050304.092144
  * [2] Howard A. Stone  “Dynamics of Drop Deformation and Breakup in Viscous Fluids”  In _Annual Review of Fluid Mechanics_ 26 Annual Reviews, 1994, pp. 65–102  DOI: 10.1146/annurev.fl.26.010194.000433
  * [3] Christophe Josserand and Sigurdur T. Thoroddsen  “Drop Impact on a Solid Surface”  In _Annual Review of Fluid Mechanics_ 48 Annual Reviews, 2016, pp. 365–391  DOI: 10.1146/annurev-fluid-122414-034401
  * [4] Shiji Lin et al.  “Impact of viscous droplets on different wettable surfaces: Impact phenomena, the maximum spreading factor, spreading time and post-impact oscillation”  In _Journal of Colloid and Interface Science_ 516 Elsevier, 2018, pp. 86–97  DOI: https://doi.org/10.1016/j.jcis.2017.12.086
  * [5] Detlef Lohse  “Fundamental Fluid Dynamics Challenges in Inkjet Printing”  In _Annual Review of Fluid Mechanics_ 54 Annual Review, 2022, pp. 349–382  DOI: 10.1146/annurev-fluid-022321-114001
  * [6] Robert D. Schroll, Christophe Josserand, Stéphane Zaleski and Wendy W. Zhang  “Impact of a Viscous Liquid Drop”  In _Phys. Rev. Lett._ 104 American Physical Society, 2010, pp. 034504  DOI: 10.1103/PhysRevLett.104.034504
  * [7] Jens Eggers, Marco A. Fontelos, Christophe Josserand and Stéphane Zaleski  “Drop dynamics after impact on a solid wall: theory and simulations”  In _Physics of Fluids_ 22 AIP Publishing, 2010, pp. 062101  DOI: 10.1063/1.3432498
  * [8] Taehun Lee and Lin Liu  “Lattice Boltzmann simulations of micron-scale drop impact on dry surfaces”  In _Journal of Computational Physics_ 229 Elsevier, 2010, pp. 8045–8063  DOI: https://doi.org/10.1016/j.jcp.2010.07.007
  * [9] Siddhartha F. Lunkad, Vivek V. Buwa and K.D.P. Nigam  “Numerical simulations of drop impact and spreading on horizontal and inclined surfaces”  In _Chemical Engineering Science_ 62 Elsevier, 2007, pp. 7214–7224  DOI: https://doi.org/10.1016/j.ces.2007.07.036
  * [10] Yu Pan and Kazuhiko Suga  “Numerical simulation of binary liquid droplet collision”  In _Physics of Fluids_ 17 AIP Publishing, 2005, pp. 082105  DOI: 10.1063/1.2009527
  * [11] N. Nikolopoulos, K.-S. Nikas and G. Bergeles  “A numerical investigation of central binary collision of droplets”  In _Computers & Fluids_ 38 Elsevier, 2009, pp. 1191–1202  DOI: https://doi.org/10.1016/j.compfluid.2008.11.007
  * [12] Ilia V. Roisman  “Inertia dominated drop collisions. II. An analytical solution of the Navier–Stokes equations for a spreading viscous film”  In _Physics of Fluids_ 21 AIP Publishing, 2009, pp. 052104  DOI: 10.1063/1.3129283
  * [13] Stéphane Popinet  “An accurate adaptive solver for surface-tension-driven interfacial flows”  In _Journal of Computational Physics_ 228 Elsevier, 2009, pp. 5838–5866  DOI: https://doi.org/10.1016/j.jcp.2009.04.042
  * [14] Stéphane Popinet  “Numerical Models of Surface Tension”  In _Annual Review of Fluid Mechanics_ 50 Annual Review, 2018, pp. 49–75  DOI: 10.1146/annurev-fluid-122316-045034
  * [15] S. Chandra and C. T. Avedisian  “On the collision of a droplet with a solid surface”  In _Proceedings of the Royal Society of London. Series A: Mathematical and Physical Sciences_ 432 The Royal Society, 1991, pp. 13–41  DOI: 10.1098/rspa.1991.0002
  * [16] Nick Laan et al.  “Maximum Diameter of Impacting Liquid Droplets”  In _Phys. Rev. Appl._ 2 American Physical Society, 2014, pp. 044018  DOI: 10.1103/PhysRevApplied.2.044018
  * [17] J. Fukai et al.  “Wetting effects on the spreading of a liquid droplet colliding with a flat surface: Experiment and modeling”  In _Physics of Fluids_ 7 AIP Publishing, 1995, pp. 236–247  DOI: 10.1063/1.868622
  * [18] Jae Bong Lee, Dominique Derome, Robert Guyer and Jan Carmeliet  “Modeling the Maximum Spreading of Liquid Droplets Impacting Wetting and Nonwetting Surfaces”  In _Langmuir_ 32 ACS Publications, 2016, pp. 1299–1308  DOI: 10.1021/acs.langmuir.5b04557
  * [19] Chijioke Ukiwe and Daniel Y. Kwok  “On the Maximum Spreading Diameter of Impacting Droplets on Well-Prepared Solid Surfaces”  In _Langmuir_ 21 ACS Publications, 2005, pp. 666–673  DOI: 10.1021/la0481288
  * [20] Kuo-Long Pan, Chung K. Law and Biao Zhou  “Experimental and mechanistic description of merging and bouncing in head-on binary droplet collision”  In _Journal of Applied Physics_ 103 AIP Publishing, 2008, pp. 064901  DOI: 10.1063/1.2841055
  * [21] J.-P. Estrade, Hervé Carentz, G. Lavergne and Y. Biscos  “Experimental investigation of dynamic binary collision of ethanol droplets – a model for droplet coalescence and bouncing”  In _International Journal of Heat and Fluid Flow_ 20 Elsevier, 1999, pp. 486–491  DOI: https://doi.org/10.1016/S0142-727X(99)00036-3
  * [22] J. Qian and C. K. Law  “Regimes of coalescence and separation in droplet collision”  In _Journal of Fluid Mechanics_ 331 Cambridge University Press, 1997, pp. 59–80  DOI: 10.1017/S0022112096003722
  * [23] J. Qian et al.  “An experimental and computational study of bouncing and deformation in droplet collision”  In _35th Aerospace Sciences Meeting and Exhibit_ American Institute of AeronauticsAstronautics, 1997  DOI: 10.2514/6.1997-129
  * [24] K. D. Willis and M. E. Orme  “Experiments on the dynamics of droplet collisions in a vacuum”  In _Experiments in Fluids_ 29 Springer Nature, 2000, pp. 347–358  DOI: 10.1007/s003489900092
  * [25] Sander Wildeman, Claas Willem Visser, Chao Sun and Detlef Lohse  “On the spreading of impacting drops”  In _Journal of Fluid Mechanics_ 805 Cambridge University Press, 2016, pp. 636–655  DOI: 10.1017/jfm.2016.584
  * [26] Šikalo et al.  “Dynamic contact angle of spreading droplets: Experiments and simulations”  In _Physics of Fluids_ 17 AIP Publishing, 2005, pp. 062103  DOI: 10.1063/1.1928828
  * [27] P. Attané, F. Girard and V. Morin  “An energy balance approach of the dynamics of drop impact on a solid surface”  In _Physics of Fluids_ 19 AIP Publishing, 2007, pp. 012101  DOI: 10.1063/1.2408495
  * [28] Song Yang, Yu Hou, Yuheng Shang and Xin Zhong  “BPNN and CNN-based AI modeling of spreading and icing pattern of a water droplet impact on a supercooled surface”  In _AIP Advances_ 12 AIP Publishing, 2022, pp. 045209  DOI: 10.1063/5.0082568
  * [29] Moussa Tembely, Damien C. Vadillo, Ali Dolatabadi and Arthur Soucemarianadin  “A Machine Learning Approach for Predicting the Maximum Spreading Factor of Droplets upon Impact on Surfaces with Various Wettabilities”  In _Processes_ 10 MDPI, 2022, pp. 1141  DOI: 10.3390/pr10061141
  * [30] Jiguo Tang et al.  “Universal Model for Predicting Maximum Spreading of Drop Impact on a Smooth Surface Developed Using Boosting Machine Learning Models”  In _Industrial & Engineering Chemistry Research_ 62 ACS Publications, 2023, pp. 15268–15277  DOI: 10.1021/acs.iecr.3c01560
  * [31] Lap Au-Yeung and Peichun Amy Tsai  “Predicting Impact Outcomes and Maximum Spreading of Drop Impact on Heated Nanostructures Using Machine Learning”  In _Langmuir_ 39 ACS Publications, 2023, pp. 18327–18341  DOI: 10.1021/acs.langmuir.3c02405
  * [32] Elham Heidari, Abolghasem Daeichian, Mohammad Amin Sobati and Salman Movahedirad  “Prediction of the droplet spreading dynamics on a solid substrate at irregular sampling intervals: Nonlinear Auto-Regressive eXogenous Artificial Neural Network approach (NARX-ANN)”  In _Chemical Engineering Research and Design_ 156 Elsevier, 2020, pp. 263–272  DOI: https://doi.org/10.1016/j.cherd.2020.01.033
  * [33] Raihan Choudhury et al.  “Uncovering the Effect of Physical Conditions and Surface Roughness on the Maximum Spreading Factor of Impinging Droplets Using a Supervised Artificial Neural Network Model”  In _Industrial & Engineering Chemistry Research_ 62 ACS Publications, 2023, pp. 21208–21221  DOI: 10.1021/acs.iecr.3c02650
  * [34] Jingzu Yee, Daichi Igarashi, Shun Miyatake and Yoshiyuki Tagawa  “Prediction of the morphological evolution of a splashing drop using an encoder–decoder”  In _Machine Learning: Science and Technology_ 4 IOP Publishing, 2023, pp. 025002  DOI: 10.1088/2632-2153/acc727
  * [35] Jingzu Yee, Akinori Yamanaka and Yoshiyuki Tagawa  “Image features of a splashing drop on a solid surface extracted using a feedforward neural network”  In _Physics of Fluids_ 34 AIP Publishing, 2022, pp. 013317  DOI: 10.1063/5.0077050
  * [36] Jingzu Yee et al.  “Correlation between morphological evolution of splashing drop and exerted impact force revealed by interpretation of explainable artificial intelligence”, 2023  arXiv: https://doi.org/10.48550/arXiv.2309.10266
  * [37] Hugo L. França and Cassio M. Oishi  “A machine learning strategy for computing interface curvature in Front-Tracking methods”  In _Journal of Computational Physics_ 450 Elsevier, 2022, pp. 110860  DOI: https://doi.org/10.1016/j.jcp.2021.110860
  * [38] Aniruddha Gaikwad et al.  “In-process monitoring and prediction of droplet quality in droplet-on-demand liquid metal jetting additive manufacturing using machine learning”  In _Journal of Intelligent Manufacturing_ 33 Springer Nature, 2022, pp. 2093–2117  DOI: 10.1007/s10845-022-01977-2
  * [39] M. Pasandideh‐Fard, Y. M. Qiao, S. Chandra and J. Mostaghimi  “Capillary effects during droplet impact on a solid surface”  In _Physics of Fluids_ 8 AIP Publishing, 1996, pp. 650–659  DOI: 10.1063/1.868850
  * [40] Nicolas G. Hadjiconstantinou  “Estimating the Maximum Splat Diameter of a Solidifying Droplet”  In _ASME International Mechanical Engineering Congress and Exposition_ Heat Transfer: Volume 1 American Society of Mechanical Engineers, 1999, pp. 311–316  DOI: 10.1115/IMECE1999-0998
  * [41] Jae Bong Lee, Dominique Derome, Ali Dolatabadi and Jan Carmeliet  “Energy Budget of Liquid Drop Impact at Maximum Spreading: Numerical Simulations and Experiments”  In _Langmuir_ 32 ACS Publications, 2016, pp. 1279–1288  DOI: 10.1021/acs.langmuir.5b03848
  * [42] Lennon Ó Náraigh and Juan Mairal  “Analysis of the spreading radius in droplet impact: The two-dimensional case”  In _Physics of Fluids_ 35 AIP Publishing, 2023, pp. 102124  DOI: 10.1063/5.0169274
  * [43] Qingzhe Liu et al.  “The role of drop shape in impact and splash”  In _Nature communications_ 12 Springer Nature, 2021, pp. 3068  DOI: 10.1038/s41467-021-23138-4
  * [44] Zhenyu Zhang and Peng Zhang  “Kinetic energy recovery and interface hysteresis of bouncing droplets after inelastic head-on collision”  In _Physics of Fluids_ 29 AIP Publishing, 2017, pp. 103306  DOI: 10.1063/1.5000547
  * [45] Karrar H. Al-Dirawi and Andrew E. Bayly  “A new model for the bouncing regime boundary in binary droplet collisions”  In _Physics of Fluids_ 31 AIP Publishing, 2019, pp. 027105  DOI: 10.1063/1.5085762
  * [46] D.V. Antonov et al.  “Energy analysis of secondary droplet atomization schemes”  In _International Communications in Heat and Mass Transfer_ 117 Elsevier, 2020, pp. 104666  DOI: https://doi.org/10.1016/j.icheatmasstransfer.2020.104666
  * [47] Antoine Riaud et al.  “Hydrodynamic constraints on the energy efficiency of droplet electricity generators”  In _Microsystems & Nanoengineering_ 7 Springer Nature, 2021, pp. 49  DOI: 10.1038/s41378-021-00269-8
  * [48] Olinka Ramírez-Soto et al.  “Lifting a sessile oil drop from a superamphiphobic surface with an impacting one”  In _Science Advances_ 34 Science Advances, 2020, pp. eaba4330  DOI: 10.1126/sciadv.aba4330
  * [49] Sepp Hochreiter and Jürgen Schmidhuber  “Long Short-Term Memory”  In _Neural Computation_ 9 MIT Press, 1997, pp. 1735–1780  DOI: 10.1162/neco.1997.9.8.1735
  * [50] A. J. Chorin and J. E. Marsden  “A mathematical introduction to fluid mechanics”  Springer-Verlag, 1979
  * [51] Cassio M. Oishi, Roney L. Thompson and Fernando P. Martins  “Transient motions of elasto-viscoplastic thixotropic materials subjected to an imposed stress field and to stress-based free-surface boundary conditions”  In _International Journal of Engineering Science_ 109 Elsevier, 2016, pp. 165–201  DOI: https://doi.org/10.1016/j.ijengsci.2016.08.004
  * [52] L. C. Jain and L. R. Medsker  In _Recurrent Neural Networks: Design and Applications_ CRC Press, Inc., 1999
  * [53] Seol-Hyun Noh  “Analysis of Gradient Vanishing of RNNs and Performance Comparison”  In _Information_ 12 MDPI, 2021, pp. 11  DOI: 10.3390/info12110442
  * [54] Ralf C. Staudemeyer and Eric Rothstein Morris  “Understanding LSTM – a tutorial into Long Short-Term Memory Recurrent Neural Network”, 2019  arXiv: https://doi.org/10.48550/arXiv.1909.09586
  * [55] “Optuna - A hyperparameter optimization framework” Accessed: 2024-03-07, https://optuna.org/
  * [56] C. M. Oishi, R. L. Thompson and F. P. Martins  “Impact of capillary drops of complex fluids on a solid surface”  In _Physics of Fluids_ 31 AIP Publishing, 2019, pp. 123109  DOI: 10.1063/1.5129640
  * [57] Hugo L. França and Cassio M. Oishi  “Numerical investigation of shear-thinning and viscoelastic binary droplet collision”  In _Journal of Non-Newtonian Fluid Mechanics_ 302 Elsevier, 2022, pp. 104750  DOI: https://doi.org/10.1016/j.jnnfm.2022.104750
  * [58] L.G. Nonato et al.  “A mass-conserving smooth method”  In _Mecánica Computacional_ 23 Asociación Argentina de Mecánica Computacional, 2004, pp. 1897–1910