Particle–droplet coalescence and jumping on superhydrophobic surfaces—A direct numerical simulations study Open Access

In the field of surface technology, preserving the surface properties becomes a challenge due to contamination caused by solid particles and microbiological agents. These contaminants deteriorate the functionality of the surfaces that are used in numerous applications, such as defrosting, deicing, and heat-exchanging condensers.^1–3 To address this problem, surface modifications have been suggested, particularly focused on achieving hydrophobicity.^4,5 The design of such surfaces was initially inspired by the water-repellent behavior that certain organisms have developed in nature. These organisms demonstrate natural self-cleaning on their surfaces, with some examples being lotus plant leaves and cicada wings.⁶ They have developed micro- and nano-scale structures to repel water and remove contaminants by droplet motion.^7–10 Studies have shown that by reproducing patterns with micro- and nano-scale structures, contact angles of 170° and above were achieved, while the sliding angle and hysteresis remained below 2°.^11,12

A significant challenge related to a successful industrial application of these bio-mimicking superhydrophobic surfaces is how to improve their durability. Typical estimated lifetimes of these low-adhesion surfaces are often short,^3,13,14 but recent studies suggest that the life-cycle of superhydrophobic surfaces can be significantly extended.^9,12 This makes the need for understanding mechanisms and devising techniques for cleaning (including self-cleaning) of such surfaces even more industrially relevant.

Multiple ways to remove particles from surfaces have been suggested in the literature, with studies explaining the dynamics of removal when an initial velocity for a liquid phase is present¹⁵ or an air flow stream carries the particles away harnessing inertia.¹⁶ Gravity often contributes as an external supporting mechanism, leading droplets to roll off the substrate, while they capture and remove particles from the surface in the process.^4,17 In addition to the mentioned active cleaning mechanisms, there is a need for devising passive mechanisms that use droplets to remove the contaminants and require no presence of an external field.

A recently suggested and very popular passive mechanism for the removal of contaminants is a somewhat counter-intuitive phenomenon of droplet coalescence and jumping on superhydrophobic surfaces.^18,19 Wisdom et al.²⁰ were one of the first to explore the self-cleaning benefits of the phenomenon, but effort was foremost focused on the physics behind the eventual jumping.^19,21–23 The interest associated with this mechanism was mainly attributed to the absence of external contribution instigating the jumping,²¹ and the phenomenon was studied extensively, both experimentally and numerically, to uncover the factors that affect its efficiency,^24,25 influence of the initial velocity,^26,27 or the existence of grooves,²⁸ which all can significantly increase the jumping velocity. The surface energy released during coalescence is translated into the kinetic energy that propels the merged droplet. The total energy gained by the droplets can vary, and it is influenced by surface modifications that promote a Cassie–Baxter wetting stage.^29–31 In our previous study,³² we emphasized the strong connection between the values for advancing and receding contact angles and the energy losses. We also studied and quantified how hysteresis and the appearance of a dynamic contact angle variation increase dissipation. We have shown that the accurate capturing of the merged droplet contact line movement and the droplets' oscillations leads to a more accurate prediction of the energy losses at the contact line. We note that most of these studies have only focused on the coalescence and jumping of droplets, without considering the presence and influence of contaminants on the jumping process. The only work that we are aware of that studied numerically the removal of a spherical particle by binary droplet coalescence and jumping³³ used the lattice-Boltzmann method to model the fluids and a separate model for the capillary forces.³⁴ In that study, the particle was already encapsulated in the droplet, a condition that is not very physically relevant.

In a recent work, Yan et al.³⁵ reported a different passive mechanism for self-cleaning, based on coalescence (and subsequent jumping) of a single droplet and a particle. The authors studied the phenomenon experimentally and performed momentum and wetting analyses to identify the effects of particle/droplet properties and wettabilities. The process starts with the liquid bridge forming during the initial droplet spreading on the particle (referred to as the spreading liquid front) and interacting with the superhydrophobic surface and eventually pushing the particle–droplet system upward with a specific jumping velocity. The key factor behind this mechanism is the wettability of the particle (especially for hydrophilic particles), which reduces the overall surface energy of the system during the droplet spreading phase. A part of this energy is converted into the kinetic energy of the particle–droplet system. In some cases, the aggregate attains a jumping velocity similar to that of jumping droplets. Moreover, the particle–droplet system shows rotation due to the size and weight mismatch of the particle and the droplet, while an oblique upward direction may appear when the merged system is lifted from the solid surface.

