2. Multiphase many-body dissipative particle dynamics simulations

Three-dimensional simulations were performed using the mDPD method (Li et al. Reference Li, Hu, Wang, Ma and Zhou2013; Xia et al. Reference Xia, Goral, Huang, Miskovic, Meakin and Deo2017; Rao et al. Reference Rao, Xia, Li, McConnell, Sutherland and Li2021), which is an extension of the traditional dissipative particle dynamics (DPD) model (Español & Warren Reference Español and Warren1995; Groot Reference Groot2004); DPD is a mesoscale simulation technique for studies of complex fluids, particularly multiphase systems, such as emulsions, suspensions and polymer blends (Ghoufi & Malfreyt Reference Ghoufi and Malfreyt2012; Lei et al. Reference Lei, Bertevas, Khoo and Phan-Thien2018; Zhao et al. Reference Zhao, Chen, Zhang and Liu2021). The relation between DPD and other coarse-grained methods and atomistic simulations has been studied and discussed by Murtola et al. (Reference Murtola, Bunker, Vattulainen, Deserno and Karttunen2009), Li et al. (Reference Li, Bian, Yang and Karniadakis2016), Chan, Li & Wenzel (Reference Chan, Li and Wenzel2023) and Español & Warren (Reference Español and Warren2017).

In DPD and mDPD models, the position ($\boldsymbol {r}_i$) and velocity ($\boldsymbol {v}_i$) of a particle $i$ with mass $m_i$ are governed by Newton's equations of motion in the form

(2.1)\begin{equation} \left.\begin{gathered} \frac{{\rm d}\boldsymbol{r}_i}{{\rm d} t} = \boldsymbol{v}_i, \\ m_i\,\frac{{\rm d}\boldsymbol{v}_i}{{\rm d} t} = \boldsymbol{F}_{i} = \sum_{j\ne i}\boldsymbol{F}_{ij}^{C} + \boldsymbol{F}_{ij}^{D} + \boldsymbol{F}_{ij}^{R}. \end{gathered}\right\} \end{equation}

The total force on particle $i$, i.e. $\boldsymbol {F}_i$, consists of three pairwise components, i.e. the conservative $\boldsymbol {F}^{C}$, dissipative $\boldsymbol {F}^{D}$, and random $\boldsymbol {F}^{R}$ forces. The latter two are identical in DPD and mDPD models, given by

(2.2)$$\begin{gather} \boldsymbol{F}_{ij}^{D} =- \gamma\,\omega_{D}(r_{ij})\,( \boldsymbol{v}_{ij} \boldsymbol{\cdot} \boldsymbol{e}_{ij})\,\boldsymbol{e}_{ij}, \end{gather}$$

(2.3)$$\begin{gather}\boldsymbol{F}_{ij}^{R} = \zeta\,\omega_{R}(r_{ij})\,({\rm d}t)^{-1/2} \xi_{ij}\,\boldsymbol{e}_{ij}, \end{gather}$$

