Droplet coalescence kinetics: Coalescence mechanisms and thermodynamic non-equilibrium effects under isothermal and non-isothermal conditions

Abstract

Droplet coalescence is a multiscale phenomenon driven by multiple thermodynamic non-equilibrium (TNE) effects. This study investigates the coalescence mechanisms and the complex interplay between various TNE effects in both isothermal and non-isothermal cases kinetically. The main findings include: (i) Coalescence initiation and cut-through mechanisms: In the non-isothermal case, latent heat release results in a temperature rise that slightly increases the surface tension gradient (driving force) near the contact point, while substantially enhancing the pressure gradient (resistance). This leads to a significantly prolonged initiation of coalescence compared to the isothermal case. The additional temperature field effects introduce greater spatial complexity to both the flow field and morphology. For the cut-through mechanism, in both cases, pressure extends the liquid-vapor interface in opposite directions, promoting liquid bridge growth. (ii) TNE effects: The latent heat-induced temperature rise in the thermal case considerably dampens the intensity of TNE effects. Before and after droplet contact, TNE quantities driven by the temperature gradient ({\mathbf{\Delta }}_{3,1}^{\ast } and {\mathbf{\Delta }}_3^{\ast }) and those driven by the velocity gradient ({\mathbf{\Delta }}_2^{\ast } and {\mathbf{\Delta }}_{4,2}^{\ast }), alternate in dominating the coalescence process. This alternating dominance results in a more complex spatiotemporal evolution of TNE effects compared to the isothermal case, where the TNE intensity is dominated by {\mathbf{\Delta }}_2^{\ast } and {\mathbf{\Delta }}_{4,2}^{\ast }. (iii) Entropy production mechanisms: In the non-isothermal case, entropy production is contributed by both {\mathbf{\Delta }}_2^{\ast } and {\mathbf{\Delta }}_{3,1}^{\ast }, with the former being the dominant contributor. The temperature field reduces the entropy production rate, but extends its duration, thereby increasing the total entropy production. The entropy production rates from {\mathbf{\Delta }}_2^{\ast } and {\mathbf{\Delta }}_{3,1}^{\ast } alternate before and after the droplet cut-through. Further investigation into effects of the initial droplet distance r_d revealed that it increases droplet cut-through time t_c, the time when total TNE reaches its maximum t_{D\mathrm{m}\mathrm{a}\mathrm{x}}, and the time when total entropy production rate reaches its maximum t_{\dot{S}\mathrm{p}\mathrm{r}\mathrm{m}\mathrm{a}\mathrm{x}} in a power-law manner: (t_c,t_{D\mathrm{m}\mathrm{a}\mathrm{x}},t_{\dot{S}\mathrm{p}\mathrm{r}\mathrm{m}\mathrm{a}\mathrm{x}})\sim {r_d}^{\alpha }. This research provides kinetic insights into droplet coalescence, offering mesoscopic support for dynamic, cross-scale regulation and multifunctional integration of coalescence processes in industrial applications.

Graphical abstract

Keywords

droplet coalescence / discrete Boltzmann method / non-equilibrium effect / entropy production

1 Introduction

Droplet coalescence is a widespread phenomenon in fields such as nature, biomedical sciences, microfluidics, engine combustion, multi-medium vehicles, and photolithography [1–8]. In the field of microfluidics, droplet-based technologies involving deformation, fragmentation, coalescence, and collision, hold great potential for applications in medical diagnostics, food science, drug development, and artificial biomimetic systems manufacturing [1, 2, 8]. Droplets, serve as multifunctional carriers for controlled reactions, mixing, and analysis, play a key role in these processes [1, 2, 9]. In combustion systems with dense spray devices, intense interactions between droplets-including atomization, collision, rupture, coalescence, and oscillation-greatly affect fuel heat and mass transfer, combustion efficiency, and combustion stability [3, 4]. In explosion dynamics, studying droplet deformation, collision, and coalescence under shock waves at various Mach numbers provides crucial insights for mitigating shockwave energy and reducing shock damage [5].

The physical mechanisms and evolutionary patterns of droplet collision, coalescence, and fragmentation have been extensively studied under isothermal conditions. Surface tension is the primary driving force for the spontaneous coalescence of static droplets, which reduces the system’s free energy by minimizing the total inter-droplet surface area. In the initial stage of static droplet coalescence, the growth of the liquid bridge radius r_b over time is directly related to the system’s viscosity [10–17]. The shape of the droplet after collision, coalescence, and fragmentation is closely linked to characteristic parameters such as the Bond number (Bo), Weber number (We), Ohnesorge number (Oh), and the initial radius ratio of the droplets [18–31]. The number of droplets, the gas-liquid density ratio of the system, and the presence of chemical reactions between droplets also influence the morphological characteristics and evolutionary patterns of droplet coalescence [32, 33].

The main difference between thermal and isothermal multiphase flow systems is the heat generation and conduction during phase separation, resulting in lower interfacial stresses and localized temperature changes due to latent heat release, which leads to new local mechanical equilibrium [34]. Under non-isothermal conditions, the variety of driving mechanisms increases significantly, and the competition and coupling between these mechanisms intensify, resulting in a more complex droplet coalescence process. Specifically, because surface tension varies sensitively with temperature, a gradient forms across the droplet surface, propelling the droplet along the temperature gradient — a process known as the Marangoni effect. This effect significantly influences coalescence dynamics [35–38]. The local temperature difference resulting from the gradient also alters the liquid’s viscosity, enhancing both heat and mass transfer, and the coupling between them. These changes modify the internal convection patterns and interface behavior of the droplet, subsequently affecting its deformation velocity and coalescence dynamics [33, 35]. Currently, the kinetic study of droplet collision, coalescence, and fragmentation under non-isothermal conditions remains limited. Non-isothermal conditions, however, are prevalent in combustion systems [3, 4, 33, 35], explosion dynamics [5], nuclear reactors [39], and microfluidics [40, 41]. For instance, in microfluidics, the thermal Marangoni effect is often used to control the droplet coalescence process [40, 41]. This highlights the importance of comparative studies on the driving mechanisms and underlying kinetic behaviors of droplet coalescence under isothermal and non-isothermal conditions.

Numerical simulations offer several advantages, such as high efficiency, low cost, strong repeatability, flexible control of variables, and the capacity to handle complex multi-scale systems. These qualities make them indispensable tools for studying droplet coalescence. Widely employed macroscopic simulation methods, based on the Navier−Stokes (NS) equations [30, 31, 42–49], and mesoscopic methods, such as the lattice Boltzmann Method (LBM) [50–65] and discrete unified gas kinetic scheme [66–68], have made significant strides in advancing our understanding of this phenomenon. Despite significant progress, several aspects of droplet coalescence still require further development.

(i) Neglect of temperature field and compressibility effects: Most studies currently rely on isothermal, incompressible models, neglecting the impact of temperature fields and compressible convection on flow behavior. As a result, these models fail to fully capture the intricate interactions between droplets in real-world multiphase flows.

(ii) Oversight of non-equilibrium behaviors and effects: While most studies focus on the evolution of conserved quantities, such as mass, momentum, and energy, alongside slow-varying variables, they often overlook the evolution of non-conserved, fast-varying variables like thermodynamic non-equilibrium (TNE) and hydrodynamic non-equilibrium (HNE) behaviors. These non-equilibrium phenomena are closely tied to small-scale structures and fast modes, such as liquid bridges, droplet boundary layers, vortices, droplet fragments, splash phenomena, and surface fluctuations. In regions near these structures, discrete effects become pronounced, challenging the reliability of the methods discussed above. Furthermore, during rapid collisions and coalescence, the onset of fast modes impedes the system’s return to thermodynamic equilibrium in a timely manner, leading to strong non-equilibrium effects that further undermine the validity of these approaches. On the mesoscopic scale, HNE and TNE arise from the deviation of the distribution function from equilibrium, highlighting system complexity. On the macroscopic scale, these non-equilibrium effects are intricately coupled with macroscopic quantities and exert a fundamental influence on the overall behavior of the system.

(iii) Limitations in data analysis and information extraction: In the droplet coalescence process, various driving mechanisms, including TNE effects, surface tension, temperature gradients, and velocity gradients, interact to produce complex, multi-scale spatiotemporal structures. To achieve a deeper understanding of these processes, a multi-perspective analysis is essential.

(iv) Lack of systematic comparison of entropy production mechanisms: Entropy production mechanisms are pivotal in describing the thermodynamic behavior during droplet coalescence, particularly with regard to energy conversion and dissipation in non-equilibrium states. Despite their importance, existing studies lack a systematic comparison of entropy production mechanisms, especially in terms of the contributions of different driving forces, such as surface tension, heat flux, and viscous stress, to the entropy production rate. This gap limits our understanding of entropy production and non-equilibrium effects during droplet coalescence and hinders the broader applicability of existing models in complex systems.

To characterize the rich and complex HNE and TNE effects during droplet merging, and to clarify the driving mechanisms and entropy production processes of droplet coalescence, we resort to the discrete Boltzmann method (DBM) based on kinetic and mean-field theories. DBM is a non-equilibrium statistical physics modeling and complex physical field analysis method based on the discrete Boltzmann equation [69, 70]. Compared to macroscopic fluid modeling, DBM provides more kinetic information, offering new physical insights for the study of complex systems; and compared to microscopic molecular dynamics, it can simulate and investigate kinetic processes over larger spacetime scales. The core advantage of DBM lies in dynamically and intuitively describing and displaying discrete/TNE effects and behaviors. It addresses two major challenges in numerical experimental research: “how to model before simulation” and “how to analyze after simulation” for complex flows, especially those with significant discrete state and effects as well as non-equilibrium state and effects [70].

