Recently,nanofluids have received enormous attention due to their potential applications in many commercial and industrial fields,e.g.,cosmetics,textiles,pharmaceutical,catalysts,and electronics. It is important to evaluate the behavior of nanofluid flow before experiments. Since the testing of nanofluid flow inside nanochannels is expensive and time-consuming,many the- oretical and numerical models are developed to evaluate and predict the behaviors and risks of nanofluid flow in various systems^{[1, 2]}. Due to the lack of systematic experimental data,it is very difficult to analyze the relative significance of the mechanism caused by nanofluids. Molecular dynamics simulation (MDS) is a powerful tool to simulate and investigate nanofluid behaviors at nanoscale channels. MDS is a deterministic method. It can track atomic trajectories of many particles by numerically solving Newton’s equation of motion for a specific interatomic potential with certain initial conditions and boundary conditions^[3]. MDS is not only a well established method for predicting the transport properties of nanofluids^{[4, 5]} but also a suitable way for observing the microscopic motion behaviors of nanoparticles at the molecular level and analyzing the enhanced heat transfer mechanism of nanofluids^{[6, 7]}.

MDS has been successfully used to calculate the properties of nanofluids^{[8, 9, 10]}. Kang et al.^[8] examined the nonequilibrium heat transfer in a copper (Cu)-argon (Ar) nanofluid by MDS. Two different methods were introduced to calculate the coupling factor between the nanoparticles and the base fluid,and the obtained coupling factors were consistent. Mohebbi^[9] performed a combined equilibrium molecular dynamics (EMD) and non-equilibrium molecular dynamics (NEMD) simulation to calculate the specific heat and thermal conductivity of the nanofluid systems. The obtained results showed that the effect of nanoparticles was not significant on increasing the thermal conductivity of nanofluids. Kondarajau et al.^[10] developed an Eulerian- Lagrangian based direct numerical simulation model to investigate the effective thermal con- ductivity of nanofluids,and suggested that it was necessary to act the particle conductivity and forces on nanoparticles. Vladkov et al.^[11] used MDS to model the thermal properties of nanofluids. The effects of the interactions among the particle mass,particle-particle,particle- fluid,and spatial distribution of the particles on the thermal conductivity were considered. The results showed that the heat transfer increased when particles were aligned in the direction of the temperature gradient.

Simple fluid problems have been extensively studied by MDS. Sofos et al.^{[12, 13]} presented the NEMD simulations of the planar Poiseuille flow of liquid Ar,and investigated the density, velocity,temperature,and strain rate across the channel whose width was in the range from 0.9 nm to 17.1 nm. Sofos et al.^[14] investigated the liquid flows in various nanochannels,and presented the effects of all physical parameters affecting the slip length at the boundary by use of MDS. The results showed that slip existed at nanoscale,and the slip was significantly affected by both the surface and the fluid characteristics. Priezjev^[15] investigated the dynamic behavior of the slip length in a fluid flow confined between atomically smooth surfaces. The results showed that the slip length could be well described by a function of a single variable that in turn depends on the in-plane structure factor,contact density,and temperature of the first fluid layer near the solid wall.

The mechanisms and microstructures of nanofluids have also been extensively studied by MDS. Li et al.^[16] investigated the molecular layering at the liquid-solid interface in a nanofluid by equilibrium MDS. It was found that an absorbed slip layer of liquid was formed at the interface between the nanoparticle and liquid. This thin layer would move with the Brownian motion of the nanoparticle. Evans et al.^[17] used a kinetic theory based analysis of heat flow in nanofluids to demonstrate the hydrodynamics effects associated with the Brownian motion. The results showed that the Brownian motion had only a minor effect on the thermal conduc- tivity of the nanofluid,while some effects,such as particle clustering,were responsible for the large thermal conductivity increases observed in experiments. Merabia et al.^[18] explored the phenomena of heat transfer in the vicinity of strongly heated nanoparticles by use of MDS of atomically realistic models or of more coarse-grained Lennard-Jones (LJ) monoatomic fluids. Physically,the possibility of sustaining very high heat fluxes was related to the strong curvature of the interface that inhibits the formation of an insulating vapor film. Lv et al.^[19] simulated an MD model of nanofluids between the flat plates under shear flow conditions. The phenomena observed during the shear flowing process were that the nanoparticles would vibrate and rotate besides the main flowing with liquid Ar and these micro-motions could strengthen the partial flowing in nanofluids. Lv et al.^[20] investigated the impact and friction model of nanofluids by MDS which consists of two Cu plates and one Cu-Ar nanofluid. The motion states and inter- action of nanoparticles at different time through the impact and friction process were obtained and the friction mechanism of nanofluids were analyzed. Cui et al.^{[21, 22]} performed the MDS on the Couette flow of nanofluids. It was found that the even-distributed liquid Ar atoms near the solid surfaces of nanoparticles seemed like a reform to the base liquid and was contributed to the heat transfer enhancement. When the shearing velocity increased,the motion of nanoparticles was strengthened significantly. The analysis of the reviewed literature shows that there still is no conclusive theory concerning the flow behavior of nanoparticle suspensions on the Poiseuille flow^[23],in which the external force drives the nanofluid.

