A Novel Pressure-Controlled Molecular Dynamics Simulation Method for Nanoscale Boiling Heat Transfer

Abstract

Pool boiling, enabling remarkable phase-change heat transfer, has elicited increasing attention due to its ubiquitous applications in solar thermal power stations. An explicit understanding of the effect of system pressure on pool boiling is required to enhance the phase-change heat transfer. Despite its wide application when exploring the potential mechanism of boiling, the molecular dynamics method still needs to be improved when discussing the working mechanism of system pressure. Therefore, in the present study, a novel molecular dynamics simulation method of nanoscale pool boiling was proposed. This method provides a way to change and control pressure during the phase-change process. Furthermore, the bubble nucleation and growth in nanoscale pool boiling are quantitatively investigated through pressure-control molecular dynamics simulations. We expect that this study will improve the present simulation method of pool boiling and provide useful insights to the physics of the process.

1. Introduction

Solar energy is available in abundance in many parts of the Earth and has zero global warming potential, unlike thermal power plants where the main heat source can contribute to climate change [1,2]. Sonawane et al., analyzing a number of published studies, found that solar desalination plants powered by solar energy have increased to 76% [3]. The highest solar energy and exergy production still using black toner is 26.9% and 27.0%, respectively, higher than that of a conventional solar still [4]. Considering the unavailability of the solar source at nighttime and intermittency during the day [5], latent heat storage that offers great storage density per unit volume is a promising solution [6]. Therefore, an explicit understanding of boiling is required to enhance the phase-change heat transfer and the efficient operation of a whole power plant [7].

Boiling is a kind of thermal phenomenon in which liquid steam latent heat transfers heat at a low surface temperature [8], which is common in the field of industrial production and thermal management [9]. The enhancement of boiling heat transfer is of great significance because it can improve, for example, the efficiency and economy of heat engines and reduce energy consumption. Therefore, it has received extensive attention in energy and environmental protection-related fields [10]. When the surface comes in contact with liquid, boiling will occur if the temperature is higher than a certain threshold, while pool boiling will occur when the heating surface is lower than the free surface of the liquid [11]. With the increased cooling capacity of facilities, the cooling fluid parameters can be better optimized, which can improve pool boiling heat transfer performance [12].

Pool boiling is a typical two-phase flow phenomenon that includes the following stages [13,14,15]. First, single-phase convective heat transfer, followed by bubble nucleation. In the process of increasing superheat, the heat transfer efficiency continues to improve, and the heat flux reaches the critical value under certain conditions. Nucleate boiling changes to film boiling, which has a negative impact on heat transfer, and then leads to minimum heat flux in a short time [16]. After a stable film is generated near the surface, the heat flux increases again [17]. However, due to the influence of the gas film, the surface temperature is very high, which may be harmful to the equipment [18]. The characteristics of the nucleate boiling process are controllable and efficient, so it can meet the requirements of industrial production well. In order to improve the performance of the cooling system, it is necessary to increase the initial superheat of nucleate boiling, which is another research hotspot in this field [19].

At present, research on nanoscale boiling is increasing, and a series of achievements have been reached [20]. Cleanroom-processed micro- and nanostructures have been exploited for systematically studying structural effects on boiling heat transfer enhancement [21]. Some scholars have found that boiling heat transfer can be enhanced in different ways, mainly by modifying the heating surface through micro/nano structures (such as microporous layers [22,23], micro-fins [15,24,25], micro-channels [23,26,27], nanowires [28], and nanoparticles [29,30]). The micro/nanostructures discussed before can achieve one or more goals, such as increasing the heat transfer coefficient or critical heat flux. However, research results show that it is generally difficult for existing technologies to achieve all the above goals at the same time [31].

The adjustment of liquid supplement and bubble coalescence can significantly affect the heat transfer performance of such structures. In the traditional experimental mode, only the macroscopic analysis of the bubble dynamic behavior can be carried out, so the microscopic mechanism of the heat transfer process cannot be clarified. The boiling process is very complex, so it is difficult to accurately observe bubble behavior and changes [32]. During boiling, heat transfer performance is easily affected by roughness and wettability, which is unfavorable to the accuracy of experimental results. Due to the complexity of the mechanism and the difficulty of observation, the bubble nucleation process cannot be accurately studied by conventional experimental techniques.