Recent studies have emphasized the importance of investigating TNE effects and behaviors. In scenarios where TNE occurs, neglecting its effects can lead to inaccurate estimations of viscosity, heat flux, and ultimately result in imprecise calculations of macroscopic quantities [71–73]. In the field of complex multiphase flows, such as phase separation [74, 75], bubble/droplet coalescence [28, 29, 76], and laser-bubble/droplet interactions [77], TNE effects offer significant advantages in characterizing the evolutionary features of scaled systems. For example, during liquid-vapor phase separation, the maximum points of the total TNE intensity and entropy production rate can be used to distinguish between the spinodal decomposition and domain growth stages [74, 75]. The non-organized momentum flux component in the coalescing direction aids in identifying stages in bubble coalescence [76]. Additionally, the intensity of non-organized momentum flux (NOMF) can provide a standard for distinguishing different stages of the droplet collision process and for differentiating between types of droplet collisions [28]. In fluid instability systems, by detecting TNE features from different perspectives, a more detailed physical structure of the interface can be provided, which can then be used for the physical identification and tracking technology design of different interfaces within the system [78–82].

It is important to note that the droplet coalescence process differs significantly under isothermal and non-isothermal conditions. This difference arises from the distinct coalescence mechanisms and the underlying TNE effects. Understanding the coalescence mechanisms under both isothermal and non-isothermal conditions is crucial for effective droplet manipulation in various scenarios. This paper compares and analyzes the differences and connections between isothermal and non-isothermal droplet coalescence processes, focusing on coalescence mechanisms, morphological features, flow field characteristics, non-equilibrium effects, and entropy production mechanism. The structure of this paper is as follows: Section 2 introduces the DBM method; Section 3 presents simulation results and an analysis of droplet coalescence kinetics under both isothermal and non-isothermal conditions; Section 4 summarizes the findings and outlines future research directions.

2 The model

The DBM is a kinetic modeling and complex physical field analysis approach developed based on discrete Boltzmann equation:

(1){\displaystyle \frac{\mathrm{\partial }{f}_{ki}}{\mathrm{\partial }t}}+{\mathit{v}}_{ki}\cdot {\displaystyle \frac{\mathrm{\partial }{f}_{ki}}{\mathrm{\partial }\mathit{r}}}=-{\displaystyle \frac{1}{\tau }}(f_{ki}-f_{ki}^{eq})+I_{ki},

where f_{ki}^{eq} is the discrete form of the local equilibrium distribution function [83], with k representing the k-th group of particle velocities {\mathit{v}}_{ki}, and i the direction of {\mathit{v}}_{ki}. The external force term, I_{ki} =-[A+\mathit{B}\cdot ({\mathit{v}}_{ki} -\mathit{u})+(C +C_q)({\mathit{v}}_{ki}-\mathit{u}{)}^2]f_{ki}^{eq}, is introduced using a mean-field approximation [71, 84] to account for the contribution of intermolecular interaction forces to the rate of change of the distribution function. Here, A=-2(C+C_q )T ensures mass conservation. \mathit{B}=\frac{1}{\rho T}\mathbf{\nabla }\cdot \left[(P-\rho T)\mathit{I}+\mathbf{\Lambda }\right] with \mathbf{\Lambda }=K\mathbf{\nabla }\rho \mathbf{\nabla }\rho -K(\rho {\mathbf{\nabla }}^2\rho +\frac{1}{2}|\mathbf{\nabla }\rho {|}^2)\mathit{I}-[\rho T\mathbf{\nabla }\rho \cdot \mathbf{\nabla }(K /T )]\mathit{I} representing the density gradient contribution to the pressure tensor, \mathit{I} is the unit tensor and K is the surface tension coefficient. The term C=\frac{1}{2\rho T^2}\{ \left(P-\rho T\right)\mathbf{\nabla } \cdot \mathit{u}+\mathbf{\Lambda }: \mathbf{\nabla }\mathit{u}+a{\rho }^2\mathbf{\nabla }\cdot \mathit{u}-K[\frac{1}{2}\mathbf{\nabla } \rho \cdot \mathbf{\nabla }\rho \mathbf{\nabla } \cdot \mathit{u}+\rho \mathbf{\nabla }\rho \cdot \mathbf{\nabla }\left(\mathbf{\nabla }\cdot \mathit{u}\right)+\mathbf{\nabla } \rho \cdot \mathbf{\nabla }\mathit{u}\cdot \mathbf{\nabla } \rho ]\} represents a partial contribution to energy. The term C_q=\frac{1}{\rho T^2}\mathbf{\nabla }\cdot (q\rho T \mathbf{\nabla }T) accounts for additional heat conduction within the system. The Carnahan-Starling equation of state (EOS) is given by P=\rho T\frac{1+\eta +{\eta }^2-{\eta }^3}{(1-\eta {)}^3}-a{\rho }^2, with \eta = \frac{b\rho }{4}, and a and b being the attraction and repulsion parameters, respectively. Here, \rho, \mathit{u}, and T represent the local density, velocity, and temperature, respectively.

As a complex physical field analysis method, DBM extracts non-equilibrium information from the system using the high-order non-conserved moment:

(2)\begin{array}{rl}{\mathbf{\Delta }}_{m,n}=& {\mathit{M}}_{m,n}(f_{ki}-f_{ki}^{eq})\\ =& \sum _{ki}{\left(\frac{1}{2}\right)}^{1-{\delta }_{m,n}}(f_{ki}-f_{ki}^{eq})({\mathit{v}}_{ki}\cdot {\mathit{v}}_{ki}{)}^{\frac{m-n}{2}}{\mathit{v}}_{ki}^n,\end{array}

where {\mathbf{\Delta }}_{m,n} describes non-equilibrium effects, including both molecular average behavior and thermal fluctuation behavior — referred to as thermo-hydrodynamic non-equilibrium (THNE) effects. When {\mathit{v}}_{ki} is replaced by {\mathit{v}}_{ki}^{\ast }={\mathit{v}}_{ki}-\mathit{u}, {\mathbf{\Delta }}_{m,n} becomes {\mathbf{\Delta }}_{m,n}^{\ast }, which exclusively represents TNE effects. The difference {\mathbf{\Delta }}_{m,n}-{\mathbf{\Delta }}_{m,n}^{\ast } represents hydrodynamic non-equilibrium (HNE) effects.

Among the TNE measures, {\mathbf{\Delta }}_2^{\ast } represents the NOMF or viscous stress, with its first-order analytical form: {\mathbf{\Delta }}_2^{\ast (1)}=-\rho T\tau \left[\mathbf{\nabla }\mathit{u}+(\mathbf{\nabla }\mathit{u}{)}^{\mathrm{T}}-\mathit{I}\mathbf{\nabla }\cdot \mathit{u}\right]. The term {\mathbf{\Delta }}_{3,1}^{\ast } represents the non-organized energy flux (NOEF) or heat flux, with its first-order analytical expression: {\mathbf{\Delta }}_{3,1}^{\ast (1)}=-2\rho T\tau \mathbf{\nabla }T. The term {\mathbf{\Delta }}_3^{\ast } represents the viscous stress flux, and {\mathbf{\Delta }}_{4,2}^{\ast } represents the flux of heat flux. To provide a comprehensive description of the system’s non-equilibrium intensity, DBM adopts the non-equilibrium intensity vector: {\mathit{S}}_{\mathrm{T}\mathrm{N}\mathrm{E}}= \{| {\mathbf{\Delta }}_2^{\ast } | ,|{\mathbf{\Delta }}_{3,1}^{\ast }|,|{\mathbf{\Delta }}_3^{\ast }|,|{\mathbf{\Delta }}_{4,2}^{\ast }|,\dots ,|{\mathbf{\Delta }}_{m,n}^{\ast }|,D^{\ast }\}, with D{}^{\ast }=\sum _{m,n}|{\mathbf{\Delta }}_{m,n}^{\ast }| the total TNE intensity. The definitions, components, and their physical meanings are displayed in the Appendix.

Taking the moments of Eq. (1) with the collision invariant vector (1, {\mathit{v}}_{ki}, v_{ki}^2/2), the generalized hydrodynamic equations (GHEs) for non-ideal fluids considering surface tension effect are obtained:

(3){\displaystyle \frac{\mathrm{\partial }\rho }{\mathrm{\partial }t}}+\mathbf{\nabla } \cdot (\rho \mathit{u})=0 ,

(4){\displaystyle \frac{\mathrm{\partial }(\rho \mathit{u})}{\mathrm{\partial }t}}+\mathbf{\nabla }\cdot (\rho \mathit{u}\mathit{u}+P \mathit{I})+\mathbf{\nabla } \cdot (\mathbf{\Lambda }+{\mathbf{\Delta }}_2^{\ast })=0,

(5)\begin{array}{rl}\frac{\mathrm{\partial }{e}_T}{\mathrm{\partial }t}+& \mathbf{\nabla }\cdot (e_T\mathit{u}+P\mathit{u})\\ +& \mathbf{\nabla }\cdot \left[(\mathbf{\Lambda }+{\mathbf{\Delta }}_2^{\ast })\cdot \mathit{u}+{\mathbf{\Delta }}_{3,1}^{\ast }+2\rho Tq\mathbf{\nabla }T\right]=0,\end{array}

with e_T=\rho T-a{\rho }^2+K{\left|\mathbf{\nabla }\rho \right|}^2/2+\rho u^2/2 the total energy density. Recovering the GHEs represents only one aspect of the physical capabilities of the DBM. In fact, the full physical description provided by DBM corresponds to the extended hydrodynamic equations [71], which encompass not only the evolution equations for conserved moments, but also for non-conserved moments closely associated with thermodynamic and hydrodynamic non-equilibrium effects. The significance of these extended equations increases notably as the discrete effects or non-equilibrium intensity becomes stronger.