In the present work,a novel nanoscale MDS is presented for the study of nanofluid behaviors on the Poiseuille flow. This study aims to model and simulate the evolution process of three- dimensional (3D) microstructures in the nanofluids under an external force with the Poiseuille flow condition. More importantly,the evolution of energy under different conditions based on the microstructures is analyzed. The MDS visualization is provided at different time steps to fully describe the kinetics of aggregation. The simulation results show that the clustering effect occurs through the nanofluid flow under an external force. The density and velocity profiles of Ar-Cu nanofluid flow and simple Ar flow under the Poiseuille condition are shown. An asymmetric ordering of the base fluid flow is observed due to the existence of nanoparticles. Furthermore,the enhanced velocity of the base fluid particles is due to the nanoparticle random and fluctuating motion. Moreover,the slip length changes due to the existence of nanoparticles. The slip length of the nanofluid flow increases compared with that of the simple flow in the nanochannel. This study is important in the cognition of nanofluids under an external driving force,which makes the nanofluids far from the equilibrium. 2 Simulation technique and details 2.1 Methodology

We perform the MDS of the Poiseuille flow in a model nanofluid comprising of four nanopar- ticles. On the basis of the molecular dynamics method,the Newtonian equation for each fluid molecule can be written as^[24]

Click Browse formula

where the subscript i represents the particle i,and 1 is the unit vector in the x-coordinate. The first term of the right-hand side of Eq. (1) is the molecular force due to any type of LJ potential between the particle i and other fluid molecules in the calculation domain. The second term is the molecular force between the particle i and all the solid wall or the nanoparticle j. The third term is the molecular force between the particle i and the nanoparticle j_np. The last term represents the external force,such as the wall shear force,the gravity force,and the electric force,which makes the fluid deviate from the equilibrium.

When a two-body potential model is applied,the Newtonian equation for each molecule of the simulation model follows this relation^[24]:

Click Browse formula

where ¢(r_ij) is the interaction potential function between particles i and j,and F_ext is the inlet external force per molecule. When a force is applied on an atom,its potential energy will change due to the movement in the applied force field.

Suitable potential functions are important factors to ensure the MDS results accurate and reliable. Empirical or semi-empirical correlations are adopted in most classic MDSs^[21]. The LJ potential is the most commonly used and simple two-body interaction potential. It is given by

Click Browse formula

where r_ij represents the intermolecular distance between atoms i and j,ε is the binding energy (depth of the potential),and σ is the molecular diameter (the distance at which the interparticle potential is zero). Both ε and σ depend on the type of the used molecules. The LJ model attempts to account for both short-range repulsive overlap forces and longer range attractive dispersive forces. The first term represents the short-range repulsive interactions preventing the overlap of the molecules,while the second term represents a dipole-induced attractive interaction. The LJ parameters for Ar and Cu are listed in Table 1^[3].

Table 1 LJ potential parameters

Table options

The parameters for Ar-Cu atoms are calculated according to the Lorentz-Berthlot mixing law^[4] given by

Click Browse formula

The schematic of the Poiseuille flow comprising four nanoparticles confined between two parallel walls is depicted in Fig. 1.

Fig. 1 Schematic of nanofluid Poiseuille flow through nanochannel

Figure options

The fluid slippage at the solid-liquid interface is an important hydrodynamic characteristic in dealing with nanoscale flows. The slip length L_s,as a measure of the fluid slippage,is defined in Navier’s slip formula^[25] as follows:

Click Browse formula

where U_s stands for the slip velocity,and (dU/dn_w) stands for the wall velocity gradient in the normal direction n. Figure 2 shows the concept of the slip length in the Poiseuille flow inside the nanochannel.

Fig. 2 Concept of slip length in Poiseuille flow

