A many-body dissipative particle dynamics with energy conservation study of droplets icing on microstructure surfaces

Droplets icing has important applications in real life. The icing process of droplets on microstructure is explored based on the MDPDE method in this study. Firstly, the correctness of the heat transfer model was verified by one-dimensional heat conduction simulation, and then the feasibility of the phase change model was verified by investigating the icing process of droplets. The influence of cold surface temperature, droplet volume and contact angle on freezing time of droplets was discussed, and it was found that the temperature of cold surfaces had a greater effect on freezing. We finally explored the influence of different microstructure surfaces on the icing of droplets, and results showed that the presence of microstructures greatly enhanced the anti-icing effect of the surface. In our research, the contact angle is a relatively large factor affecting the icing of droplets. In addition, it was discovered that the droplet had the strongest ability to delay freezing on the surface of triangle microstructures with a contact angle of 157.1°.

impact process. Additionally, during the recoil of droplets, the radius of cylindrical surface has a greater influence on the diffusion coefficient of droplet in the circumferential direction, while a smaller influence on the axial diffusion coefficient. Among the freezing process of water droplets, changes in the temperature and radius of cylindrical surfaces did not cause significant influences on the shape of ice droplets. Next, Wang Y et al. [7] observed the freezing of droplets with different contact angles on surfaces of polished silicon wafers and micro-structured silicon wafers, and captured the initiation of ice nucleation of droplets. This investigation found that the free energy barrier of nucleation on graded surfaces was enhanced, which obviously delayed the nucleation of ice. The unique micro-micro-scale topological structure determined the dominant effect of wetting characteristics on ice nucleation. Pan Y et al. [8] studied the effect of surface wettability on the impact of droplets on surface and freezing. The results revealed that droplets spread less on hydrophobic surfaces and the oscillation time is longer. Because droplets contact the surface less, the heat transfer rate is lower, resulting in a delay in icing, and the ice time is significantly prolonged. On inclined surfaces, droplets show the process of sliding and stretching. Because of the increasement of contact area, the total ice time decreases at higher inclination angles. Yutaka Y et al. [9] demonstrated the freezing delay of pre-cooled water droplets on surfaces with four different wetting characteristics, and measured that the freezing delay time enhanced significantly with the increase of surface hydrophobicity. This is because the limited actual contact area between solid and liquid phases limited the possibility of heterogeneous ice nucleation. Reduction of the actual contact area of solid-liquid phase restricts the formation of ice nuclei. When the actual contact area of surface is constant, a surface structure with a smaller radius of curvature is required to further extend the freezing delay time. In addition, Ding B et al. [10] used high-speed cameras to visualize the dynamic behavior of water droplets on superhydrophobic surfaces with different inclination angles and subcooling degrees. It was found that increasing the surface inclination could promote the rebound process of droplets. With gradual decrease of surface temperature, the droplet undergoes complete rebound, partial rebound and no rebound on the surface. The boundary temperatures of three modes are − 31.25°C and − 33.75°C, respectively. Ju J et al. [11] explored the impact and freezing process of hot water droplets on ice surfaces. The water droplet temperature and ice surface temperature were parameterized. Experimental results show that when the temperature of ice surfaces is the same, the increased temperature of water droplets leads to an increase in the maximum diffusion factor and a slight increase in the freezing time. What's more, ice surface temperature has a small effect on the linear contact velocity, and it was also found that the temperature of water droplets has a significant effect on thickness of the thin layer during the expansion process.
On the other hand, in numerical simulation, Hindmarsh JP et al. [12] determined a simple heat balance model. The numerical model established can predict the temperature transition and freezing time of droplets. Subsequently, Chaudhary G et al. [13] conducted a thermal simulation by numerically solving the heat conduction equation based on enthalpy. The numerical calculation results of the droplet freezing are compared with experimental data, and results show that droplets freezing time is in good agreement with experimental data. Then Zhang X et al. [14] proposed a model to predict the icing of droplets. This model assumed that the freezing front was a spherical surface and took the effects of gravity and supercooling into account. Simulation results showed that the position of three-phase line, droplet volume and expansion rate were all in great agreement with experimental results. The geometric angle between the gas-liquid interface and horizontal plane tended to decrease most of the time during the solidification process, but suddenly increased at the end of solidification, forming a sharp tip with a sharply curved apex. Zhang X et al. [15] also used the volume of fluid (VOF) multiphase model and the solidification/melting phase transition model to establish a numerical model considering the influence of subcooling on the physical properties of droplet and the influence of dynamic contact angle on the contact line motion to simulate the impact-freezing behavior of droplets. During the impact freezing process of supercooled droplets on a cold hydrophobic surface, three different forms of complete rebound, partial rebound and complete adhesion were found, and a unified rebound and adhesive impact freezing behavior morphology was proposed. Additionally, Bodaghkhani A et al. [16] studied the freezing time of droplets with different wettability in horizontal and oblique directions. Using the twodimensional VOF method, results showed that the lower the surface temperature of droplets, the smaller the static contact angle and the larger the inclination angle, the freezing speed of droplets will be faster.
Despite there having been quantities of studies on droplets icing, there are still unsolved problems. This paper aims to study the icing of droplets on microstructures based on the energy-conserved many-body dissipative dynamics method (MDPDE). Three different microstructure surfaces are proposed, and the icing process of liquid droplets on a flat surface and three types of microstructure surfaces is simulated.