where $\boldsymbol {e}_{ij}$ is a unit vector, $\omega _{D}$ and $\omega _{R}$ are weight functions for the dissipative and random forces, and $\xi _{ij}$ is a pairwise conserved Gaussian random variable with zero mean and second moment $\langle \xi _{ij}(t)\,\xi _{kl} (t') \rangle = ( \delta _{ik} \delta _{jl} + \delta _{il} \delta _{jk})\,\delta (t-t')$, where $\delta _{ij}$ is the Kronecker delta, and $\delta (t-t')$ is the Dirac delta function. Together, the dissipative and random forces constitute a momentum-conserving Langevin-type thermostat. The weight functions and the constants $\gamma$ and $\zeta$ are related via fluctuation–dissipation relations first derived by Español & Warren (Reference Español and Warren1995),

(2.4)$$\begin{gather} \omega_{D} = (\omega_{R})^2, \end{gather}$$

(2.5)$$\begin{gather}\zeta = \sqrt{2\gamma k_{B}T}, \end{gather}$$

in which $k_{B}$ is the Boltzmann constant, and $T$ is the temperature. This relation guarantees the canonical distribution (Español & Warren Reference Español and Warren1995) for fluid systems in thermal equilibrium. The functional form of the weight function is not specified, but the most common choice (also used here) is

(2.6)\begin{equation} \omega_{D}(r_{ij}) = \left\{\begin{array}{@{}ll} (1 - r_{ij}/r_d)^s & \mbox{for}\ r_{ij} \le r_{d}, \\[2pt] 0 & \mbox{for}\ r_{ij}>r_{d}. \end{array}\right. \end{equation}

Here, $s=1.0$ is used, and $r_{d}$ defines a cutoff distance for the dissipative and random forces.

Although the above equations are the same for both DPD and mDPD, they differ in their conservative forces. Here, we use the form introduced by Warren (Reference Warren2001, Reference Warren2003):

(2.7)\begin{equation} F^{C}_{ij} = A\,\omega_{C}(r_{ij})\,\boldsymbol{e}_{ij} + B (\rho_i + \rho_j) \omega_{B}\,\boldsymbol{e}_{ij}. \end{equation}

The functional form of both weight functions $\omega _C$ and $\omega _B$ is the same as $\omega _{D}$ in (2.6), but with different cutoff distances $r_c$ and $r_b$. The first term in (2.7) is the standard expression for the conservative force in DPD, and the second is the multi-body term. The constants $A$ and $B$ are chosen such that $A <0$ for attractive interactions, and $B>0$ for repulsive interactions; note that in conventional DPD, $A >0$ and $B=0$. The key component is the weighted local density

(2.8)\begin{equation} \rho_i = \sum_{j \ne i} \omega_\rho(r_{ij}). \end{equation}

There are several ways to choose the weight function (Zhao et al. Reference Zhao, Chen, Zhang and Liu2021), and here, the normalized Lucy kernel (Lucy Reference Lucy1977) in three dimensions is used:

(2.9)\begin{equation} \omega_\rho(r_{ij}) = \frac{105}{16 {\rm \pi}r^3_{{c}\rho}} \left(1+ \frac{3 r_{ij}}{r_{{c}\rho}}\right) \left(1- \frac{r_{ij}}{r_{{c}\rho}}\right)^3, \end{equation}

with a cutoff distance $r_{c\rho }$ beyond which the weight function $\omega _\rho$ becomes zero.

2.1. Simulation parameters and system set-up

To simulate molten CMAS, the parameter mapping of Koneru et al. (Reference Koneru, Flatau, Li, Bravo, Murugan, Ghoshal and Karniadakis2022) was used together with the open-source code LAMMPS (Thompson et al. Reference Thompson2022). In brief, the properties of molten CMAS at approximately 1260 $^\circ$C, based on the experimental data from Naraparaju et al. (Reference Naraparaju, Gomez Chavez, Niemeyer, Hess, Song, Dingwell, Lokachari, Ramana and Schulz2019), Bansal & Choi (Reference Bansal and Choi2014) and Wiesner et al. (Reference Wiesner, Vempati and Bansal2016), were used. In physical units, density was 2690 kg m$^{-3}$, surface tension 0.46 N m$^{-1}$, and viscosity 3.6 Pa s. Using the density and surface tension to estimate the capillary length ($\kappa = (\sigma /(\rho g))^{1/2}$) gives 4.18 mm. The droplets in the simulations (details below) had linear sizes shorter than the capillary length, hence gravity was omitted. In terms of physical units, time was $6.297\times 10^{-6}$ s, length was $17.017\times 10^{-6}$ m, and mass was $1.964 \times 10^{-8}$ kg.

Using the above values, droplets of initial radii $d = 8, 9 , 10, 11, 12$ in mDPD units, corresponding to $R_{0}= 0.136, 0.153, 0.17, 0.187, 0.204$ mm, respectively, were used in the simulations; all of them are smaller than the capillary length. The time step was 0.002 (mDPD units) corresponding to 12.59 ns. In addition, $k_{B}T =1$, $r_c =1.0$, $r_b=r_{{c}\rho }=0.75$, $r_d=1.45$, $\gamma =20$ and $B=25$. The attraction parameter $A$ in (2.7) has to be set for the interactions between the liquid particles ($A_{ll}$), and the liquid and solid particles ($A_{ls}$). The former was set to $A_{ll}=-40$, and $A_{ls}$ was chosen based on simulations that provided the desired $\theta _{eq}$, thus allowing for controlled variation of $\theta _{eq}$ (see figure 2).

Figure 2.

The equilibrium contact angles $\theta _{eq}$ for the different attraction parameters between the liquid and solid particles ($A_{ls}$; see (2.7)). It is worth noting that the data for this figure have been extracted from Koneru et al. (Reference Koneru, Flatau, Li, Bravo, Murugan, Ghoshal and Karniadakis2022).

The initial configuration of the droplet and the solid wall was generated from a random distribution of equilibrated particles with number density $\rho =6.74$. This amounts to approximately 60 660 particles in the wall, and depending on the initial radius, anywhere between 14 456 and 48 786 particles in the droplet. Periodic boundary conditions are imposed along the lateral directions, and a fixed, non-periodic boundary condition is imposed along the wall-normal direction. Since mDPD is a particle-based method, the spreading radius and the dynamic contact angle are approximated using surface-fitting techniques. First, the outermost surface of the droplet is identified based on the local number density, i.e. particles with $\rho \in [0.45, 0.6]$. The liquid particles closest to the wall are fitted to a circle of radius $r$, i.e. the spreading radius. On the other hand, a sphere with the centroid of the droplet as the centre is fitted to the surface particles to compute the contact angle. The contact angle is defined as the angle between the tangent at the triple point (liquid–solid–gas interface) and the horizontal wall. The wall in these simulations is made up of randomly distributed particles to eliminate density and temperature fluctuations at the surface. Following Li et al. (Reference Li, Bian, Tang and Karniadakis2018), the root mean square height ($R_q$) of the surface scales linearly with $1/\sqrt {N_w}$, where $N_w=(2{\rm \pi} r_{cw}^3/3) \rho _w$ is the number of neighbouring particles. In this work, $R_q$ comes out to be around 0.0708 mDPD units, or $1.2\,\mathrm {\mu }$m.

As the CMAS droplet spreads on the substrate, it loses its initial spherical shape and begins to wet the surface as depicted in figure 3, forming a liquid film between the droplet and the substrate. Understanding how droplets behave on surfaces is important for a wide range of applications, including in industrial processes, microfluidics, propulsion materials, and the design of self-cleaning surfaces (Pitois & François Reference Pitois and François1999; Chen et al. Reference Chen, Bonaccurso, Deng and Zhang2016; Hassan et al. Reference Hassan, Yilbas, Al-Sharafi and Al-Qahtani2019; Jain et al. Reference Jain, Lemoine, Chaussonnet, Flatau, Bravo, Ghoshal, Walock and Murugan2021; Nieto et al. Reference Nieto, Agrawal, Bravo, Hofmeister-Mock, Pepi and Ghoshal2021).

Figure 3.

(a) Illustration of the spreading behaviour of a CMAS droplet on a hydrophilic surface at different times. The droplet with initial size $R_0$ spreads on the surface with radius $r(t)$ and contact angle $\theta (t)$. (b) A series of snapshots from a simulation of a droplet with initial size $R_0 = 0.136$ mm and equilibrium contact angle $\theta _{eq} = 93.4^\circ$.

3. Simulation results

The size of a droplet changes over time. By tracking the changes, we can gain insight into the physical processes involved in spreading. The time evolution of the droplet radius ($r(t)$) is shown in figure 4. The log-log plots show the effect of the initial drop size $R_0$ and equilibrium contact angles $\theta _{eq}$ on the radius $r(t)$. Figure 4(a) displays $r(t)$ for initial drop sizes $R_0 = 0.136, 0.153, 0.17, 0.183, 0.204$ mm and equilibrium contact angle $\theta _{eq}=54.6^\circ$. Similarly, figure 4(b) shows the spreading radius for different equilibrium contact angles ($\theta _{eq}=93.4^\circ$, $85.6^\circ$, $77.9^\circ$, $70.1^\circ$, $62.4^\circ$, $54.6^\circ$, $45.3^\circ$ and $39.1^\circ$) with initial drop size $R_0=0.136$ mm.

Figure 4.

The impact of $\theta _{eq}$ and $R_0$ on the droplet radii as a function of time for various (a) initial drop sizes with equilibrium contact angle $\theta _{eq} = 54.6^\circ$ corresponding to $A_{ls} = 30.0$, and (b) equilibrium contact angles (corresponding to $A_{ls} = -25.0, -25.8, -27.0, -28.0, -29.0, -30.0, -31.4, -32.2$) and initial drop size $R_{0}=0.136$ mm.

Eddi, Winkels & Snoeijer (Reference Eddi, Winkels and Snoeijer2013) used high-speed imaging with time resolution covering six decades to study the spreading of water–glycerine mixtures on glass surfaces. By varying the amount of glycerine, they were able to vary the viscosity over the range 0.0115–1.120 Pa s. They observed two regimes, the first one for early times with $\alpha$ changing continuously as a function of time from $\alpha \approx 0.8$ to $\alpha \approx 0.5$. This was followed by a sudden change to the second regime in which $\alpha$ settled to $0.1 < \alpha < 0.2$. As pointed out by Eddi et al. (Reference Eddi, Winkels and Snoeijer2013), the second regime agrees with Tanner's law (Tanner Reference Tanner1979). All of their systems displayed complete wetting.

Based on $r(t)$ of the CMAS drops and $\alpha$ shown in figure 4, it can be observed that $r(t)$ (and based on the power law, $\alpha$) depends on both the initial drop size and the equilibrium contact angle. The values of $\alpha$ for some simulation datasets are plotted over time in figure 5. The plot shows the behaviour of $\alpha$ for different initial drop sizes and equilibrium contact angles $\theta _{eq}=62.4^\circ$ and $85.6^\circ$.

Figure 5.

The value of $\alpha$, calculated using (1.2), varies for different initial radii and fixed equilibrium contact angles (a) $\theta _{eq}=85.6^\circ$ and (b) $\theta _{eq}=62.4^\circ$. The figure illustrates that $\alpha$ is influenced by both the initial drop size $R_0$ and the equilibrium contact angle $\theta _{eq}$.

Inspired by the experimental results of Eddi et al. (Reference Eddi, Winkels and Snoeijer2013), the simulations of Koneru et al. (Reference Koneru, Flatau, Li, Bravo, Murugan, Ghoshal and Karniadakis2022) and the current simulations, we propose a simple sigmoid type dependence for $\alpha$:

(3.1) \begin{equation} \frac{{\rm d}\ln (r)}{{\rm d}\ln (t)} = \alpha(t, R_0, \theta_{eq}):= \eta \Big[\frac{1}{1+\exp(\beta (\tau-\ln(t)))}-1\Big]\!. \end{equation}

The two constant values of $\alpha$ discussed in the above references are the two extrema of the sigmoid curve, given that the transition between the two regimes occurs at $\ln (t_{transition})=\tau$.

The parameters of (3.1) are discovered by PINNs, and their dependence on $R_0$ and $\theta _{eq}$ is then expressed using symbolic regression. The general steps in the discovery of the droplet spreading equation and the extraction of the unknown parameters $\eta (R_0, \theta _{eq})$, $\beta (R_0, \theta _{eq})$ and $\tau ( R_0, \theta _{eq})$ are shown in figure 6 and can be summarized as follows.

1. (i) Data collection: for this study, data are collected by conducting three-dimensional simulations using the mDPD method in LAMMPS with varying initial drop sizes $R_0$ and equilibrium contact angles $\theta _{eq}$.

2. (ii) PINNs: the input of the network is time $t$, and the output of the network is the spreading radii $\tilde {r}(t)$. The physics-informed part of the PINNs is encapsulated in designing the loss function. In this study, the ‘goodness’ of the fit is measured by (a) deviation from trained data together with (b) deviation of network predictions and those from ordinary differential equation (ODE) (3.1) solutions. This optimization process reveals the the unknown parameters $\eta ( R_0, \theta _{eq})$, $\beta ( R_0, \theta _{eq})$ and $\tau ( R_0, \theta _{eq})$.

3. (iii) Data interpolation: after PINNs are trained, we used their predictions together with two additional multilayer perceptron neural networks to fill the sparse parameter space. This step helps our next goal, which is relating the ODE parameters to $R_0$ and $\theta _{eq}$ without performing three-dimensional simulations.

4. (iv) Symbolic regression: discovering a mathematical expression for each unknown parameter $\eta (R_0, \theta _{eq})$, $\beta (R_0, \theta _{eq})$ and $\tau ( R_0, \theta _{eq})$ of (3.1).

5. (v) In order to quantify the uncertainty associated with our predictions, we utilize B-PINNs and leverage the insights gained from the PINNs prediction and symbolic regression, with a specific emphasis on the known value of $\eta$. By employing B-PINNs, we can effectively ascertain the values of two specific parameters, $\beta$ and $\tau$, which in turn enable us to quantify the uncertainty in our predictions.

Figure 6.

The process of utilizing PINNs to extract three unknown parameters of the ODE (3.1), using three-dimensional mDPD simulation data. First, a neural network is trained using simulation data, where the input is time $t$ and the output is spreading radii $\tilde {r}(t)$. This neural network comprises four layers with three neurons, and is trained for 12 000 epochs. Subsequently, the predicted $\tilde {r}(t)$ is used to satisfy (3.1) in the physics-informed part. The loss function for this process consists of two parts: data matching and residual. By optimizing the loss function, the values of $\eta (R_0, \theta _{eq})$, $\beta (R_0, \theta _{eq})$ and $\tau (R_0,\theta _{eq})$ are determined for each set of $R_0$ and $\theta _{eq}$. After predicting the unknown parameters using PINNs, two additional neural networks, denoted as $NN_\beta$ and $NN_\tau$, are trained using these parameters to generate values for the unknown parameters at points where data are not available. The outputs of these networks, together with the outputs of the PINNs, are then fed through a symbolic regression model to discover a mathematical expression for the discovered parameter.

4. Physics-informed neural networks

In this section, we delve into the process of uncovering the parameters of the proposed ODE using the available simulation data.

In recent years, the concept of approximating the behaviours of complex nonlinear dynamical systems with experimental or numerically simulated data has gained significant traction, driven by the exponential growth in available data. The conventional approach for model fitting involves nonlinear least squares, which relies on iterative optimization algorithms (Srinivasan & Mason Reference Srinivasan and Mason1986). These algorithms can be computationally expensive and sensitive to initial conditions, sparsity, and high dimensionality of the data. To circumvent these issues, Brunton et al. (Reference Brunton, Proctor and Kutz2016) proposed a framework named SINDy (sparse identification of nonlinear dynamical system), which selects a class of algebraic functions along with an arbitrary order of derivatives. By performing sparse regression, which assumes that the contributions of many of these functions are negligible, SINDy can accurately discover ODEs and PDEs. On the other hand, equation learner (EQL; Martius & Lampert Reference Martius and Lampert2016) leverages the ability to use algebraic functions as neural network units, which brings interpretability to the model fitting process. Both SINDy and EQL share a reliance on predefined basis functions, constraining their ability to represent complex functional forms beyond those that are predefined. Therefore, applying these methods in the present study is not feasible. These coefficients not only are the coefficients of the basis functions but also cannot be utilized to define a library of functions, primarily due to the structure of the proposed ODE.

Physics-informed neural networks form another promising class of approaches that leverages the flexibility and scalability of deep neural networks to discover governing equations, incorporating physical laws (PDEs, ODEs, integro-differential equations (IDEs), etc.) and data (initial conditions plus boundary conditions) into the loss function (Raissi et al. Reference Raissi, Perdikaris and Karniadakis2019; Chen et al. Reference Chen, Lu, Karniadakis and Dal Negro2020; Mao, Jagtap & Karniadakis Reference Mao, Jagtap and Karniadakis2020; Mishra & Molinaro Reference Mishra and Molinaro2020; Shukla et al. Reference Shukla, Di Leoni, Blackshire, Sparkman and Karniadakis2020; Karniadakis et al. Reference Karniadakis, Kevrekidis, Lu, Perdikaris, Wang and Yang2021). In contrast, EQL and SINDy focus on deriving equations from data without enforcing the governing laws, lacking a guarantee of adherence to these principles. In addition, PINNs leverage automatic differentiation for computing the exact derivatives of differential equations, while SINDy, considered as an optimization method, relies on numerical differentiation for gradient approximations. The PINNs represent a more adaptable approach, capable of adjusting to complex systems by learning directly from data, even when a predefined mathematical model is completely or partially absent.

The method employed to find unknown parameters in this study exhibits potential for broader applications, such as different fluid properties, reduced-order PDEs, and scenarios in which only a part of the equation is absent (grey box learning framework; Kiyani et al. Reference Kiyani, Shukla, Karniadakis and Karttunen2023). While, potentially, the parameters of the proposed ODE could be discovered using other aforementioned approaches for the purpose of this paper, the scalability and extensibility to PDEs or more complex ODEs, and insensitivity to initial conditions in the fitting process compared to other methods, suggest that PINNs are the most suitable choice for this task.

4.1. Discovering parameters of the ODE

As discussed in § 3, our study aims to identify the values of the parameters $\eta (R_0,\theta _{eq})$, $\beta (R_0,\theta _{eq})$ and $\tau (R_0, \theta _{eq})$ in the ODE (3.1). The general steps of the framework are shown in figure 6(I). PINNs take time $t$ as an input to predict $\tilde {r}(t)$ at each time.

For each unique combination of a contact angle and initial drop size, the data were divided randomly into training, validation and test sets at ratio 80 : 10 : 10. A neural network with four layers of three neurons each was trained on the training data using the JAX and FLAX Python libraries on the M1 Apple chip. The training process consisted of 12 000 epochs, and took approximately 5 minutes to complete for each case.

The predicted values for the time evolution of the radii $\tilde {r}(t)$ should satisfy the data and the physics-informed step, i.e. meet the requirements of the ODE (3.1). The two-component loss function, designed to meet the requirements, consists of ${Loss_{data}}$ and ${Loss_{ODE}}$, given by

(4.1)\begin{equation} \left.\begin{gathered} {Loss_{data}}=\frac{1}{N_{k}} \sum_{i=1}^{i=N_{k}} | (r_{i}(t) - \tilde{r}_{i}(t)) | ^{2}, \\ {Loss_{ODE}} = \frac{1}{N_{r}} \sum_{i=1}^{i=N_{r}} \left | \frac{{\rm d}\ln (\tilde{r}) }{{\rm d}\ln (t)} - \eta \left(\frac{1}{1+\exp(-\beta(\ln(t)-\tau))}-1\right) \right | ^{2}, \end{gathered}\right\} \end{equation}

where $\tilde {r}(t)$ and $r(t)$ stand for the radii from the prediction and simulation, respectively. Here, $N_{k}$ is the number of training points, and $N_{r}$ is the number of residual points.

Figures 7, 8 and 9 illustrate the results obtained by utilizing PINNs to discover the parameters of the ODE (3.1), which describes the dynamics of the radii of the CMAS drops.

Figure 7.

Comparison of the time evolution of the droplet radii: mDPD simulations (symbols), ODE model (3.1) (solid lines) and PINN predictions (dashed lines) for $\theta _{eq}=\{39.1^\circ, 62.4^\circ,93.4^\circ \}$ and $R_0=\{0.136, 0.153, 0.187,0.204\}$ mm parameter sets.

Figure 8.

(a,b,c) The evolution of parameters $\eta (R_0,\theta _{eq})$, $\beta (R_0,\theta _{eq})$ and $\tau (R_0,\theta _{eq})$, respectively, over multiple epochs. These plots demonstrate that the parameters gradually converge to a stable state after $12\,000$ epochs. (d) The traces of the loss function for the PINNs framework. The learning curves demonstrate the decreasing trend of the loss functions, indicating that they converge to a stable point for all initial drop sizes and $\theta _{eq}$.

Figure 9.

The values of (a) $\eta$, (b) $\beta$ and (c) $\tau$ obtained through PINNs. These values exhibit varying behaviours depending on the initial radius $R_{0}$ and equilibrium contact angle $\theta _{eq}$. The horizontal axes display the equilibrium contact angles $\theta _{eq}$. The vertical axes of all plots represent the values of $\eta$, $\beta$ and $\tau$. Here, $\eta$ remains nearly constant within a small range of values between $-0.325$ and $-0.200$, and $\beta$ and $\tau$ change within ranges from $1.0$ to $5.0$, and $6.5$ to $8.0$, respectively.

Figure 7 shows comparisons of the simulation data ($r(t)$), the prediction ($\tilde {r}(t)$) and the solution of the ODE (3.1). The figure demonstrates a remarkable degree of agreement between simulations, PINNs and our ODE model. The first three panels in figure 8 show the convergence of $\eta (R_0,\theta _{eq})$, $\beta (R_0,\theta _{eq})$, and $\tau (R_0,\theta _{eq})$ parameters during training. One can conclude that, all three parameters stabilize roughly after 12 000 epochs. In the rightmost panel of figure 8, the loss function (4.1) history is plotted against the training epochs, stabilizing around $10^{-4}$. This indicates successful training of the PINNs model. The results shown in figures 7 and 8 demonstrate the capability of our proposed framework to accurately predict the spreading radius of CMAS across the different initial radii and equilibrium contact angles.

Figures 9(a,b,c) show the values of $\eta (R_0,\theta _{eq})$, $\beta (R_0, \theta _{eq})$ and $\tau (R_0,\theta _{eq})$, respectively, obtained by PINNs for different initial radii $R_{0}$ ranging from $0.136$ to $0.204$ mm, and equilibrium contact angles $\theta _{eq}$ ranging from $39.1^\circ$ to $93.4^\circ$. The results show that $\eta$ changes within a small window between $-0.325$ and $-0.200$ for all $R_{0}$ and $\theta _{eq}$. However, the changes in $\beta$ (between $1$ and $5$) and $\eta$ (between $6.0$ and $8.0$) are significant, indicating that these parameters depend strongly on $R_{0}$ and $\theta _{eq}$.

4.2. Generate more samples of feasible radii and contact angles

As discussed earlier, the parameters in our ODE model (3.1) are functions of the initial radius and the equilibrium contact angle. Using PINNs, we were able to find the values for those parameters. To find a closed-form relation between the parameters $R_0$ and $\theta _{eq}$, more data than the rather small current set are needed. Performing three-dimensional mDPD simulations is, however, computationally expensive. In this subsection, we train two additional neural networks to capture the nonlinear relation between the ODE parameters found by PINNs, and the variables $R_0$ and $\theta _{eq}$. Then we will use these trained networks to fill our sparse parameter space to perform symbolic regression in the next section.

Specifically, we generate values for $R_{0}$ in the interval $[0.136, 0.204]\ \mathrm {mm}$ and $\theta _{eq}$ in the range $[40^\circ, 95^\circ ]$, as shown in figure 6(II).

Two fully connected networks $NN_{\beta }$ and $NN_{\tau }$ consist of eight dense layers with $256/256/256/128/64/32/16/8$ neurons. These networks are trained using an Adam optimizer with learning rate $10^{-2}$ for a total of 4000 epochs.

Discovered parameters from $NN_{\tau }$ and $NN_{\beta }$ are visualized in figures 10(a,b), respectively, where the parameter values obtained from PINNs are denoted by green and red circles, indicating the training data for $NN_{\beta }$ and $NN_{\tau }$, respectively. Furthermore, the parameter values generated by $NN_\beta$ and $NN_\tau$ are depicted as light orange and light blue dots. Visually, it is evident that these dots have filled the gaps between parameters that were absent in our LAMMPS dataset.

Figure 10.

Predictions of the parameters (a) $\tau$ and (b) $\beta$ using the trained neural networks $NN_{\beta }$ and $NN_{\tau }$. The horizontal axes show $R_0$ and $\theta _{eq}$. The green and red circles correspond to the obtained values of $\tau$ and $\beta$ using the PINNs that were used to train $NN_{\beta }$ and $NN_{\tau }$. Additionally, the orange and blue dots represent the predicted values for grid interpolations from $R_0=0.136\ {\rm mm}$ to $R_0=0.204\ {\rm mm}$, and $\theta _{eq}=40^\circ$ to $\theta _{eq}=95^\circ$.

5. Symbolic regression

In this section, we use symbolic regression to find the explicit relation between the ODE parameters, and the initial radii and equilibrium contact angles. Symbolic regression is a technique used in empirical modelling to discover mathematical expressions or symbolic formulas that best fit a given dataset (Billard & Diday Reference Billard and Diday2002). The process of symbolic regression involves searching a space of mathematical expressions to find the equation that best fits the data. The search is typically guided by a fitness function that measures the goodness of fit between the equation and the data. The fitness function is optimized using various techniques, such as genetic algorithms, gradient descent, or other optimization algorithms. The equations discovered by symbolic regression are expressed in terms of familiar (i.e. more common) mathematical functions and variables, which can be understood and interpreted easily by humans.

In this study, we used the Python library gplearn for symbolic regression (Stephens Reference Stephens2016). As discussed in § 4.2, in order to have more accurate formulation for the ODE parameters before using symbolic regression, we trained two networks, $NN_{\beta }$ and $NN_{\tau }$, using the discovered values $\beta (R_0,\theta _{eq})$ and $\tau (R_0,\theta _{eq})$, enabling us to predict the parameter values for grid interpolations where no corresponding data points were available. The predicted parameters from both PINNs and the $NN_{\beta }$ and $NN_{\tau }$ networks are fed through the symbolic regression model to discover a mathematical formulation for each parameter. For this purpose, $\theta _{eq}$ and $R_0$ are fed as inputs, and $\eta$, $\beta$ and $\tau$ are the outputs. We set the population size to $5000$, and evolve $20$ generations until the error is close to $1\,\%$. Since the equation consists of basic operations such as addition, subtraction, multiplication and division, we do not require any custom functions.

The following results of symbolic regression can be substituted in (3.1):

(5.1)\begin{equation} \left.\begin{gathered} \eta =-0.255, \\ \beta = 0.283 + 0.27 \left(\frac{\theta_{eq}}{d}\right)\!, \\ \tau = 6.13 \left(\frac{d}{\theta_{eq}} + 1\right)\!, \end{gathered}\right\} \end{equation}