This mentioned work³⁵ indeed provided significant information related to this novel passive self-cleaning mechanism, but certain details can be uncovered solely by numerical simulations that fully resolve the contact-line physics and capillary wave dynamics taking place during droplet–particle coalescence and jumping. A fundamental challenge, but also one of the major contributions of such numerical studies, is to establish the connection between adhesion of the particle to the droplet and inertia of the liquid bridge at the early spreading stage. In the mentioned numerical study,³³ the system was initiated with the particle partially wetted in equilibrium or fully encapsulated in a droplet. This leads us to conclude that there is great relevance in investigating the wetting of a particle surface from a fully dry state, and in that way identifying the driving factors for spreading dynamics, estimating the energy losses during spreading of the droplet on the particle, and studying the interaction of the oscillating droplet with the substrate. The present study provides qualitative and quantitative information about these characteristics and their influence on the jumping process by fully resolving the relevant physics and investigating different particle, droplet, or wetting properties. By creating a relevant parameter space (mirrored by the cases we will run), we will identify and quantify the influence of these parameters to the overall efficiency of the proposed mechanism.

The novelty and complexity of the task are highlighted by the absence of similar numerical studies. Numerous challenges exist in studying numerically the dynamics of complex interfaces, movement of a three-phase contact line, and wetting of a (highly) curved boundary such as that of the particle surface. The most important ones comprise the lack of knowledge related to the adhesion physics, particularly in the dynamic spreading during the liquid bridge formation and the limited (and insufficiently validated) models formulated to describe capillarity adhesion in Eulerian frameworks. The capillary forces on the particle from the droplet result from a line tension force acting at the three-phase contact line. The total of these forces combines the effects of the van der Waals bonds, electrostatic forces, and capillarity. The direction and magnitude of the resulting force are described and modeled for static applications such as colloids in rest,^36,37 captured contaminants,³⁸ or elastic solids.³⁹ However, we expect that the contact-line dynamics has a strong effect on the adhesive capillary force acting on the particle. So, any numerical framework that aspires to study this self-cleaning mechanism needs to successfully deal with three fundamental factors: (i) the contact line movement on the superhydrophobic surface, (ii) the spreading of the droplet on the curved particle surface, and (iii) the very small temporal and spatial scales determined by the fluid properties and droplet and/or particle sizes. In addition, it is not to be forgotten that the movement of the particle in space, as a result of the acting forces, creates an additional complexity since the numerical framework needs to account for the volume of the solid, which will be populated by the fluid and vice versa. To achieve these goals, we combine in this work an immersed boundary method (IBM)⁴⁰ with a sharp interface representation in an Eulerian framework.

Hence, this study addresses the aforementioned challenges and introduces a comprehensive and robust numerical method able to describe all stages of the particle–droplet coalescence and jumping process. The method comprises a combined volume of fluid (VOF)–IBM approach that is for the first time, to the best of the authors' knowledge, used for this type of problems. We also note that other studies used VOF-IBM solvers previously, for example, for a water flooding process,⁴¹ droplets and wetting cylinders,^42,43 droplet–particle interactions,⁴⁴ and wetting of non-spherical particles.⁴⁵ Specifically, we are particularly interested in accurately capturing the initial spreading phase for which only limited information can be extracted through experimental studies. We will improve the understanding of that phase of the process by looking in great detail at the liquid spreading on a curved surface, the partial contact of the droplet with the substrate, and the eventual launch of the particle–droplet system into the air. Furthermore, a comprehensive study will be conducted regarding the variations in the kinetic and surface energies of the system, along with an analysis of the total energy budgets for systematically varied droplet, particle, and surface properties. Our results will highlight, among other things, the key mechanisms of total viscous dissipation before and after particle–droplet departure.

The paper is organized as follows: in Sec. II, we describe in detail our numerical framework and how we configure our simulations and calculate fundamental geometrical and dynamical properties of the system (Sec. III). A thorough validation process is presented in Sec. IV. In Sec. V, we present and discuss our results and go through individual effects of the chosen geometrical and wetting properties of the system. Finally, in Sec. VI, we summarize our findings and implications of our study and present suggestions for possible directions of future work.

To counter issues stemming from the nature of co-located grids, we use the balanced-force method^52,53 in the momentum equation [Eq. (1)]. This method provides improved values for the velocities at the faces by accounting contributions from the forces and the pressure gradient. Crucially, it also includes the surface tension body force for the cells near the interface and helps to balance the variations in the velocity field due to the phase change. Moreover, it contains a time derivative in the velocity interpolation. This stable method for the face velocities reduces spurious currents and avoids oscillations caused by pressure instabilities. The correct implementation of this method has been considered key in achieving perfect advection of the sharp interface on top of curved surfaces. Our framework is able to capture the interface even for the smallest length scales that appeared in our study.

To account for the wetting of a particle, the contact angle is imposed at the particle surface in the same manner as we did for droplets in our previous works.^51,54 In summary, a Neumann boundary condition is implemented for the volume fraction in the vicinity of the immersed boundary. A correction for ghost cells near the contact line is used where the volume fraction is extrapolated to the ghost cell from the main field, in the tangential direction of the interface. We consider important to separate advancing and receding behaviors. We have tested here using both a constant contact angle for each wetting motion, and to let the framework compute the angles from a dynamic contact angle model, dependent on the contact line velocity.^55–57 A contact angle is always imposed with reference to the normal vector at the surface of the solid body. Additionally, the Navier-slip boundary condition is used for the superhydrophobic substrate, to account for roughness effects and to deal with the stress singularity problem.⁵⁴ 