Molecular dynamics (MD) simulation method is currently attracting attention in the heat transfer mechanism research [33,34]. Its advantage is that it can accurately analyze the microscopic behavior of molecular motion, thus providing support for the study of heat transfer enhancement [35,36]. Therefore, it is widely concerned at present. Some scholars have studied the nanoscale boiling and bubble nucleation process based on MD technology, and found that it can meet the application requirements well [37]. The MD technology can accurately control the size, wettability, and roughness of liquid surface structure, and improve the reference value of the results obtained. It can analyze the interface evolution in the nucleation process from a microscopic perspective, which provides support for the in-depth exploration of the nucleation mechanism [38]. Some scholars have studied the nucleate boiling phenomenon caused by the solid wall under heated liquid, and discussed the correlation between the surface wettability and explosive boiling when the liquid is heated. The results show that boiling intensity increases with increased wettability. The MD simulation results also found that the energy transfer rate efficiency would be greatly improved under the conditions of strong wettable surface.

Nevertheless, many problems still exist that need to be solved to clarify MD simulations. The boiling pressure affects the phase-change heat transfer performance. Bubbles are expected to exhibit various dynamic behaviors during different environmental pressures [39]. However, on reviewing the literature, the pressure effects has not been discussed in previous studies. Specifically, all the cited works that were found employed MD simulations that apply the fixed volume method, where the change of the fluid temperature induced by the heat transfer results in a change of the system pressure. Therefore, a novel pressure-controlled MD simulation method for nanoscale boiling heat transfer was proposed in the present paper. The main aim of the proposed method is to provide a way to change and control the pressure during the nanoscale phase-change process. Moreover, the bubble nucleation and growth in nanoscale pool boiling can also be quantitatively studied based on the proposed numerical method.

2. Computational Methods

2.1. MD Simulation System

All simulations were performed using the software package for the large-scale atomic/molecular massively parallel simulator (LAMMPS). Figure 1a below shows the initial configuration of this system, which includes a 16 nm (x) × 5 nm (y) × 60 nm (z) cuboid. The working gas is set as argon (Ar), which is used as the simulation medium in this experiment. The solid wall is made of copper, which has good thermal conductivity and is easy to process, and has been widely used in many engineering fields and experimental research. The bottom solid wall is composed of three layers of Cu atoms, which are arranged in a face-centered cubic (FCC) lattice structure with a lattice constant of 0.35 nm. Liquid Ar atoms are initially arranged on the surface of ordinary copper in the form of FCC, and the lattice constant is 0.58 nm. The thickness of liquid Ar film is 10 nm, and the top wall consists of two layers of copper atoms. In order to facilitate simulation analysis, periodic boundary conditions are set in the vertical and horizontal directions of the simulation box.

2.2. Simulation Method and Procedure

In this paper, the interactions between argon–argon (Ar–Ar) atoms, copper–copper (Cu–Cu) atoms and argon–copper (Ar–Cu) atoms are expressed as follows [40]:

$$E_{\mathrm{Ar}–\mathrm{Ar},\mathrm{Cu}–\mathrm{Cu},\mathrm{Ar}–\mathrm{Cu}}=4\epsilon \left[{\left(\frac{\sigma }{r}\right)}^{12}-{\left(\frac{\sigma }{r}\right)}^6\right],r<r_c$$

(1)

where ε represents energy, σ represents distance; r is inter-atomic distance, and r_c is the cut-off radius. The interactions were determined using Equations (2) and (3) in terms of the combined Lorentz-Berthelot rule:

$${\epsilon }_{\mathrm{Ar}–\mathrm{Cu}}=\alpha \sqrt{{\epsilon }_{\mathrm{Ar}–\mathrm{Ar}}{\epsilon }_{\mathrm{Cu}–\mathrm{Cu}}}$$

