Molecular Dynamics Simulation Of Smooth And Rough Channels Using Short Range Attractive Potential

Introduction

Fluid behavior at the nano–scale cannot be predicted by the well–known Navier–Stokes equations since the continuum assumption–upon which these equations are based–becomes no longer valid. Molecular Dynamics (MD) is a powerful approach to study the flows at the nano-level and it is widely used due to the progress in the computation capabilities and improvement in computer hardware and software over the past few decades. The intermolecular potential that represents the energy distribution and particle-particle interaction is one of the main ingredients of the MD simulations. Fluid interactions are usually described by Lennard-Jones “LJ” potential [1-MD Simulation, haile] and it’s very successful in simulating many fluid flow problems. LJ potential depends on the particle-particle separation as well the cut-off distance since the potential itself is not finite. The Short-Range Attractive Potential “SHRAT” was introduced by Hess and Kroger [2,3 Elastic and Plastic Behavior of Model Solids]. SHRAT potential–Like the LJ Potential–describes both the repulsive and attractive interactions between molecules. Unlike the LJ potential, the SHRAT potential is finite, and its first and second derivatives are zero at the cutoff distance and there is no need to correct any quantity that the potential or its derivatives take part in its calculation. Besides, the SHRAT potential has a cutoff distance of 1. 5σ, whereas LJ potential usually has a cutoff distance of 2. 5σ –3. 2σ which means less computation efforts. The phase diagram of SHRAT potential was constructed by Bandow B, et al, [4].

SHRAT Potential was successfully used in shockwave simulation [5,6,7]. Sock–waves were formed by accelerating piston into a gas channel, shock-wave propagation, reflection, transmission and attenuation by a solid wall with various porous structures were monitored. The results indicate that placing a nano–porous material layer in front of the target wall reduced the stress magnitude and the energy deposited inside the solid by about 30 percent, while at the same time substantially decreasing the loading rate. The nano–channel flow characteristics mainly depend on: channel height, the fluid-solid particles interaction strength near the solid boundary and the structure of the rigid wall [noorian]. The Nano-scale velocity driven flow problem was studied by many researchers using LJ potential. Thompson and Troian [8] studied the slip on solid–fluid interface and the found that slip length is dependent on the shear rate. Barisik and Beskok [9] have studied the near wall interaction and found that the fluid transport characteristics, such as momentum and energy, significantly deviate from predictions of kinetic theory they also found out that the wall force field penetration depth is significant parameter for gas flows in nano-channels in addition to the Knudsen and Mach numbers. Soong,et al [10], have investigated the wall-fluid interaction parameters on the nano-channel hydrodynamics. They found that wall density, lattice plane, flow orientation, and LJ interaction energy have a significant effect on the nanochannel flow characteristics. Ghorbanian and Beskok [11] have simulated liquid argon flows in gold nano-channels; they have varied the interaction strength between argon and gold molecules to investigate the effects of slip–velocity on the slip-heating term and the resulting temperature profiles. In nano–channels, surface roughness is a very important parameter that would affect the flow properties, when compared to channel dimensions–channel height in particular.

At atomic level, walls are composed by individual particles with finite roughness and even smooth walls will always have some kind of wall roughness. Lin, et al [12] have simulated a velocity driven nano-channel with rough wall, they found that the wall roughness parameters affect flow characteristics. Noorian, et al [13] have investigated the effects of the checker wall roughness geometry as well as the attraction energy on the flow characteristics of a pressure-driven nano–channel, two types of roughness shapes were used; Cube-shaped and Sphere-shaped. Ziarani et al [14] carried out a two-dimensional simulation of a nano- channel with rough wall, roughness achieved by distributing a solid rectangle on the wet surfaces of the solid walls, with specified amplitude and periodicity. They found out that slip depends on topology and number of roughening elements as well as the Knudson number. Moreover, they reported that although the adsorption wall has impact on the slip in a smooth channel, it doesn’t influence the slip in rough channel. Esmaeilzadeh et al [15], studied the velocity driven liquid Argon flow in smooth and rough channels, results showed that flow velocity and slip lengths are highly dependent on the flow parameters including the shear rate, the channel heights, solid- fluid interaction strength as well as surface roughness.