Washino et al.⁵⁹ introduced the continuum capillary force method (CCF) to directly estimate the force exerted on a solid surface due to the presence of an interface of the fluid that spreads on that surface. It follows an approach similar to the CSF model in VOF,⁴⁹ which involves integrating the volume fraction gradient across a cell to locate the interface location and direction. With this technique to account for the existence of two interfaces, the method is able to estimate the location of the contact line from their intersection. In CCF, a field $χ$ denoting the solid volume fraction is generated, imposing a smooth transition between phases, in order to compute the Eulerian gradient of the solid volume fraction. The contact line is then obtained from the volume integral of two Dirac functions,^59,60 corresponding to the liquid–gas interface location with tangential direction to the solid and the solid–fluid boundary normal. Both are integrated over the volume of a cell to compute the line force in that cell.

We also needed to develop a method to estimate a possible instantaneous force using the analytical expression and considering the contact line to be fully circular, with a constant contact angle of $θeq,p$⁠. To compute the magnitude and direction of this analytical force, the mass center of the droplet is initially computed by our numerical framework by volume averaging. The unit vector $d̂$ between the two centers is selected to represent the direction of the force. The magnitude is obtained by computing the spreading angle $a$⁠, and subsequently the angle $ψ$⁠. To estimate the spreading angle, the contact line position is located in a plane created by $d̂$ and the z direction unit vector $eẑ$⁠. This was performed by evaluating in the investigated plane the liquid volume fraction $a$ in a certain radius $Rp$ from the particle center.

It is important to note that we have observed potential problems when using the CCF model. In particular, when the contact line is experiencing highly dynamic or volatile behavior, and with abrupt retractions and advancements of the contact line on the particle surface, the force $FCCF$ [Eq. (15)] was unable to capture the anticipated increased attractive forces. To assess the impact of that observation, we have chosen to compare the CCF model (⁠ $FCCF$⁠) with the analytical expression given in Eq. (17).

By looking at the directions of the two forces, a rather similar behavior was observed, with the CCF force expectedly showcasing a more volatile change, while the analytical formulation capturing the trend well. On the other hand, the magnitudes of the forces showed opposite behaviors with the analytical force attaining local maxima in the same instants that the force from the CCF model showcased local minima. The forces were normalized by the analytical capillary force for an ideal contact area between droplet and particle $Spl*$ and subsequently also compared with the magnitude of the total hydrodynamic forces acting on the particle [as computed with IBM from Eq. (8)]. The physical explanation for the different trends in the magnitude is at the moment unclear, and a parallel study is undergoing to understand this behavior. Nonetheless, it became apparent to us that the total integrated Laplace pressure force coincided with the weighted average of these two forces, performing exceptionally during all stages of the particle–droplet coalescence and jumping processes. This, and a successful validation against experimental results (shown below), made us comfortable to proceed with our investigation.

For the particle–fluid coupling that requires exchanging of information for the stresses to the rigid-body solver and the updated location and velocity of the particle to the fluid solver, the Newmark time-stepping method⁶¹ is selected.

We describe here how we formulate our simulation cases and obtain fundamental geometrical and dynamical properties of the droplet–particle system. In essence, the numerical setup is based on the experiments performed by Yan et al.³⁵ In that study, a spherical particle was positioned on top of a superhydrophobic surface with a droplet, in contact with the substrate, exhibiting spontaneous spreading over the particle. While that study explored particle–droplet coalescence and jumping using three types of particles of various densities and sizes, our investigation focuses on a single particle type that, in our belief, shows the most interesting features of this complex multiphase process. Therefore, we have chosen here to work with the stainless steel particle, which has the highest density of the three particle types in the experiments and is of sufficient size, which facilitates obtaining fine details from the available experimental images.

In all the simulations, the particle has a radius of $Rp=80$ μm and is meshed by a volume tetrahedron mesh. In the immersed boundary method (IBM), the surface triangulation is extracted and connected to the Eulerian mesh. The IBM can handle arbitrary shaped objects. We vary the droplet size and define the size ratio $kR$ as the radius of the particle to the radius of the droplet (⁠ $kR=Rp/Rd$⁠). The domain size and the minimum cell size are selected based on the droplet size, since the curvature of the droplet dictates the capillary-inertial scales of the problem. The size of the computational domain should be large enough to avoid any boundary effects on the droplet–particle system. To be specific, the domain dimensions are set to $10Rd×8Rd×6Rd$⁠, with the smallest grid size $Δx$ in the vicinity of the interface ranging from 20 to 40 cells per $Rd$⁠. Additionally, the time step has been adjusted to maintain a constant Courant number $CFL=||v|| Δt/Δx$ for each simulation. The velocity $v$ in the system is scaled by the capillary-inertial scaling $UCI=σRd/ρl$⁠. The simulation time t is normalized for the different setups to $τ$ by the capillary-inertial timescale (⁠ $τ=t/τCI$⁠), where $τCI$ is defined as $τCI=ρl/(σRd 3)$⁠.