Figure options

Dimensionless physical quantities (LJ reduced units) are used for convenience and to facilitate the scale-up of the system. Table 2 lists the units used in our MDS,where k_B is the Boltzmann constant. Hereafter,all quantities will be given in terms of LJ reduced units.

Table 2 Units used in MDS

Table options

For computational efficiency,we select a cutoff distance r_c = 2.5σ. Then,the LJ potential will be modified such that both energy and force are equal to zero at the cutoff distance. The detailed simulations on the cut-off effect can be referred to Ref. ^[26]. All the atoms in the simulation model are in a face-center-cubic (FCC) lattice. The solid nanoparticles are formed by carving spheres out of an FCC lattice of atoms.

In this study,the large-scale atomic/molecular massively parallel simulator (LAMMPS)^[27], which is an open-source software written in C⁺⁺ and developed at Sandia National Laboratories,is used for the MDS. To mimic solid Cu particles,the mass of atoms forming nanoparticles is 1.590 53 times larger than that of the fluid atom. The edge length of the corresponding periodic cubic simulation box is L_x = 10 nm,L_y = 8 nm,and L_z = 2.6 nm. In our model,four spherical Cu nanoparticles are suspended in the Ar fluid including 3 419 atoms. There are 361 atoms in the nanoparticles. The nanoparticle diameter r_np is 1 nm,and the total amount of the atoms in the model is 5 220. Liquid Ar is chosen to be the base liquid for its stability and simple potential function,to avoid dealing with unnecessary complex interactions which will cost a quantity of time computing^[21]. Two Cu walls,which are fixed,lie in the xz-plane,and are separated in the y-direction with the distance of 6 nm. The periodic boundary conditions are imposed in the flow direction,i.e.,the x- and z-directions. The nanofluid flows in the x-direction due to the external force applied on the system. Canonical ensemble (NVT) is chosen for the MDS,and the temperature of the simulation system is fixed at T = 86 K (T * = 0.71K),which is controlled by the Nose-Hoover thermostat. The starting orientation of each molecule is random,and the starting speed of each molecule is random,according to the Maxwell-Boltzmann distribution^[3]. 3 Results and discussion

Prior to the Poiseuille flow simulations,the whole system is equilibrated at constant volume and temperature. The simulation time step of 0.001τ,where τ is the characteristic time,and the Verlet integration algorithm^[3] are employed in all the simulations. A total of 2×10⁵ time steps are requisite for the accomplishment of the macroscopic data. The entire simulation model is shown in Fig. 3,where the atoms are in the relative equilibrium situations.

Following relaxation,all flow particles are given the inlet driving force in the flow direction along the x-axis. The simulation time step at the force applying process is the same as that at the relaxation process.

To confirm the simulation results on the Poiseuille flow condition,the system energy or enthalpy must be converged^[21]. The system enthalpy is defined and calculated as follows^[21]:

Click Browse formula

where e_total is the total energy of the system,which is the sum of the kinetic energy and the potential energy,and is given by the following expression^[21]:

Click Browse formula

where m_j is the mass,and v_j is the velocity of the jth particle. Figure 4 presents the evolution of the system enthalpy during the equilibrium process.

The initial velocity estimations at the beginning of simulation are not genuine. By updating the quantities during the progress in time,the system will be confined by the simulation con- ditions,and reaches an equilibrated state. Therefore,it can be seen that the system enthalpy increases at first,and then keeps stable with time. The time for the overall energy stabilization is about 5×10⁵ simulation time steps. 3.1 Microstructure of nanofluid flow

The microstructure of the Ar-Cu nanofluid flow inside the nanochannel is discussed under several driving forces. To explore the effect of the external force on the nanofluid flow,three magnitudes of external force,i.e.,F_ext = 0.1PN,0.2PN,0.3PN,are examined. The values of the external force are chosen optional based on the fact that the system enthalpy can be converged. In these situations,the modeled system and the obtained results are not genuine. The driving force applied on the fluid particles in the x-direction,generating a streaming velocity,is similar to the experimental process in which the fluid is forced by gravity or any pressure driven force. In fact,this parameter (force magnitude) can control the velocity and acceleration of the flowing particles and consequently the enthalpy (or energy) of system.

Fig. 3 Atomic positions of nanofluid parti- cles through nanochannel

Figure options

Fig. 4 Convergence of system enthalpy at equilibrium process

Figure options

