Hybrid AI-Analytical Modeling of Droplet Dynamics on Inclined Heterogeneous Surfaces

1. Introduction

The motion of droplets is a phenomenon that most people encounter in their everyday lives, either in the form of vapor condensation on a cold window, water spilling on leaves during the watering of plants, etc. It is also a phenomenon that is relevant to many applications; from microfluidic lab-on-a-chip devices that can process a large number of microscopic chemical or biological samples [1], to coating and 3D printing [2]. For these reasons, there is a great need to understand and control droplet motion in these and many other applications [see 3, for a review].

From a modeling perspective, the most intricate aspect of a moving droplet is the description of the multiscale physical processes affecting the contact line dynamics, which originate from the molecular interactions at the nanoscale. The contact line is defined as the line on the periphery of a droplet, separating the solid, liquid and gas phases. The associated local contact angle is the angle formed between the solid-liquid and the liquid-gas interfaces and is typically dictated by the chemistry of the surface and may exhibit spatial variations and features due to unavoidable randomness on a natural surface or deterministic by surface design. A solid surface is considered hydrophilic when a fully equilibrated droplet assumes angles that are smaller than ${90}^{\circ }$, and hydrophobic when these are larger than ${90}^{\circ }$. When the droplet is in motion, what we observe macroscopically is that the apparent contact angle is different from the local one, causing the displacement of the contact line. In surfaces with no hysteresis effects, the droplet reaches equilibrium when both the local and apparent contact angles along the contact line match.

The focus of this study is the modeling of droplet motion on inclined, chemically heterogeneous surfaces. This specific setting is very relevant to water harvesting, with numerous examples in nature and relevant bio-inspired applications [4]. Droplets are typically observed to remain stationary on the inclined surface up to a critical inclination angle, beyond which the droplet moves downhill as it is able to overcome the surface heterogeneities be they natural or artificial. For larger inclination angles, instabilities may develop, leading to the formation of cusps and even droplet breakup [5,6]. The impact of structured surface heterogeneities has also attracted interest in recent years, especially with the advent of specialized surface fabrication technologies. Several experimental studies focused on uncovering the influence of heterogeneous features on the critical inclination angle and the subsequent dynamic sliding behaviour [7,8,9,10,11,12]. Going even further, Varagnolo et al. [13] studied the influence of the shape of the hydrophobic regions, suggesting practical criteria for surface designing in order to tune the static and dynamic behavior of droplets. In more “exotic” applications, researchers were able to make droplets run uphill either with a heterogeneity gradient opposing the gravity field [14], or with vibrations [15]. Moreover, approaches based on mathematical modelling and analyses were also utilized to study idealized geometries, typically adopting simplifying assumptions [16,17] which limited their applicability in more realistic scenarios.

Thus, numerical modelling emerges as an attractive alternative approach to go past the limitations of analytical results and mitigate the need for specialized equipment for accurate experiments. Numerical studies, however, also adopt different approaches. For example, lubrication-type models that are derived as approximations of the Navier–Stokes equations and apply for small contact angles and negligible inertia, have been employed to study droplets on inclined surfaces, with or without chemical heterogeneities [18,19,20,21]. Moreover, Afkhami et al. [22] conducted direct numerical simulations (DNS) to study forced dewetting, comparing their results against Cox’s contact line velocity model [23]. Relevant computational fluid dynamics (CFD) studies emerged, studying various wetting hydrodynamic settings for realistic apparent contact angles [24,25,26,27,28]. Other numerical approaches include Lattice–Boltzmann methods [29,30,31], phase-field methods [32,33], and to a lesser extend molecular dynamics [34,35], and smooth particle hydrodynamics [36]. Even though these numerical approaches offer a reliable alternative to costly experiments, they are typically associated with a high computational cost and significant runtime.

The astounding uptake of data-driven methods in science and engineering has showcased the potential to accelerate numerical solutions in different disciplines. Specifically for fluid dynamics [37], common examples of enhancing the solution efficiency through data-driven methods include turbulence modelling through Reynolds averaged Navier–Stokes and large eddy simulations [38,39] and super-resolution, where coarse simulations are processed by data-driven models to infer the solution to higher resolution grids [40,41,42]. Another promising approach is the construction of hybrid surrogate models, which combine data-driven and reduced-order models. Characteristic examples of this approach are two works that developed surrogate models for the prediction of the kinematics of spherical particles [43] and bubbles [44], as well as the recent study by the present authors for the dynamics of thin droplets on heterogeneous surfaces without gravity [45]. With this kind of hybrid modelling, the data-driven component of the model learns to correct and extend the applicability of an analytically derived reduced-order model using the full simulation data as a reference solution. This approach demonstrated improved agreement with the reference solution compared to the base reduced-order model, expanding its regime of applicability to tackle more complex scenarios, and also exhibited better generalizability compared to fully data-driven models.

Specifically, the aforementioned study by Demou and Savva was a proof-of-concept study in which a hybrid analytical/data-driven modelling approach was utilized for the first time in the context of wetting hydrodynamics, which was based on training datasets that were efficiently generated for lubrication-type models [45]. In the limit of thin droplets, the asymptotic theory of Lacey [46] applies, which provides a reasonably accurate estimate for the normal component of the contact line velocity. Although higher-order corrections to this result have been obtained via analytical means [47], this study demonstrated that the data-driven route is also applicable, especially for cases where improving upon a simplified low-order model becomes a formidable challenge. Within this approach, the hybrid model for the contact line velocity consists of the superposition of Lacey’s model and its data-driven counterpart, which was then evolved using the method of lines with a standard ordinary differential equation (ODE) solver. The model was shown to considerably shorten simulation times compared to a full numerical simulation, and to be able to capture the dynamics accurately even for heterogeneity profiles that were very different from the heterogeneity profiles seen during training.