where $d$ is the initial size of the droplet in mDPD units.

Figure 11 shows the history of $\alpha (t, R_0, \theta _{eq})$ using the the ODE (3.1) with parameters from (5.1). The conversion between $d$ in mDPD units and $R_0$ in physical unit is $R_{0} = d\times 1.701\times 10^{-2}$ mm.

Figure 11.

The behaviour of $\alpha$ from the right-hand side of (3.1) with parameters from (5.1): (a) different contact angles with fixed initial radius $R_0=0.136$ mm; (b) varying initial radii with fixed contact angle $\theta _{eq}=77.9^\circ$.

Figure 12(a) depicts the values of $\alpha$ for initial drop size $R_0 = 0.127$ mm and contact angles $\theta _{eq}=93.4^\circ$ and $87.2^\circ$. Figure 12(b)compares the solution of (3.1), using the discovered parameters (5.1), with the simulation data. The figure demonstrates the agreement between the ODE solution and the actual simulation results for this particular, unseen dataset. It is important to note that this particular drop size lies outside the training interval for initial drop sizes, $[0.136, 0.204]$ mm.

Figure 12.

(a) The behaviour of the parameter $\alpha$ using (3.1) for $R_{0} = 0.127$ mm, which falls outside the range of the initial drop sizes used for training the networks. (b) The simulation data and the solution obtained from solving the ODE (3.1) with parameters from symbolic regression (5.1).

