In this work non- equilibrium molecular dynamics (MD) simulations were performed to invest igate the liquid argon flow passing a stationary carbon nanotube and to esti mate kinetic energy of the flow and flow forces exerted on the nanotube. Effects of different factors were also investigated . CFD simulations of flow around a circular cylinder were also performed to compare the results with those of MD around the carbon nanotube. In molecular dynamics simulation, the fluid forces on the body are calculated directly from the summation of the molecular forces exerted on solid atoms by fluid atoms. In this work, this method is used to investigate the ef fects of different flow parameters such as flow velocity, flow temperature on the flow behavior over the carbon nanotube , atom positioning around the nanotube and forces exerted on it. The simulation is 3D , and the computational domain consists of a maximu m of 33,700 liquid argon atoms as the fluid and 240 atoms of carbon which represent the solid nanotube. A single -walled carbon nanotube is simulated as a rigid body of fixed carbon atoms , and both argon -argon and argon- carbon interactions are modeled by the standard Lennard- Jonnes potential function . Flow is driven by rescaling fluid particles velocities at the inlet region with a length of 3% of the domain length in that direction every 50-time steps. Based on the information given, a parallel code was developed, and investigations have been conducted using this code. Results show that among all the parameters which were investigated in this paper, the flow velocity has the most significant effects on the exerte d forces and ki netic energy on the carbon nanotube and the drag force is increased with an increment of the flow inlet velocity in all cases. Lift fluctuates around zero in all cases of the stationary carbon nanotube. Inlet flow velocities and temperature s have also changed atoms positioning in the vicinity of the nanotube and so atomic density in the nearby bins which is used to further clarify the difference in the drag force. Keywords : Molecular Dynamics, Flow Forces, Stationary carbon nanotube, Kinetic Energy . 1. Introduction For many years, many researchers have show n interest in studying carbon nanotubes and related phenomena in detail . Recently, many efforts have been devoted to studying of adsorption of many gas flows by carbon nanotubes [1-3]. Micro -scale fluid flow structure in microchannels or over cylinder s, the physical interactions between them have become subject of many investigations recently [4, 5] . An important subject in these studies was to investigate the behavior of flow and also the exerted force between that and the cylinder [6-15]. In recent years, the study of fluid flow in nanoscale has become a subject of interest , especially the direct interactions of flow and carbon nanotubes , which behave like cylindrical fibers in many aspects, with respect to a high number of applications for that and also the unusual properties of them [8, 16 -18]. Numerical analyses are significant ly important to study the hydrodynamic effects of the flow on different solid bodies to unravel the underlying physics behind the phenomena and it can be really expensive and time- consuming using experimental tests [19-21]. In 2 addition, microfluidic lab -on-chip devices have gained popularity over the recent years to study blood related diseases, food industry and cli nical purposes [22 -24]. A recent review on numerical methods available for such studies is presented in [25] . Moreover, due to the lack of experimental information in Nano -scales, the roles of numerical methods are even more significant in this regard [26] . At length scales less than ten molecular diameters, the continuum theory breaks d own in both gas and liquid flow phases and, thus, the common continuum methods for simulation of fluid flows are not applicable , and the particle- based approaches must be used [27]. A prevalent and popular one among them is molecular dynamics (MD), which is especially suitable for problems with length scales of 10 nm and below [28] . One of the first and outstanding studies for understanding th e behavior of nanoscale flows passing an obstacle was do ne by [26] , where they conducted a 2D molecular dynamics simulation on nanoscal e fluid flowing over a circular -shape obstacle and compared their results with those obtained from other particle -based methods and proved MD to be accurate for studying micro/nano flows and their structures . [27] performed a non- equilibrium molecular dynamics (NEMD) simulation of uniform water flow over a collection of ca rbon nanotubes and investigated drag forces exerted on the nanotube s in detail for different sizes of carbon nanotubes and flow velocities between 50 -200 m/s. They concluded that the inclusion of the mobility of carbon atoms in the nanotube would not affect the hydrodynamic properties of a system of water flow over a carbon nanotube. [29] simulated 2D fluid flow over two side by side molecular cylinders using the molecular dynamics method. They investigated the effects of different flow parameters on flow structure around the body and the interacting hydrodynamic forces on it . In this work, they used two different approaches for calculating the hydrodynamic forces; the balance of flow momentum using momentum conservation law (continuum approach) and direct summation of molecular forces (direct method) . In the continuum approach, they used the velocity profile behind the cylinders to estimate the drag forces from the momentum change of flow passing the cylinders without presenting any details of how to predict the velocity profile around the cylinder . They concluded that the momen tum balance method predicted the drag force with 30 percent error compare to the direct method . [30] introduced a novel indirect model to calculate the force components in molecular dynamics simulations, in which they implemented a modified averagi ng method similar to those of macro -scale to calculate the average values of forces exerted on a nanotube. [7] performed a non-equilibrium molecular dynamic simulation (NEMD) and investigated uniform liquid argon fl ow over a fixed carbon nanotube . Direct summation of molecular forces was used to investigate drag forces exerted on the carbon nanotube for a w ide range of inlet flow velocities and domain sizes for a single - wall and multi -wall carbon nanotube . They discussed details about forces exerted on the nanotube and compared their results with those of finite element analysis and empirical equations. Based on their work, it was concluded that at very low flow inlet velocities, MD results for drag coefficients are a lot larger than those of FE analysis and empirical equations and showed that continuum methods c ould not be implemented for simulation of nano -flows . [31] presented details about monitoring and controlling the pressure -field indirectly by monitoring the atoms’ interactions based on the potential functions in each bin. Nowadays, nanotechnology devices are becoming more and more practic al and with respect to these developments in technology, application of carbon nanotubes has also become larger and larger and therefore, more efforts have been devoted to understanding these nanotubes better to find the possible application in the future. Targeted drug delivery via nanotubes significantly increases the efficacy of the drug, as well as decreases the toxic effects of the drug on healthy/normal cells [32-37]. Microfluidic devices have been developed to facilitate the sorting process for healthy and unhealthy cells in blood Many researchers devoted their efforts to investigate the possibility of using vari ous nanomaterials in cancer treatment process using molecular dynamics method [38-41]. As an example of these studies, [17] fully investigated the possibility of using carbon nanotubes as lamivudine drugs carriers instead of the recent common carriers and concluded that in terms of binding energy, these drugs are preferred to be delivered in capsules of the carbon nanotube. [18] investigated the effects of nanotube rotation on fluid flow behavior around the nanotube. However, there is still a need for in detail understanding of thermodynamics changes in fluid around a stationary nanotube with respect to different parameters. [42] investigated the possibility of building a fluidic gas -driven motor based on using a multi -wall carbon nanotube in a fluid flow and concluded that the density and flow rates play significant roles in the performance of the specific nano -motor. 3 In spite of numerous applicat ions for carbon nanotubes in variety of industries , such as drug delivery , biotechnology, nanomotors, nanotube oscillators and etc., few investigations have been focused on numerical simulations of nanotube and flow interactions and structure of flow around that in nano - scale and also the hydrodynamic forces exerting on it and possible factors playing important roles in this behavior. Thus, in this paper, structure , and behavior of nano -scale liquid argon flow passing over a stationary carbon nanotube has been fully investigated with the objectives of comprehensively studying the effects of various flow parameters on carbon nanotube drag force, atom positioning around it and the average kinetic energy of the flow. E ffects of flow inlet velocity and flow initial temperature on the exerted forces on the carbon nanotube and kinetic energy have been studied in detail . This research provides useful information about effective parameters and how they affect hydrodynamic forces and make more improvements in the a pplications above of carbon nanotubes in industries. 2. Simulation In this paper , all variables and results in the MD section are presented in reduced molecular dynamics units unless otherwise stated. Reduced molecular units used in this paper include ; the molecular mass (m) and diameter ( σ) of liquid argon are used as units of mass and length, respectively and t he ε, which is the strength of interaction for Lennard- Jones Potential, is considered to be the unit of energy. All these units are presented in Table 1 below. Table 1. Molecular units Quantity Symbol Equivalent in SI Length σ 3.40×10 -10 (m) Mass m 6.625×10 -26 (kg) Energy ε 1.657×10 -21 (J) Other par ameters based on molecular units Time ( )1/221/m τσ ε−= 2.15×10- 12 (s) Force 1Fεσ−≡ 4.873×10 -12 (N) Velocity 1U στ−= 1.58×102 (m.s-1) 2.1. Simulation domain specification In this work , the liquid flow of argon over a stationary carbon nanotube is simulated. Carbon nanotube , which is used in this work , is a zigzag (12,0) type with diameter and length of 2.77 and 6.26 molecular units , respectively. The carbon nanotube was defined by providing coordinates of each atom in a separate input file. This st ructure presents a simple carbon nanotube and also allows the results to be compared with those of Tang and Advani [7], which are used for verification of the in -house code. 4 The simulation domain is a cubic box with dimensions of Lz=50 and Ly=60 units in z and y direction s, respectively , and Lx=6.26 units in x-direction , which is equal to the length of the nanotube. Schematic of the simulation domain is illustrated in Figure 1. As can be seen in the figure, the system contains two types of atoms; first , gray and yellow argon atoms representing liquid argon flow and purple carbon atoms forming the carbon nanotube in the domain. The carbon nanotube is assumed to be a rigid body , which means there is no interaction between the carbon atoms , based on the conclusions of [27], that the hydrodynamic properties of water flow over a carbon nanotube are not significantly affected by the mobility of the nanotube itself . The periodic boundary condition is applied to simulation domain bound aries in all three directions and helps the domain to represent a bigger scale area. Figure 1. Schematic of the three -dimensional solution domain 2.2. Simulation details In each case, argon atoms are distributed in the channel in a lattice form of FCC with the bulk density of 0.8 units . Components of the initial velocity of fluid atoms are assigned based on Maxwell – Boltzmann distribution , which results in a zero-value resultant velocity . Initial thermal velocities of atoms are assigned according to the ir initial temperature and initial flow is generated by adding a velocity of 0U in y-direction to all argon atoms . To maintain the flow, velocities of atoms located in the 0.03Ly of the inlet region are reset every 50 time steps. In addition to maintaining the flow, this technique removes excess heat from the system [43] and there is no need to use a thermostat or other techniques for this purpose . Nanotube assumed to be fixed in a specific location in all sections of this work and flow over a fixed nanotube is fully analyzed . 2.3. MD simulation All the atoms in the system interact according to Lennard- Jones (LJ) 12- 6 potential function, however with different parameters for argon -argon and carbon- argon interactions. 12 6 4ij ij ijVrrσσε  = −    , (1) Where ijVis the LJ potential, ε is the strength of interaction , σis the molecular diameter of argon and ijris the distance between a pair of atoms. Based on LJ potential, the interaction force between a pair of atoms is calculated by the following formula: 14 8 248 1 2ij ij ij ijFrrrεσ σ σ  = −    , (2) where ijris the position vector from atom i to atom j and cris the cut -off distance , beyond which the interaction forces are neglected . The interaction forces between argon- argon and carbon- argon atoms are computed using the LJ parameters pr ovided in Table 2 . Table 2. LJ potential parameters [7] U0 5 Argon -Argon Carbon -Argon ε 1.6567 × 10-21 J 1.9646 × 10-21 J σ 3.4 × 10-10 m 3.573 × 10-10 m A cut -off distance of 2.5σis used for both types of interactions and the neighbor list method is used to calculate interactions among atoms. The e quation of motion is integrated using a Verlet scheme with a time step of 0.001t τ ∆= , equivalent to 152.15 10 s−× . With these assumptions , a parallel MD code is developed using atom decomposition approach to reduce computational time. The simulation is performed for 58 10× time steps corresponding to 1.7 ns. 3. Results and discussion In this section, the effects of inlet induced velocity and flow temp erature on the forces exerted on the nanotube have been studied . To evaluate the drag and lift forces on the nanotube, in all of the cases , simulations are performed t o reach an equilibrium time 𝑡𝑡 𝑒𝑒𝑒𝑒 and then conti nued to a desired time t f to get enough data for extracting forces and other macroscopic properties. Some of the methods reported in the literature to detect system equilibration are to record instantaneous values of the energies and pressure during this p eriod [44], or to monitor the positional disorder and velocity distribution of at oms [45]. In this wor k, the trend of kinetic energy is used to detect the equilibrium time. In all cases, equilibrium was reached approximately after 200-time units , however, a more conservative value of 800-time units was assumed for equilibrium time in all cases. 3.1. Force evaluation In MD simulations , exerted forces on solid bodies can be calculated by different methods [46] . The most common method is to calculate them from the summation of forces exerted on that body by each of the fluid particles directly. In this paper, the same approach is used to evaluate drag, lift, and lateral forces of carbon nanotube in y, z, and x -directions, respectively. Instantaneous forces of drag, lift, and lateral fluctuate with time about their mean values in all cases . In Figure 2 part a , instantaneous values of drag and lateral forces, for instance, are plotted versus time for flow around a nanotube with inlet velocity of U 0=1.0 units for simulation of domain 2, which is describe d in Table 3 in the following. In z direction fluctuations are around zero as expected, however , in the flow direction (y) , the instantaneous forces fluctuate around a higher value , which is called the mean drag force . These values are not physical ly the one s, which will be achieved by macr oscopic methods and so, the mean force value s should be computed. Mean forces are averaged values of instantaneous forces over a period of time until equilibrium is reached . In this study mean value s of forces are calculated using the following formulation. This is the same method used by [7] . 1(t) eq eqt ave inst t ttFFN−= − ∑ , (3) where eqttN− is the number of time steps from equilibrium time ( 𝒕𝒕𝒆𝒆𝒆𝒆) to time ( 𝒕𝒕𝒇𝒇) used for data averaging. In these cases, according to the conve rgence tren d of kinetic energy , the value of 800 and 1000 are considered for 𝒕𝒕𝒆𝒆𝒆𝒆 and 𝒕𝒕𝒇𝒇, respectively and conservatively . 6 Figure 2, Typical force components on nanotube : (a) instantaneous values of forces (b) mean values of forces presented in [18, 46] 3.2. Verification To Verify this in -house parallel -code , many simulations have been done for different domain sizes defined by [7] and given in Table 3. Note that the width of the domain is not changed and is equal to the length of carbon nanotube ( Lx=6.265) for all cases and it is not defined in the table . In all the cases, molecular density is 0.8 units and the ratio of L y/Lz is approximat ely constant. Table 3. Different domain sizes used in this work Size Ly Lz Number of Atoms 1 30 25 3741 1.33 40 33 6584 1.67 50 42 10475 2 60 50 14964 2.33 70 58 20251 2.67 80 67 26736 3 90 75 33669 The results obtained for mean drag force s are compared with those of [7] as shown in Figure 3. Figure 3. Comparison of drag forces obtained from present simulation and Tang and Advani [7] 7 As can be seen in Figure 3, there is a good agreement between the results of the two simula tions . Results show that the maximum error betwe en these results was less than 4 % and thus, it can be concluded that the accuracy of the in-house code is verified . Also , the same simulation was done with LAMMPS software. Drag Coefficient obta ined from all simulation s for domains 1, 2 and 3 and the flow inlet velocity of 1.0 are presented compared with each other in Table 4. Table 4. Comparison of drag coefficient of present work, LAMMPS simulation and (7) Domain Size 1.0 2.0 3.0 Present Work 25.55 16.52 14.09 LAMMPS 27.12 17.03 15.21 Tang and Advani 24.95 16.45 14.38 Max D ifference 8.7% 3.53% 7.95% The difference is less than 10% for the simulations. 3.3. Effects of inlet velocity on drag force After verifying the code, the effects of other parameters such as inlet velocity can be studied. For this study , the flow of liquid argon over the nanotube in a domain of size 2 is simulated for 14 different inlet velocities in the range of 0.25 to 3.5 units. The effects of flow inlet velocity on drag force and the drag coefficient of the carbon nanotube for these 14 cases have been studied, and the results are plotted versus inlet velocity in Figure 4. As it can be seen in Figure 4 part a, drag forc e on nanotube is approximately proportional to the inlet velocity. Figure 4 part b shows the change of drag coefficient with inlet velocity. Since the square of inlet velocity appears in the denominator of drag coefficient definition, such decreasing trend in drag coefficient with inlet velocity is expected . Note that the drag coefficient is defined as ( )2 0/ 0.5d avgC Drag U A ρ = − , where 2xA rL= × and ρis the bulk density of the flow, which is stated before and remains the same in all simulations o f this work . Note that the results in Figure 4 were also presented in our previous work [18] . However , it is included here for a better understanding of the findings in this paper. (a) (b) Figure 4. Drag force and drag coefficient of stationary nanotube, a) Drag force versus inlet velocity, b) Drag coefficient versus inlet velocity also presented with different details in [18] 3.4. Effects of t emperature on kinetic energy and d rag force Comparing the kinetic energy over time is an alternative method of ensuring the equilibrium. This method is not used in this work, however, Figure 5 presents an example of kinetic energy convergence and how changing the flow temperature from 95 to 120 kelvin affects that for inlet flow velocity of 1.0 . This figure shows that given the range selected for the convergence, the kinetic energy is converged as well. In addition, by increasing the flow temperature the kinetic energy is only shifted up and increased 8 and the pattern stayed the same. Note that this behavior was observed for all other cases of different flow inlet velocity and was independent of the velocity. Figure 5. Effects of inlet flow temperature on the average flow kinetic energy for different inlet velocities For this study , the flow of liquid argon over the nanotube in a domain of size 2 and with i nlet velocities of 0.5, 1.0,1.5, 2.0 and 2.5 is simulated for flow temperature of 95, 105 and 120K . Note that argon can naturally be liquid in a range of temperature vary ing from a bout 95 to 120 K. First, the effects of inlet flow temperature are considered on the average flow kinetic energy for different inlet velocities. The results are presented in Figure 6 as follows. As can be seen from the figure, increasing flow temperature caused the average kinetic energy of the flow to increase for all different cases. Increasing temperature increases the molecular activity of the fluid particles and thu s, the vibrational motion of the particles. Taking the average of this quantity, the macroscopic kinetic energy of the flow is also expected to increase with increasing the temperature. In addition, these effects were all the same in different inlet velocities because the average particle energy induced was the same in all the cases. Figure 6. Effects of inlet flow temperature on the average flow kinetic energy for different inlet velocities After considering the effects o f flow temperature on the kinetic energy, one might expect the drag force to increase by increasing the temperature in all cases. Drag forces are calculated and plotted versus flow 9 temperature in this section and in Figure 7 part (a) to part (e). As it can be seen in these figures , in low flow velocities , the drag force on nanotu be is decreased as flow temperature increased , however, in higher flow velocities , drag is increased as the flow temperature increased , which is not exactly according to the kinetic energy plots. The rate, at which the drag was reduced was also decreased as the flow inlet velocity was increased. There appears to be a reason for this constant behavior change with flow inlet velocities. To justify the re sults, this behavior originates from the physical basis of the viscosity coefficient of the flow and that is the reason for different behaviors of these plots. The viscosity of a liquid flow decreases with increasing the flow temperature . In liquid flow , fluid viscosity causes by cohesion forces. These cohesion forces decrease with an increment in flow temperature and this causes a viscosity coefficient to decrease. In lower inlet velocities, the effects of the viscosity coefficient are dominant compared to the kinetic energy of the flow, however, in higher flow inlet velocities, the effects of the kinetic energy become stronger and therefore , more dominant compare d to viscosity effects. This comparison can also be seen by having a deeper look over the graphs in Figure 6 again. As the inlet velocity increases, the differen ce between the kinetic energ ies becomes larger , which can be interpreted from the larger gap between two subsequent graphs . All in all, i t should be noted that change of drag force due to temperature change is smaller in comparison with the effect of inlet velocity on it, but even small changes are important in some applications, where the accuracy of the results are of significant importance such as in drug delivery applications. Moreover, it is important to the exact phy sics and possible reasons behind it, which might enable us to create better ideas for the future applications. Note that these behaviors have been seen for all other domain sizes and are independent of that. Based on the numerical formulations of method, i t is highly expected that the same kind of behavio r can be seen for other liquids as well. (a) (b) (c) (d) 10 (e) Figure 7. Non-dimensional drag force of stationary nanotube, a) versus flow temperatu re at inlet velocities of a) 0.5, b) 1.0, c) 1.5, d) 2.0 and e) 2.5 To quantify the effects of temperature change on the drag force exerted on the nanotube, the results for two different flow inlet t emperatures of 95 and 120 K are presented in Table 5 below. As it can be seen from the data in the table, the temperature can change the drag force exerted on the nanotube by a maximum of about 6 percent, which can be an important difference , given the highly sensitive applications of such problems [7, 32] . Table 5. Comparison of non-dimensional drag force by increasing the flow temperature Inlet Velocity Inlet Temperature (K) 0.5 1.0 1.5 2.0 2.5 95 58.14 117.6 2 162.74 209.89 256.50 120 54.72 113.56 165.72 215.43 264.94 Difference Percentage -5.88% -3.45% 1.83% 2.64% 3.29% To show the se effects of flow temperature on the drag force in a clear way , there is a need to analyze the argon atom positions around the nanotube after the equilibrium is reached . Figure 8 and Figure 9 show the fluid atom positions in the vicinity of the carbon nanotube for flow inlet temperatures of 95 and 120K and for two different inlet velocities of 1.0 and 2.0 respectively to be able to compare this difference in effects of temperature on the drag force. Figure 8 shows the positions of argon atoms around the carbon nanotube after the equilibrium is reached for flow inlet velocity of 1.0 and two different inlet temperature s. As can be seen, by increasing the temperature , the kinetic energy of the atoms were increased and they are located closer to the carbon nanotube behind i t. This increased the density of argon atoms behind the carbon nanotube and they exert a higher force component in the opposite direction of the flow. This negative force will decrease the magnitude of the average drag force on the nanotube. That can justi fy the decreasing behavior of the drag force with increasing temperature for lower inlet velocities . Note that the reason for this phenomenon was justified based on atoms viscosity and kinetic energy in the previous section. 11 (a) (b) Figure 8. Comparison between argon -atom positioning in the vicinity of carbon nanotube which is exposed to inlet flow velocity of 1.0 for a flow inlet temperature of a) 95K and b) 120K Figure 9Figure 8 shows the positions of argon atoms around the carbon nanotube after the equilibrium is reached for flow inlet velocity of 2.0 and two different inlet temperature s. As discussed in the previous sections of this work, by increasing inlet velocity, the kinetic energy will become the dominant factor in controlling the drag force exerted on the nanotube. This can be seen from Figure 9 below. According to these figures, by increasing the temperature, the atoms have more internal energy to come closer to the nanotube in the front region of the carbon nanotube and thus , the liquid density was increased in that region , which increases the drag force component exerted on the nanotube . Note that the density was also increased in the region behind the nanotube , however , this increase is more significant in the front region as there are a lot more argon atom in front of the nanotube than there are behind it and the overall average drag force was increased with increasing the temperature of the flow. (a) (b) Figure 9. Comparison between argon -atom positioning in the vicinity of carbon nanotube which is exposed to inlet flow velocity of 2.0 for a flow inlet temperature of a) 95K and b) 120K 4. Conclusion In this work , liquid argon flow over a stationary carbon nanotube was simulated using molecular dynamics method. The flow was driven by rescaling particle velocities at the inlet layer . The d rag force on the nanotube is calculated using the direct summation of forces exerted from each atom in the domain on the nanotube . It was found that drag force increases as the flow inlet velocit y increases, as it would be expected from the macroscopic results. Also, in low flow velocities of liquid flow , drag force exerted on nanotube is decreased as flow temperatur e increased, however, in higher flow velocities , drag is increased as the flow temperature increased. Flow inlet velocity and temperature also change the atom positioning in the 12 vicinity of the carbon nanotube , and it is a significant factor in changing the force components exerted on the nanotube. Results for lift force and lateral force show that instantaneous lift and lateral force quantities fluctuate around zero as expected from continuum results and thus, the mean forces in these two directions are about zero. It should be noted that the effects of flow temperature on drag force are not significant in comparison with effects of flow inlet velocity . These results are based on a simulation of a simple molecular system of monoatomic liquid argon flow . It is firmly expected that these conclusions can also be extended to the flow of other fluid system s such as the system of polyatomic molecule s as it is proposed in [7]. It also should be noted that t he structure of carbon nanotube r emained the same in all of the simulations of this work and it is expected that the results of drag force s for different carbon nanotube structures show the same behavior s and just the quantities would be different due to different nanotube sizes and thicknesses of the nanotubes . 5.