Both types of non-equilibrium effects, NOMF and NOEF contribute to the entropy production rate,

(6){\dot{S}}_{\mathrm{p}\mathrm{r}}= {\dot{S}}_{\mathrm{N}\mathrm{O}\mathrm{M}\mathrm{F}}+ {\dot{S}}_{\mathrm{N}\mathrm{O}\mathrm{E}\mathrm{F}},

where

(7){\dot{S}}_{\mathrm{N}\mathrm{O}\mathrm{M}\mathrm{F}}= {\displaystyle \int \left(-{\displaystyle \frac{1}{T}}{\mathbf{\Delta }}_2^{\ast }:\mathbf{\nabla }\mathit{u}\right)\mathrm{d}\mathit{r},}

and

(8){\dot{S}}_{\mathrm{N}\mathrm{O}\mathrm{E}\mathrm{F}}= {\displaystyle \int \left({\mathbf{\Delta }}_{3,1}^{\ast }+2\rho Tq\mathbf{\nabla }T\right)\cdot \mathbf{\nabla }\left({\displaystyle \frac{1}{T}}\right)\mathrm{d}\mathit{r}.}

3 Simulations and results

Droplet coalescence is a multiscale phenomenon driven by the interplay of multiple physical mechanisms, particularly the competing TNE effects. This study systematically explores these mechanisms by examining: (i) slowly varying macroscopic variables, (ii) interface dynamics, (iii) rapidly varying non-equilibrium quantities, and (iv) entropy production rates under both isothermal and non-isothermal conditions.

All parameters in the manuscript and simulations are dimensionless. The reference variables can be chosen as the reference density {\rho }_{\mathrm{r}}, reference temperature T_{\mathrm{r}}, and the global size of the domain L_{\mathrm{r}}. The choice of the reference quantity is arbitrary; for instance, the critical density {\rho }_c and the critical temperature T_c can be selected as reference variables, calculated from the non-ideal equation of state. The relationship between the parameters in DBM and the physical ones are as follows:

(9)\hat{\rho }={\displaystyle \frac{\rho }{{\rho }_{\mathrm{r}}}}, \hat{T}= {\displaystyle \frac{T}{T_{\mathrm{r}}}}, {\hat{r}}_{\alpha }={\displaystyle \frac{r_{\alpha }}{L_{\mathrm{r}}}}, {\hat{u}}_{\alpha }={\displaystyle \frac{u_{\alpha }}{\sqrt{R{T}_{\mathrm{r}}}}},

(10)\left(\hat{t},\hat{\tau }\right)={\displaystyle \frac{\left(t,\tau \right)}{L_{\mathrm{r}}/\sqrt{R{T}_{\mathrm{r}}}}},{\hat{c}}_p={\displaystyle \frac{c_p}{R}},\hat{\sigma }={\displaystyle \frac{\sigma }{{\rho }_{\mathrm{r}}R{T}_{\mathrm{r}}{L}_{\mathrm{r}}}},

(11)\hat{\mu }={\displaystyle \frac{\mu }{{\rho }_{\mathrm{r}}{L}_{\mathrm{r}}\sqrt{R{T}_{\mathrm{r}}}}},\hat{\kappa }= {\displaystyle \frac{\kappa }{{\rho }_{\mathrm{r}}{L}_{\mathrm{r}}R\sqrt{R{T}_{\mathrm{r}}}}},