To increase computational efficiency, adaptive mesh refinement (AMR) has been employed to maintain the desired refinement near the interface and at the location of the contact line. This allows using coarser resolutions in the regions distant from these critical areas.

The physical properties were defined for a temperature state of 20 °C with the density of air being $ρg=1.204$ kg/m³ and that of the liquid $ρl=998.21$ kg/m³. The dynamic viscosity of the liquid is set to $μl=100.16×10−5$ Pa s and for the air to $μg=1.813×10−5$ Pa s. Finally, the surface tension between air and water has been specified as $σ=7.286×10−2$ N/m.

The wetting properties of the substrate are first taken in accordance with the experiment. The advancing angle in a static test was recorded as $θadv,w=170.3°$ and the receding one as $θrec,w=167.7°$⁠. Under these conditions, our framework imposes these angles by first establishing from the contact line velocity whether we deal with an advancing or a receding motion. The particle, on the other hand, had a wider range in the measured static contact angles, with the values being between 35° and $76°$⁠. Therefore, we have chosen to start with a conservative approach and decided to have an equilibrium contact angle of $θeq,p=55°$ and a static hysteresis of $Δθp=5°$⁠. The advancing and receding angles provided to the system are, thus, $θadv,p=57.5°$ and $θrec,p=52.5°$⁠. In Table I, we summarize the changes that will be introduced as additional simulation cases and as compared to the discussed main configuration (case 1). We vary the droplet size, the wall and particle contact angles, the particle density and, finally, the droplet viscosity. The goal with having those nine additional cases has been to identify and quantify the relative influence of geometrical and wetting droplet, particle, and surface properties on the efficiency of the coalescence and jumping process. We particularly focus on the balances of different components of energy (both the energy beneficial for jumping and that representing losses). Such parameter-related information will be of great relevance when self-cleaning surfaces are to be developed.

According to the literature, the spreading factor $β$ remains constant in a logarithmic scale during the initial spreading phase, which suggests that the spreading dynamics is governed by the velocity of the contact line $UCI$⁠.^63–65 To obtain the spreading factor, a fitting of the arch length is performed in line with the following functional dependence on time $lcon,norm∼cτβ$⁠. Studies in this field^62,64 suggest that the spreading factor may vary from 0.25 to 0.66, with Yan et al.³⁵ estimating it being between 0.40 and 0.70 in the particle–droplet coalescence and jumping case.

The next step in our study is to provide evidence of the temporal and spatial convergence of our numerical framework, and we will start with the former. Within the context of the droplet–particle interaction, a numerical simulation is required to capture the transient spreading of the droplet on the particle. To this end, we selected three time step sizes and compared the resulting spreading dynamics. Throughout each of the simulations, the selected time step remained constant. We choose the initial time step on the basis that the Courant number (CFL) did not exceed 1, peaking slightly above 0.9. The time step was subsequently reduced twice, resulting in the simulations where the CFL remained below 0.5 and 0.25, respectively.

In Fig. 2, we measure the spreading of the droplet on the particle by calculating the arch length $lcon$ and comparing it to the ideal spreading length $lcon*$⁠. The study compares the spreading behavior in a logarithmic scale as a function of normalized time $τ$⁠. In agreement with prior studies on solid surfaces,^62,63,66 the spreading follows an exponential growth rate and is governed by inertia. We computed the spreading factor to be 0.366, with all the cases obtaining almost the identical behavior of the spreading length over time. It is in good agreement with the previous experimental results in Yan et al.,³⁵ regarding both the spreading time and the spreading factor. The jumping took place at the normalized time of $τ=3.25$ for all the cases, and it is indicated by an asterisk in the graph. The case with the largest time step with a maximum CFL of a bit more than 0.9 required a higher number of inner SIMPLEC iterations to converge to the residual criterion. For all the cases in this work that will serve for a parametric analysis, we have chosen the median of the three time-steps. The decision is based on our desire to maintain CFL values below 0.5 for the entirety of the simulation, which ensured a continuous behavior of the convective model (CICSAM), enhanced interface sharpness, and reduced the number of inner iterations for convergence.