6. Bayesian physics-informed neural network results

Bayesian physics-informed neural networks integrate the traditional PINN framework with Bayesian neural networks (BNNs; Bykov et al. Reference Bykov, Höhne, Creosteanu, Müller, Klauschen, Nakajima and Kloft2021) to enable quantification of uncertainty in predictions (Yang et al. Reference Yang, Meng and Karniadakis2021). This framework combines the advantages of BNNs (Bishop Reference Bishop1997) and PINNs to address both forward and inverse nonlinear problems. By choosing a prior over the ODE and network parameters, and defining a likelihood function, one can find posterior distributions using Bayes’ theorem. The B-PINNs offer a robust approach for handling problems containing uncorrelated noise, and they provide aleatoric and epistemic uncertainty quantification on the parameters of neural networks and ODEs. The BNN component of the prior adopts Bayesian principles by assigning probability distributions to the weights and biases of the neural network. To account for noise in the data, we add noise to the likelihood function. By applying Bayes’ rule, we can estimate the posterior distribution of the model and the ODE parameters. This estimation process enables the propagation of uncertainty from the observed data to the predictions made by the model. We write (3.1) as

(6.1a)$$\begin{gather} \mathcal{N}_t(r; \boldsymbol{\lambda}) = f(t),\quad t \in \mathbb{R}^+, \end{gather}$$

