Effect of closed-loop nanochannels on the onset of explosive boiling: a molecular dynamics simulation study

As a primary boiling mode, explosive boiling has shown a promising future in many applications and received much research attention. The topology of the solid surface in contact with the liquid, particularly nanostructured surfaces, significantly affects the onset time of explosive boiling of a liquid nanofilm. Most studies investigated explosive boiling on non-closed-loop (parallel) nanochannel surfaces. Here, for the first time, explosive boiling in a closed-loop nanochannel was studied by the non-equilibrium molecular dynamics simulation method. Explosive boiling of liquid argon nanofilm on solid copper surfaces with different topologies, including an ideally smooth, a non-closed-loop, and a closed-loop nanochannel, was simulated. The results showed that, compared with the ideally smooth surface, the onset time of explosive boiling decreased for the non-closed-loop and closed-loop nanochannel surfaces. However, it turned out that compared to the non-closed-loop nanochannel, using the closed-loop nanochannel has an adverse effect on heat flux and the onset time of explosive boiling.


Introduction
There is significant scientific and practical value in investigating the phase transition processes of ultrathin liquid films on a solid surface.The surface temperature determines whether the phase transition occurs during normal or explosive boiling.Explosive boiling occurs when a liquid is heated above its saturation temperature but below its critical thermodynamic temperature; this causes a sharp pressure increase and the rapid ejection of liquid near the superheated surface, giving the phenomenon that looks like an explosion [1].
The utilization of explosive boiling can serve as an effective way to enhance the efficiency of heat transfer on heated surfaces.Therefore, enhancing the comprehension of this phenomenon is an imperative and crucial subject matter to improve thermal energy conversion efficiency.
Although nearly all solid surfaces display some degree of molecular-level roughness, creating multiple distinct nanopatterns on a flat surface is feasible by advancing nano-manipulation technology [2].Although explosive boiling of ultrathin liquid films has been experimentally investigated in recent years [3][4], a profound comprehension of physical phenomena based on macroscopic experiments is not fully understood owing to the intricacy of physical phenomena at the nanoscale, the limits on the time scale, and extremely high heat fluxes and temperatures.Therefore, this phenomenon must still be better understood, especially in nanoscale.Recent developments in hardware and software have made it possible to use computational methods to probe explosive boiling at the nanoscale.Molecular dynamics simulation (MDS) method has proven to be an effective tool for investigating the dynamics of explosive boiling on nanostructured surfaces in recent years.
Table 1 summarizes and tabulates the results of several MDS research on using parallel (non-closedloop) nanochannel surfaces in explosive boiling.Results have shown that using parallel nanochannels can decrease the onset time of explosive boiling [5][6][7][8] and increase heat flux [7][8].Therefore, the MDS studies proved that non-closed-loop nanochannel surfaces could drastically enhance the explosive boiling performance.To the best of the authors' knowledge, the effect of closedloop nanochannel surfaces on explosive boiling has not been investigated within available literature and remains to be unveiled.However, using closed-loop nanochannels could alter the absorption of heat flux by the liquid nanofilm.Hence, they should also be of interest from a scientific viewpoint.The main novelty of this study is the investigation of the effect of closed-loop nanochannels, together with the parallel nanochannels and the ideally smooth surfaces, on explosive boiling behavior of liquid nanofilm using the MDS method.