A grid resolution investigation was also conducted to determine the minimum acceptable cell size $Δx$⁠, which allows physics-based interface movement and contact line spreading within the particle–droplet–substrate system. We examined three configurations for the cell length at the interface, corresponding to the grid sizes 20, 30, and 40 times smaller than the droplet diameter. The droplet and particle sizes align here with the video images of all stages of the process, provided by Yan et al.,³⁵ featuring a particle size of $Rp=80$ μm, the droplet size of $Rd=120$ μm, and having the particle–droplet size ratio of $kR=0.67$⁠. We have deliberately chosen this setup (with the particle having a relatively high curvature) in order to observe the scales dominated by the droplet pressure jump and the capillary-inertial regime. To this end, the curvature resolution for the particle corresponds to 13.3, 20, and 26.7 cells per particle radius.

To compare the three grid configurations, we examined the calculated different energy components and the obtained trajectories of centers of both particle and the droplet throughout the process. In Fig. 3, we show the energy in the z direction, as well as the translational and oscillatory energies for the three different grid sizes, having 20, 30, and 40 cells $/Rd$⁠. The moment of jumping of the particle–droplet system for each of the grid sizes is indicated with an asterisk. In addition, we tracked the centers of mass of both the particle and the droplet and plotted their locations in the x–z mid-plane of the system (shown in Fig. 4). We note that the center of the droplet was determined from the volume-averaged center of the liquid volume. The results demonstrate that jumping occurs in all three cases and showcase the converging upward motion for the 30 and 40 cell cases, in contrast to the 20 cells $/Rd$ case.

Figure 5 compares the forces acting on the particle. The hydrodynamic forces acting on the surface of the particle (mainly attributed to the higher inner pressure of the droplet on the wetted area, $Spl$⁠) are compared to the total contact line forces resulting from capillarity. For the three different grid resolutions, a very consistent trend is observed, even though the particle response to these forces exhibits a less converged trajectory. We note that the slight force mismatches between the cases with 30 and 40 cells per radius in Fig. 5 (blue and black line, respectively) will yield a minor inaccuracy for the particle location as the force integrates in time. Therefore, we have chosen the following cornerstones of the grid study, as seen in Fig. 3: the jumping time, the jumping velocity, and the upward energy (which divided by the available energy gives the total efficiency of the process), but also the behavior of the system that is clearly observed in the oscillational part of the kinetic energy. Of course, the displacement, which is a double temporal integration of the acceleration due to the forces acting on the particle, will highlight any small mismatches. We, however, argue that with obtaining so similar velocities for both x and z directions the grid validation satisfies the requirements to be able to carry out an extensive investigation of the particle jumping properties. As a result of this analysis, a configuration with 30 cells $/Rd$ is selected for moving forward in our parametric investigation as presented in Sec. V.

Finally, Fig. 6 shows a detailed qualitative comparison between our simulations and a sequence of video snapshots from the experiments by Yan et al.³⁵ of all stages of the coalescence and jumping process. Overall, we argue that the simulation accurately captures the spreading phase, the propagation of the capillary wave toward the particle, and the entire process of the apparent attraction between the particle and the droplet. We further observe that the oscillations of the droplet, resulting from the movement of the capillary wave, are accurately reproduced. As a result, the developed liquid bridge hits the surface causing a reaction force that is a main mechanism behind the jumping of both the particle and droplet. We acknowledge that the simulated departure shape from the substrate exhibits minimal differences when compared to the experimental images, with the oscillations of the droplet still fully agreeing throughout the process in both types of studies. It is crucial to observe that the moment of departure takes place at the expected instant and that the normalized jumping velocity was $vjump*=0.10$⁠, while Yan et al.³⁵ estimated it 0.075. After detachment, the droplet rotates with regard to the particle position (otherwise explained as sliding on the solid surface), with a similar rotational speed as in the experimental images. Finally, a pronounced disagreement between the two cases is observed during the later stages of the jumping process. We hypothesize that the simulation may incorporate a time-integrated error in the capillary forces, elevating the particle–droplet aggregate with a seemingly higher velocity. In summary, the presented validation makes us comfortable to continue using the developed numerical framework for an analysis of the relative importance of the individual particle and droplet properties on the entire process.

We will present and analyze here the results of our simulations by looking at the influence of a chosen system parameter on the outcome of the process. We start by looking at the overall energy budget in the system and continue our analysis with investigating the influence of the chosen particle, droplet, and substrate properties on the identified energy components. These properties can be divided into two groups, and this is how our analysis will be presented: (i) wetting properties of the substrate and the particle and (ii) the droplet and particle physical properties. The former group comprises the degree of hydrophobicity of the substrate, contact angles at the particle surface (including the influence of hysteresis) and, finally, the effect of implementing a dynamic contact angle model at both the substrate and the particle surface. In the latter group, we look at the influence of changing the droplet size, the particle density (having a significantly lighter particle), and the fluid viscosity.