Nucleation theory
The nucleation of droplets is affected by many factors and is a complicated process. According to whether nucleation is spontaneous, it can be divided into homogeneous nucleation and heterogeneous nucleation [17][18][19]. According to classical nucleation theory [20], the nucleation rate J is: The kinetic parameter K satisfies: K 0 is a constant and is related to the concentration of nucleation sites, k B is the Boltzmann constant, T is the absolute temperature, and ΔG * is the work required for nucleation, which is expressed as: Where n * is the number of molecules or atoms in the nucleus, Φ * is the effective excess energy on the surface of nucleus, and Δμ is the chemical potential difference between the liquid phase and crystalline phase, which is expressed as [21]: P sw and P si are the saturated vapor pressures of water and ice, respectively.
In the case of non-uniform nucleation, the energy barrier for nucleation is lower than that of uniform nucleation. Therefore, the nucleation efficiency parameter ϕ is introduced to describe the influence of inhomogeneity [22]: Where θ is the contact angle between the ice and the substrate. Relative to homogeneous nucleation, heterogeneous nucleation requires a lower energy barrier. In most cases, droplet icing starts under heterogeneous nucleation.

Numerical method
MDPDE is a particle-based numerical method. Because of its ascendency in capturing micron-level mesoscopic features such as droplet splash during droplet collisions, this method has great advantages in simulating droplet collisions. So far, there have been quantities of related studies [23][24][25][26][27][28]. Many-body dissipative dynamics method (MDPD) is based on the original dissipative dynamics method (DPD) [29,30], and by modifying the original conservative force term that only contains repulsive force into a multi-body potential energy function, it solves the problem that a large number of gas-phase particles are required for the traditional DPD method to simulate gas-liquid two-phase flow. The temperature of the original DPD method is constant and there is no heat transfer. After adding an energy equation to the original equation of state to calculate the internal energy of each particle, the temperature simulation is achieved, and the energyconserved dissipative dynamics method (DPDE) is obtained. Combining the above two models, the MDPDE method is obtained, and the critical theoretical formulae are as follows: Where F C ij , F D ij and F R ij are the conservative force, dissipative force and random force of the system respectively. A and B are the attraction and repulsion coefficients between particles, ω c (r ij ) and ω d (r ij ) are their respective weight functions. Also, γ and ξ are the coefficients of dissipation force and random force, respectively, and they satisfy the thermodynamic fluctuation-dissipation theorem. θ ξ is a white noise term with zero mean unit variance.
The local density function is: For solid and fluid particles, according to the study of Wang L et al. [25], the solidliquid interaction function is: Energy-related parts: Where C v is the specific heat capacity, q V h ij , q C h ij and q R h ij are the internal energy change items caused by viscosity, collision and random thermal motion, respectively. Among them: When phase change occurs, the consumption of latent heat needs to be considered. The relationship between temperature, latent heat and internal energy is as follows: Where L and ε represent latent heat and internal energy, respectively.

Verification example one: one-dimensional heat conduction problem
Based on the above numerical model, the one-dimensional heat conduction problem was first verified. Relevant simulation parameters are shown in Table 1, and the symbol with l represents liquid phase, while s represents solid phase. As shown in Fig. 1a, a rectangular two-dimensional computational domain is established. Periodic boundary  condition is used in the x-direction, and the lower boundary in the y-direction fixed temperature is 0.8T R . Above the solid wall is a fluid area composed of liquid particles with a number density of 8.6 and temperature of 1.2T R . 7980 fluid particles are consumed and the total number of solid wall particles is 1353. Figure 1b is the temperature distribution in the calculation domain when t = 10t R . It can be seen that as time progresses, the temperature of the cold wall gradually spreads to the liquid. Figure 2 shows the temperature change with time in the y direction of the calculation domain. It can be seen that the numerical solution obtained based on the MDPDE model is in great consistency with the theoretical value. We count the errors of different time steps and find that they are all within the range of 1%, which verifies the correctness of the heat transfer model.

Verification example two: the droplets on the cold wall freeze
The initial temperature of droplets is 0°C in all simulation situations in our study. When the droplet is placed on a cold surface, it begins to nucleate when its highest temperature drops to T n , and then the entire droplet begins to condense.
Next, the icing process of droplets was simulated and compared with the experimental results of Zhang X et al. [31]. Results are shown in Fig. 3, and the same label of experiment and simulation corresponds to the same proportion of entire icing time. The size of droplets in our simulation is 20ul, the temperature of cold wall surfaces is maintained at − 18.4°C, and the contact angle of the droplet is 78°. In the experiment, the droplet experienced 11.7 s from the end of recalescence to complete freezing. The simulation result was 14 s, and the error (19.7%) between the two was within a reasonable range. Corresponding experimental and simulated icing processes are shown in Fig. 3. The sharp angle obtained by simulation is 148°, and the experiment is 140°. The error between the two is only 5.7%, which verifies the accuracy of the numerical results.