(6.1b)$$\begin{gather}\mathcal{I} (r, \boldsymbol{\lambda}) = r_0, \quad t=0, \end{gather}$$

where $\boldsymbol {\lambda } = [\eta,\beta, \tau ]^\textrm {T}$ is a vector of the parameters of the ODE (3.1), $\mathcal {N}_t$ is a general differential operator, $f(t)$ is the forcing term, and $\mathcal {I}$ is the initial condition. This problem is an inverse problem; $\boldsymbol {\lambda }$ is inferred from the data with estimates on aleatoric and epistemic uncertainties. The likelihoods of simulation data and ODE parameters are given as

(6.2) \begin{equation} \left.\begin{gathered} P(\mathcal{D} \mid \boldsymbol{\theta}, \boldsymbol{\lambda}) = P(\mathcal{D}_{r} \mid \boldsymbol{\theta})\, P(\mathcal{D}_f \mid \boldsymbol{\theta}, \boldsymbol{\lambda})\, P(\mathcal{D}_{\mathcal{I}} \mid \boldsymbol{\theta}, \boldsymbol{\lambda}),\quad \text{where} \\ P(\mathcal{D}_{r} \mid \boldsymbol{\theta}, \boldsymbol{\lambda}) =\prod_{i=1}^{N_{r}} \frac{1}{\sqrt{2 {\rm \pi} \sigma_r^{(i)^2}}} \exp \Big[-\frac{({r}(\boldsymbol{t}_r^{(i)}; \boldsymbol{\theta}, \boldsymbol{\lambda})-\bar{r}^{(i)})^2}{2 \sigma_r^{(i)^2}}\Big]\!, \\ P(\mathcal{D}_f \mid \boldsymbol{\theta}, \boldsymbol{\lambda})=\prod_{i=1}^{N_f} \frac{1}{\sqrt{2 {\rm \pi} \sigma_f^{(i)^2}}} \exp \Big[-\frac{(\,f(\boldsymbol{t}_f^{(i)} ; \boldsymbol{\theta}, \boldsymbol{\lambda})-\bar{f}^{(i)})^2}{2 \sigma_f^{(i)^2}}\Big]\!, \\ P(\mathcal{D}_\mathcal{I} \mid \boldsymbol{\theta}, \boldsymbol{\lambda})= \prod_{i=1}^{N_\mathcal{I}} \frac{1}{\sqrt{2 {\rm \pi} \sigma_\mathcal{I}^{(i)^2}}} \exp \Big[-\frac{(\mathcal{I}(\boldsymbol{t}_i^{(i)} ; \boldsymbol{\theta}, \boldsymbol{\lambda})- \bar{\mathcal{I}}^{(i)})^2}{2 \sigma_{\mathcal{I}}^{(i)^2}}\Big]\!, \end{gathered}\right\} \end{equation}