Detailed information related to energy budget is essential when assessing the overall efficiency of the particle–droplet coalescence and jumping process. Such information is readily available from our numerical simulations. The available energy of the system stems from the excess surface energy before the droplet initiates spreading on the particle. The surface energy in the particle–droplet contact area minimizes the total energy of the system since our particle is hydrophilic. On the contrary, the surface energy for the area of the droplet in contact with the superhydrophobic substrate increases the total surface energy. When the droplet detaches, the total surface energy of the system is dependent on the particle–droplet area and the droplet–gas interface area. In Fig. 7, the available interface energy is presented for case 1 for all the interfaces in the system. The energy starts from the ideal interface energy measured for starting at rest up to the point when the system jumped and reached an ideal contact area between the particle and the droplet, as computed for the equilibrium contact angle. Moreover, the total kinetic energy of the system (computed from the velocity at each cell of the droplet and the kinetic energy of the particle) is presented. The mentioned plots in the figure are accompanied by those of the total translational kinetic energy computed from the separate translational kinetic energies of the two bodies, the kinetic energy of the two-body system in the vertical z direction, and the rotational kinetic energy of the system. We observe that the peak kinetic energy of the system reaches roughly $38%$ of the available energy, while the final gained momentum in the z direction corresponds to approximately $4%$ of the initial energy. Additionally, the rotation of the system gains energy in the period before detachment, and then, it decreases gradually. When the total kinetic energy of the system is subtracted from the total available energy in the system, we are able to estimate the fraction of the initial energy that is lost either by forms dissipation or by excess contact with the substrate. The latter can indeed be recovered by the system in the form of the kinetic energy, and therefore, the designation of this energy as lost is not necessarily appropriate. Still, we prefer using this term as compared to calling it the dissipated energy. We distinctly observe that the detachment of the droplet from the substrate causes a sudden decrease in the rate of energy being dissipated. After detachment, the energy is lost mostly due to viscosity in the liquid phase or due to contact line movement on the particle surface. Table II gives a summary of the results from the performed simulation cases (following the configurations of Table I), giving the jumping velocity, the total dissipated energy, and the efficiency of the process. An in-depth analysis of these results is included in Secs. V B and V C.

We analyze in this section the effects of different both particle and substrate wetting properties on the overall efficiency of the coalescence and jumping process. In all the cases, we will look at the behavior of different energy components, resulting from variations of the degree of substrate hydrophobicity and the value of the particle surface contact angle hysteresis.

We look here at the effect of the degree of the substrate hydrophobicity on the particle–droplet coalescence and jumping. Superhydrophobic surfaces with high water repellency typically exhibit contact angles of 170° and above, while minimizing pinning with a hysteresis of 3° or less. We compare the behavior of the high-performance substrate (case 1), with a second one having a reduced equilibrium contact angle and an increased hysteresis that corresponds to the configuration of case 4 (see Table II). The advancing and receding contact angles selected were 164° and 158°, respectively, presenting a hysteresis of 6°. Our selection is motivated by the minimum possible contact angle that is needed for this case, considering similar values for a classical droplet–droplet coalescence and jumping, that would result in decreased efficiency or even non-jumping. The comparison of the velocities of the two systems is depicted by the evolution of the kinetic energy along the vertical z direction, as shown in Fig. 8. A decrease in that energy is observed for the case with the less superhydrophobic substrate. Additionally, our analysis includes the total translational kinetic energy (adding the rotational energy of the particle) and the oscillating kinetic energy, calculated as the difference between the total kinetic energy accumulated in the system and the translational kinetic energy. Even though there is a minimal difference in the oscillational behavior for the droplet, the translational energy is higher for the high-performance substrate, while the particle–droplet system jumps from the surface at an earlier stage (as indicated by an asterisk). These results suggest that further investigation to obtain the minimum contact angle for jumping scenarios requires significant effort involving a wide parameter space, but that it also represents an intriguing topic for future investigations. The parameter space would also be strongly influenced by the particle–droplet size ratio as our previous experimental work suggested.³⁵ 

We know that motion in the system is initiated as the droplet starts spreading on the particle. A decrease in the contact angle with the particle increases the available wetting area on the particle surface. To examine this effect, we simulated two additional cases, with decreased (case 5) and increased (case 6) particle contact angles. From the analysis of the different kinetic energies in the system (Fig. 9), we see that for case 5 with the increased contact area, a higher vertical velocity upon jumping is exhibited. Additionally, a slight delay in reaction is recorded until the maximum spreading and jumping. The data, seen also in Table II, suggest that the non-normalized values for the total kinetic energy increase when the contact angle is reduced. This is expected due to the greater forces and the available wetting area. However, for the normalized kinetic energies presented in Fig. 9, we observe a slight decrease in the total translational (and rotational for the particle) kinetic energy for $θeq=45°$⁠. This effect is attributed to the increase in the available surface energy when the equilibrium contact angle on the particle decreases. Therefore, the rate of transformed energy into the total kinetic energy decreases when the particle has a larger contact area with the droplet. We consider this to be the result of the higher dissipation near the contact line, as both the area and the contact line increase in size. In contrast, the upward momentum of the system becomes higher, so even if a fraction of the energy is to be absorbed into the kinetic energy, there is still enough excess surface energy and a greater contact area with the substrate, causing the particle–droplet system to jump with more energy. We present in Fig. 10 the increase in the maximum contact area between the particle and the droplet due to the decreasing contact angles.