(2)

$${\sigma }_{\mathrm{Ar}–\mathrm{Cu}}=\frac{{\sigma }_{\mathrm{Cu}–\mathrm{Cu}}+{\sigma }_{\mathrm{Ar}–\mathrm{Ar}}}{2}$$

(3)

The copper–argon interaction strength can be adjusted based on the energy parameter of α in the present research. Energy parameters used in the present study are detailed in

Table 1. Potential energy parameters.

Inter–Particle Interaction	Potential Energy Parameter
	ε/eV	σ/Å
Ar–Ar	0.0104	3.405
Ar–Cubottom	Variable	2.87
Ar–Cuup	0.02	2.87
Cubottom–Cubottom, Cuup–Cuup	0.406	2.338
Cubottom–Cuup, Cubottom–Cuup	0	0

.

In the simulation research, the large-scale atom/molecule parallel simulator is applied, and the atomic parameters are updated based on the velocity Verlet algorithm. In order to meet the temperature control requirements, the Langevin thermostat is applied. The phase-change heat transfer process can be divided into equilibrium stage and wall temperature-change stage, which are simulated separately. First, the temperature of the whole system is controlled at 80 K and kept in balance when it lasts for 0.5 ns. Then, the temperature of the bottom wall is heated appropriately to reach 145 K. The microcanonical ensemble (NVE) is adopted and the simulation lasts for 5 ns. The data is acquired simultaneously in these two stages. The data are output every 1000 time-steps, and the atom trajectories are visualized by the open-source visualization software OVITO [41].

The absorbed heat flux q can be determined as follows [42]:

$$q=\frac{1}{V}\left[{\displaystyle \sum _{i=1}^N\left(v_i\cdot e_i+\frac{1}{2}{\displaystyle \sum _{j=1,j\ne i}^N{r}_{ij}\left(F_{ij}\cdot v_i\right)}\right)}\right]$$

(4)

where V is the volume of liquid computation domain, N denotes the number of liquid argon atoms, v_i and e_i represent the thermal velocity and the total energy (including kinetic energy and potential energy) of atom i, respectively, and F_ij refers to the force acting on atom i from atom j. To constrict the noise influence associated with oscillations in the original data, the exhibited value of q in the present paper is the average of the data every 10,000 timesteps.

2.3. Pressure Control Method

Usually, in the MD simulation of pool boiling, the size of the simulation box remains unchanged; however, during the boiling process, the liquid will undergo a gas–liquid phase change process under the heating of the wall, and the system pressure will change if the simulation area remains unchanged, which will influence the whole boiling process. Therefore, a systemic pressure-control method is proposed herein to keep the pressure in the simulated area constant in the boiling process. Specifically, a control wall is set above the simulation area, as shown in Figure 1. In the simulation process, firstly, the force on the solid wall in z-direction is set to 0, and then an external force is applied to each solid particle on the control wall. The external force value is calculated by a given system pressure, as shown in Figure 1b. The applied external force in -z-direction is balanced with the pressure applied by the system to the pressure control wall in z-direction, to keep the whole system pressure constant and avoid the influence of pressure change on the bubble dynamics.

Considering the wall-fix method can only simulate the atmospheric pressure, 1 atm was set in the pressure-control method. The force applied on every single atom of the pressure-control wall F_atom is calculated by

$$F_{\mathrm{atom}}=\frac{F_{\mathrm{total}}}{N_{\mathrm{up}}}=\frac{p\times S_{xy}}{N_{\mathrm{up}}}$$

(5)

where F_total is the total force applied on the pressure-control wall, N_up is the number of atoms on the pressure-control wall, p denotes the pressure set in the simulation, and S_xy is the cross-sectional area in the x-y plane. Moreover, the wall-fixing method was also calculated and compared setting and freezing an upper wall to the effect of phase-change on the system pressure. Different lengths of the simulation box in the z-direction were set and computed.

2.4. Quantitative Description of Bubble Dynamic Behavior