where $D=D_r \cup D_f \cup D_{\mathcal {I}}$, with $\mathcal {D}_r=\{({\ln t}_r^{(i)}, {\ln r}^{(i)})\}_{i=1}^{N_r}$, $\mathcal {D}_f=\{(t_f^{(i)}, f^{(i)})\}_{i=1}^{N_f}$, $\mathcal {D}_\mathcal {I}= \{(t_\mathcal {I}^{(i)}, \mathcal {I}^{(i)})\}_{i=1}^{N_\mathcal {I}}$ are scattered noisy measurements. The joint posterior of $[\boldsymbol {\theta }, \boldsymbol {\lambda }]$ is given as

(6.3)\begin{align} P(\boldsymbol{\theta}, \boldsymbol{\lambda} \mid \mathcal{D}) &= \frac{P(\mathcal{D} \mid \boldsymbol{\theta}, \boldsymbol{\lambda})\, P(\boldsymbol{\theta}, \boldsymbol{\lambda})}{P(\mathcal{D})} \nonumber\\ &\approx P(\mathcal{D} \mid \boldsymbol{\theta}, \boldsymbol{\lambda})\, P(\boldsymbol{\theta}, \boldsymbol{\lambda}) \nonumber\\ &=P(\mathcal{D} \mid \boldsymbol{\theta}, \boldsymbol{\lambda})\, P(\boldsymbol{\theta})\,P(\boldsymbol{\lambda}). \end{align}