The present work expands upon this modelling framework to describe the dynamics of droplets that are driven by the interplay of chemical heterogeneities and the gravitational force, using (i) full DNS datasets and (ii) a reduced-order model which is based on the theory of Cox [23] which does not explicitly account for the effects of gravity, so that they are entirely accounted for by the data-driven model. In Section 2, the mathematical and numerical framework used to generate the datasets and train the data-driven model is described. This is followed by the presentation of the results of this study, including the testing of the hybrid surrogate model in different settings in Section 3. Finally, Section 4 concludes this study with a summary of the key results.

2. Governing Equations and Numerical Methods

The goal of this study is to train data-driven models to predict the droplet dynamics for a system that is schematically shown in Figure 1. More specifically, a liquid droplet with density ${\widehat{\rho }}_1$ and viscosity ${\widehat{\mu }}_1$, which is surrounded by a gas phase with corresponding properties ${\widehat{\rho }}_2$ and ${\widehat{\mu }}_2$ (hats denote dimensional quantities), slides under the influence of gravitational acceleration $\widehat{g}$ over a chemically heterogeneous surface with inclination angle ${\alpha }_i$. The surface tension $\widehat{\sigma }$ of the two-fluid interface, and the chemical character of the substrate, as expressed through the local contact angle ${\theta }_{*}$, influence the dynamic behaviour of the droplet and its shape. Considering the case set-up, the relevant dimensionless groups are the Reynolds number $\mathrm{Re}={\widehat{U}}_{ref}{\widehat{L}}_{ref}{\widehat{\rho }}_1/{\widehat{\mu }}_1$, characterizing the relative influence of inertial over viscous effects, the capillary number $\mathrm{Ca}={\widehat{\mu }}_1{\widehat{U}}_{ref}/\widehat{\sigma }$, characterizing the relative influence of viscous over capillary effects, and the Bond number, $\mathrm{Bo}={\widehat{\rho }}_1\widehat{g}{\widehat{L}}_{ref}^2/\widehat{\sigma }$, for capturing the relative importance of gravity to surface tension. In these dimensionless groups, ${\widehat{L}}_{ref}$ is the characteristic length scale taken to be the domain size, and ${\widehat{U}}_{ref}={(\widehat{\sigma }/\left({\widehat{\rho }}_1{\widehat{L}}_{ref}\right))}^{1/2}$ is the characteristic velocity scale.

To generate the training dataset, a number of DNS cases are run and the evolution of the contact line is extracted under different conditions. In all cases, an initially hemispherical droplet of radius 0.125 was placed on a solid surface at the bottom of a $1\times 1\times 1$ box, with open boundaries everywhere else. The density and viscosity ratios between the two fluids were set to 0.1, and Reynolds and Capillary numbers to $\mathrm{Re}=0.1$, $\mathrm{Ca}=10$. The surface inclination angle and the Bond number were varied within ${\alpha }_i\in [0,{80}^{\circ }]$, and $\mathrm{Bo}\in [0,5]$ respectively. Furthermore, the initial position of the droplet and the heterogeneity profile on the inclined substrate were also varied, with more details provided in Section 3. The following paragraphs provide a description of the mathematical and numerical framework used to generate the datasets and train the data-driven models.

2.1. Governing Equations for DNS

The governing equations for the wetting scenarios considered are the two-phase Navier–Stokes equations, including the effects of surface tension and gravity. The two different fluids (resembling a liquid and a gas) are differentiated with a colour function $C(\mathit{x},t)$ [48]. Considering that volume ${\Omega }_1\left(t\right)$ is occupied by fluid (1) and ${\Omega }_2\left(t\right)$ by fluid (2), with boundaries $\partial {\Omega }_1\left(t\right)$ and $\partial {\Omega }_2\left(t\right)$, respectively, the colour function is defined as