where r_{\alpha } indicates the space coordinate in the \alpha direction, c_p is the specific heat at constant pressure, \mu and \kappa are the viscosity and heat conductivity coefficients, respectively, \sigma is the surface tension, computed as \sigma = K{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{(\frac{\mathrm{\partial }\rho }{\mathrm{\partial }{r}_{\alpha }})}^2\mathrm{d}{r}_{\alpha }. Dimensionless variables are denoted by a “\wedge” sign on the left side. When the parameters from the non-ideal EOS and the physical property parameters of the real fluid are provided, it is feasible to recover the actual physical quantities from the numerical results [75]. In the manuscript, all quantities are dimensionless, and the symbol “\wedge” is dropped for simplicity.

To conduct simulations, the particle velocity space is discretized using the D2V33 discrete velocity set [83], while spatial derivatives are discretized using a 16th-order fast Fourier transform (FFT) scheme [34] and temporal derivative is handled using the second-order Runge-Kutta finite difference method. The validity of the physical model and numerical methods has been confirmed in Ref. [76]. The computational grid is set to N_x\times N_y=256\times 256, with spatial steps of \mathrm{\Delta }x= \mathrm{\Delta }y=1/128. The parameters a and b in EOS are set to a=2.0 and b=0.4, which determine the critical point at T_c= 1.88657, {\rho }_c=1.3044, and P_c= 0.8832.

Two initially stationary droplets of identical size are positioned as

(12)\begin{array}{rl}\rho (x,y)=& {\rho }_{\mathrm{L}}-\frac{({\rho }_{\mathrm{L}}-{\rho }_{\mathrm{V}})}{2}\tanh \left[\frac{\sqrt{{(x-x_{0\mathrm{l}})}^2+{(y-y_{0\mathrm{l}})}^2}-r_0}{0.5w}\right]\\ & -\frac{({\rho }_{\mathrm{L}}-{\rho }_{\mathrm{V}})}{2}\tanh \left[\frac{\sqrt{{(x-x_{0\mathrm{r}})}^2+{(y-y_{0\mathrm{r}})}^2}-r_0}{0.5w}\right].\end{array}

Here, {\rho }_{\mathrm{L}} = 2.0658 and {\rho }_{\mathrm{V}} = 0.6894 are the densities of the liquid and vapor phases at a temperature of T=1.8. The total width of the boundary layer is 2w=12\mathrm{\Delta }x. The initial radius of both droplets is the same, denoted as r_0. The center coordinates of the left and right droplets are ( x_{0\mathrm{l}},y_{0\mathrm{l}} ) and ( x_{0\mathrm{r}},y_{0\mathrm{r}} ), respectively. The distance between the two centers is \left|x_{0\mathrm{l}}-x_{0\mathrm{r}}\right|= 2r_0+ r_d. For all cases, y_{0\mathrm{l}}=y_{0\mathrm{r}} and r_d=\omega \mathrm{\Delta }x, where \omega is a real constant. The relaxation time used in the simulation is \tau =0.0025, with Pr=0.5. Unless otherwise specified, K=0.0001.

3.1 Differences in HNE effects and coalescence mechanisms

Fig.1 shows the spatial distribution of density at eight characteristic moments with similar morphology under isothermal and non-isothermal conditions. The droplet coalescence process is significantly slower in the non-isothermal case compared to the isothermal case. Specifically, under non-isothermal condition, the initiation of the coalescence process is delayed, cut-through occurs more slowly, the rapid coalescence phase lasts longer, and exhibits noticeable damped oscillations.

(i) Slower initiation and cut-through processes: The coalescence initiation duration under non-isothermal condition is much longer than that under isothermal condition. For the isothermal case, the cut-through time of two static droplets is t_c=1.8 [see Fig.1A(b)], while for the non-isothermal case, it is t_c=5.6 [see Fig.1B(b)]. The time difference in both cases is \mathrm{\Delta }{t}_c=3.8. The differences in cut-through time can be more obviously reflected in the evolution of morphological factor. As shown in Fig.2, in both cases, at t=t_c, the morphological factor L_d=L_y/L_x and L_y begin to increase, where L_x and L_y are the long and short semi-axes of the new droplet, respectively.

(ii) Longer duration of the rapid coalescence process: As shown in Fig.1 and Fig.2, under isothermal condition, the time at which the new droplet first forms an unstable spherical shape, corresponding to the morphological factor L_d first equaling 1, is t_{l1} =7.0. Then the time difference between cut-through and the first formation of an unstable spherical droplet is \mathrm{\Delta }{t}_{l1c}=t_{l1} -t_c=5.2. Under non-isothermal condition, the time for the new droplet to first forms an unstable spherical shape is t_{l1}=15.2, and the time difference between cut-through and the first formation of an unstable spherical droplet is \mathrm{\Delta } t_{l1c}=t_{l1}-t_c=9.6. The time difference for the first formation of an unstable spherical droplet in both cases is \mathrm{\Delta }{t}_{l1} =8.2. Under non-isothermal condition, the growth rate of the morphological factor L_d of the new droplet is smaller than under isothermal condition, because the slopes of both the long semi-axis L_x and the short semi-axis L_y of the droplet are smaller under non-isothermal condition (see Fig.2).

(iii) Greater oscillation amplitude and longer oscillation period: After coalescence, the newly formed larger droplet often exhibits a highly non-circular shape. Since the system is not initially at its minimum surface energy, surface tension drives the droplet toward restoring a circular configuration. During this relaxation process, due to fluid inertia, the droplet typically does not reach equilibrium in a single step. Instead, it undergoes damped oscillations around the spherical shape until its kinetic energy is fully dissipated. The timescale of these oscillations is governed jointly by viscous dissipation, heat transfer, and surface tension effects [18]. As shown in Fig.2(a), under isothermal condition, the time at which L_d reaches its maximum is t_{l\mathrm{m}\mathrm{a}\mathrm{x}}=9.6, whereas under non-isothermal condition, t_{l\mathrm{m}\mathrm{a}\mathrm{x}}=20.0, with a difference of \mathrm{\Delta } t_{l\mathrm{m}\mathrm{a}\mathrm{x}}=10.4. After the first oscillation period, under isothermal condition, the oscillation amplitude of L_x and L_y is almost zero, as shown in Fig.2(b). In contrast, under non-isothermal condition, small oscillations persist. Furthermore, the maximum value of L_d under non-isothermal condition is greater than that under isothermal condition. This can be attributed to the more complex stress and flow field distributions near the vapor−liquid interface under non-isothermal condition, which prevent the liquid−vapor interface from stabilizing easily [85].

To explain the above differences, Fig.3 presents the spatial distribution along the horizontal midline of various quantities: temperature T, latent heat of phase change \mathrm{\Delta } H=H_t-H_0 (with enthalpy H=ϵ+P /\rho, the internal energy density ϵ= \rho T-a {\rho }^2+K{\left|\mathbf{\nabla }\rho \right|}^2/2, and H_t and H_0 being the enthalpies at time t and initial state), velocity in the coalescence direction u_x, density \rho, pressure P, the xx-component of the total pressure tensor \mathbf{\Gamma }=P\mathit{I}+\mathbf{\Lambda } (denoted as {\mathrm{\Gamma }}_{xx}), and the xx-component of the pressure tensor contribution from the macroscopic gradient \mathbf{\Lambda } (denoted as {\mathrm{\Lambda }}_{xx}).

In the isothermal case, it is assumed that the system is in contact with a sufficiently large heat bath, resulting in instantaneous heat exchange with the environment [34]. This assumption leads to no temperature gradient within the system. In contrast, the non-isothermal case incorporates the temperature field’s inhomogeneity, introducing additional HNE and TNE effects, such as the temperature-induced Marangoni effect and heat conduction. These effects significantly alter the system’s dynamics and flow field morphology.

As shown in Fig.3(a), the release of latent heat in the non-isothermal case creates a localized high-temperature region in the liquid bridge between the droplets, leading to a prolonged local increase in pressure P. In the isothermal case, the pressure peak at the center of the liquid bridge oscillates but gradually decreases after t\ge 0.2, reducing by approximately 11.5\mathrm{\%} at t_c=1.8. In contrast, the pressure peak in the non-isothermal case remains nearly constant during the period 0.2\le t\le t_c=5.6.

The pressure gradient in the non-isothermal case is larger, resulting in stronger resistance to droplet coalescence. Although the gradient of surface tension \mathbf{\nabla }\cdot \mathbf{\Lambda } is larger in the non-isothermal case [with larger positive peaks and smaller negative troughs, as shown by the curves on the right axis of Fig.3(d) and (f)], the gradient of the total pressure tensor \mathbf{\nabla }\cdot \mathbf{\Gamma } is greater in the isothermal case [as shown by the curves on the left axes of Fig.3(d) and (f)]. This explains why the coalescence initiation time is significantly prolonged in the non-isothermal case.

As shown in Fig.3(b), at the same moment and state, the flow velocities inside and outside the liquid−vapor boundary layer in the non-isothermal case are consistently much smaller than those in the isothermal case. Even at t=t_c=5.6, the coalescence velocity is still smaller compared to t=t_c=1.8 in the isothermal case.

Interestingly, for t\le t_c, while the pressure gradient decreases in both cases, the manner in which it decreases is distinctly different. In the isothermal case, the pressure gradient decreases by maintaining the trough value approximately constant, while the pressure peak in the liquid bridge region gradually decreases. In contrast, in the non-isothermal case, the pressure peak remains high, and the trough value gradually increases. This behavior indicates that, in the droplet coalescence process, the pressure gradient influences the expansion of the liquid−vapor interface in two opposite ways in the isothermal and non-isothermal cases. This distinction is more clearly observed in the velocity field, as shown in Fig.4.

Fig.4 illustrates the spatial distribution of the flow field at eight characteristic moments with similar morphology, under both isothermal and non-isothermal conditions. The blue solid lines represent the liquid-vapor interface during coalescing. A comparative analysis of the flow field distribution in Fig.4 reveals significant differences in the velocity vortices between the two conditions.

(i) Under isothermal condition, when t\le 3.8, the velocity in the vapor phase along the vertical direction of coalescence (the y-direction) converges toward the saddle point in the saddle region between the two droplets. No vortices are formed at the liquid-vapor interface in this region. In contrast, under non-isothermal condition, when t\le 10.8, the velocity in the vapor phase outside the saddle point diverges along the vertical direction of coalescence. This divergence creates vortices with opposite rotational directions at the liquid−vapor interface in the saddle region. The non-isothermal case results in a complex flow pattern near the interface, with the Marangoni effect significantly enhancing the flow outside the interface.

(ii) Under isothermal condition, due to the faster release of surface energy, the peak of the velocity in merging direction is higher compared to the non-isothermal case. As shown in Fig.5, the coalescence speed under isothermal condition starts faster and reaches its peak earlier. Specifically, the time when the coalescence speed reaches its maximum for the isothermal case is t_{u\mathrm{m}\mathrm{a}\mathrm{x}}=3.8, while for the non-isothermal case, t_{u\mathrm{m}\mathrm{a}\mathrm{x}}=10.8. Furthermore, the peak value {\overline{u}}_{\mathrm{m}\mathrm{a}\mathrm{x}} is larger in the isothermal case. However, under isothermal condition, the coalescence speed decays more rapidly, particularly during the damping oscillation phase (t\ge t_{l1}). After the first oscillation cycle, the coalescence speed under isothermal condition almost decays to zero, while under non-isothermal condition, the coalescence speed continues to oscillate for several cycles before gradually decaying to zero.

3.2 Differences in TNE effects

The upper section discusses the differences in HNE between isothermal and non-isothermal cases during droplet coalescence. Below, we analyze the more fundamental TNE differences, as TNE is both the source of HNE.

Fig.6 shows the spatial distribution of the three components of the NOMF, {\mathbf{\Delta }}_2^{\ast }, at four characteristic moments under both isothermal and non-isothermal conditions. As a TNE characteristic, the spatiotemporal evolution of the flow field in both cases governs the distribution and evolution of {\mathbf{\Delta }}_2^{\ast }, since {\mathbf{\Delta }}_2^{\ast } primarily driven by the velocity gradient. As a TNE effect, {\mathbf{\Delta }}_2^{\ast } influences the flow field evolution at faster time scales and smaller spatial scales.

The key differences between the two cases are as follows: (i) In the isothermal case, the spatial distributions of {\mathrm{\Delta }}_{2xx}^{\ast } and {\mathrm{\Delta }}_{2yy}^{\ast } are not completely opposite due to the different directions of the velocity vortices. In contrast, these distributions are fully antisymmetric in the non-isothermal case. (ii) In the non-isothermal case, the spatial distribution of {\mathrm{\Delta }}_{2xy}^{\ast } is more complex due to the increased variety of non-equilibrium driving forces. {\mathrm{\Delta }}_{2xy}^{\ast } in the non-isothermal case evolves from an interior−exterior antisymmetric double quadrupole, while in the isothermal case, it evolves from an antisymmetric single quadrupole structure. The coupling of multiple TNE effects in the non-isothermal system makes the spatiotemporal evolution of non-equilibrium behaviors more complex.

Fig.7 illustrates the time evolution of several non-equilibrium intensities, including the total non-equilibrium intensity {\overline{D}}^{\ast }, NOMF intensity {\overline{D}}_2^{\ast }, NOEF intensity {\overline{D}}_{3,1}^{\ast }, strength of flux of NOMF {\overline{D}}_3^{\ast }, and strength of flux of NOEF {\overline{D}}_{4,2}^{\ast } in both isothermal and non-isothermal systems. The main differences between the two cases are as follows:

(i) Duration and peak of TNE effects: In the isothermal case, non-equilibrium effects last for a shorter time but have higher peak values. For example, all intensities reach their maximum at approximately t=2.2 for the isothermal case, which is earlier than in the non-isothermal case, where they peak at t=8.6. Additionally, the peak value of {\overline{D}}^{\ast } is about twice as large in the isothermal case compared to the non-isothermal case. Similar to the non-isothermal case [29], in the isothermal condition, the total TNE intensity {\overline{D}}^{\ast } reaches its maximum value at t=2.2, when the boundary length decays the fastest and the merging acceleration is the greatest; it reaches its first local trough at t=3.8, when the merging velocity is the highest. That is, the characteristic points on the evolution curve of {\overline{D}}^{\ast } can be used to delineate the various stages of the droplet coalescence process.

(ii) Damping oscillation: The non-isothermal case exhibits more pronounced damping oscillations in {\overline{D}}_2^{\ast }, {\overline{D}}_{4,2}^{\ast }, and {\overline{D}}^{\ast } compared to the isothermal case, with larger amplitudes, longer periods, and extended durations.

(iii) Interplay of various TNE effects: In the isothermal case, there is no temperature gradient, so {\overline{D}}_3^{\ast } and {\overline{D}}_{3,1}^{\ast } tend to zero, and the total non-equilibrium intensity is mainly driven by the velocity gradients ({\overline{D}}_2^{\ast } and {\overline{D}}_{4,2}^{\ast }). In the non-isothermal case, when t\le 5.0, the total non-equilibrium intensity is dominated by {\overline{D}}_3^{\ast } and {\overline{D}}_{3,1}^{\ast } due to the temperature gradient. After t=5.0, the velocity gradients increase, and {\overline{D}}_2^{\ast } and {\overline{D}}_{4,2}^{\ast } take over as the dominant contributors. The alternation between the dominance of non-equilibrium behaviors driven by temperature gradients and velocity gradients corresponds to the moment when the two droplets begin to make contact.

This is because, in the early stage of the relaxation process, when the droplets approach each other, the velocity gradient is negligible, and the temperature gradient due to latent heat of phase change governs the system’s evolution. After the droplets make contact, the relative velocity between particles increases, causing the velocity gradient to dominate. As shown in Fig.5, at t=5.0, the droplets begin to make contact but have not yet merged, marking the moment when the boundary length rapidly decreases. In the later stage of coalescence, {\overline{D}}_{3,1}^{\ast } and {\overline{D}}_3^{\ast } decay smoothly to zero, while {\overline{D}}_2^{\ast } and {\overline{D}}_{4,2}^{\ast } exhibit damped oscillations and also decay to zero as the droplet shape evolves. In summary, the increasing number of driving factors in the non-isothermal case adds the interplay and complexity to the system, as evidenced by the alternating dominance of behaviors driven by temperature gradients ({\mathbf{\Delta }}_3^{\ast }, {\mathbf{\Delta }}_{3,1}^{\ast }) and velocity gradients ({\mathbf{\Delta }}_2^{\ast }, {\mathbf{\Delta }}_{4,2}^{\ast }).

(v) Higher-order TNE effects in the isothermal case: Theoretically, the NOEF contains contributions from various orders of non-equilibrium effects, as expressed by {\mathbf{\Delta }}_{3,1}^{\ast } ={\mathbf{\Delta }}_{3,1}^{\ast (1)}+{\mathbf{\Delta }}_{3,1}^{\ast (2)}+\cdots, where {\mathbf{\Delta }}_{3,1}^{\ast }= {\mathbf{\Delta }}_{3,1}^{\ast (1)} in the first-order model. According to theoretical derivation, {\mathbf{\Delta }}_{3,1}^{\ast (1)} only includes contributions from the temperature gradient, whereas {\mathbf{\Delta }}_{3,1}^{\ast (2)} incorporates both temperature and velocity gradient contributions [71]. Thus, in the isothermal case, despite the absence of a temperature gradient, {\overline{D}}_{3,1}^{\ast } and {\overline{D}}_3^{\ast } are not exactly zero. This phenomenon reflects the complex cross-coupling of multiple physical quantities in constitutive relations, which is one manifestation of the multiscale system’s complexity. This also indicates that the current DBM framework primarily focuses on first-order non-equilibrium effects, but naturally captures partial second-order non-equilibrium behaviors without requiring additional modifications.

Overall, as a phenomenon, non-equilibrium fundamentally embodies system complexity by revealing the nonlinear coupling among various physical mechanisms. As a driving effect, non-equilibrium significantly influences system evolution through its intrinsic, dominant, and integrative characteristics. Non-equilibrium quantities thus serve as sensitive probes for identifying competing mechanisms, providing earlier, more sensitive, and more precise detection of dominant processes and their transitions compared to conventional macroscopic variables.

3.3 Differences in entropy production mechanisms

Fig.8 illustrates the spatial distribution of the total entropy production rate {\dot{S}}_{\mathrm{p}\mathrm{r}} at four characteristic times, showing similar morphology for both isothermal and non-isothermal conditions. For t\le t_c, {\dot{S}}_{\mathrm{p}\mathrm{r}} is mainly concentrated in the quasi-contact region between the two droplets in both two cases, as shown in Fig.8. The difference lies in the fact that {\dot{S}}_{\mathrm{p}\mathrm{r}} in the isothermal case is more concentrated, with its peak about three times higher than in the non-isothermal case. During the rapid coalescence stage (t_c\le t\le t_{l1}), {\dot{S}}_{\mathrm{p}\mathrm{r}} in the isothermal case is primarily concentrated on the outer side of the saddle region of the liquid-vapor boundary layer between the two droplets. In contrast, in the non-isothermal case, {\dot{S}}_{\mathrm{p}\mathrm{r}} is mainly concentrated in the liquid bridge region between the two droplets.

The reasons are as follows. In the isothermal case, {\dot{S}}_{\mathrm{p}\mathrm{r}} is solely influenced by NOMF where the velocity gradient is concentrated corresponds to the main distribution of entropy production. In the non-isothermal case, {\dot{S}}_{\mathrm{p}\mathrm{r}} is influenced by both NOMF and NOEF. The inhomogeneous temperature field caused by latent heat broadens the entropy production region, while simultaneously suppressing local flow velocity, viscous stress, and the resulting entropy production.

Fig.9 presents the entropy production rate evolution curves over time for both isothermal and non-isothermal conditions. A detailed comparison of the total entropy production rate {\overline{\dot{S}}}_{\mathrm{p}\mathrm{r}} shows distinct behaviors between the two scenarios. In the isothermal case, {\overline{\dot{S}}}_{\mathrm{p}\mathrm{r}} reaches its peak at t=2.2, with a higher peak value compared to the non-isothermal case. However, it decays rapidly, displaying a relatively smooth, unimodal structure with small fluctuations and stabilizing after t\ge 10. In contrast, the non-isothermal case reveals a more intricate evolution pattern. Although the initial peak is lower than in the isothermal case, the system experiences significant fluctuations around t=7.2 and t=24.4, leading to a multi-peak structure. This behavior suggests that, in the non-isothermal case, the disturbances induced by the temperature gradient alongside the velocity gradient contribute to extended and prolonged non-equilibrium dissipation processes. In the isothermal case, the isothermal entropy production curve can also be used to identify various stages of droplet coalescence: the peak point (t=2.2) and the first local trough point (t=3.8) on the evolution curve correspond to the moments of maximum merging acceleration and maximum merging velocity, respectively.

Comparing the entropy production rates of the NOMF, {\overline{\dot{S}}}_{\mathrm{N}\mathrm{O}\mathrm{M}\mathrm{F}}=\sum \rho {\dot{S}}_{\mathrm{N}\mathrm{O}\mathrm{M}\mathrm{F}}/\sum \rho, and that of the NOEF, {\overline{\dot{S}}}_{\mathrm{N}\mathrm{O}\mathrm{E}\mathrm{F}}=\sum \rho {\dot{S}}_{\mathrm{N}\mathrm{O}\mathrm{E}\mathrm{F}}/\sum \rho, in the non-isothermal case, shown in Fig.9, it is can be observed that for t\le 3.8, the total entropy production is primarily driven by NOEF, with {\overline{\dot{S}}}_{\mathrm{N}\mathrm{O}\mathrm{E}\mathrm{F}} gradually decreasing. Between , both {\overline{\dot{S}}}_{\mathrm{N}\mathrm{O}\mathrm{M}\mathrm{F}} and {\overline{\dot{S}}}_{\mathrm{N}\mathrm{O}\mathrm{E}\mathrm{F}} increase gradually, with similar amplitudes and growth rates. This suggests that during this period, the two entropy production mechanisms are competing and growing together. After t=5.6, when the droplets fully cut-through, {\overline{\dot{S}}}_{\mathrm{N}\mathrm{O}\mathrm{M}\mathrm{F}} surpasses {\overline{\dot{S}}}_{\mathrm{N}\mathrm{O}\mathrm{E}\mathrm{F}} and begins to dominate the system’s entropy production.

In summary, akin to the system’s TNE effects, the two entropy production mechanisms associated with temperature and velocity gradients alternate in their dominance throughout the droplet coalescence process. During the damping oscillation phase, {\overline{\dot{S}}}_{\mathrm{N}\mathrm{O}\mathrm{M}\mathrm{F}}, also exhibits noticeable damping oscillations.

Fig.10 displays the total entropy production of the system in both isothermal and non-isothermal conditions. As mentioned previously, in the isothermal case, entropy production is primarily due to NOMF; in the non-isothermal case, entropy production is contributed by both NOMF and NOEF, with the former being the dominant contributor. Compared to the isothermal case, the non-isothermal system exhibits a larger total entropy production due to the contribution of NOEF, though the entropy increase from NOEF is smaller, with S_{\mathrm{N}\mathrm{O}\mathrm{E}\mathrm{F}} being approximately 1/3 of S_{\mathrm{N}\mathrm{O}\mathrm{M}\mathrm{F}}. While the peak value of {\overline{\dot{S}}}_{\mathrm{N}\mathrm{O}\mathrm{M}\mathrm{F}} in the isothermal case is about 2.6 times higher than in the non-isothermal case, the range of {\dot{S}}_{\mathrm{N}\mathrm{O}\mathrm{M}\mathrm{F}} distribution is wider in the non-isothermal case, and its duration is longer. As a result, the entropy production from NOMF in the non-isothermal case, S_{\mathrm{N}\mathrm{O}\mathrm{M}\mathrm{F}}, is slightly greater than in the isothermal case.

3.4 Effects of initial distance on the coalescence process

The droplet coalescence process is influenced by various factors, with the initial distance between droplets playing a crucial role in coalescence velocity, mode, and final morphology. The initial separation between droplets determines the intensity of their interaction during coalescence, the formation and expansion of the liquid bridge, and the changes in surface tension and pressure gradients. For instance, in applications such as spraying and atomization, the initial distance between droplets significantly affects atomization efficiency, particle size, and spray uniformity. Optimizing this distance can improve atomization performance and enhance the effectiveness of processes like combustion, fertilizer spraying, and spray-based drug delivery [3, 4, 86].

Fig.11 shows the evolution curves of the droplet’s morphology factor L_d over time for different initial distances r_d between the droplets, under both isothermal and non-isothermal conditions. As shown, the evolution trends of L_d are similar for both isothermal and non-isothermal cases, except for the significant difference in the cut-through time t_c. In the isothermal case, t_{l\mathrm{m}\mathrm{a}\mathrm{x}}-t_c\approx 7.8, while in the non-isothermal case, t_{l\mathrm{m}\mathrm{a}\mathrm{x}}-t_c\approx 14.4. This difference arises because after the cut-through, the rapid coalescence process is primarily controlled by surface tension. When droplet scale, surface tension coefficient, and viscosity remain constant, the droplet deformation, flow field distribution, and energy and momentum transfer are similar.

Fig.12 shows how the initial distance r_d between two droplets affects the cut-through time t_c. In both isothermal and non-isothermal cases, the cut-through time t_c increases with the initial droplet distance r_d, following a power-law relationship. For identical parameters, the exponent in the non-isothermal case is smaller than that in the isothermal case. This behavior is due to the similarity in the system’s dynamic. In the non-isothermal case, the temperature gradient induces energy dissipation mechanisms, such as the Marangoni effect and heat conduction, which delay the droplet coalescence stage.

When r_d is large, the contact mode between droplets changes in the isothermal case, reducing the effect of surface tension and significantly prolonging the droplet approach process, which sharply increases the cut-through time t_c. The functional relationships between t_c and r_d in the isothermal case are

t_c=( 11.5r_d{)}^{13.92} -0.093(K={10}^{-4})

and

t_c=( 9.3r_d{)}^{11.35} -0.043(K=1.5 \times {10}^{-4}).

The functional relationships between t_c and r_d in the non-isothermal case are

t_c=( 13.3r_d{)}^{7.02} +0.96(K={10}^{-4})

and

t_c=( 12.4r_d{)}^{6.13} +0.83(K=1.5 \times {10}^{-4}).

The power exponent decreases as the surface tension coefficient K increases. In other words, a larger surface tension coefficient K results in a smaller crossing time t_c. As shown in Fig.13, t_c and K follow a negative power function relationship. This can be attributed to two factors: First, surface tension is the primary driving force for droplet coalescence, regardless of whether the case is isothermal or non-isothermal. A larger surface tension coefficient K increases the surface tension, accelerating droplet approach and shortening the crossing time. Second, a larger surface tension coefficient increases the liquid-vapor interface width w=\sqrt{2K/(T_c/T-1)} [87], effectively reducing r_d and shortening the crossing time.

In the isothermal case shown in Fig.13, the specific functional relationships are

t_c=1.77\times {10}^{-7}{K}^{-1.62}+0.085 (r_d=10\mathrm{\Delta } x),

t_c=3.03\times {10}^{-12}{K}^{-2.93}+0.19 (r_d=12\mathrm{\Delta } x),

and

t_c=2.68\times {10}^{-21}{K}^{-5.58}( r_d=15 \mathrm{\Delta }x).

In the non-isothermal case shown in Fig.13, the specific functional relationships are

t_c=4.94\times {10}^{-8}{K}^{-1.9}+0.69 (r_d=10\mathrm{\Delta } x),

t_c=1.06\times {10}^{-9}{K}^{-2.4}+1.36 (r_d=12\mathrm{\Delta } x),

and

t_c=1.66\times {10}^{-12}{K}^{-3.27}+2.86 (r_d=15\mathrm{\Delta } x).

It should be noted that, in the isothermal case with K=0.0001, when r_d=0.117, the cut-through time t_c approaches infinity, indicating that the droplets no longer merge. In fact, each set of parameters defines a limiting initial merging distance, r_{d\mathrm{m}\mathrm{a}\mathrm{x}}. For the same parameters, the limiting initial merging distance, r_{d\mathrm{m}\mathrm{a}\mathrm{x}}, is larger in the non-isothermal case than in the isothermal case. This is because, in the non-isothermal case, the temperature gradient-induced Marangoni effect and thermal convection increase the complexity of the flow field around the droplets, expanding the droplet coalescence range.

Fig.14(a) and (b) present the time evolution curves of the total non-equilibrium intensity, {\overline{D}}^{\ast }, for the isothermal and non-isothermal cases, respectively. As shown in Fig.14, the time evolution of {\overline{D}}^{\ast }, exhibit similarity for different r_d. In both cases, a power-law relationship exists between the time, t_{D\mathrm{m}\mathrm{a}\mathrm{x}}, at which the total non-equilibrium intensity reaches its maximum, and r_d. Moreover, the power exponent is greater in the isothermal case than in the non-isothermal case. The specific relationships between t_{D\mathrm{m}\mathrm{a}\mathrm{x}} and r_d are: t_{D\mathrm{m}\mathrm{a}\mathrm{x}}=3.46 \times{10}^{14} r_d^{13.63} +0.21 (isothermal case) and t_{D\mathrm{m}\mathrm{a}\mathrm{x}}=2.45\times{10}^7{r}_d^{6.5} +3.33 (non-isothermal case).

Fig.15(a) and (b) illustrate the time evolution of the total entropy production rate, {\overline{\dot{S}}}_{\mathrm{p}\mathrm{r}}, in the isothermal and non-isothermal cases, respectively. The initial distance, r_d, between the droplets has minimal impact on the time evolution and values of {\overline{\dot{S}}}_{\mathrm{p}\mathrm{r}}, except that a larger r_d delays the time at which {\overline{\dot{S}}}_{\mathrm{p}\mathrm{r}} reaches its maximum. Both isothermal and non-isothermal cases also exhibit power-law relationships between t_{\dot{S}\mathrm{p}\mathrm{r}\mathrm{m}\mathrm{a}\mathrm{x}} (when {\overline{\dot{S}}}_{\mathrm{p}\mathrm{r}} reaches its maximum) and r_d, with the exponent in the isothermal case larger than that in the non-isothermal case. Additionally, r_d has little influence on the amounts of the three types of entropy production (S_{\mathrm{N}\mathrm{O}\mathrm{M}\mathrm{F}}, S_{\mathrm{N}\mathrm{O}\mathrm{E}\mathrm{F}}, S_{\mathrm{p}\mathrm{r}}), as shown in Fig.16. The relationships between t_{\dot{S}\mathrm{p}\mathrm{r}\mathrm{m}\mathrm{a}\mathrm{x}} and r_d in Fig.15(c) also follow power-laws: t_{\dot{S}\mathrm{p}\mathrm{r}\mathrm{m}\mathrm{a}\mathrm{x}}=5.99\times {10}^{14}{r}_d^{13.89}+0.08 (isothermal) and t_{\dot{S}\mathrm{p}\mathrm{r}\mathrm{m}\mathrm{a}\mathrm{x}}=2.63\times {10}^7{r}_d^{6.54}+2.03 (non-isothermal).

4 Conclusions and outlooks

This paper presents a detailed analysis of the coalescence process of two initially static droplets under both isothermal and non-isothermal conditions, focusing on the impact of the initial distance between the droplets. Using the DBM method, the study explores several key aspects. The main findings include:

(i) Coalescence dynamics: In both cases, surface tension drives droplet coalescence, while the resistance arises from an isotropic pressure gradient. However, in the non-isothermal case, the temperature rise due to phase transition latent heat significantly increases the isotropic pressure gradient, though it slightly increases the surface tension gradient near the droplet contact point. This results in a smaller total pressure tensor gradient leading to the deferment of droplet coalescence in the non-isothermal case.

(ii) Flow field complexity: The non-isothermal case exhibits a significantly more complex flow field due to the effects of thermal Marangoni forces and heat conduction, which contribute to greater complexity in the evolution of droplet morphology. The spatial distributions of various non-equilibrium driving forces also differ markedly between the two cases.

(iii) Thermodynamic non-equilibrium: In the isothermal case, non-equilibrium behaviors are primarily driven by velocity gradients, with a stronger non-equilibrium intensity. In the non-isothermal case, however, non-equilibrium behaviors alternate between temperature gradients and velocity gradients as the droplets begin to merge, leading to more complex oscillations in non-equilibrium and entropy production characteristics. The shift in the dominant factor reflects the complexity of the system, which arises from the combined influence of multiple driving factors.

(iv) Entropy production: In the isothermal case, entropy production is solely due to NOMF, while in the non-isothermal case, both NOMF and NOEF contribute, with NOMF being the dominant contributor. The contributions from NOMF and NOEF alternate before and after droplet cut-through.

(v) Effect of initial distance: Increasing the initial distance between the droplets delays the coalescence time (t_c) but does not significantly alter the system’s morphological, non-equilibrium, or entropy production characteristics. The relationships between key times (e.g., t_c, t_{D\mathrm{m}\mathrm{a}\mathrm{x}}, t_{{\dot{S}}_{\mathrm{p}\mathrm{r}\mathrm{m}\mathrm{a}\mathrm{x}}}) and r_d follow a power-law, with larger exponents in the isothermal case.

This study demonstrates that the dynamics of droplet coalescence are heavily influenced by thermodynamic non-equilibrium effects, particularly in the non-isothermal case. Non-equilibrium measures can provides an earlier and more sensitive means of detecting the dominant processes associated with different mechanisms and their transitions. The initial distance between droplets primarily impacts the timing of coalescence, with minimal effect on other characteristics, and the system exhibits similar behavior in its evolution.

In future work, the following three research directions are proposed to expand and deepen the current study:

(i) Extension beyond the local equilibrium assumption: The current external force are introduced base on the equilibrium distribution function. This assumption requires that the phase separation process proceeds much slower than the relaxation of TNE effects, i.e., t_{\mathrm{p}\mathrm{s}}\gg \tau, where t_{\mathrm{p}\mathrm{s}} denotes the characteristic time of phase separation. However, as the phase separation process accelerates and the degree of discreteness or nonequilibrium increases, the validity of the local equilibrium assumption becomes questionable. In particular, when the equation of state becomes invalid or unknown, the current modeling approach becomes inadequate. Under such conditions, it becomes necessary to explicitly incorporate intermolecular potential functions to accurately describe molecular interactions and recover the correct thermodynamic behavior [88, 89].

(ii) Influence of solid walls and shear flows: In microfluidic systems, the presence of solid walls and external shear flows has a significant impact on the droplet coalescence process. Specifically, factors such as the contact angle, wall temperature, and shear strength play critical roles in determining the dynamics of droplet deformation, coalescence, and stabilization. Future investigations should incorporate wetting boundary conditions and shear effects to better capture the complex physical phenomena occurring in confined or structured environments.

(iii) Droplet deformation and breakup under shock waves: Future research will focus on the deformation, fragmentation, cavitation, and jetting processes of droplets under shock waves [63, 90, 91] — phenomena with far-reaching applications. In areas like spray combustion, engine injection, and aerospace, the shape and breakup of droplets are crucial in determining combustion efficiency, thrust, and dynamic performance. From a physical standpoint, this process is marked by increasingly significant thermodynamic non-equilibrium effects, making it an essential domain where kinetic methods are necessary for precise modeling and analysis.

References

[1] T. S. Kaminski and P. Garstecki, Controlled droplet microfluidic systems for multistep chemical and biological assays, Chem. Soc. Rev. 46(20), 6210 (2017)

[2] T. Moragues, D. Arguijo, T. Beneyton, C. Modavi, K. Simutis, A. R. Abate, J. C. Baret, A. J. deMello, D. Densmore, and A. D. Griffiths, Droplet-based microfluidics, Nat. Rev. Methods Primers 3(1), 32 (2023)

[3] H. Li, N. Riefler, T. Wriedt, and L. Mädler, Reference data set for three-dimensional measurements of double droplet combustion of p-xylene, Proc. Combust. Inst. 38(2), 3151 (2021)

[4] M. P. Sitte and E. Mastorakos, Large Eddy simulation of a spray jet flame using doubly conditional moment closure, Combust. Flame 199, 309 (2019)

[5] J. Zhao, Q. Li, L. Zhang, S. Lin, and L. Jiang, Experimental study on mitigation effects of water mist on blast wave, Explosion and Shock Waves 43(10), 105401 (2023)

[6] L. Yang, Y. Wei, and Y. Guo, Water impact of deformable spheres with vertical and oblique entries, Phys. Fluids 35(3), 033320 (2023)

[7] L. Yang, Y. Wei, W. Xia, G. Xia, and Z. Rao, A comparative study on water entry and water walking events of elastic spheres, Phys. Fluids 35(8), 085111 (2023)

[8] R. Chen, S. Zhou, L. Zhu, L. Zhu, W. Yan, and H. Yu, A new criterion of coalescence-induced microbubble detachment in three-dimensional microfluidic channel, Phys. Fluids 33(4), 043320 (2021)

[9] F. Shan, Z. Chai, B. Shi, J. Xiao, and C. Wang, Directional transport of a droplet on biomimetic ratchet structure, Phys. Fluids 35(11), 113314 (2023)

[10] J. D. Paulsen, J. C. Burton, S. R. Nagel, S. Appathurai, M. T. Harris, and O. A. Basaran, The inexorable resistance of inertia determines the initial regime of drop coalescence, Proc. Natl. Acad. Sci. USA 109(18), 6857 (2012)

[11] J. D. Paulsen, Approach and coalescence of liquid drops in air, Phys. Rev. E 88(6), 063010 (2013)

[12] J. D. Paulsen, R. Carmigniani, A. Kannan, J. C. Burton, and S. R. Nagel, Coalescence of bubbles and drops in an outer fluid, Nat. Commun. 5(1), 3182 (2014)

[13] J. Eggers, J. R. Lister, and H. A. Stone, Coalescence of liquid drops, J. Fluid Mech. 401, 293 (1999)

[14] C. R. Anthony, M. T. Harris, and O. A. Basaran, Initial regime of drop coalescence, Phys. Rev. Fluids 5(3), 033608 (2020)

[15] J. C. Burton and P. Taborek, Role of dimensionality and axisymmetry in fluid pinch-off and coalescence, Phys. Rev. Lett. 98(22), 224502 (2007)

[16] V. S. Sivasankar, S. A. Etha, D. R. Hines, and S. Das, Coalescence of microscopic polymeric drops: Effect of drop impact velocities, Langmuir 37(45), 13512 (2021)

[17] G. K. Seevaratnam, Y. Suo, E. Ramé, L. M. Walker, and S. Garoff, Dynamic wetting of shear thinning fluids, Phys. Fluids 19(1), 012103 (2007)

[18] J. Qian and C. K. Law, Regimes of coalescence and separation in droplet collision, J. Fluid Mech. 331, 59 (1997)

[19] J. R. Adam, N. R. Lindblad, and C. D. Hendricks, The collision, coalescence, and disruption of water droplets, J. Appl. Phys. 39(11), 5173 (1968)

[20] N. Ashgriz and J. Y. Poo, Coalescence and separation in binary collisions of liquid drops, J. Fluid Mech. 221, 183 (1990)

[21] Y. J. Jiang, A. Umemura, and C. K. Law, An experimental investigation on the collision behaviour of hydrocarbon droplets, J. Fluid Mech. 234, 171 (1992)

[22] M. Orme, Experiments on droplet collisions, bounce, coalescence and disruption, Pror. Energy Combust. Sci. 23(1), 65 (1997)

[23] J. P. Estrade, H. Carentz, G. Lavergne, and Y. Biscos, Experimental investigation of dynamic binary collision of ethanol droplets – a model for droplet coalescence and bouncing, Int. J. Heat Fluid Flow 20(5), 486 (1999)

[24] K. Pan, C. K. Law, and B. Zhou, Experimental and mechanistic description of merging and bouncing in head-on binary droplet collision, J. Appl. Phys. 103(6), 064901 (2008)

[25] K. Pan, P. Chou, and Y. Tseng, Binary droplet collision at high Weber number, Phys. Rev. E 80(3), 036301 (2009)

[26] G . Brenn, Droplet collision, in: Handbook of Atomization and Sprays: Theory and Applications, edited by N. Ashgriz, Springer US, Boston, MA, pp 157–181 (2011)

[27] C. L. Tang, P. Zhang, and C. K. Law, Bouncing, coalescence, and separation in head-on collision of unequal-size droplets, Phys. Fluids 24(2), 022101 (2012)

[28] Y. D. Zhang, A. G. Xu, J. J. Qiu, H. T. Wei, and Z. H. Wei, Kinetic modeling of multiphase flow based on simplified Enskog equation, Front. Phys. (Beijing) 15(6), 62503 (2020)

[29] G. Sun, Y. Gan, A. Xu, and Q. Shi, Droplet coalescence kinetics: Thermodynamic non-equilibrium effects and entropy production mechanism, Phys. Fluids 36(3), 032109 (2024)

[30] C. He and P. Zhang, Nonaxisymmetric flow characteristics in head-on collision of spinning droplets, Phys. Rev. Fluids 5(11), 113601 (2020)

[31] C. He, X. Xia, and P. Zhang, Vortex-dynamical implications of nonmonotonic viscous dissipation of off-center droplet bouncing, Phys. Fluids 32(3), 032004 (2020)

[32] A. Tiribocchi, A. Montessori, M. Lauricella, F. Bonaccorso, S. Succi, S. Aime, M. Milani, and D. A. Weitz, The vortex-driven dynamics of droplets within droplets, Nat. Commun. 12(1), 82 (2021)

[33] C. He, Z. He, and P. Zhang, Droplet collision of hypergolic propellants, Droplet 3(2), e116 (2024)

[34] Y. Gan, A. Xu, G. Zhang, Y. Li, and H. Li, Phase separation in thermal systems: A lattice Boltzmann study and morpho-logical characterization, Phys. Rev. E 84(4), 046715 (2011)

[35] K. Sun, F. Jia, P. Zhang, L. Shu, and T. Wang, Marangoni effect in bipropellant droplet mixing during hypergolic ignition, Phys. Rev. Appl. 15(3), 034076 (2021)

[36] A. Vrij and J. Th. G. Overbeek, Rupture of thin liquid films due to spontaneous fluctuations in thickness, J. Am. Chem. Soc. 90(12), 3074 (1968)

[37] S. Perumanath, M. K. Borg, M. V. Chubynsky, J. E. Sprittles, and J. M. Reese, Droplet coalescence is initiated by thermal motion, Phys. Rev. Lett. 122(10), 104501 (2019)

[38] M. Geri, B. Keshavarz, G. H. McKinley, and J. W. M. Bush, Thermal delay of drop coalescence, J. Fluid Mech. 833, R3 (2017)

[39] X. Liang, V. Kumar, F. Ahmadi, and Y. Zhu, Manipulation of droplets and bubbles for thermal applications, Droplet 1(2), 80 (2022)

[40] S. S. Dixit, H. Kim, A. Vasilyev, A. Eid, and G. W. Faris, Light-driven formation and rupture of droplet bilayers, Langmuir 26(9), 6193 (2010)

[41] A. Davanlou and R. Kumar, Thermally induced collision of droplets in an immiscible outer fluid, Sci. Rep. 5(1), 9531 (2015)

[42] J. C. Chrispell, R. Cortez, D. B. Khismatullin, and L. J. Fauci, Shape oscillations of a droplet in an Oldroyd-B fluid, Physica D 240(20), 1593 (2011)

[43] X. Xia, C. He, and P. Zhang, Universality in the viscous-to-inertial coalescence of liquid droplets, Proc. Natl. Acad. Sci. USA 116(47), 23467 (2019)

[44] C. He, X. Xia, and P. Zhang, Non-monotonic viscous dissipation of bouncing droplets undergoing off-center collision, Phys. Fluids 31(5), 052004 (2019)

[45] S. Ray, Y. Chi, P. Zhang, and S. Cheng, Head-on collision of unequal-size droplets on a wetting surface, Phys. Fluids 35(2), 022114 (2023)

[46] P. Zhang and C. K. Law, An analysis of head-on droplet collision with large deformation in gaseous medium, Phys. Fluids 23(4), 042102 (2011)

[47] Z. Hu, F. Chu, X. Wu, S. Ding, and Y. Lin, Bidirectional transport of split droplets, Phys. Rev. Appl. 18(4), 044057 (2022)

[48] S. Gao, Z. Yuan, and X. Wu, Coalescence-induced jumping of in-plane moving droplets: Effects of initial velocity and sideslip angle, Chem. Eng. Sci. 265, 118247 (2023)

[49] S. Gao and X. Wu, Numerical investigation on coalescence-induced jumping of centripetal moving droplets, Langmuir 38(41), 12674 (2022)

[50] J. Zhang, H. Liu, and X. Zhang, Modeling the deformation of a surfactant-covered droplet under the combined influence of electric field and shear flow, Phys. Fluids 33(4), 042109 (2021)

[51] S. Li, H. Liu, J. Zhang, F. Jiang, and G. Xi, Modeling of three- phase displacement in three-dimensional irregular geometries using a lattice Boltzmann method, Phys. Fluids 33(12), 122108 (2021)

[52] Y. Yu, D. Liang, and H. Liu, Lattice Boltzmann simulation of immiscible three-phase flows with contact-line dynamics, Phys. Rev. E 99(1), 013308 (2019)

[53] Y. Yu, H. Liu, D. Liang, and Y. Zhang, A versatile lattice Boltzmann model for immiscible ternary fluid flows, Phys. Fluids 31(1), 012108 (2019)

[54] H. Liu and Y. Zhang, Phase-field modeling droplet dynamics with soluble surfactants, J. Comput. Phys. 229(24), 9166 (2010)

[55] A. Mazloomi Moqaddam, S. S. Chikatamarla, and I. V. Karlin, Simulation of droplets collisions using two-phase entropic lattice Boltzmann method, J. Stat. Phys. 161(6), 1420 (2015)

[56] A. Zarghami and H. E. A. Van den Akker, Thermohydrodynamics of an evaporating droplet studied using a multi-phase lattice Boltzmann method, Phys. Rev. E 95(4), 043310 (2017)

[57] X. An, B. Dong, Y. J. Zhang, and X. Zhou, Discussion on interface deformation and liquid breakup mechanism in vapor-liquid two-phase flow, Chin. Phys. B 32(9), 094702 (2023)

[58] H. Liang, H. Liu, Z. Chai, and B. Shi, Lattice Boltzmann method for contact-line motion of binary fluids with high density ratio, Phys. Rev. E 99(6), 063306 (2019)

[59] X. Xu, F. Wang, Z. Qin, and B. Wen, Electrowetting lattice Boltzmann method for micro- and nano-droplet manipulations, Phys. Rev. E 107(4), 045305 (2023)

[60] B. Huang, H. Liang, and J. Xu, Lattice Boltzmann simulation of binary three-dimensional droplet coalescence in a confined shear flow, Phys. Fluids 34(3), 032101 (2022)

[61] K. H. Luo, L. Fei, and G. Wang, A unified lattice Boltzmann model and application to multiphase flows, Philos. Trans. Royal Soc. A 379(2208), 20200397 (2021)

[62] G. Wang, L. Fei, and K. H. Luo, Unified lattice Boltzmann method with improved schemes for multiphase flow simulation: Application to droplet dynamics under realistic conditions, Phys. Rev. E 105(4), 045314 (2022)

[63] Y. Yang, J. Tu, M. Shan, Z. Zhang, C. Chen, and H. Li, Acoustic cavitation dynamics of bubble clusters near solid wall: A multiphase lattice Boltzmann approach, Ultrason. Sonochem. 114, 107261 (2025)

[64] Y. Wei, Y. Li, Z. Wang, H. Yang, Z. Zhu, Y. Qian, and K. H. Luo, Small-scale fluctuation and scaling law of mixing in three-dimensional rotating turbulent Rayleigh-Taylor instability, Phys. Rev. E 105(1), 015103 (2022)

[65] Z. Wang, Y. Wei, and Y. Qian, A novel thermal lattice Boltzmann model with heat source and its application in incompressible flow, Appl. Math. Comput 427, 127167 (2022)

[66] Z. Yang, S. Liu, C. Zhuo, and C. Zhong, Free-energy-based discrete unified gas kinetic scheme for van der Waals fluid, Entropy (Basel) 24(9), 1202 (2022)

[67] W. Zeng, C. Zhang, and Z. Guo, Well-balanced discrete unified gas-kinetic scheme for two-phase systems, Phys. Fluids 34(5), 052111 (2022)

[68] S. Tao, L. Wang, Q. He, J. Chen, and J. Luo, A sharp interface immersed boundary-discrete unified gas kinetic scheme for fluid-solid flows with heat transfer, Int. Commun. Heat Mass Transf. 139, 106424 (2022)

[69] A. Xu,Y. Zhang, Complex Media Kinetics, Science Press, Beijing, 2022

[70] A. Xu, D. Zhang, and Y. Gan, Advances in the kinetics of heat and mass transfer in near-continuous complex flows, Front. Phys. (Beijing) 19(4), 42500 (2024)

[71] Y. Gan, A. Xu, H. Lai, W. Li, G. Sun, and S. Succi, Discrete Boltzmann multi-scale modelling of non-equilibrium mul-tiphase flows, J. Fluid Mech. 951, A8 (2022)

[72] Y. Zhang, A. Xu, G. Zhang, Z. Chen, and P. Wang, Discrete Boltzmann method for non-equilibrium flows: Based on Shakhov model, Comput. Phys. Commun. 238, 50 (2019)

[73] C. Lin and K. H. Luo, Discrete Boltzmann modeling of unsteady reactive flows with nonequilibrium effects, Phys. Rev. E 99(1), 012142 (2019)

[74] Y. Gan, A. Xu, G. Zhang, and S. Succi, Discrete Boltzmann modeling of multiphase flows: Hydrodynamic and thermodynamic non-equilibrium effects, Soft Matter 11(26), 5336 (2015)

[75] Y. Zhang, A. Xu, G. Zhang, Y. Gan, Z. Chen, and S. Succi, Entropy production in thermal phase separation: A kinetic-theory approach, Soft Matter 15(10), 2245 (2019)

[76] G. Sun, Y. Gan, A. Xu, Y. Zhang, and Q. Shi, Thermodynamic nonequilibrium effects in bubble coalescence: A discrete Boltzmann study, Phys. Rev. E 106(3), 035101 (2022)

[77] D. Zhang, A. Xu, Y. Gan, Y. Zhang, J. Song, and Y. Li, Viscous effects on morphological and thermodynamic non-equilibrium characterizations of shock-bubble interaction, Phys. Fluids 35(10), 106113 (2023)

[78] Y. Gan, A. Xu, G. Zhang, C. Lin, and Z. Liu, Nonequilibrium and morphological characterizations of Kelvin-Helmholtz instability in compressible flows, Front. Phys. (Beijing) 14(4), 43602 (2019)

[79] H. Lai, A. Xu, G. Zhang, Y. Gan, Y. Ying, and S. Succi, Nonequilibrium thermohydrodynamic effects on the Rayleigh−Taylor instability in compressible flows, Phys. Rev. E 94(2), 023106 (2016)

[80] J. Chen,A. Xu,Y. Zhang,D. Chen,Z. Chen, Kinetics of Rayleigh−Taylor instability in van der Waals fluid: The influence of compressibility, arXiv: 2024)