Jian-Fei Xie and Bing-Yang Cao [16] investigated gas flow in a permeable wall channel with rectangular nano-pores, a density jump on the permeable surface existed and hence the density distribution across the permeable wall is discontinuous. Besides, energy was consumed in the pores due to the presence of nano–vortices inside the pores that dissipate energy which in-turn significantly affect the slip lengths and boundary friction coefficient on the permeable surface. Priezjev [17] studied shear driven flow in rough channels and investigated the effect of surface roughness on slip behavior in thin liquid films: he concluded that there is a gradual transition in the functional dependence of the slip length on shear rate for simple fluids by varying the strength of the wall-fluid interaction. Besides, thermally smooth walls can be achieved by increasing the elastic stiffness of the solid wall particles. Rahmatipour H, et al,[18] studied the flow behavior in smooth and rough nano-channels with oscillating wall, the slip length in smooth channels depends on wall velocity amplitude, frequency and oscillation period. Slip was reduced by introducing roughing elements to the lower wall of the channel, two types were used: rectangular and triangular roughing elements. Besides, slipis inversely proportional to the fluid -solid interaction strength.

Mathematical Model

The SHRAT potential is used to simulate the Fluid-Fluid, Fluid-Solid and Solid-Solid interaction, it takes the formWhere ε, σ are length and energy scales, respectively, while r is the separation distance between two molecules/atoms. The first derivative of the potential USHRAT(r) yields the force/acceleration, while subsequent integrations can provide the velocity and the position. Here, the cutoff is rather short ranged and smooth, such that not only the potential but also its first and second derivatives vanish at the cutoff distance. The solid molecules are fixed to their equilibrium position with elastic springs and they interact according to the simple harmonic potentialWhere k is the spring stiffness and r0 is the equilibrium “initial” position of the wall molecules, r is the instantaneous position of the solid molecule. The interaction force between a pair of molecules () can be evaluated form the 1st derivative of the potential for fluid particles for Solid particles.

The stress tensor is composed of a kinetic part and a potential part The pressure, however, is the average sum of the diagonal components of the stress tensor Channel set up and reference valuesThe physical computational domain is shown in Fig. 1 is a channel of square cross section where fluid molecules are confined between two parallel walls, the size of the domain is Lx × Ly × Lz= 17. 46σ × 25. 82σ × 11. 11σ, Bottom solid wall thickness; L1 = 4. 76σ, Top solid wall thickness; L2=3. 17 σ and flow region thickness; L3 = 17. 88σ. Where σ is characteristic length in the SHRAT potential. Molecular dynamics simulations codes are written in such way that all quantities are unitless. The reduced unit system normalize all quantities based on the fundamental quantities: m: mass of one atom, Length: σ, energy: ε. For Argon molecules, σ = 0. 34 nm, m = 6. 69e-26 kg, ε =Tref*kB, where kB is Boltzmann constant, kB = 1. 38054E-23, Tref is the reference temperature –for SHRAT Potential Tref = 160 K. The derived quantities can be summarized in table 1 Derived quantity Reference quantity Reference value Time 1. 8711 pS Density 1702 kg/m3 Temperature 160 K Pressure 56. 2 MPa Velocity 181. 7 m/sTable 1. Derived and reference quantities for SHRAT potentialSolution method and boundary conditionsA Fortran code was developed to solve general MD problems including velocity and pressure driven flows. Periodic Boundary Conditions are applied in the flow direction as well as the transverse direction only. Gear’s fifth order Predictor-Corrector algorithm is used to solve the equation of motion, F=ma. Algorithm is built as follows: first, predict positions and up tofifth-order derivatives at time t+dt via Taylor expansions. Force is then computed at time t+dt from the derivative of the pair potential.

Finally, correct the predicted values of the first-step approximations by a scaled term which represents the discrepancy between the second derivative obtained from Taylor series expansion and the acceleration calculated explicitly from force. Mass of each fluid molecule is set to one unit where that of solid molecule is set to 100 units and spring stiffness is also set to 100 reduced force units. Local quantities such as density, velocity, temperature and stress can be computed by subdividing the computation domain into 2D cells or “sampling bins” in the X-Y plane, which fill the computational domain without gaps nor over-laps. There are 48 cells with width of Lx and thickness Ly/48. In each cell, one accumulates the value of the field of interest, given by a local function of atomic positions and velocities,System equilibrationThe fluid molecules are initially located at face-centered cubic lattice sites within the fluid region of the domain. The time step is set 0. 005τ and a cut–off radius of rc= 1. 5σ is then equilibrated at a reduced temperature T = 0. 825 for 20,000 time–steps by rescaling the atomic velocities.