The total number of the time steps for the following simulation is 2 × 10⁶. Figures 5−7 illustrate the microstructures of the nanofluid flow at different time steps under different external forces.

In the simulation process,nanoparticles move quickly in the flow direction,accompanying with the random motion. Based on the observed microstructures in Figs. 5−7,it is obvious that Cu nanoparticles aggregate into chain-like structures first,and then make the cluster flow in the middle of the channel under an externally applied force. The interaction potential and external force are considered to be the main culprits behind the agglomeration of nanoparticles. When the interaction potential between the nanoparticles is much stronger than that of the base fluid and external force potential,the clustering effect will occur. By applying external force on the nanofluid flow,when two nanoparticles are near each other and attach together, they will assemble to be one bigger nanoparticle. Once the bigger nanoparticle approaches other nanoparticles,a new bigger nanoparticle will form. At the end,all four nanoparticles

Fig. 5 Trajectory of nanoparticles and liquid atoms on Poiseuille flow with Fext = 0.1 PN

Figure options

Fig. 6 Trajectory of nanoparticles and liquid atoms on Poiseuille flow with Fext = 0.2 PN

Figure options

Fig. 7 Trajectory of nanoparticles and liquid atoms on Poiseuille flow with Fext = 0.3 PN

Figure options

join together and assemble into a cluster-like structure. This cluster flows in the flow direc- tion,especially in the middle of the nanochannel due to the large velocity at the core of the nanochannel. As seen from Figs. 5−7,the nanoparticles aggregate at 1 840 ps,1 860 ps,and 2 500 ps when F_ext = 0.1PN,at 1 560 ps,1 570 ps,and 2 155 ps when F_ext = 0.2PN,and at 1 290 ps,1 387 ps,and 1 650 ps when F_ext = 0.3PN. Therefore,the augmentation in the external force results in the reduction of agglomeration instant with the intention that the final cluster is created at 2 500 ps,2 155 ps,and 1 650 ps for the external forces of 0.1PN,0.2PN, and 0.3PN,respectively.

Clustering is a problem that is inevitable and must be solved before nanofluids can be considered for long-term practical uses. Although the increase in the thermal conductivity will increase the efficiency of a system where nanofluids are used,the life of the system may be decreased over time if the particles begin to form clusters. From the molecular dynamics viewpoint,the reasons behind this phenomenon are due to the strong interaction potential of the solid nanoparticles besides the the external driving force on the nanofluid. The driving force leads to the increases in the velocities of the nanoparticles,especially in the middle of the nanochannel. The strong interaction energy results in the aggregation of nanoparticles when they reach near each other.

It is obvious that the aggregation of nanoparticles plays an important role in their en- vironmental risks by influencing their transport,fate,bioavailability,and biological effects. Understanding the fundamental principles underlying the aggregation process of nanoparti- cles and quantitatively describing this process are essential for characterizing the behavior of nanoparticles and further quantifying the risk. 3.2 System enthalpy

The thermodynamic data are presented during the simulation time steps. The reduced enthalpy of the nanofluid flow is exhibited in Fig. 8 for

Fig. 8 Development of reduced enthalpy throughout simulation time steps with different Fext

Figure options

As exposed in Fig. 8,nanoparticle aggregation makes the system become far from equilib- rium. Therefore,the system enthalpy reduces at every aggregation instant,and finally attain a fixed extent until the next aggregation. It can be seen from Fig. 8 that the enhancement of the external force causes extra reduction in the system enthalpy so that the equilibrium state attains at a law period of time under F_ext = 0.3PN during 2×10⁶ simulation time steps. This result has proven that the nanofluids are hydrodynamically instable. The reason is that by the agglomeration of the nanoparticles or clusters,the interaction surfaces of the nanoparticles be- tween each other and with the base fluid particles decreases,and thus the interaction potential of the system reduces and the total energy and enthalpy of the system drop off. The energy keeps dropping until all nanoparticles are attached. 3.3 Density,velocity profiles,and slip length

The thermodynamic statistics is recorded after 2 × 10⁶ simulation time steps. To compute the density and velocity distributions,the simulated system is segregated into many segments along the y-direction,each of these segments is called “bin”. The density and velocity profiles are obtained by statistical averaging the velocities of each atom in each bin. 3.3.1 Density profiles

Figure 9 shows the spreading of fluid particles for the nanochannel comprising single Ar fluid without nanoparticles and nanochannel including Ar-Cu nanofluid at F_ext = 0.1PN.

