Simulation of Boiling Process in a Nanochannel with 700000 Argon Particles by Molecular Dynamics

In this paper, the boiling flow inside a nanochannel with 700000 argon particle has been simulated by molecular dynamics (MD) simulation. This approach has been employed to analyse the superheated flow and its heat transfer pattern as well. For all simulations, an external thrust force varying from 1 PN to 12 PN is exerted on inlet nanoparticles along the channel to have the forced annular boiling flow. Computations reveal that increase of saturation condition and superheated region have significant impacts on the liquid-vapor interface. Furthermore, because of the major influence of surface tension throughout the nanochannel, the x-velocity of liquid film and vapor core has no considerable fluctuations and stay smooth. Results show the behaviors of what significantly similar with the author’s previous published simulation with 700000 solid particles.


Introduction
Nowadays, the majority of macroscale and nanoscale heat and mass transfer research is focused on boiling and phase change flows issues.Improvement of thermal energy efficiency of engineering processes is rooted in more vivid and clear understanding of these phenomena.Experimental measurements can considered as efficient means for determining the flow structure, analysis the boiling flow and describing the heat transfer characteristics in micro and nanochannels.Besides, huge progresses in computer science and computational techniques make an adequate possibility with more resolution as well as time consuming method to solve such programmes.
The extremely powerful technique of molecular dynamics simulation involves solving the classical many-body problem in contexts relevant to the study of matter at the atomic level.Since there is no alternative approach capable of handling this extremely broad range of problems at the required level of detail, molecular dynamics methods have proved indispensable in both pure and applied research [1].Molecular dynamics simulation is a powerful tool to capture the microscopic view of the nanoscale phenomena.Previously, molecular dynamics simulations have been carried out by a number of researchers to understand the molecular behavior of the phase transition from liquid to vapor.It has been demonstrated that molecular dynamics simulation is an efficient tool to investigate these complicated phenomena.The molecular dynamics study of these effects could provide information that is difficult to be obtained by experiments, and thereby improve our understanding in the technological processes such as flow in microporous media, and transport of biological and chemical materials through nanodevices or microdevices.Furthermore, Newton's law can be simulated and averaging of the molecular system provides some thermodynamical parameters including pressure velocity and temperature [1].
Most of the researchers' focus was on the heterogeneous phase transition of liquid from a molecular level flat solid surface.Phase transition from a solid surface can occur in a normal boiling mode or in an explosive boiling mode depending on the surface superheat temperature [2].The complete evaporation of a three-dimensional submicron droplet under subcritical conditions has been modeled by Long et al [3] using molecular dynamics.The two-phase system consisted of 2048 argon atoms was first allowed to relax to equilibrium, then the droplet was evaporated by increasing the temperature of the vapor phase atoms at the boundaries of the system until only the vapor phase remained.The computed evaporation rate agrees with that predicted by the Knudsen aerosol theory.
The thermal conductivity of solid thin films in the direction perpendicular to the film plane has been carried out by Lukes et al. [4] by the molecular dynamics computational technique.They also delineated the conditions necessary for meaningful thermal conductivity calculations and offered recommendations for efficient simulations as well.Their research demonstrates the potential of molecular dynamics for ascertaining microscale thermophysical properties in complex structures.Yi et al. [5] performed molecular dynamics simulations of the vaporization phenomenon of an ultra-thin layer (2 nm) of liquid argon on a platinum surface.The results reveal trends that agree with our knowledge of vaporization of a similar macroscopic system.For example, for the high surface temperature the vaporization process is reminiscent of the Leidenfrost phenomenon and after the formation of a vapor region between the surface and the liquid mass, the latter deforms and tends to approximately acquire a spherical "droplet" shape, as one would have expected from macroscopic considerations.Wang et al. [6] investigated the net evaporation process by performing non-equilibrium molecular dynamics (NEMD) simulations.They found the evaporation coefficients of the nonequilibrium molecular dynamics agree with those from the transition state theory and some data of equilibrium molecular dynamics calculations in literature within 15%.Toghraie and Azimian [7] investigated Molecular dynamics simulation of annular flow boiling of 3200 Argon atoms with the modified Lennard-Jones potential function.They found that saturation condition and superheat degree have great influences on the liquid-vapor interface.Also, the results show that due to the relatively strong influence of the surface tension in small channels, the interface between the liquid film and the vapor core is fairly smooth, and the mean velocity along the stream-wise direction does not change anymore.Morshed et al [8] investigated the boiling phenomena of thin liquid films adsorbed on a nanostructured solid surface.Their results indicated nanostructures play a significant role in both cases: Argon responds very quickly for the nanostructured surface, the transition from liquid to vapor becomes more gradual, and the evaporation rate increases with the nanoposts height.
In the present study, molecular dynamics simulations have been employed to investigate the boiling flow of simple Lennard-Jones fluids subject to an inlet driving force in a nanochannel.The main purpose of this research is to scrutinize the influence of modified Lennard-Jones potential function on the temperature, heat flux, velocity and surface tension profiles in order to comprehend the boiling flow behavior in nanochannel geometry.As far we know, no one has reported any investigation for these phenomena.