At the end of the equilibration phase, the gas molecules fill their region randomly and uniformly on average, while the solid atoms have small thermal fluctuations about their original lattice positions. This equilibrated configuration provides the initial condition for upper wall movement. Velocity driven flow formationAfter the equilibration is achieved, additional run of 30,000 time–steps is made, the upper rigid wall molecules are accelerated to the desired velocity during the first 2,000 time–steps and they kept moving at this speed till the end of simulation, desired results in each bin are averaged over 10,000 times steps. Flow in Rough channelsA rough channel was constructed by inserting cylindrical poles to the inward side of the bottom wall “xz-plane”. Nine identical staggered poles were introduced to the geometry, each has a diameter “d” and depth “Lp” as shown in figure 3. Roughness of the channel depends on the number of poles and their size “depth and diameter” and it can be better represented by the roughness number η. The roughness number is the ratio of the particles that can be placed between the poles “ghost particles” to the maximum number of particles that can be placed in a slab of size “Lx×Lz×Lp” at the same spacing. In other word, the number of ghost particles that can be placed in the voids to achieve a solid slab to the number of particles that can be present in a solid slab of size “Lx×Lz×Lp” at that particular spacing. The number of ghost particles is proportional to the void space within the rough wall and the density number is affected by the voids. In order to achieve the accurate density number of the flow, the number of fluid particles was adjusted to account for the change in the flow volume due to the presence of voids. The flow volume is the volume of the main flow channel plus the volume of the voids.

The voids volume was calculated by multiplying the volume of rough region “Lx×Lz×Lp” by the roughness number η. Results and discussions1. Code validationTwo runs were made for two identical channels with density number of 0. 81 and Temperature of 132 K, one using SHRAT potential with εwf = 0. 5 and the other using LJ potential with εwf = 0. 6 which is identical to that of Thompson and Troian []. Figure 2 shows that alignment between the two velocity curves is satisfactory. 2. Flow in smooth channel Four channels of different initial density numbers of [0. 65, 0. 70, 0. 75, 0. 80] were simulated with upper wall moving at velocity of 1. 0 σ/τ. The Density, velocity and pressure profiles are shown in figures 4,5 and 6 respectively. The density distribution– shown in figure 4– is not uniform due to strong fluid–wall interactions; were solid particles tend to attract the fluid particles and fluid experience large density oscillation near the solid boundaries in response to the oscillatory motion of solid particles. However bulk density in middle section of the channel is less – on average – but close to the initial density number.

The effect of relatively slow wall velocity of 1. 0 σ/τ is not significant were little variation is noted when density profiles are compared to those of the stationary walls case. Oscillations close to the stationary wall were observed in smaller area compared to those close to the moving wall mainly due to the disturbance caused by the wall motion. Oscillation lengths from stationary wall are almost the same for all density number cases, however, Oscillation lengths from moving wall are proportional to the density number. The velocity distribution is illustrated in figure 5, all velocity profiles are almost the same and density has little effect on slip length and the slope of the velocity profile–i. e. viscosity. Figure 6 shows the effect of density on the fluid pressure inside the channel. The stress tensor and consequently pressure is dependent on the particle–particle interaction and particle velocity. The higher the density, the higher the interaction and consequently the pressure, since the particles are denser at the wall, stresses and pressure become higher. Small wall velocity “1. 0 σ/τ” has no significant effect on the resulting pressure distribution. Figure 7 shows the effect of fluid–solid SHRAT interaction strength on the velocity. The simulation code was validated by comparing the resulting velocity profiles with that of Thompson and Troian [] and the best alignment achieved at εwf = 0. 5, it clear that for higher values the of εwf slip increases at the moving wall and for lower value slip increases at the stationary wall.

Flow in rough channel

Three rough channels with roughness numbers of [ 0. 18, 0. 27, 0. 47] were simulated at initial density number of 0. 8 and compared to that of the smooth channel flow. The smooth upper wall was moving at speed of 1. 0 σ/τ. Profiles of density, velocity and pressure inside the channel are plotted in figures 8,9 and 10 respectively. As illustrated in figure 8, density at channel walls is much higher than the bulk density inside the channel due to fluid–wall strong interactions.

Furthermore, the channel roughness increases the solid wall surface area and hence more fluid particles are trapped between the solid poles on the stationary wall. The surface area of the solid structure increases with the increase of the roughness number as long as solid poles don’t intersect or overlap. The velocity profiles are shown in figure 8. The velocity profiles are linear in all cases, there is a slight decrease in velocity slope with the increase in roughness number. Figure 10 shows the pressure distribution in smooth and rough channels. There is a very high pressure build up near the stationary wall with the increase of roughness number. The pressure is highly dependent on the particle–particle interaction, with the increase of stationary solid wall surface area due to the increase of the roughness number, more fluid particles are present near the stationary wall and pressure becomes higher.