[81] F. Chen, A. Xu, J. Song, Y. Gan, Y. Zhang, and N. Guan, Surface tension effects on Rayleigh−Taylor instability in nonideal fluids: A multiple-relaxation-time discrete Boltzmann study, Sci. China Phys. Mech. Astron. 67(12), 124611 (2024)

[82] J. Song, A. Xu, L. Miao, F. Chen, Z. Liu, L. Wang, N. Wang, and X. Hou, Plasma kinetics: Discrete Boltzmann modeling and Richtmyer−Meshkov instability, Phys. Fluids 36(1), 016107 (2024)

[83] M. Watari and M. Tsutahara, Two-dimensional thermal model of the finite-difference lattice Boltzmann method with high spatial isotropy, Phys. Rev. E 67(3), 036306 (2003)

[84] G. Gonnella, A. Lamura, and V. Sofonea, Lattice Boltzmann simulation of thermal nonideal fluids, Phys. Rev. E 76(3), 036703 (2007)

[85] H. Chen, Z. Hong, and D. Zang, New insights into suspended drops: When soft matter meets acoustic levitation, Droplet 3(1), e95 (2024)

[86] S. Xue, X. Xi, Z. Lan, R. Wen, and X. Ma, Longitudinal drift behaviors and spatial transport efficiency for spraying pesticide droplets, Int. J. Heat Mass Transf. 177, 121516 (2021)

[87] A. J. Wagner and C. M. Pooley, Interface width and bulk stability: Requirements for the simulation of deeply quenched liquid−gas systems, Phys. Rev. E 76(4), 045702 (2007)

[88] H. Liang, B. C. Shi, Z. L. Guo, and Z. H. Chai, Phase-field-based multiple-relaxation-time lattice Boltzmann model for incompressible multiphase flows, Phys. Rev. E 89(5), 053320 (2014)

[89] H. Liang, J. Xu, J. Chen, H. Wang, Z. Chai, and B. Shi, Phase-field-based lattice Boltzmann modeling of large-density-ratio two-phase flows, Phys. Rev. E 97(3), 033309 (2018)

[90] Y. Yang, M. Shan, Y. Zhang, H. Li, X. Kan, and Q. Han, Numerical investigation of acoustic cavitation behavior and cavitation-induced thermal effects using lattice Boltzmann method, Comput. Fluids 279, 106325 (2024)

[91] M. Shan, Y. Zha, Y. Yang, C. Yang, C. Yin, and Q. Han, Morphological characteristics and cleaning effects of collapsing cavitation bubble in fractal cracks, Phys. Fluids 36(6), 063337 (2024)

RIGHTS & PERMISSIONS

Higher Education Press