Next, we investigate the influence of different degrees of hysteresis on the particle surface, again focusing on the energy balance for the system. Emphasis is placed on quantifying the percentage of the energy that is lost due to dissipation when the droplet is experiencing a higher degree of hysteresis with the particle. Hysteresis, being defined as the difference between the angles observed in advancing and receding motions of the contact line, can be represented in the framework in two ways. The first one is by prescribing a constant degree of hysteresis or, alternatively, by varying the exact values of the contact angle imposed for the two different motions. The second way to estimate the influence of particle contact angles is by obtaining the effective contact angle from a dynamic contact-angle model. In our study, we utilize the Kistler dynamic model,^55–57 where the value for the contact angle is dependent on the contact line velocity. The model requires the advancing and receding contact angles as an input, as they are observed from a static experimental configuration. For a setup with a static hysteresis of 5°, equivalent to the low hysteresis configuration for the benchmark case (case 1), the effective hysteresis for the simulation that corresponds to case 8 is measured to be 15°. Therefore, the described case provides a suitable comparison with a setup that imposes a constant hysteresis of 15° and using the quasi-static model (case 7).

The kinetic energy in the vertical direction is decreased for the cases with a more pronounced hysteresis, as depicted in Fig. 11. Furthermore, the total translational kinetic energy in the system has also decreased, both for the higher hysteresis with the static angles case and the one with the dynamic contact-angle model. The percentage of the decrease is consistent for the two simulations, indicating that having a higher hysteresis, mostly due to a greater advancing angle, results in increased losses in the system after jumping. The jumping in the non-benchmark cases takes place moments after that in the benchmark case (case 1), as indicated by an asterisk in the respective cases. Case 7 with static hysteresis jumps after the benchmark case, with the Kistler case 8 detaching slightly afterward. Finally, the energy from the droplet oscillations appears to have slightly decreased for the dynamic model case, without affecting the kinetic energy saved in the system, while the static hysteresis case of 15° turns nearly identical to the benchmark case (case 1). Our analysis suggests that for the simulations that aim to more accurately estimate dissipation stemming from the particle surface hysteresis, precise values are required for the advancing and receding contact angles. Such values may be recorded in a dynamic test. For the cases with high effective hysteresis, using a dynamic contact-angle model can be a better option, especially if exact values of the angles cannot be provided. In this work, the benchmark case (case 1) has been designed having a low degree of hysteresis, since the information on the angles in the experimental tests varied significantly. Therefore, the effect of hysteresis is otherwise neglected in further investigations.

We will here investigate again the variation of different energy components, but this time as a function of changing the droplet and particle physical properties. We start with the effect of changing the droplet size and continue the analysis with working with a lighter particle and having a less viscous liquid phase in the system.

Experiments of Yan et al.³⁵ have shown that different types of particles required different particle–droplet size ratios to achieve efficient jumping and attain higher jumping velocities. For the stainless steel particle, we are working with in our simulations, and this ratio is $kr≃0.6$⁠, with the maximum normalized jumping velocity of 0.09.³⁵ We have chosen to test a larger droplet with the radius of 240 μm for case 2 (from the 120 μm radius of the benchmark case 1) and a smaller with an 85 μm radius defining case 3 in our study. Such a change influences the characteristic capillary-inertial velocity and time scales of the system. From the process efficiency table (Table II), it is observed that the normalized jumping velocity for the benchmark case (case 1) is 0.11, with its values being reduced to 0.8 and 0.5 for the larger (case 2) and the smaller droplet size (case 3), respectively. The kinetic energy in the z direction, as shown in Fig. 12, demonstrates a notable increase in the velocity in the early stage of coalescing for the larger droplet case. We attribute this effect to the stronger capillary forces attracting the particle toward the droplet with an energy component in the z direction, as well as due to the particle–droplet size mismatch. Since in case 2, the Laplace pressure decreases, the effect of the attracting force becomes stronger during the spreading phase. This energy is primarily transferred to the particle, with the droplet displacing at a later stage. Additionally, the increased inertia forces are evident from the oscillatory kinetic energy of the system, as the oscillations in the droplet are not dampened for the duration of the simulation. The total kinetic energy of the system displays a volatile behavior, with dissimilar time scales for maximum spreading and jumping. Furthermore, the particle–droplet system experiences in this case only a slight elevation and retains a minimal upward velocity. On the other hand, the smaller droplet (case 3) retains the same trends for jumping as the benchmark case (case 1), albeit with a significantly reduced jumping velocity. The oscillations are now dampened earlier due to the higher dissipation present in that configuration. Notably, the efficiency of the translational kinetic energy decreases for both the smaller and the larger droplet, aligning with the experimental findings from Yan et al.³⁵ 