To sample the parameters from the the posterior probability distribution defined by (6.3), we utilized the Hamiltonian Monte Carlo approach (Radivojević & Akhmatskaya Reference Radivojević and Akhmatskaya2020), which is an efficient Markov Chain Monte Carlo (MCMC) method (Brooks Reference Brooks1998). For a detailed description of the method, see e.g. Neal (Reference Neal2011, Reference Neal2012) and Graves (Reference Graves2011). To sample the posterior probability distribution, however, variational inference (Blei, Kucukelbir & McAuliffe Reference Blei, Kucukelbir and McAuliffe2017) could also be used. In variational inference, the posterior density of the unknown parameter vector is approximated by another parametrized density function, which is restricted to a smaller family of distributions (Yang et al. Reference Yang, Meng and Karniadakis2021). To compute the uncertainty in the ODE parameters by using B-PINN, a noise of $5\,\%$ was added to the original dataset. The noise was sampled from a normal distribution with mean $0$ and standard deviation ${\pm }1$.

Here, the neural network model architecture comprises of two hidden layers, each containing 50 neurons. The network takes time $t$ as the input, and generates a droplet radius $r(t)$ as the output. Additionally, we include a total of 2000 burn-in samples.

The computational expense of B-PINNs compared to traditional neural networks arises primarily from the iterative nature of Bayesian inference and the need to sample from the posterior distribution. The B-PINNs involve iterative Bayesian inference, where the posterior distribution is updated iteratively based on observed data. This iterative process requires multiple iterations to converge to a stable solution, leading to increased computational cost compared to non-iterative methods. Moreover, B-PINNs employ sampling-based algorithms such as MCMC or variational inference to estimate the posterior distribution of the model parameters. These algorithms generate multiple samples from the posterior distribution, which are used to approximate uncertainty and infer calibrated parameters.