Simulation method
Four simulation cases, including one liquid-vapor argon coexistence system (simulation Case I) and three explosive boiling simulation cases with distinct surface topologies (an ideally smooth, a nonclosed-loop, and a closed-loop nanochannel), were studied to examine the model validation and impact of nanostructured surfaces on explosive boiling, respectively.The details of the simulation procedure, which are presented here, were similar for all but one.Simulation details for the exception, simulation Case I, are presented in section 3.1.
All simulations were carried out using the open-source MDS code LAMMPS (large-scale atomic/molecular massively parallel simulator) package [9].Due to the lack of visualization, the Open Visualization Tool (OVITO) was used for the post-processing analysis [10].
The interactions between argon (Ar) and copper (Cu) atoms were modeled by the Lennard-Jones 12-6 function potential (U ij ) as follows [11]: where r ij denotes the distance between interacting atoms (i and j), and ε ij and σ ij represent the energy parameter (the depth of the potential well) and the length parameter (the specified distance at which the interatomic potential is zero), respectively.Table 2 lists the energy and length parameters obtained from Refs.[12,13], the Lorentz combination rule [14], and a modified version of the Berthelot rule proposed by Din and Michaelides [15].In order to guarantee the calculation's precision and accelerate the simulations, the cut-off radius was adjusted at 9×10 -10 m.Moreover, to maintain the position of Cu atoms throughout the simulations, an individual spring force was applied to each atom; hence, the Cu atoms could vibrate around their initial lattice position.A force of magnitude -Kr was used at every single step, where K is the spring constant, and r is the distance between the atom's current location and its original location.The spring constant is closely related to Young's Modules; therefore, it is possible to estimate it as [16]: where E is Young's Modules (128.20 GPa at 300 K [17]), and d is its corresponding lattice constant (3.61×10 -10 m) of copper.The above data and correlation give the following value for the spring constant: K= 46.34 J/m 2 .

Atom pairs σ (m) ε (J)
Cu-Cu 1.93×10 -10 32.80×10 -21 Ar-Ar 3.40×10 -10 1.67×10 -21 Ar-Cu 2.67×10 -10 1.04×10 -21 As presented in Figure 1, three explosive boiling simulation cases in this work employed a rectangular simulation cell consisting of three different regions: a solid copper surface (with different topologies), a liquid argon region, and a vapor argon region.The solid copper surface was created of copper (in face-centered cubic (FCC) with a lattice constant of 3.61×10 -10 m which corresponds to the copper density of 8936.78 kg/m 3 ) and set at the bottom of the simulation boxes so that its 〈0 0 1〉 surface in touch with the liquid argon region.All the solid copper surfaces had the same dimensions in X-and Y-directions (shown in Figure 1).In order to better understand the influences of different topologies, a non-closed loop nanochannel (the surface P) and a closed-loop nanochannel (the surface C) with the same dimensions in nanowall width and height (3.61×10 -10 m and 9.04×10 -10 m, respectively) were constructed.The nanowalls were formed in the upper layer of the ideally smooth surface (the surface F).As can be seen in Figure 1, for all simulation cases, to create a more realistic substrate, the solid copper surfaces had three different regions from the bottom to the top, namely, the fixed region, the phantom region, and the conduction region.The first monolayer at the bottom of the surface (the fixed region) was kept stagnant to avoid any migration, penetration, and thermal deformation of the solid surface.The next four monolayers (the phantom region), whose temperature was regulated by a thermostat, acted as the heat source for generating heat flux.The remaining monolayers (the conduction region) were set as real atoms from which the liquid argon was heated, thus triggering explosive boiling process.The nanowall and ideally smooth surfaces are marked with different colours (light and dark green, respectively) to clearly show the topologies.The liquid argon region of the initial thickness of 40.01×10 -10 m, containing 2476 atoms, was placed over the solid copper surfaces.The liquid argon density was 1396.20 kg/m 3 at 87.18 K [18].The space above the liquid argon region was filled by 152 vapor argon atoms (with an initial thickness of 601.71×10 -10 m) with a density of 5.70 kg/m 3 at 87.18 K [18] to simulate a real phase change system.Considering this, the size of the vapor space in the Z-direction was selected to be large (much greater than the liquid argon film thickness) so that the top boundary had almost no influence on explosive boiling process.In addition, the periodic boundary conditions were applied in the X-and Y-directions, so an infinite plane of explosive boiling is simulated, and the non-periodic fixed boundary condition was applied in the Zdirection.This configuration complies with actual experiment procedures.
The simulation process was divided into four steps: energy minimization, equilibrium preparation, equilibrium relaxation, and non-equilibrium explosive boiling.Step I (energy minimization): at the beginning of each simulation, after the initial configuration and the interaction force field were prepared, using the conjugate gradient (CG) algorithm, the energy of the whole simulation system was minimized.The stopping tolerances for energy and force were taken as 1.60×10 -24 J and 1.60×10 -9 N, respectively.
Step Ⅱ (equilibrium preparation): in this step, the uniform temperature of the simulation systems was set at 87.18 K, which is the boiling point of liquid argon at 1 bar [18] under the Langevin thermostat for 2×10 -9 s.During this step, all simulation systems reached their equilibrium state.
Step III (equilibrium relaxation): in this step, the simulation systems were run for another 2×10 -9 s while the thermostat was changed into NVE (i.e., constant atom number, domain volume, and energy) ensemble for the fluid argon domain; however, the temperatures of the solid copper surfaces were still fixed at 87.18 K by the Langevin thermostat.
Step IV (non-equilibrium explosive boiling): in this step, the phantom monolayer atoms on solid copper surfaces were abruptly elevated to a higher temperature (250.00K), simulating ultrafast laser heating via Nose-Hoover thermostat, while the rest of atoms (Ar atoms as well as Cu atoms in the fixed and conduction monolayers) were free to interact with the NVE ensemble, and simulations were run for 1×10 -9 s.
Usually, the onset of explosive boiling is defined as the first appearance of a complete layer of vapour film with a specific thickness (20.00×10 -10 m in this work) over the solid surface.For the density studies, the simulation systems were divided into equal bins with a thickness of 2.00×10 -10 m in the direction of the phase change (Z-direction).The heat flux between the solid copper surface to Ar atoms was calculated by [19]: where A represents the cross-sectional heat transfer area between solid copper surfaces and Ar atoms, projected onto the X-Y plane (the nominal surface area), and E f stands for the total energy of whole Ar atoms.