To quantify the changes in bubble shape and size with time in the boiling process, a bubble-related statistical method was used. Statistical calculations are divided into three steps. First, the calculation file is exported, as shown in Figure 2a; the area is calculated and divided into cubic bins along the three principal directions using a custom MathWorks MATLAB script. Second, the number of argon molecules in each bin is counted, and the bin is marked with the number of molecules less than a certain value (Figure 2b). Third, the processed file is exported, the marked bin is read in OVITO, and the bubble shape and volume are calculated (Figure 2c).

3. Results and Discussion

3.1. Visualization of the System Equilibrium Stage

Four different kinds of boiling systems were constructed and simulated in the present study, which can be seen in Figure 3. The x-z plane view is exhibited to observe the interface evolution process. The pressure-control method was applied in the pressure-control case, shown in Figure 3a. The initial height of the simulation box H in case (a) equals 433 Å. The other three cases apply the wall-fix method, which means the simulation box remains unchanged during the simulation process. The heights of the simulation box H in cases (b), (c), and (d) are 216 Å, 433 Å, and 650 Å, respectively. There are two equilibrium processes in the pressure-control case. The liquid is controlled at 80 K to reach an equilibrium state. Under the given system pressure, atoms of the upper control wall are subject to external forces F_atom, which can be calculated by Equation (5). Therefore, to maintain the set system pressure, the upper control wall is observed to move in the -z-direction in Figure 3a. As the equilibrium process continues, the upper wall reaches a certain position and remains stationary. Finally, both the liquid and the system box reach the equilibrium state at 2500 ps. In contrast, there is only one equilibrium process in the wall-fixing cases. The box heights parameter is kept constant and only the liquid shifts to the equilibrium state.

3.2. Bubble Dynamic Behaviors of the Phase-Change Stage

Under the condition of entering equilibrium state, the temperature of the bottom wall increases from 80 K to 160 K to simulate the heating process. Thus, the liquid is heated and the phase-change phenomena can be observed. Representative snapshots of the four cases during the phase-change process are chosen and illustrated in Figure 4. This stage starts at 2500 ps, which is the time marking completion of the first equilibrium stage.

The simulation box height keeps increasing with time during the phase-change process, as shown in Figure 4a. The heights of the other three cases remain constant, as shown in Figure 4b–d. These results are consistent with the analysis in Section 2.3. The system pressure increases with the continuous occurrence of the vapor–liquid phase change process. The upper control wall moves upwards as a result. The model can also be verified by combining the results of Figure 3a and Figure 4a. Specifically, the initial height of the simulation box is 433 Å, and the system pressure is set as one bar. As the system moves through the equilibrium stage, the upper wall moves downwards because the argon atoms are in the liquid state under 80 K and there is no argon vapor at this stage. In contrast, when the argon atoms absorb heat from the contact wall, the phase-change process occurs and the volume of the vapor increases. As a result, the upper control wall moves upwards. The upper wall moves along with the system pressure changes, which verifies the pressure control method.

The bubble dynamic behavior was also investigated. As illustrated in Figure 4a, a suitable bubble nucleus is produced on the wall surface at about 4450 ps. This indicates that the tin is 1950 ps. According to the experimental results, small vesicular embryos also appeared at the base center. However, their size is small and does not reach the critical bubble core size, so its stability is poor and it will disappear after a period of time. In this process, the stable bubble nucleus continues to grow upward and laterally, finally merging to form a large bubble and generate a stable vapor film at 4950 ps.

The effects of the simulation box height on bubble dynamic behavior were also examined. Figure 4b shows bubble nucleation at 4100 ps; however, the box height is relatively small for the vapor-liquid phase change. No room remained for subsequent bubble growth. Both cases c and d were observed to have undergone complete bubble nucleation, bubble growth, and bubble detachment on the heated surface, as marked in Figure 4. Nevertheless, different droplet dynamic behaviors were observed. Bubble nucleation was observed at 4800 ps and 4900 ps in cases c and d, respectively, which were later than the time when they were observed in the pressure-control case. The bubble detachment time also varied in different cases (at 4950 ps, 5350 ps, and 5600 ps). The height of the simulation box of the wall-fix method affected the bubble dynamic behavior. Furthermore, it is very necessary to research the mechanisms of the interface evolution process, which will be discussed later in this paper.