Sampling from the posterior distribution can be computationally expensive, particularly for high-dimensional parameter spaces or complex physics models. Furthermore, B-PINNs often require running multiple forward simulations of the physics-based model for different parameter samples. Each simulation represents a potential configuration of the model parameters. Since physics-based simulations can be computationally intensive, conducting multiple simulations significantly increases the computational cost of training B-PINNs. Achieving a high acceptance rate for posterior samples, especially for high-dimensional data, demands running a large number of simulations. This further adds to the computational complexity.

Due to the computational expense associated with B-PINNs, we opt to use the method selectively for a few cases only. Utilizing the insights gained from PINNs prediction and symbolic regression, specifically the known value of $\eta = -0.255$, we can leverage the power of B-PINNs to uncover and ascertain the values of the parameters $\beta$ and $\tau$. Figure 13 showcases the comparison between the mean values of the radii $\tilde {r}(t)$ predicted by B-PINNs represented by solid lines, the corresponding standard deviations denoted by highlighted regions, and the simulation data used for training presented as stars, while the test data are indicated by coloured circles. The horizontal axes represent time, while the vertical axes depict the spreading radii $\tilde {r}(t)$. This comparative analysis is conducted for two distinct initial drop sizes, namely $R_0 = 0.136$ and $0.170$ mm, considering various equilibrium contact angles.

Figure 13.

The mean and uncertainty (mean $\pm$2 standard deviations) of B-PINN predictions of the spreading radii history are given as solid lines and shaded regions, respectively. The test simulation data are depicted by solid circles, and training data are indicated by stars. This analysis is carried out for two different initial drop sizes, namely $R_0 = 0.136$ and $0.170$ mm, for three equilibrium contact angles.

Figure 14 illustrates the mean values of the parameters $\beta$ and $\tau$ obtained using B-PINNs along with their corresponding standard deviations. The solid lines represent the average values of the discovered parameters, while the highlighted regions indicate the standard deviations. The parameters discovered by PINNs are represented by the dashed lines. In figures 14(a,c), the vertical axes represent the values of $\beta$, while figures 14(b,d) display the values of $\tau$. The results are presented for two initial drop sizes, $R_0=0.136$ and $0.17$ mm. From the plots, it can be observed that the parameter $\beta$ exhibits a range of values between $1.0$ and $3.0$. On the other hand, the parameter $\tau$ fluctuates within the range $5.0$ to $7.0$. These ranges provide insight into the variability and uncertainty associated with the estimated values of $\beta$ and $\tau$ obtained through the B-PINN methodology.

Figure 14.

Comparison between B-PINNs and PINNs discovered parameters for a range of equilibrium contact angles and two initial radii. The mean values (solid lines) and the standard deviations (mean values $\pm$2 standard deviations, shaded region) of (a,c) $\beta$ and (b,d) $\tau$. The dashed lines represent the parameters discovered by PINNs.

By comparing figures 9 and 14, it becomes evident that the discovered parameters $\beta$ and $\tau$ using PINNs of B-PINNs frameworks exhibit remarkable similarity. This striking similarity reinforces the efficacy and capability of our models in identifying accurately the parameters of the ODE (3.1). The close alignment between the discovered parameters in both figures demonstrates the robustness and reliability of our models. It highlights their ability to capture effectively the underlying dynamics and characteristics of the spreading behaviourof CMAS, leading to accurate parameter estimation. This consistency and agreement between the PINN and B-PINN results provide further validation of the power and effectiveness of our modelling approaches in uncovering the true values of the parameters $\beta$ and $\tau$ in the ODE. Additionally, (3.1) with anticipated parameters obtained using B-PINNs is solved using Odeint. The results are presented in figure 15, which provides a comparison between the simulated spreading radii $r(t)$ (circles) and the solution of the ODE (solid lines). This comparison is conducted for initial drop sizes $R_0 = 0.136$ and $0.17$ mm, considering different contact angles.

Figure 15.

Comparison between the ODE solution with parameters found by B-PINNs (solid lines), and the simulation radii (circles). Two initial drop sizes, $R_0 = 0.136$ and $0.170$ mm, and three equilibrium contact angles are shown.