The particle density significantly influences the particle response to the forces acting on it, since the density proportionally changes the momentum of the particle. We test here the behavior of a lighter particle with the density of 2530 kg/m³ (case 9), corresponding to the soda-lime glass particle of the experimental study by Yan et al.³⁵ We keep the remaining parameters (contact angles and size) consistent with the benchmark case (case 1). Figure 13 displays the normalized kinetic energies of the system. The results suggest that the contact angles and the droplet size are indeed the primary factors influencing the oscillations in the droplet. The rationale behind this observation is that we could not observe significant differences between cases 1 and 9. However, less total translational energy is observed in the system in the latter case, clearly due to the decreased mass of the particle. The main conclusion from these results is that the peak upward velocity for the particle–droplet system increases significantly, achieving more than a doubled efficiency for the jumping process compared with the benchmark case (exceeding $8%$⁠, see case 9 in Table II).

Having noted variations in the dissipation with changes in the droplet size and capillary-inertial scales, we investigated how the change of the liquid-phase viscosity affects the particle–droplet coalescence and jumping. We now define case 10 by reducing this viscosity to one tenth of that in the benchmark case (case 1), which corresponds to a decrease in the characteristic Ohnesorge number from 0.0107 to 0.001 07. Although the maximum spreading area remains constant, oscillations in the contact line persist longer. Figure 14 compares the kinetic energies of the system. As expected, the wave-like behavior of the translational kinetic energy persists with a higher amplitude and less energy losses over time compared to those in case 1. This outcome clearly suggests that there is less viscous dissipation in the later stage of the process, where the droplet detaches from the substrate. The vertical velocity is also increased throughout the simulation, indicating a minimized viscous dissipation at the particle surface and near the contact line. Moreover, the higher inertia in the system causes the energy stored by oscillations from the droplet to remain constant during the airborne stage of the simulation.

We have performed a comprehensive multiphase direct numerical study (using a combined VOF—immersed boundary method) of coalescence and subsequent jumping from a superhydrophobic surface of a particle–single droplet system. Our simulations reveal fundamental dynamics of all relevant stages of that passive self-cleaning mechanism for the removal of contaminants from surfaces in both nature and technological applications: (i) initial spreading of the droplet on the particle (and development of a liquid bridge), (ii) oscillation of the droplet normal to the surface and the impact of the liquid bridge to the surface, and (iii) departure of the particle–droplet system from the surface. A thorough validation demonstrated temporal and spatial convergence of the developed numerical framework. A fully resolved three-dimensional simulation was compared with the previously published experimental images,³⁵ with great agreement in capturing the overall droplet–particle behavior and good conformity related to the removal of the particle from the substrate due to the capillary forces acting at the three-phase contact line. We have implemented a model for the estimation of capillary forces from a continuum perspective⁵⁹ and introduced our suggestions on how to reduce oscillations of the adhesion forces in a multiphase direct numerical simulations (DNS) framework.

We have systematically investigated particle-, droplet- and surface wetting, and geometrical properties and their influence on the coalescence and jumping processes. In the chosen parameter space, we have quantified energy efficiencies and dissipation, information that is not straightforward to elucidate from experiments. We have shown that pronounced superhydrophobicity enhances the efficiency of the jumping process, while the average contact angle at the particle surface is of great significance for the available energy and the dynamics of the process. We demonstrate an increase in dissipation when both higher degrees of hysteresis at the particle surface are present and when the viscosity of the droplet is increased. The simulations prove that the size of the droplet is an important parameter, with an increasing dissipation when the size of the droplet tends to that of the particle.

Our particle properties (wettability, size, and density) are representative of typical deposits on many outdoor surfaces in nature and industrial applications. The findings of this study thus provide guidelines for the design of functional surfaces for various technological applications.

The authors would like to acknowledge the Swedish Research Council for the financial support of this project (Vetenskapsrådet, No. 2019-04969). In addition, we would like to acknowledge that the handling of data and other computations were enabled by resources provided by the National Academic Infrastructure for Supercomputing in Sweden (NAISS), partially funded by the Swedish Research Council through Grant Agreement No. 2022-06725.

The authors have no conflicts to disclose.

K. Konstantinidis: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Methodology (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). J. Göhl: Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Visualization (supporting); Writing – review & editing (equal). A. Mark: Formal analysis (supporting); Funding acquisition (supporting); Investigation (equal); Methodology (equal); Resources (equal); Software (lead); Supervision (equal); Validation (equal); Visualization (supporting); Writing – review & editing (equal). X. Yan: Formal analysis (equal); Investigation (equal); Resources (equal); Validation (equal); Visualization (supporting); Writing – review & editing (equal). N. Miljkovic: Formal analysis (equal); Investigation (equal); Resources (equal); Validation (equal); Visualization (supporting); Writing – review & editing (equal). S. Sasic: Conceptualization (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon request.