It is known that the density ordering of single Ar fluid in the Poiseuille flow condition is symmetric^{[12, 13]},as confirmed in Fig. 9(a). The distribution of fluid particles near the walls is high due to the wall particle effect,and it has a uniform state in the middle part of the nanochannel. As seen in Fig. 9,the distribution of the base fluid particles changes by utilizing the nanoparticles in the nanochannel on the Poiseuille flow condition. Therefore,the nanofluid Poiseuille flow in the nanochannel has an asymmetric density distribution. Nanoparticles move randomly and push the base fluid particles into different sides of the channel. Because of large mass of Cu particles relative to the Ar particles,the base fluid particles are influenced by the nanoparticles. It is noted that the arrangement of the base fluid particles changes at different time steps since the size of the agglomerated nanoparticles differs during the simulation run time. At the end time step,i.e.,time step= 2 × 10⁶,the formed cluster of nanoparticles moves around the center part of the nanochannel due to the large velocity of the channel core. Therefore,the base fluid ordering reduces at the middle part of the channel.

Fig. 9 Ordering of base fluid particles at different time steps at Fext = 0.1 PN for single Ar fluid and Ar-Cu nanofluid

Figure options

For extracting the fluid velocity profile in the nanochannel,the simulation area is divided into 24 bins along the y-direction. Three bins at each side of the nanochannel designating the lower wall,the upper wall,and the fluid domain are assigned,whose numbers are in the range from 4 to 21. Figure 10 displays the evolution of the velocity distribution on the nanochannel Poiseuille flow with Ar-Cu nanofluid at F_ext = 0.1PN.

The parabolic velocity profile for the Poiseuille flow in the nanochannel is in well agreement with previous studies^{[12, 13]}. Once force is put on the system,the fluid velocity will vary until all nanoparticles stick together and the system energy will become at the minimum degree,and then the velocity will remain fixed.

Fig. 10 Evolution of base fluid velocity for Ar-Cu nanofluid at different time steps at Fext = 0.1 PN

Figure options

The final velocity distribution of the base fluid particles obtained with the MD predictions is also depicted in Fig. 11 for the system comprising single Ar fluid and the system including Ar-Cu nanofluid. As perceived from Fig. 11,the augmentation of the external driving force results in an increase in the velocity due to the reinforcement of the fluid particles momentum^{[12, 13]}. As shown in Figs. 11(a)−11(c),the base fluid particle velocity of the Ar-Cu nanofluid is larger than that of the single Ar fluid. Adding the same external driving force to the Poiseuille nanofluid flow causes collisions of heavy Cu nanoparticles with the base fluid particles,and increases their velocities. The results show that the nanochannel including nanofluids has larger velocity distribution relative to that of single Ar fluid.

Fig. 11 Final velocity distribution of base fluid particles for single Ar fluid and Ar-Cu nanofluid with different Fext

Figure options

Table 3 Variations of slip length in nanochannel poiseuille flow with variation in external force for single Ar fluid and Ar-Cu nanofluid

Table options

Based on the obtained results,the slip length increases by means of the nanoparticles in nanochannel fluid flow. Adding nanoparticles into the single Ar flow drags the fluid particles into the main flow region. Therefore,the sticking of fluid particles near the wall reduces,and then the slip length increases. Derived from the outcome of this study,the discrepancy between the experimental and the theoretical results may be due to the difference in the slip length, which can be corrected by introducing the observed slip into the experimental parameters to reduce the deviation in the results. 4 Conclusions

In this work,the microstructure and behavior of the Ar-Cu nanofluid flow exposed to an external driving force are studied by the MDS. The system energy is investigated for the mi- crostructure evolution. The simulation results indicate that the Cu nanoparticles in the Ar fluid aggregate into clusters like microstructures under driving force. The strong potential in- teraction of the Cu nanofluids contributes to the increase in the attractive force between the nanoparticles. The steady external force rearranges the structure of the nanoparticle cluster into the pattern moving along the middle part of the nanochannel due to the large velocity. Moreover,the nanofluids under external steady force result in more aggregated structures and consequently lead to clustering effects. The velocity distribution in the Poiseuille flow with nanoparticles is in good agreement with the result obtained by the analytical second-order polynomial prediction. The fluid slips the wall-fluid interface investigated by the MDS. Based on the obtained results,addition of nanoparticles into base fluid can increase the slip length and reduce multilayer-sticking phenomena. The obtained results of the present research are useful for the future experiments studied on the nanofluids flow at nanoscale.