3.3. Quantitative Study of Bubble Dynamics Behavior

In Section 3.2, the boiling heat transfer over the copper surface is exhibited intuitively in terms of representative snapshots. The snapshots in Figure 4 show that the simulation box height of the wall-fix method affected bubble dynamic behavior, which was also different from the results of the pressure-control case. To reveal the potential mechanisms of this process, more quantitative studies of bubble dynamic behavior are examined herein.

The bubble volume was determined to quantify the bubble nucleation process. The growth curves of bubbles on surfaces with and without pressure control are presented in Figure 5a,b. The bubble growth process includes two stages based on its volumetric variations: the nucleation and growth stages. In the first stage, the size and volume of bubbles are increasing. When heat is introduced, the stable bubble nucleus begins to enter the rapid growth stage; under this condition, the inflection point of the bubble growth curve corresponds to t_in. When the solid–liquid interaction intensity is low, the nucleation process is slow. After that, the bubbles begin to grow rapidly and maintain a certain stability.

The liquid near the bottom wall is heated by coming in contact with the substrate, which is closely related to energy transfer. Therefore, the liquid here is of great significance in the study of the boiling heat transfer mechanism. Figure 5c shows the change trend of liquid film temperature near the substrate during the time increasing process (z = 20–40 Å). The curves include three stages, marked in Figure 5a: Stage I refers to the equilibrium stage of the simulation. The temperature near the heated wall is balanced at 80 K in this stage. The bottom wall is heated at 2500 ps and the phase-change process begins. The liquid temperature increased quickly in the initial stage. Then in Stage II, horizontal fluctuation appeared. After the bubble grew and detached from the heated wall, there was no heat transfer between the near wall region and the bottom wall. The temperature fluctuated horizontally, and this period of activity is referred to as Stage III.

To compare the heat transfer of the four cases, the absorbed heat flux q in the near wall region of liquid argon was determined. The heat flux curves with and without pressure control are illustrated in Figure 5d. In the initial period, the heat flux is basically kept constant during this process. When the liquid is continuously heated, a group of liquid atoms absorb more heat, and bubble nucleation will occur when a certain critical value is reached. After that, the bubble volume increases continuously, and its interface with the liquid needs to absorb energy to evaporate, which will greatly increase the heat flux. It also shows that nucleate pool boiling can promote heat transfer on the solid surface, and it has the advantages of high efficiency and safety.

3.4. Mechanisms of Nucleate Boiling with Pressure-Control

3.4.1. Illustrations of Nucleate Boiling from the Perspective of Energy

To quantify the bubble nucleation of the pressure-control case, the energy was computed. The simulation domain was divided into bins with a size of 0.5 Å (x) × 0.5 Å (y) × 0.5 Å (z). The number density, E_k and E_p of argon atoms in each cell were determined. The mean total energy E of atoms can be calculated by adding E_k and E_p. According to the above discussion results, on the macro scale, the energy barrier of liquid gas phase transition belongs to latent heat, which can be overcome to produce phase transition after continuous heating and absorption of sufficient heat. According to the same idea, the macro energy theory and the micro kinetic energy related theoretical tools are jointly applied; when the average E_k of a group is higher than the corresponding E_p, bubble nucleation will occur. It is noteworthy that the E_p is less than zero, which reflects the restriction of interactions between atoms, and the value of E_k and E_p is affected.

The pressure profile change at the initial time of the non-equilibrium stage is shown in Figure 6 below. Analysis of this figure shows that, at the initial nucleation time, and the moment of bubble detachment, the initial positions of the solid-liquid interface are z = 0 Å. In the nonequilibrium phase, the density and energy in the corresponding liquid region basically remain uniform, with few fluctuations. Previous studies [22] have found that a small number of liquid atoms are always crystallized on the surface of hydrophilic substrate due to the interaction of solid and liquid atoms. In this case, the density of liquid atoms near the substrate is improved. In the analysis of boiling heat transfer mechanism, such liquid atoms are the focus of research.