Simulation Procedure
The molecular dynamics simulation is employed to describe the forced annular boiling flow in a nanochcannel with 700000 argon atoms placed and simulated in cuboids of dimension where σ is the length scale parameter and it is included a solid platinum wall, liquid argon layers, and an argon vapor.Each atom is spaced 1.15σ from its neighboring atoms in all three directions.Figure 1 illustrates a schematic of MD domain simulation.According to this figure, argon molecules are divided into three regions: the vapor molecules are moving along the center of the channel and the liquid molecules comprise two thin annular films along the channel walls sandwich the vapor phase portion.At the bottom and top surfaces, wall boundary conditions are considered.Also, the system has doubly periodic boundary conditions in the x and y directions.The velocity of an atom crossing the x and y boundaries remains constant.An external thrust force varying from 1 PN to 12 PN is exerted on inlet nanoparticles confined between two triple-layer solid Platinum walls.

Figure 1. Schematics of MD domain simulations
Newton's second law as in equation ( 1) in all MD simulations should be solved for all fluid molecules.Before MD simulation, an integration method for the mentioned equation, the potential model and the simulation cell are computed.
The exerting force of each molecule can be described as the gradient of the potential function of a molecule . Lennard-Jones particles with the MD simulations lead to modified Lennard-Jones potential function in the NVT ensemble where N, V and T mean the number of particles in system, system's volume and temperature respectively, [9].Potential function presented by Stoddard and Ford [10] can be used to modeling the relations between fluid particles separated by a distance rij as follow: where the distance between the centres of i th and j th spherical particles is defined as , ε is a parameter corresponding to the energy and is called the maximum potential depth and rc=4.5σ is the cut-off radius.Thompson and Robbins [11] understood that decreasing rc to less than 2σ made considerable changes in the potential, boundary conditions and other parameters.Therefore, argon atom is used as the LJ fluid with the following potential parameters in the present work: presented by [12][13][14][15][16][17] are used as the solid wall.
The influence of solid walls on the velocity, surface tension and temperature distribution can be simulated by a bottom wall held at a superheated temperature, i.e.Tsat+ΔT and an upper wall held at a saturation temperature Tsat and the degree of superheat ΔT.The superheat degree not only makes a constant temperature gradient between the bottom and top walls in the z direction but also provides a driving force for heat flux.
Tsuruta and Nagayama [17] presented a few useful relations for modeling of the interaction between solid wall and liquid phase of flow as follow: where σsf and εsf can be computed according to the Kong mixing rule [18] as follow: 6 Verlet [19] presented the following equations to compute the particle's trajectories: The time step δt=5 fs has been employed in all simulations.Gaussian distribution based on the specific temperature can be used as the initial condition for each molecule as the following formula: NEMD simulation is employed when the system reaches equilibrium state after 500000 time steps.At least 2000000 time steps are used for sample data gathering for all simulations.Additionally, the velocity rescaling Vnew=λVold is first updated from the force acting on the particles and λ is defined as λ=(Ttarget/Tglobal) 1/2 where Ttarget is the target temperature and Tglobal is the global temperature of the fluid particles in the system.The instantaneous temperature is defined as a fluctuation of kinetic energy per particle as follow: where n(z) is the number of atoms at height z.
In the z direction, the system is confined by two stationary walls and periodic conditions are employed along the x-and y-directions as well.Based on the results of Alexander and Garcia [20], the thermal wall boundary condition is applied for particles near the walls.They showed when the liquid particle is near the wall surface, the repulsive force will push the liquid particle away from the wall surface, but some liquid particles may collide with the wall surface.Then the bottom wall at the temperature of Tsat+ΔT and all velocity components are reset to a biased Maxwellian distribution presented by [20]: In above equations kB is the Boltzmann constant, mi is the mass of each molecule, 3 ψ is a uniformly distributed number in range of (0,1), 1 ψ and 2 ψ are Gaussian-distributed random numbers with zero mean and unit variance.Also, when a particle collides with the upper wall at the temperature of Tsat , all three components of the velocities are reset to a biased Maxwellian distribution [20]: When the system reaches the steady state, the flow velocity, temperature, surface tension and density profiles are computed as a function of channel height.The program was run for 1500000 time steps for time averaging.In order to compute density, surface tension, liquid film velocity, vapor core velocity and temperature profiles, the domain is divided into Nbin,z=512 rectangular bins along the z-direction, each bin with height Δz=Lz/Nbin,z, where Lz is the box length along the z-direction.Also, the domain is divided into Nbin,x=12 bins along the x-direction, each bin with length Δx=Lx/Nbin,x, where Lx is the box length along the x-direction.Hence, the volume of each bin is , where Ly is the box length along the y-direction.For a large number of particles, more bins are needed.For adding molecules in each layer, Hn(zi) is defined as: where the subscript j represents the j th time step.Hence, the local average density in the n th slab from time step Istart to Iend is [21]: where m is the mass of each molecule and xi and zi is the coordinate of the mid-point of the n th slab.When sampling begins, the molecule number in each slab is accounted for and the local density can be considered.
At the atomic scale, the surface tension can be expressed as the integrated imbalance of normal and tangential pressure near the interface.The following expressions are used to calculate the two components of the pressure tensor presented by Haile [22], ( ) ( ) where the n(z) is the value of density at the height z, the i<j notation indicates a summation over all distinct pairs i and j without counting any pair twice (i.e. as ij and ji), B k is the Bolzmann's constant, denotes an ensemble average taken over the duration of the simulation for which thermal equilibrium exists.In terms of atomic position and liquid-vapor potential ϕ the expression for surface tension is presented by Haile [22], [ ]