Icing on a flat surface
Based on the above verification results, the feasibility of the MDPDE method to simulate the icing process of droplets is proved. In this study, by means of changing the volume, contact angle of droplets and temperature of solid walls, the change of freezing time of droplets with them was explored. The icing process of droplets starts from the end of recalescence to complete freezing. The specific values of the droplet volume, contact angle and wall temperature are shown in Fig. 5 and Table 2. Figure 6 shows the result of freezing time of 10ul droplets with distinct contact angles on cold surfaces with different temperatures. It can be proved that the freezing time enhances with the increase of the droplet contact angle and solid wall temperature. When the solid wall temperature is higher, the variation in freezing time between droplets with diverse contact angles is more apparent. From Fig. 6a and b, it can be seen that under identical conditions, the larger the volume of droplets, the longer the corresponding freezing time. In the case of an equivalent contact angle, temperature of different cold wall makes a great divergence in freezing time. In summary, the contact angle has a relatively significant influence on droplets freezing time.

Freezing on the surface with microstructures
Eventually, the icing of droplets on surfaces with microstructures is studied. As shown in Fig. 7, three different types of microstructures are considered, namely square,  6 Freezing time of (a) 10ul (b) 20ul droplets on cold surfaces with different contact angles varies with wall temperature semicircle and triangular. The height and bottom width of these three patterns are kept identical respectively. Droplet volume in this simulation includes two cases: 20ul and 40ul. The temperature of cold wall surfaces is kept at − 15.9°C. The icing process of droplets starts from the end of recalescence to the end of droplets' complete freezing. Figure 8 is the change of freezing time of droplets on different cold surfaces with parameter A (the force between solid and liquid). Contact angles of droplets on plane and three types of microstructure surfaces corresponding to different parameter A are shown in Tables 3 and 4, and the maximum θ corresponding to each A is bolded. It can be seen that with the increase of the plane contact angle, the freezing time on four cold surfaces all improves, and the increase of droplet volume will lead to the enhancement of freezing time. It can be concluded from Fig. 8 that the presence of microstructures greatly enhances the anti-icing effect of a surface. For each solid-liquid force A, a maximum freezing time exists. According to the results, surfaces with the largest contact angle correspond to the longest freezing time in each case, which means the foremost determinant of freezing time of different microstructure surfaces is the contact  angle. In general, for a fixed solid-liquid force A, the larger the contact angle, the better the corresponding delayed icing effect, which is consistent with the experimental results of Liu Y et al. [32].
Taking the icing process of 40ul droplets on four types of surfaces under solid-liquid force A = 48 as an example, the temperature distribution inside the droplets during the icing process is shown. It can be seen from Fig. 9 that after adding the microstructure, the heat transfer between the droplet and the wall becomes slower, and the droplet on the plane completely solidifies into an ice phase at 18.9 s. On the other three types of microstructures, the liquid droplets were completely solidified in about 30 s. Next, further analysis of the results in each case is carried out.
Count the longest freezing time among the four surfaces corresponding to each fixed solid-liquid force A, and the result is shown in Fig. 10, where 1-6 on the abscissa correspond to contact angles θ = 157.1°, 135.1°, 134.9°, 126.7°, 114.6°, 102.6°(20ul), θ= 144.1°, 133.6°, 124.0°, 118.2°, 99.8°, 98.5°(40ul), respectively. Obviously, it can be known that under the premise of a single variable, either the contact angle or volume of droplets increases, the freezing time improves significantly. The greater the hydrophobicity, the more beneficial it is to prevent freezing. In addition, the droplet has the strongest ability to delay icing on triangular microstructure surfaces with a contact angle of 157.1°.

Conclusion
This paper explored the influence factors of droplets icing on a flat surface and the effect of the type of micro-structured surface on droplets icing. The results demonstrated the following: (1) The freezing time enhances with the increase of the droplet contact angle, droplet volume and solid wall temperature. The contact angle has a relatively significant influence on droplets freezing time; (2) The presence of microstructures obviously improves the anti-icing effect of a surface. The droplets freezing time of different microstructures mainly depends on the contact angle. Freezing time of the same surface increases with enhancement of the contact angle. Under the same solid-liquid force A, the contact angle of droplets on a microstructure surface is larger, so its effect of delaying freezing is better than the flat surface; (3) The droplet has the strongest capability to delay icing on triangle microstructure surfaces with a contact angle of 157.1°.
Based on the research in this paper, the hydrophobicity of a surface should be ensured as much as possible, and the ice-repellency of a cold surface can be increased by  adding microstructures on the plane. This research has crucial guiding significance for actual anti-icing problems.