While time increases, the liquid film in the near wall area continuously absorbs heat. Some of the heat energy is converted into atomic kinetic energy, which will increase the average kinetic energy of the liquid atoms in this area, and some of the heat is converted into atomic potential energy, which will lead to the decrease of E_p, thus reducing the potential barrier of the horizontal potential energy. According to the above results, it can be judged that the average kinetic energy of liquid atoms in the near wall area will increase significantly during this process, which will lead to nucleating boiling phenomenon as shown in Figure 6d. The comparative analysis also shows that the bubble shape in the area with total energy greater than 0 eV is very similar to that in Figure 4 and Figure 6a.

When the bubble nucleus appears in the near wall liquid, the potential energy barrier at the interface decreases significantly. In the process of continuous heat transfer, the kinetic energy of the liquid atoms at the interface of the bubble core continues to increase, along with the number of atoms overcoming the weak barrier, which also promotes the formation and growth of bubbles. After that, the bubbles enter a rapid growth stage and maintain certain thermodynamic stability. The heat transfer efficiency has a significant effect on the bubble growth rate. According to the above discussion, the mechanism of bubble nucleation and growth is closely related to energy, so it can be explained from an energy perspective.

3.4.2. Separate Energy Analyses in the Near-Wall Region

The atoms between the solid and liquid walls have a strong effect, which makes the atoms in them unable to evaporate due to a large force. Under the influence of this factor, the liquid layer is always confined to the substrate surface during the heating process. This also reflects the fact that solid-like liquid atoms do not play a role or participate in the bubble nucleation process. In this way, when analyzing the nucleation energy, the liquid film in the wall area (z = 0–20 Å) can be divided into two unrelated areas: a solid-like liquid layer (from z = 0–5 Å), and a bulk-liquid layer (z = 15–20 Å). See Figure 7a below for relevant information. Based on the analysis of the collected images, it is estimated that the thickness of such a liquid layer is 5 Å. The kinetic energy and total energy of these two regions as a function of time were further calculated separately and presented in Figure 7b–d.

The body liquid layer near the wall is the key object of this analysis. Theoretical analysis shows that the heat transfer here is mainly related to the collision between solid and liquid atoms, and there is a positive correlation between the heat exchange efficiency and the collision probability. When the density of liquid atoms is large, the corresponding collision frequency also increases. There is a positive correlation between atomic density near substrate surface and collision frequency. According to the above discussion, it can be concluded that the solid liquid layer plays an important role in this interface heat transfer, which needs to be analyzed emphatically. The contrast analysis shows that the density of the solid liquid layer is higher, as shown in Figure 6a.

When the temperature of the bottom wall increases to the target value, the absorbed energy is mainly converted into atomic kinetic energy, which can increase the temperature of the liquid. However, during this process, the energy of the solid liquid layer has basically not changed. The liquid layer near the wall absorbs the most heat energy. It should also be noted that in various heating stages, the heat absorbed by the bulk liquid layer is mainly consumed to overcome the energy barrier, which also indicates that the near wall liquid layer should be the focus when analyzing the nucleation mechanism.

4. Conclusions

In summary, in this study, the molecular dynamics simulation method for nanoscale boiling heat transfer was modified by inducing a pressure-control wall. The pressure-control wall can move upwards or downwards according to the pressure of the simulation domain. Thus, the system pressure during the nanoscale phase-change process is able to be set and maintained at a given value. Moreover, the bubble nucleation and growth in nanoscale pool boiling were quantitatively investigated based on the pressure-control method. Mechanisms behind the nanoscale nucleate boiling were elaborated from the aspect of energies. This research is expected to improve current methods of pool boiling simulation. Moreover, the effects of the system pressure on bubble behavior and heat transfer of the nanoscale boiling can also be explored based on the proposed pressure-control method, which will be examined in future research.