Model Validation: Simulation Case I
In order to determine the reliability of computational accuracies, simulation model validation for a liquid-vapor argon coexistence system (the simulation Case I) was carried out in this work.Since density is the common evaluation property for validation of MDS method, this property was calculated and compared with the NIST-data [18], a famous fluid properties database.The methodology, results, and discussion are presented here.A liquid-vapor argon coexistence system, as shown in Figure 2, was developed to validate the simulation model for two-phase equilibrium.The liquid film size was 40.00×40.00×40.00(×10 -10 m) 3 (consisting of 1347 Ar atoms) with a 40.00×40.00×400.06(×10 -10 m) 3  long vapor domain (consisting of 55 Ar atoms) on each side.Then it was set at 87.18 K (the boiling point of argon at 1 bar [18]) through NVT ensemble (constant atom number, constant domain volume, and temperature canonical ensemble) for 2.5×10 -9 s.In order to extract the mean value of density distributions along the Z-direction, the system was divided into 210 bins with a height of 4.00×10 -10 m. Figure 3 demonstrates a good agreement between the liquid and vapour densities from the simulation and the experimental data.Figure 5 shows the unsmoothed and smoothed heat flux of the fluid argon domain over different solid copper surfaces.As can be seen, for all simulation cases, the heat fluxes reached the maximum as the solid temperature grew up continuously under the thermostat.Then they reduce to a relatively low value.The transition from the solid-liquid to the solid-vapor interface occurs during this time.When the vapor argon layer appears between the liquid argon cluster and the solid copper surface, it significantly increases the interfacial resistance to the heat flux (the Kapitza thermal resistance) because the Kapitza thermal resistance of solid-vapor is much higher than the solid-liquid one.Finally, heat fluxes reach almost zero after the departure of the liquid argon cluster and its vapor-liquid tail from the surface.Moreover, Figure 5 indicates that the non-closed-loop and closed-loop nanochannel surfaces show a higher value for the maximum heat flux for the maximum heat flux than the ideally smooth surface.This result is also compatible with the onset time of explosive boiling (Figure 4).The presence of nanochannels increases the contact area between the liquid argon nanofilms and the solid copper surface.Therefore, more liquid Ar atoms will be at the solid-liquid interface, which in turn excels in heat transfer from the solid surface to the fluid argon domain, decreasing the onset time of explosive boiling.Furthermore, the potential energy interaction between liquid Ar atoms and copper atoms is crucial in heat transfer at the nanoscale.The potential energy between the atoms of argon and copper is most concentrated close to the surface of the copper, especially in the nanowalls.Consequently, using nanochannels enhances the interaction of potential energy between the atoms of liquid Ar and the atoms of the copper surface.Thus, it features a significant thermal vibrational coupling, which leads to greater heat flux and an earlier onset time of boiling.
The simulation Cases P and C have an equal wetted surface area; however, interestingly the simulation Case C has a lower heat flux and longer onset time.Therefore, the onset time of explosive boiling cannot be determined solely by the value of the interaction potential energy, and there must be other reasons.Herein, this unexpected behavior is attributed to the space for the motion of Ar atoms (movement space).Therefore, for further analysis, the motion of Ar atoms should be studied via their line trajectories plot.
Figure 6 indicates the snapshots of line trajectories of liquid Ar atoms between the nanowalls for Cases P and C after 4×10 -10 s of beginning the simulation and before the onset of explosive boiling.Inspection of Figure 6 manifests that using closed-loop lead to a significant decrease in the mobility of Ar atoms.In other words, the movement of Ar atoms for simulation Case P is much better than the simulation Cases C. Therefore, other than the interaction potential energy, the effect of the movement space on the onset of explosive boiling is significant.When the solid-liquid interactions are very large, and on the other hand, the movement space is very small, the liquid Ar atoms cannot break their high potential energy barrier in a short time and move toward the liquid nanofilm bulk; therefore, accounts for the relatively long time required for the onset of explosive boiling.Therefore, the onset time of explosive boiling is delayed.Moreover, the lower movement of liquid atoms, the less they collide, so heat transfer occurs at a lower flux.