Results
The influence of the saturation temperature on the density distribution according to the NEMD simulation has been illustrated in Figure 2. The highest magnitude of density is found adjacent to the solid walls of channel where a thin film of non-evaporative fluid is.On the contrary, the lowest value of density is observed at the centre of channel where the effects of solid walls as well as saturation temperature are insignificant.As can be seen, because of significant solid-liquid interactions, the highest and lowest values of density in the snapshots are related to saturation temperature 108 K and 145 K respectively.There is a continuous transition for density from the liquid to the vapor phase under all simulation conditions in a far distance from solid walls of channel.The behaviors of density distribution and interface structure are in good agreement with the previous study [23].With regards to the nanoscale of the space between the walls, the impact of the solid wall adsorption on the concentration of the liquid particles is very remarkable.As mentioned before, the bottom wall temperature and the top wall temperature can be considered as superheated and saturation temperatures respectively.Figure 3 illustrates the variation of density on the cross section area of the channel with regards to different superheated and saturation temperatures.It can be easily found that the bottom wall temperature has no significant impact on density distribution in particular at the centre of the channel where the stationary walls have negligible effect on flow structure and physical parameters as well.Therefore, density value and its distribution are independent of driving force and bottom wall temperature.The effect of different superheat degrees on temperature profile on the cross section area of the channel is illustrated in Figure 5 where the magnitude of temperature is computed from equation ( 12) and saturation temperature is considered a constant value of 108 K for all simulations.As can seen, the temperature magnitude reaches a value higher than saturation temperature of 108 K due to the superheat degree.On the plus side, the fluctuating domain of temperature becomes wider when superheat temperature rises.The temperature jump adjacent to the solid wall causes a rise of discontinues across the liquid-solid wall interface.Since, the local temperature of this kind of problem is described as a combination of velocity of atoms and transferred momentum between them, therefore, the motion of each atom can be influenced by other neighboring atoms significantly.Moreover, each atom is affected by other neighboring atoms of liquid more than neighboring atoms of walls.This behavior makes a sudden change on kinetic energy distribution at wall-liquid interface.The variation of x-velocity components of mean liquid film and mean vapor core for different external force has been illustrated in Figures 6 and 7, respectively.As it can be seen, the significant impact of surface tension throughout a nanochannel, the x-velocity of liquid film and vapor core has not considerable fluctuations and stay smooth.This trend is in good agreement with the majority of earlier works such as [9,24].
According to Figure 8, since ΔT is the driving force for heat flux, a higher ΔT results in a larger heat flux value, which is still in accordance with Fourier's law of heat transfer, qualitatively.These figures show that a higher external force results in a larger heat flux value.It could be explained that an increase in driving force will increase the fluid mobility and as a result, the heat dissipation will rise.This phenomenon could increase the heat flux in the simulation domain.Hence, it is concluded that the heat flux heavily depends on the degree of superheat and external force.On the other hand, the heat flux also depends on the boundary conditions.

Conclusions
Molecular dynamics (MD) simulations have been carried out to investigate the boiling flow in a nanochannel with 700000 particles in this research.The computational results reveal that the bottom wall temperature has no important impact on density distribution in particular at the centre of the channel where the solid walls have small effect on flow pattern and physical parameters too.Furthermore, the surface tension distribution varies remarkably in far distance from solid walls or at the centre of channel where the liquid-vapor interface zone is dominant.At higher temperatures and at a higher degree of superheat, the surface tension drops by increasing the temperature as well.
Additionally, as it can be seen, the significant influences of surface tension throughout a nanochannel and the interface between the liquid film and the vapor core, the x-velocity has not considerable fluctuations and stay smooth along the streamwise.
. Furthermore, three layers of face centered cubic (FCC) surface of Platinum molecules with parameters

Figure 2 .
Figure 2. Influence of saturation temperature on the density profiles

Figure 3 .
Figure 3. Influence of bottom wall temperature on the density profiles

Figure 4 .
Figure 4. Local surface tension at different superheat degrees at saturation temperature of 96 K

PreprintsFigure 5 .
Figure 5. Temperature profiles at different superheated temperature for the bottom wall and at saturation temperature of 108 K for the top wall

Figure 6 .Figure 7 .
Figure 6.Variation of the x-velocity component of mean liquid film.

Figure 8 .
Figure 8. Variation of heat flux versus bottom wall temperature for different external force The number i th and j th of particles