$$C(\mathit{x},t)=\left\{\begin{array}{c}1\mathrm{if}\mathit{x}\in {\Omega }_1\left(t\right),\hfill \\ 0\mathrm{if}\mathit{x}\in {\Omega }_2\left(t\right),\hfill \end{array}\right.$$

(1)

with the interface between the two fluids denoted by $S\left(t\right)=\partial {\Omega }_1\left(t\right)\bigcap \partial {\Omega }_2\left(t\right)$, which is the part of the boundary that is common to both regions. The physical fluid properties in the whole domain can be obtained as a weighted average between the fluid properties of each phase,

$$\psi (\mathit{x},t)=C(\mathit{x},t){\psi }_1+(1-C(\mathit{x},t)){\psi }_2,$$

(2)

where $\psi$ is the value of a relevant property in the entire domain and subscripts “1” and “2” indicate the property values in each phase.

Within this setting, the governing equations are recast in dimensionless form as,

$$\begin{array}{cc}\hfill \nabla ·\mathit{u}=& 0,\hfill \end{array}$$

(3a)

$$\begin{array}{cc}\hfill \frac{\partial C}{\partial t}+\nabla ·\left(C\mathit{u}\right)=& 0,\hfill \\ \hfill \frac{\partial \rho \mathit{u}}{\partial t}+\nabla ·\left(\rho \mathit{u}\mathit{u}\right)=& -\nabla p+\frac{1}{\mathrm{Re}}\nabla ·\left[\mu \left(\nabla \mathit{u}+{\left(\nabla \mathit{u}\right)}^T\right)\right]\hfill \end{array}$$

(3b)

$$\begin{array}{c}\hfill +\frac{1}{\mathrm{Re}\mathrm{Ca}}\kappa \delta (\mathit{x}-{\mathit{x}}_s)\widehat{\mathit{n}}+\mathrm{Bo}\rho sin\left({\alpha }_i\right).\end{array}$$

(3c)

The density and viscosity are made dimensionless using the corresponding values of fluid (1), while the pressure scale is taken as ${\widehat{p}}_{ref}={\widehat{\rho }}_1{\widehat{U}}_{ref}^2$. The two-fluid interface is characterised by a unit normal vector $\widehat{\mathit{n}}$, a delta function centered on the interface located at ${\mathit{x}}_s$, $\delta (\mathit{x}-{\mathit{x}}_s)$, and a local surface curvature $\kappa$.

The open source code Basilisk ([49,50,51]) is utilized to solve equations (3), subject to appropriate boundary and initial conditions (described in Section 3). Basilisk specializes in solving partial differential equations on adaptive Cartesian meshes. The dynamics of the two-fluid interface (Equation (3b)) is captured by a conservative, non-diffusive geometric volume-of-fluid method. The chemical heterogeneities on the inclined substrate are expressed through the local contact angle, which is implemented numerically using height functions [52,53]. The discretization of viscous forces follows an implicit formulation, while the surface tension forces are modelled using the continuous surface force method [54]. The resulting Poisson equation for the pressure is solved with a build-in iterative multigrid solver. To further increase the efficiency of the solution, Basilisk provides the option of using adaptive mesh refinement, based on customized criteria. In the context of the present study, it is desired to maintain a higher spatial resolution in the vicinity of the two-fluid interface, compared to less active regions. This leads to obtaining accurate solutions on a discretized domain with significantly less grid nodes. The code was extensively validated in a variety of multiphase flow problems (see http://basilisk.fr/) and is widely used in contact line studies (see, e.g., [27,55,56]).

2.2. Reduced-Order Model

In the limit of negligible inertia and small Ca, Cox developed an asymptotic model for the velocity normal to the contact line ${\widehat{u}}_{cl}$ [23], in the form

$${\widehat{u}}_{cl}=\frac{\widehat{\sigma }}{\widehat{{\mu }_1}}\left(\frac{G\left({\vartheta }_{*}\right)-G\left(\vartheta \right)}{ln\left({\displaystyle \frac{\lambda }{c_0}}\right)+{\displaystyle \frac{Q_o}{g\left(\vartheta \right)}}-{\displaystyle \frac{Q_i}{g\left({\vartheta }_{*}\right)}}}\right)$$

(4a)

where $\vartheta$ and ${\vartheta }_{*}$ are the apparent and local contact angles, respectively,

$$G\left(\theta \right)={\int }_0^{\theta }\frac{dx}{g\left(x\right)},$$

(4b)

and

$$g\left(x\right)=\frac{2sin\left(x\right)\left[q^2\left(x^2-{sin}^{2}x\right)+2q\left(x(\pi -x)+{sin}^{2}x\right)+\left({(\pi -x)}^2-{sin}^{2}x\right)\right]}{q\left(x^2-{sin}^{2}x\right)\left(\pi -x+cosxsinx\right)+\left({(\pi -x)}^2-{sin}^{2}x\right)\left(x-cosxsinx\right)}.$$

(4c)

In Equation (4a), $c_0$ is the azimuthally averaged droplet radius, and $q={\mu }_2/{\mu }_1$ in Equation (4c) is the viscosity ratio between the two fluids. In this formulation, $\vartheta$ is derived from the solution to the Young–Laplace equation for specified contact line shape.

In Cox’s approach, $u_{cl}$ was based on an asymtptotic expansion into the vicinity of the contact line by assuming a wedge-like geometry, with the details of the dynamics in the bulk of the droplet and the dynamics of the entire flow field in the vicinity of the contact line being encoded by the parameters $Q_o$ and $Q_i$, respectively. These, in addition to the parameter $\lambda$, which corresponds to the slip length, were fitted to the DNS data to best approximate the output of the simulations. To that end, four preliminary simulations were carried out, describing the spreading of an initially hemispherical droplet on homogeneous substrates with local contact angles of ${\vartheta }_{*}={45}^{\circ }$, ${75}^{\circ }$, ${105}^{\circ }$, and ${135}^{\circ }$, which span a broad spectrum of local contact angles. The simulations were conducted for 3 dimensionless time units, within which the droplets spread on the substrate, but are still far from equilibrium so that contact line velocities do not become vanishingly small. For carrying out the fitting of DNS data with equation (4a), an explicit expression for $\vartheta$ is used (see Equation (A13)), and contact line snapshots were collected at intervals of 0.1 dimensionless time units. Using least-squares fitting for the contact line velocity, the values obtained were $\lambda =0.220$, $Q_o=1.94$ and $Q_i=0.537$, with a mean relative deviation ${\sum }_{i=1}^N\left|\left(u_i^{cox}-u_i^{dns}\right)/u_i^{dns}\right|/N=11.1\%$, where $N=120$ is the total number of contact line snapshots considered. These parameters were held constant for all tests presented in this study.

The only difficulty in using Equation (4) for chemically heterogeneous substrates is the calculation of the apparent contact angles $\theta$. To overcome this bottleneck, an approximation to the contact angle model is obtained in Appendix A, see Equation (A16). This model is based on a perturbation expansion to the Young–Laplace equation, assuming a spherical and a circular contact line at leading order, and it naturally deviates from the expected solution as the contact line deformation becomes larger (see Figure A3 for the mean absolute error estimation). Nonetheless, due to the parameter ranges considered in the present study (described in Section 3), contact line deformations are generally weak and the approximate contact angle model yields sufficiently accurate results to be used for the training of the data-driven models.

2.3. Data-Driven Method

The aim of the data-driven modeling procedure is to approximate a mapping $G_{\mathsf{\Theta }}$ between an input and an output function space, i.e. $\mathcal{O}=G_{\mathsf{\Theta }}\left[\mathcal{I}\right]$, where $\mathcal{I}$ is an instance of the input function space, $\mathcal{O}$ is an instance of the output function space, and $\mathsf{\Theta }$ are the model parameters. Kovachki et al. [57] presented a formulation based on neural operators, approximating this mapping as,

$$G_{\mathsf{\Theta }}=\mathcal{Q}\circ \sigma \left({\mathcal{W}}_L+{\mathcal{K}}_L\right)\circ \cdots \circ \sigma \left({\mathcal{W}}_2+{\mathcal{K}}_2\right)\circ \sigma \left({\mathcal{W}}_1+{\mathcal{K}}_1\right)\circ \mathcal{P}.$$

(5)

Within this approach, the input is first lifted to a higher parameter space with operator $\mathcal{P}$. Then, a local linear operator ${\mathcal{W}}_1$ and an integral kernel operator ${\mathcal{K}}_1$ are applied, followed by the activation function $\sigma$. This last operation is applied L times with different ${\mathcal{W}}_n$ and ${\mathcal{K}}_n$ operators, before operator $\mathcal{Q}$ projects the output back to the solution space. Choosing Fourier transforms to construct the integral kernel operators, leads to the Fourier neural operator (FNO) architecture [58], where each $\sigma \left({\mathcal{W}}_n+{\mathcal{K}}_n\right)$ operator called a Fourier layer, and

$${\mathcal{K}}_n\tilde{\mathit{I}}={\mathcal{F}}^{-1}\left[{\mathcal{R}}_n\mathcal{F}\left[\tilde{\mathit{I}}\right]\right].$$

(6)

In the above expression, operator $\mathcal{F}$ transforms the input $\tilde{\mathit{I}}$ to the Fourier space, where the weights of the Fourier modes ${\mathcal{R}}_n$ are learned during training. Then, the inverse Fourier transform ${\mathcal{F}}^{-1}$ transforms the result back to the physical space. Overall, the model parameters $\mathsf{\Theta }$ consist of the elements of $\mathcal{Q}$, $\mathcal{P}$, ${\mathcal{R}}_i$ and ${\mathcal{W}}_i$, for $i=1,\dots ,L$. Figure 2 shows a schematic representation of the FNO architecture, and the interested reader is referred to the contribution by Li et al. [58] for more details.

The role of FNO-based models in the present work is correcting Cox’s reduced-order model presented in Section 2.2 in order to approximate the DNS solution. As previously mentioned, this approach is an extension of a previous work by Demou and Savva [45] that developed a data-driven workflow for wetting hydrodynamics in the limit of very small apparent contact angles, without gravity effects. In the referred study, the data-driven model provided higher order corrections to the reduced-order model. These corrections were approximately independent from the chemical heterogeneities on the substrate, leading to data-driven models were able to generalize well even for cases with significantly different heterogeneity profiles. A different approach is required in the present study because the droplets migrate mainly due to the gravity force, but their movement is also affected by the presence of the chemical heterogeneities on the substrate. As the effects of gravity are not incorporated in Cox’s reduced-order model, the developed data-driven model will focus on correcting Cox’s predictions with respect to the gravity induced translation motion. Therefore, the trained model will provide the correction to the first Fourier harmonic of the velocity normal to the contact line (i.e. translational droplet motion) and will be combined with Cox’s reduced-order model for the prediction of the higher Fourier harmonics (i.e. droplet spreading).

More specifically, the inputs to the data-driven model are the physical space representation of the first Fourier harmonics of the apparent contact angle, ${\mathcal{F}}^{-1}\left[{\mathcal{F}}_1\left[\vartheta \right]\right]\left(\varphi \right)$, the local contact angle, ${\mathcal{F}}^{-1}\left[{\mathcal{F}}_1\left[{\vartheta }_{*}\right]\right]\left(\varphi \right)$, which are essentially low-pass filters of $\vartheta \left(\varphi \right)$ and ${\vartheta }_{*}\left(\varphi \right)$, and the gravity component which is approximately normal to the contact line, $\widehat{g}sin{\alpha }_{i}sin\varphi$. All three inputs are functions of the azimuthal angle, $\varphi$, that parametrizes the contact line in a moving coordinate system centred at the centroid of the wetted area. The output of the model is expected to approximate the first Fourier harmonic of the difference between the velocities normal to the contact line as predicted by the DNS and the reduced-order model, $u_{\mathrm{cl}}^{\mathrm{AI}}\approx {\mathcal{F}}^{-1}\left[{\mathcal{F}}_1[u_{\mathrm{cl}}^{\mathrm{DNS}}-u_{\mathrm{cl}}^{\mathrm{COX}}]\right]\left(\varphi \right)$. The model is trained with $N_{\mathrm{tot}}=N_{\mathrm{train}}+N_{\mathrm{test}}$ contact line snapshots from the DNS cases, therefore the training error is defined as,

$$E_{\mathrm{train}}=\frac{1}{N_{\mathrm{train}}}\sum _{n=1}^{N_{\mathrm{train}}}\frac{{∥{\mathcal{O}}_{\mathrm{AI}}^n-{\mathcal{O}}_{\mathrm{ref}}^n∥}_2}{{∥{\mathcal{O}}_{\mathrm{ref}}^n∥}_2},$$

(7)

where, ${\mathcal{O}}_{\mathrm{AI}}$ is the output of the data-driven model and ${\mathcal{O}}_{\mathrm{ref}}$ is the corresponding reference solution. The testing error, $E_{\mathrm{test}}$, is similarly estimated using $N_{\mathrm{test}}$ samples. In all snapshots considered in the training and testing datasets, the DNS solution and Cox reduced-order model prediction never exactly match. Therefore, ${\mathcal{O}}_{\mathrm{ref}}$, which quantifies their difference, never vanishes and the summands in Equation (7) never diverge. This is even true for DNS snapshots for which the droplets are at equilibrium, since, due to the estimation errors in the contact angle, ${\widehat{u}}_{cl}$ does not vanish. Once the data-driven model is successfully trained, the data-driven velocity predictions are added to the predictions of Cox’s reduced-order model to form the hybrid contact line velocity model, $u_{\mathrm{cl}}^{\mathrm{h}}=u_{\mathrm{cl}}^{\mathrm{COX}}+u_{\mathrm{cl}}^{\mathrm{AI}}$. This hybrid model is then used to evolve the solution in time using standard time integration schemes that are typically adopted for the numerical solution of ODEs.

Many variations of the training procedure described above were tested in order to identify the best available approach to model droplet dynamics on inclined heterogeneous surfaces. These preliminary training approaches included (i) different forms of the loss function, (ii) different activation functions, (iii) retaining more Fourier harmonics in the input and output of the AI model, (iv) including the contact line shape and the gravity component normal to the inclined substrate as input to the AI model, among other attempts. The training procedure described above emerged as the best approach that balances both accuracy and robustness in predicting the future state of a moving droplet in the wide range of tests presented in Section 3.

2.4. Error Measure Based on the Fréchet Distance

The accuracy of the developed model must be assessed on its ability to predict the contact line position at a given time, a measure that is not provided by Equation (7) as it is based on the difference between the predicted and reference contact line velocities. Therefore, a more relevant error measure is utilized, defined at any specific time instance of a simulation as,

$$E_{\mathrm{Fr}}\left(t\right)=d_F\left({\mathbf{c}}_{\mathrm{AI}}\left(t\right),{\mathbf{c}}_{\mathrm{DNS}}\left(t\right)\right)\sqrt{\frac{\pi }{4A_{\mathrm{DNS}}\left(t\right)}}.$$

(8)

Here, ${\mathbf{c}}_{\mathrm{AI}}\left(t\right)$ and ${\mathbf{c}}_{\mathrm{DNS}}\left(t\right)$ are the contact line positions as predicted by the AI and the DNS at time t respectively, and $d_F({\mathbf{c}}_{\mathrm{AI}}\left(t\right),{\mathbf{c}}_{\mathrm{DNS}}\left(t\right))$ is the Fréchet distance, characterising the similarity between two curves [59]. This similarity measure takes into account both the location and ordering of the points along the curves. By normalizing the Fréchet distance with a characteristic diameter (derived from the area $A_{DNS}$ of the wetted region), this error measure becomes equal to 100% for two externally tangent circles of the same radius. In the following sections, the progress of the training and testing of the data-driven models is quantified using $E_{\mathrm{train}}$ and $E_{\mathrm{test}}$ (defined in Equation (7)), while the model performance in accurately predicting the contact line position in representative demonstration cases is quantified using $E_{\mathrm{Fr}}$ at specific times.

3. Results

3.1. Modeling Parameters

All DNS cases consider a single droplet moving on a chemically heterogeneous inclined surface, in a setup that is illustrated in Figure 1. As mentioned in Section 2, the computational domain is a $1\times 1\times 1$ box, with the solid surface at the bottom and open boundaries everywhere else. An initially hemispherical droplet of radius 0.125 was placed with the center of its contact line at some random location $(x_c,y_c)$ towards the middle upper half of the domain, where $x_c\in [0.45,0.55]$ and $y_c\in [0.2,0.25]$. The gravity field was given by $\mathbf{g}=g(sin{\alpha }_i\widehat{\mathbf{y}}-cos{\alpha }_i\widehat{\mathbf{z}})$, where $\widehat{\mathbf{y}}$ and $\widehat{\mathbf{z}}$ are the unit vectors along the y- and z-directions, and ${\alpha }_i\in [0,{80}^{\circ }]$. The density and viscosity ratios between the two fluids are set to 0.1. Considering the scales introduced in Section 2.1, the adopted values for the dimensionless groups are $\mathrm{Re}=0.1$, $\mathrm{Ca}=10$ and $\mathrm{Bo}\in [0,5]$, where the Bond number varies randomly to produce different values for g. These dimensionless groups were selected instead of the arguably more physical choice of employing droplet-based dimensionless groups using, e.g., the initial droplet radius and the contact line velocity as the length and velocity scales, leading to Re_d, Ca_d, and Bo_d, because the contact line velocity is not available a priori. Scrutinizing the simulations, a typical contact line velocity is in the order of 0.01, leading to droplet-based dimensionless groups in the order of ${\mathrm{Re}}_{\mathrm{d}}\sim {10}^{-3}$, ${\mathrm{Ca}}_{\mathrm{d}}\sim 0.1$ and ${\mathrm{Bo}}_{\mathrm{d}}\lesssim 0.1$. Considering the small Ca_d and Bo_d values, it is expected that the droplets will assume weakly perturbed spherical cap shapes as they slide on the inclined surface, such that the approximate contact angle model in Equation (A16) can provide reliable estimates.

The bottom surface assumes heterogeneity profiles generated by the multiparameter function

$$\mathsf{\Theta }(x,y)=\pi /2+p_{1}tanh\left[p_{2}cos\left(p_3\left(\tilde{x}sin{p}_5+\tilde{y}cos{p}_5\right)\right)cos\left(p_4\tilde{x}\right)\right],$$

(9)

where $\tilde{x}=xcos{p}_6-ysin{p}_6$ and $\tilde{y}=xsin{p}_6+ycos{p}_6$. Parameters $p_1,\dots ,p_6$ are uniformly distributed random numbers within the ranges,

$$p_1\in [0,\pi /36],p_2\in [-5,5],p_3,p_4\in [0,20],p_5\in [-\pi /2,+\pi /2],p_6\in [-\pi ,+\pi ].$$

A similar function was also adopted in our previous study [45], validating that the rich distribution of heterogeneity profiles enhances the diversity of the training samples that will be subsequently used for the data-driven model training.

As mentioned in Section 2.1, the simulations were carried out with Basilisk, using a Courant–Friedrichs–Lewy (CFL) criterion of CFL=0.8, and an adaptive Cartesian grid that is dynamically modified in the visinity of the two-fluid interface. More specifically, the grid is locally refined or coarsened in an octree structure, by monitoring the discretisation error of the indicator function C, as estimated by wavelet transforms. The present study adopted maximum and minimum grid spacing values of $\Delta x_{max}=2^{-5}$ and $\Delta x_{min}=2^{-8}$ respectively. To verify the accuracy of the solutions with such a grid, Figure 3 shows a comparison against a higher resolution grid with $\Delta x_{max}=2^{-6}$ and $\Delta x_{min}=2^{-9}$. Even though the higher resolution simulations used 4 times more CPU cores, the wall time was approximately 6 times longer compared to the lower resolution simulations. Nonetheless, the final solutions predicted with the two different grids are in satisfactory agreement, with an error based on the Fréchet distance of 4.8% after 60 dimensionless time units. This small deviation highlights the fact that, with $\Delta x_{max}=2^{-5}$ and $\Delta x_{min}=2^{-8}$ the available computational resources can be efficiently utilised to generate a significant sample size, without a considerable adverse effect on the solution accuracy.

Overall, 160 DNS cases were carried out for an average duration of 50 dimensionless time units. Cases that either reached the domain boundaries or got pinned on the heterogeneity patterns were ended prematurely, while cases that still exhibited droplet movement far from the boundaries were run beyond 50 time units. These simulations provided over 80,000 snapshots, taken at intervals of 0.1 dimensionless time units, to be used as training and testing samples. Each snapshot was post-processed to extract the contact lines and parametrise them in terms of 128 uniformly distributed points along the azimuthal direction $\varphi$. These contact line profiles were subsequently used to calculate: (i) the velocities normal to the contact line $u_{cl}^{DNS}(\varphi ,t)$ using finite differences, (ii) the local contact angles along the contact line from Equation (9), i.e. ${\vartheta }_{*}(\varphi ,t)=\mathsf{\Theta }(x_{cl}(\varphi ,t),y_{cl}(\varphi ,t))$, and (iii) the apparent contact angles $\vartheta (\varphi ,t)$ from Equation (A16). Overall, each sample (corresponding to a snapshot in time) comprised of

• input data: $[{\vartheta }_{*}\left(\varphi \right),\vartheta \left(\varphi \right),gsin{\alpha }_{i}sin\varphi ]$, where the last term corresponds to the component of the acceleration of gravity on the $xy$-plane which is approximately normal to the contact line, and
• target data: $u_{cl}^{DNS}\left(\varphi \right)-u_{cl}^{COX}\left(\varphi \right)$

In total, 64,000 samples were used for training (i.e. 80% of the total sample size) and the remaining 16,000 for testing. The training lasted for 500 epochs, adopting a Rectified Linear Unit (ReLU) as the activation function, a batch size of 20 samples, and an initial learning rate of ${10}^{-3}$, which was halved every 50 epochs. Taking advantage of the use of Fourier transforms, the trained model was structured to retain only the first Fourier harmonic which can be interpreted as providing corrections to Cox’s reduced order model only for the translation motion of the droplet, as opposed to local spreading which is expected to be captured by Cox’s model fairly accurately. The remaining FNO hyperparameters, i.e. the width w of the channel space after the application of the lifting operator $\mathcal{P}$, and the number of Fourier layers $FL$ (see Figure 2), were tuned after a series of preliminary experiments, considering values $w=\{32,64,128,256\}$ and $FL=\{1,2,4,8\}$. Table 1 presents the testing errors $E_{\mathrm{test}}$ for all these preliminary trained models, revealing that the values $w=64$ and $FL=8$, lead to the best performing model, with 173,121 learning parameters. Adopting the best performing set of hyperparameters, Figure 4 shows the evolution of the training and testing errors as a function of the training epochs, considering two different sample sizes. The final testing error for the smaller sample size is 16.4%, approximately 1.5% larger than the larger sample size. This observation suggests that further increasing the sample size will only have a small impact on the model accuracy. Admittedly, such testing errors may appear relatively large in the context of data-driven surrogate models. Nonetheless, the role of this data-driven model is to correcting an existing reduced-order model which, in some cases, provides predictions that are close to the target solution, with very small correction margins and, therefore, large relative errors.

Figure 5, shows the velocity contributions from the reduced-order model and the data-driven model, as well as a comparison between the overall model prediction and the reference DNS solutions for two cases, considering the same Bond number but different inclination angles. For reference, the droplet trajectories for these two cases and heterogeneity profile are depicted in Figure 8(a) and (c). As clearly shown, in these two cases the reduced-order model cannot predict the large-scale translational motion of the droplet, as explained in Section 2.3. This observation provided the motivation to combine the reduced-order model, which describes the spreading motion relatively accurately, with an AI model that focuses on the translational motion and does not interfere with spreading. As expected, the reference contact line velocities in the large inclination angle case are significantly larger than the small inclination angle case, requiring the AI contribution to bridge the gap between reduced-order model and the DNS solution. In both cases, the predictions of the overall model follow the DNS solution closely, which can be used in a method-of-lines time integration procedure to predict the droplet position in time. Within this framework, the whole numerical integration takes a few seconds on a typical laptop computer, as opposed to days of runtime on high-performance computing resources for the corresponding DNS run, amounting to a speedup of at least five orders of magnitude.

3.2. Tests

To assess the accuracy of the hybrid surrogate model in different wetting scenarios, this section presents a series of test cases adopting heterogeneities that were not encountered during training. First, Figure 6 shows comparisons between the hybrid surrogate model predictions and the corresponding DNS solutions for different heterogeneity profiles, considering ${\alpha }_i={60}^{\circ }$ and $\mathrm{Bo}=5$. Figure 6(a) uses a homogeneous profile with $\mathsf{\Theta }(x,y)={80}^{\circ }$, while panels (b) and (c) of the same figure assume profiles that are defined by specific parameter choices in Equation (9). For the homogeneous case, Cox’s reduced-order model plays a minor role in the contact line velocity prediction since the contact line quickly assumes a constant shape and the droplet descends with constant speed throughout the simulation. Hence, to a large extend, the movement of the droplet is described by the AI component of the hybrid model. The hybrid surrogate model results in a 2.9% Fréchet-based error, which attests to the capability of the AI component to provide a reliable estimation of the translation motion of the droplet, at least on simple homogeneous substrates. In Figure 6(b), the droplet experiences pinning and depinning states as it moves over the horizontal striped features, altering the contact line shape and the apparent contact angles, therefore the importance of Cox’s reduced-order model contribution is elevated. Still, in this scenario the predictions of the hybrid surrogate model exhibit an error as low as 4.8% compared to the reference DNS solution. Finally, Figure 6(c) exhibits a case where the droplet is quickly pinned due to the configuration of the chemical heterogeneities, in which case the hybrid surrogate model excellently reproduces the dynamics. It is worth noting that in all cases presented in Figure 6, the trajectory of the centroid of the wetted area is accurately reproduced by the hybrid surrogate model.

Figure 6. Comparisons between the predictions of the hybrid surrogate model (red lines) and DNS solutions (semi-transparent blue regions) at the final time step of the simulations $t_f=60$. The initial conditions are also shown, centered approximately around $(0.5,0.2)$. The trajectories of the centroid of the wetted area are shown with black and red-dashed lines for the DNS and the hybrid surrogate model solutions respectively. The heterogeneity profile in (a) assumes a uniform value of ${80}^{\circ }$. The heterogeneity profiles in all other figures are described by Equation (9), with specific parameters, $\{p_1,\cdots ,p_6\}$= (b) $\{\pi /36,5,10,0,0,0\}$, and (c) $\{\pi /36,5,10,10,0,0\}$. In all cases ${\alpha }_i={60}^{\circ }$ and $\mathrm{Bo}=5$. The errors based on the Fréchet distance are (a) 2.9%,(b) 4.8%, and (c) 1.1%. The background follows the same shading convention as Figure 3.

Figure 6. Comparisons between the predictions of the hybrid surrogate model (red lines) and DNS solutions (semi-transparent blue regions) at the final time step of the simulations $t_f=60$. The initial conditions are also shown, centered approximately around $(0.5,0.2)$. The trajectories of the centroid of the wetted area are shown with black and red-dashed lines for the DNS and the hybrid surrogate model solutions respectively. The heterogeneity profile in (a) assumes a uniform value of ${80}^{\circ }$. The heterogeneity profiles in all other figures are described by Equation (9), with specific parameters, $\{p_1,\cdots ,p_6\}$= (b) $\{\pi /36,5,10,0,0,0\}$, and (c) $\{\pi /36,5,10,10,0,0\}$. In all cases ${\alpha }_i={60}^{\circ }$ and $\mathrm{Bo}=5$. The errors based on the Fréchet distance are (a) 2.9%,(b) 4.8%, and (c) 1.1%. The background follows the same shading convention as Figure 3.

Figure 7 shows various droplet snapshots in time, for a wetting scenario with a more intricate dynamical behavior, considering the same inclination angle and Bond number as before. The presence of diagonal striped features influences the droplet descent on the inclined substrate, where the centroid of the wetted area follows a curved trajectory. Initially, the droplet shifts to the left until the contact line meets a more hydrophobic stripe (darker shaded) at the left side of the bubble, somewhere between $t=10$ and 20. Afterwards the droplet moves diagonally, following the stripes configuration until $t=40-50$, when the right part of the droplet lies completely on a more hydrophobic stripe. At that point, the droplet movement aligns with the gravity force until it becomes diagonal once again. This complicated movement, where the droplet is significantly influenced by both the gravity force and the chemical heterogeneities, is accurately captured by the hybrid surrogate model, exhibiting an error of only 4.3%, compared to DNS solution at the end of the simulation.

Figure 7. Similar to Figure 6, with different panels corresponding to different time instances. The heterogeneity profile adopted is described by Equation (9), with specific parameters, $\{p_1,\cdots ,p_6\}=\{\pi /36,2,0,10,0,-\pi /6\}$. The error based on the Fréchet distance at the final time is 4.3%.

Figure 7. Similar to Figure 6, with different panels corresponding to different time instances. The heterogeneity profile adopted is described by Equation (9), with specific parameters, $\{p_1,\cdots ,p_6\}=\{\pi /36,2,0,10,0,-\pi /6\}$. The error based on the Fréchet distance at the final time is 4.3%.

To quantify the sensitivity of the hybrid surrogate model on the characteristics of the gravity force, Figure 8 shows comparisons against DNS solutions for different values of the inclination angle and the Bond number. As expected, keeping a constant Bond number value of $\mathrm{Bo}=5$ and increasing the inclination angle (panels (a)-(c)) increases the distance travelled by the droplet within the 80 time units considered for the simulation. Similar to Figure 7, the centroid follows a curved trajectory, guiding the droplet towards the left side, in accordance to the structure of the chemical heterogeneities. Quantifying the agreement between predicted and reference solutions, the errors based on the Fréchet distance increase from 8.2% at ${\alpha }_i={40}^{\circ }$ to 16.6% at ${\alpha }_i={60}^{\circ }$ and remain relatively unaffected to 14.1% at ${\alpha }_i={80}^{\circ }$. Furthermore, Figure 8(d)-(f) shows cases considering a fixed inclination angle ${\alpha }_i={60}^{\circ }$ and a varying Bond number. For the lowest Bond number, the droplet barely moves and soon gets pinned on the heterogeneity features. By increasing the Bond number, the droplet is able to overcome the heterogeneity barriers and induce a downhill sliding motion, with a corresponding increase of the Fréchet-based error, noting a maximum error of 25.6% for the largest Bond number considered. A major reason behind this significant error is the that $\mathrm{Bo}=8$ is outside the Bond number ranges used to train this model and, therefore, the data-driven model is trying to extrapolate beyond the training distribution. Another source of error is that the increasing Bond number can be thought of as an equivalent increase in g, something that influences the three-dimensional droplet shape and has an impact on the assumption of a weakly perturbed spherical cap used to develop the approximate contact angle model in Equation (A16). Overall, the errors of the hybrid surrogate model increase with increasing inclination angle and Bond numbers, while the model exhibits a satisfactory performance as long as the adopted parameters are within the ranges considered during training.

Figure 8. Similar to Figure 6, with different panels corresponding to different values of the inclination angle ${\alpha }_i$ and Bond number Bo. The comparison is documented at $t=80$, except panel (f) where $t=40$ is used. The heterogeneity profile adopted is described by Equation (9), with specific parameters, $\{p_1,\cdots ,p_6\}=\{\pi /72,2,5,10,\pi /4,\pi /3\}$. The errors based on the Fréchet distance are (a) 8.2%,(b) 16.6%, (c) 14.1%, (d) 5.1%, (e) 7.6%, and (f) 25.6%.

Figure 8. Similar to Figure 6, with different panels corresponding to different values of the inclination angle ${\alpha }_i$ and Bond number Bo. The comparison is documented at $t=80$, except panel (f) where $t=40$ is used. The heterogeneity profile adopted is described by Equation (9), with specific parameters, $\{p_1,\cdots ,p_6\}=\{\pi /72,2,5,10,\pi /4,\pi /3\}$. The errors based on the Fréchet distance are (a) 8.2%,(b) 16.6%, (c) 14.1%, (d) 5.1%, (e) 7.6%, and (f) 25.6%.

4. Conclusions

In this study, a new approach was proposed to model droplet dynamics on inclined, chemically heterogeneous surfaces, a type of process that, up to this point, could only be simulated through DNS. The proposed approach is a hybrid surrogate model for the contact line velocity, comprising of both analytically derived reduced-order model and data-driven components. The reduced-order model by Cox, coupled to a developed approximate model for the apparent contact angle, is able to model the effects of the chemical heterogeneities, but it does not include the effects of gravity. To complement the reduced-order model, the data-driven component was trained to capture the effects of gravity by modeling the large-scale translation motion of the droplet, using DNS datasets. In this way, the hybrid surrogate model is able to extend the applicability of the reduced-order model by Cox to gravity-driven droplet scenarios. As the hybrid surrogate model predicts the contact line velocities, the predictions are integrated in time (similar to an ODE solver) to simulate the droplet movement much more efficiently compared to the corresponding DNS.

To assess the performance of the hybrid surrogate model a number of different wetting scenarios were considered. The model was found to accurately capture the droplet movement on different heterogeneity profiles that are within the distribution used during the training phase. Nonetheless, the hybrid model exhibited increased errors as the inclination angle and the Bond number increased, but a reasonable agreement was observed as long as the model did not extrapolated to scenarios that were outside the training distribution.

This approach can be further refined to improve its accuracy and generalizability in the near future. For example, the approximate model for the calculation of the apparent contact angle can be replaced by the full calculation of the Young–Laplace equation (e.g. by using Surface Evolver [60]), but doing so is likely to add a substantial computational overhead. In a similar manner, we may separately improve the current contact angle approximation with an another data-driven model in order to better match the Young–Laplace solution. The improved hybrid surrogate model can then be extended to other settings of specific engineering interest. As an example, incorporating phase change dynamics in the form of evaporation and condensation, can prove beneficial in a range of industrial applications such as desalination [61], water harvesting [62], biomedical applications [63], printing, coating and cooling [64], among others.