Conclusion
This study investigates the effect of using a closed-loop nanochannel on the nanoscale explosive bubbling process.By comparing the explosive boiling process on substrates with different topologies, including an ideally smooth surface, a non-closed-loop nanochannel, and a closed-loop nanochannel, it was found that the onset time of explosive boiling decreased for the non-closed-loop and closed-loop nanochannel surfaces when compared to the ideally smooth surface.In contrast to the non-closed-loop nanochannel, however, the use of the closed-loop nanochannel has a negative effect on heat flux and the onset time of explosive boiling.Based on the abovementioned results and analysis, closed-loop nanochannel surfaces should be avoided in the design of nanostructured surfaces.

Figure 1 .
Figure 1.Snapshot of the initial configuration of explosive boiling simulation cases.Purple and yellow represent Ar and Cu atoms, respectively.Brown, red, and light/dark green represent the fixed, phantom, and conduction region of the solid copper surfaces, respectively.

Figure 2 .
Figure 2. Snapshot of the initial configuration of the simulation Case I.

Figure 3 .
Figure 3.Comparison between density profile along the Z-direction (at 2 and 2.5×10 -9 s) for the simulation Case I and the experimental data.The experimental data are from Ref. [18].3.2.Effects of surface topology: Simulation Cases F, P, and C Figure 4 illustrates the snapshots and onset time of explosive boiling for liquid argon nanofilms on solid copper surfaces with different topologies.As can be seen, adding nanochannels could expedite the explosive boiling occurrence.Moreover, the non-closed-loop nanochannel surface (the simulation Case P) shows a lower onset time than the closed-loop nanochannel surface (simulation Case C).

Figure 4 .
Figure 4. Snapshots and onset time of explosive boiling on different surfaces.

Figure 5 .
Figure 5. heat flux for different surface.

Figure 6 .
Figure 6.Projections of selected Ar atoms and their trajectory motion lines for simulation Cases P and C.

Table 1 .
Literature overview of MDS studies of explosive boiling on parallel nanochannel surfaces.