Local stress tensor calculation by the Method-of-Plane in microscopic systems with macroscopic flow: a formulation based on the velocity distribution function

In this work, we showed a calculation method of local stress tensor applicable to non-equilibrium MD systems based on the Method of Plane (MoP). From the relation between the macroscopic velocity distribution function and the microscopic molecular passage across a fixed control plane, we derived a method to calculate the basic properties of the macroscopic momentum conservation law including the density, the velocity, the momentum flux, the interaction and kinetic terms of the stress tensor defined on a surface with a finite area. Any component of the streaming velocity can be obtained on a control surface, which enables the separation of the kinetic momentum flux into the advection and stress terms in the framework of MoP. We verified the present method through the extraction of the density, velocity and stress distributions in a quasi-1D steady-state Couette flow system and in a quasi-2D steady-state system with a moving contact line. In our method, as opposed to volume average method, the density, mass and momentum fluxes are defined on a surface, which is essential to be consistent with the mass and momentum conservation laws in dynamic systems.


I. INTRODUCTION
With the increasing interest in microfluidic devices and nanotechnologies, molecular dynamics (MD) simulations have become a powerful computational tool to examine the fluid behavior for small scale systems.For the understanding in terms of flow fields, the microscopic motion of individual molecules must be averaged, and the stress tensor plays a key role in such macroscopic flow.Within the framework of fluid mechanics, the stress tensor is determined from the velocity fields through the constitutive equation typically including the viscosity, and the local acceleration of the fluid is given by the gradient of the stress tensor as well as the external field to satisfy the momentum conservation.On the other hand, the molecular motion is governed by the intermolecular interaction, and the stress tensor should be defined through the average of the molecular motion and interaction.
More concretely, the macroscopic equation of continuity is given with the density ρ(x, t) and velocity vector u(x, t) both as functions of position x and time t by where x k and u k are the k-direction components of the position and the velocity vector, respectively.The Einstein notation is used with a dummy index k for the vectors.The second term of the LHS is called the advection or streaming term.Equation (1) comes from the mass conservation satisfied for an arbitrary volume V in an enclosing surface S, where n k is the k-direction component of the outward unit normal vector n with respect to the infinitesimal surface element dS.This integral form in Eq. ( 2) means that the mass in V can be changed only by the mass flux passing the surface S. By applying Gauss' divergence theorem for the RHS as it follows for Eq. ( 2) that for an arbitrary V , which corresponds to Eq. ( 1).
Similarly, the macroscopic momentum equation, the Navier-Stokes equation, is given by where the fluid stress tensor component τ kl expresses the stress in the l-direction exerted on a surface element with an outward normal in the k-direction, and F l denotes the external force per mass.This is also derived from the momentum conservation for an arbitrary volume V enclosed by S: meaning that the total momentum in V can be changed by the momentum flux passing the surface S as well as the impulse due to the stress exerted on the surface S and external force exerted on the volume V .Specifically note that the advection ρu l u k and stress τ kl in the 1st and 2nd terms of the RHS, respectively, are separated.
In contrast to the above-mentioned macroscopic feature, microscopic molecules have their own velocities, and their local average and variance correspond to the macroscopic velocity u and temperature under the assumption of local equilibrium in space and time, respectively.Hence, an instantaneous microscopic mass or momentum transfer across an arbitrary surface S exists due to the passage of the constituent molecules even in macroscopically static system without mean flow with u = 0.In the kinetic theory of gases based on the Boltzmann equation with respect to the velocity distribution function (VDF) of the constituent molecules, the macroscopic velocity is subtracted from the microscopic molecular velocity upon the definition of the stress.Before the establishment of the MD method, Irving and Kirkwood 1 (IK) put forward the connection between the macroscopic conservation laws and microscopic molecular motion governed by the inter-molecular interaction through the statistical mechanical theory based on the distribution function in the phase space, and derived an expression of the local pointwise stress comprised of kinetic and interaction parts including Taylor series expansion of differences in delta functions to express the microscopic particle feature.
After the introduction of numerical MD simulations, 2,3 the calculation of average stress in homogeneous bulk systems, equivalent to the average bulk pressure with its sign inverted, was enabled based on the virial theorem, which indeed corresponds to the average of the IK form integrated in space and time.From this bulk stress, the viscosity as a transport coefficient can also be obtained based on the Green-Kubo relation from the time-fluctuation of the off-diagonal stress component averaged in space. 4,5Regarding the local stress implemented for MD simulations with a discrete time-step, Tsai 6 proposed a pragmatic scheme to calculate the averaged stress defined on a flat plane in quasi-1D planer systems, and Thompson et al. 7 extended the approach toward a spherical curved surface for the analysis of the surface tension of a spherical droplet.In both cases, all the momentum flux and intermolecular force across the plane or the sphere, which divide the computational domain, were summed up during the time integration in macroscopically static systems.This type of stress definition is usually called the Method of Plane (MoP), 8 or Hardy 9 stress, where the momentum conservation law with the above-mentioned MoP for quasi-1D systems was proved to be an exact consequences of Newton's laws, 9 i.e., the MoP meets Gauss' divergence theorem for any control volume (CV) surrounded by enclosing surface(s) irrespective of whether the local system in the control volume is homogeneous or not.[23] The momentum conservation is satisfied for the whole system if the VA is properly summed up; however, special care is needed to consider the momentum conservation for local CVs because the VA originally was not designed to examine local momentum conservation to be satisfied through the link to Gauss' divergence theorem. 9,14,18This feature is similar to the atomic stress, 24 for instance provided as stress/atom command in LAMMPS package, 25 often used to simply visualize the stress field.
Going back to the momentum conservation in Eq. ( 6), for systems with a non-zero local flow, i.e., macroscopically dynamic systems with u = 0, the local macroscopic velocity u must be defined on a surface S enclosing a control volume V in MD systems so that the stress calculated in MD simulations may be consistent with the macroscopic momentum Eq. ( 5).
In other words, the momentum transfer due to the microscopic molecular passage across a control surface should be separated into advection and stress contributions to examine the local flow from a macroscopic point of view.
In this paper, we show a calculation method of the MoP-based local stress tensor applicable to non-equilibrium molecular dynamics (NEMD) systems.We provide the formulation for systems consisting of single-component mono-atomic fluid molecules for simplicity while the present framework is also applicable to systems of multi-component or poly-atomic fluid molecules.For the derivation, we introduced the VDF to clearly give the average of physical properties defined on a fixed control plane as we shall see later.With this procedure, we provide not only the expression of the kinetic term of stress tensor but also those of the components in the advection term, i.e., the density and macroscopic velocity as an extension of the MoP.To check its validity, we performed test calculation in two systems: 1) a quasi-1D Couette flow system and 2) a quasi-2D system with liquid-solid-vapor contact lines, both consisted of a Lennard-Jones fluid between parallel solid walls moving in the opposite directions tangential to the walls.In the first system, we compared the density and velocity distributions obtained by the present method and the VA, and we calculated the distributions of the stress components and advection term.Furthermore, we showed that the same velocity distribution was obtained on bin faces with different normal directions, which is essential to determine the advection term.In the second system, the density, velocity and stress distributions are calculated in the complex flow with liquid-vapor interfaces and contact lines.

II. THEORY
We show the derivation of the stress averaged on a finite bin face in a Cartesian coordinate system for single-component mono-atomic fluid in the following for simplicity.Note that the Einstein notation with dummy indices used in Sec.I is not applied hereafter.The fluid stress tensor component τ kl which expresses the stress in the l-direction exerted on a surface element with an outward normal in the k-direction, is given by the kinetic term τ kin kl and the inter-molecular interaction term τ int kl as In the standard MoP for equilibrium MD systems without mean flow consisting of singlecomponent mono-atomic fluid molecules, 8,10 the kinetic term τ kin kl in Eq. ( 7) on a bin face of area S k with its normal vector pointing to the k-th Cartesian direction is calculated by where m i and v i l denote the mass and l-component of the velocity vector v i of fluid particle i, respectively.We also denote the bin face by S k hereafter.The angular brackets denote the ensemble average, and the summation across S k i∈fluid,δt is taken for every fluid particle i passing through S k within a time interval of δt, which is equal to the time increment for the numerical integration.Considering that we deal with a single-component fluid molecules of an identical mass m, we substitute m i with m hereafter.A sign function to ±1 is multiplied to the momentum transfer mv i l across S k to evaluate the kinetic effect on the stress depending on the passing direction.Note that in static equilibrium systems, i.e., systems without macroscopic local mean flow, the advection term is zero in the whole system.
On the other hand, the intermolecular interaction term τ int kl in Eq. ( 7) in the case of simple two body potential is calculated by where r ij k and F ij l denote the k-component of the relative position vector r ij ≡ x j − x i and the l-component of the force vector F ij on particle j at position x j from particle i at position x i , respectively.The summation across S k (i,j)∈fluid is taken for all line segments of the inter-particle interaction between x i and x j which cross S k .A sign function is multiplied for this interaction term to evaluate the force effect depending on the force direction.Note that technically the fluid-solid interaction can also be included as i-j pair in the summation across S k (i,j) in Eq. ( 9), but only the fluid-fluid interaction was taken into account as the fluid stress, and fluid-solid contribution was considered as an external force field. 10,26,27Also note that for multi-component systems or systems with poly-atomic molecules, difficulties mainly arise to treat the interaction force between different kind of molecules or the constraint force 28 of the polyatomic molecules, where the interaction forces should be properly implemented into the stress calculation to satisfy the conservation laws. 29 extend the standard MoP to steady-state NEMD systems with a non-zero macroscopic mean local flow, the mean velocity should be properly subtracted from the kinetic term τ kin kl in Eq. ( 8) so that the macroscopic momentum flux as the advection term due to the mean velocity u may be included not in the stress term but in the advection term within the macroscopic description of the momentum conservation, i.e., in the Navier-Stokes equation ( 5).In the following, we provide a general framework to connect a microscopic variable ξ i of particles and a macroscopic field value ξ(x, t) averaged on S k under non-zero mean velocity based on the local VDF in the Cartesian xyz-coordinate system.
At first, we define the VDF f (x, v, t) for the mass with a velocity v = (v x , v y , v z ) at position x = (x, y, z) at time t, which gives the local density ρ(x, t) by where we rewrite Then a microscopic variable ξ i per mass of particle i can be related to a corresponding macroscopic field variable ξ(x, t) as The RHS denotes the integral weighted with VDF considering an oblique pillar of a base area S k and a height |v k |δt with its central axis parallel to v , which is typically assumed upon the derivation of the equilibrium pressure in the kinetic theory of gases.
With the limit δt → 0, and by rewriting the average of f (x, v, t) and ξ(x, t) in the oblique pillar by f (S k , v, t) and ξ(S k , t), respectively, the integral with respect to x k in the RHS of Eq. ( 11) writes and it follows for Eq. ( 11) that lim δt→0 crossing S k i∈fluid,δt Hence, by dividing both sides by S k δt, is derived as a basic equation for the connection between the macroscopic field variable ξ(S k , t) and microscopic variable ξ i which belongs to the constituent particle i upon crossing S k .Now, we proceed to the expressions of the macroscopic field variables averaged on S k .By substituting ξ(S k , t) and ξ i in Eq. ( 14) with 1 |v k | and 1 , respectively, and using Eq. ( 10), it follows where v i k denotes the velocity component in the k-direction of particle i.Similarly, regarding the macroscopic mass flux ρu l given by substituting ξ(S k , t) and ξ i in Eq. ( 14) with v l |v k | and From Eqs. ( 15) and ( 17), the macroscopic velocity u l results in Finally, to write the kinetic contribution of the stress τ kin kl , we use the expression in the kinetic theory of gases given by By expanding Eq. ( 19), it follows Hence, by substituting ξ(S k , t) and ξ i in Eq. ( 14) with , respectively, it follows where Eq. ( 17) is used in the final equality.Note that the second term in the rightmost-HS can be obtained by using Eqs.( 17) and (18), which correspond to the advection term in the macroscopic momentum conservation in the Navier-Stokes equation (5).By subtracting ρu l u k (S k , t) from the rightmost-HS and leftmost-HS of Eq. ( 21), it follows meaning that the microscopic total momentum transfer in the RHS corresponds to the stress minus the advection term in the LHS.Technically, the summation in the RHS of Eq. ( 23) is calculated during the MD simulation, and as the post process, the stress τ kin kl (S k , t) is obtained by adding the advection term ρu l u k to the total microscopic momentum transfer as where the advection term is calculated by the dividing ρu l (S k , t) • ρu k (S k , t) by the density ρ(S k , t) as in Eq. ( 22): all obtained also as the post process.
The relation between the macroscopic variables in Eqs. ( 1) and ( 5) corresponding microscopic expressions are summarized in TABLE I.
In practice, within the framework of MD, δt (→ 0) must be replaced by a small non-zero time step of ∆t for the numerical integration.Upon this procedure without this limit, we have to assume the following: 1) the change of the distribution function f (x, v, t) within Eq. ( 15) Eq. ( 17) Eq. ( 9) Eq. ( 23) ρu l u k -Eq.( 22) 24) the distance range of |v|δt is negligibly small, and 2) the values of v i k and v i l upon 'crossing' should be properly evaluated based on the position update procedure of particles depending on the time integration scheme.For the velocity Verlet method, which is applied in the numerical test in Sec.III, we adopted v i ≡ x i (t+∆t)−x i (t) ∆t using the positions x i (t) and x i (t + ∆t) of fluid particle i at time t and t + ∆t before and after crossing the bin face to avoid the discrepancy of the mass flux by the MoP calculation and by the position update.
Note that Eq. ( 23) without the limit δt → 0 is the same as the RHS of Eq. ( 8), which simply sums up the momentum transfer across the bin face S k with a sign function Hence, if one locates a control volume with a closed surface consisting of the MoP bin faces, then the momentum conservation is strictly satisfied with Eq. ( 23).Different choices are indeed possible to determine the advection term ρu l u k (S k , t) in Eq. ( 22) to separate the stress τ kin kl (S k , t) from τ kin kl (S k , t) − ρu l u k (S k , t) by Eq. ( 24), and this may sound that the definition of τ kin kl (S k , t) is not unique.However; by setting l = k in Eq. ( 17), the surface normal mass flux is evaluated as the simple sum of the mass passage with a sign function , and this strictly satisfies the mass conservation in Eq. ( 2), meaning that one can choose a unique definition of ρu l u k (S k , t) that simultaneously satisfies the macroscopic mass and momentum conservation.
Another point to be noted is that the final forms in Eqs. ( 15), ( 17) and ( 23) are formally equivalent to the MoP expressions by Daivis, Travis, and Todd 30 , which were derived for a quasi-1-dimensional flow through the expressions of the time derivative of the fluxes in a control volume with the Fourier transform, and were in principle applicable for the average on an infinite plane under a periodic boundary condition.On the other hand, our non-fluxbased derivation with a definition of physical properties averaged on a face through the VDF enables the calculation of physical properties on a finite area.In addition, taking advantage of this non-flux-based feature, one can calculate, for instance, the velocity component u l on a bin face S k (l = k) tangential to the velocity component by Eq. ( 18).This point will be discussed more in detail with a quasi-1D Couette-type flow in Sec.III A as an example.

III. NUMERICAL TEST
The extended MoP was tested through the calculation of the density, macroscopic mean velocity and stress distributions in two systems with a Lennard-Jones (LJ) fluid: a quasi-1D Couette-type flow and a quasi-2D shear flow with solid-liquid-vapor contact lines.Note that both systems are in steady state and we applied time average instead of ensemble average.

A. Quasi-1D Couette-type flow
Figure 1 (a) shows the MD simulation system of a quasi-1D Couette-type flow, where the basic setup is a standard one similar to our previous study. 31,32The two parallel solid walls were fcc crystals and every pair of the nearest neighbors in the walls was bound through a harmonic potential Φ h (r) = k 2 (r − r eq ) 2 , with r being the interparticle distance, r eq = 0.277 nm, and k = 46.8N/m.Interactions between fluid particles and between fluid and solid particles were modeled by a 12-6 LJ potential , where r ij was the distance between the particles i and j, while ǫ and σ denoted the LJ energy and length parameters, respectively.This LJ interaction was truncated at a cutoff distance of r c = 3.5σ and quadratic functions were added so that the potential and interaction force smoothly vanished at r c . 22We used the following parameters for fluid-fluid (ff) and solid-fluid (sf) interactions: σ ff = 0.340 nm, ǫ ff = 1.67 × 10 −21 J, σ sf = 0.345 nm, ǫ sf = 0.646 × 10 −21 J.The atomic masses of fluid and solid particles were m f = 39.95 u and m s = 195.1 u, respectively.Finally, the equations of motion were integrated using the velocity-Verlet algorithm, with a time step ∆t of 5 fs.
The periodic boundary condition was set in the x and y-directions, and 4000 LJ particles were confined between two parallel solid walls consisting of the fcc crystal located on the bottom and top sides of the calculation cell, which directed (001) and (00−1) planes normal to the z-direction.Both had eight layers so that the possible minimum distance between the fluid particle and the solid particle in the outmost layer was longer than the cutoff distance.
The relative positions of the solid particles in the outmost layers of each base crystal were fixed and the temperature of those in the second outermost layers was controlled at a control temperature of 100 K by using the standard Langevin thermostat. 33The system was first equilibrated for 10 ns using the top wall as a piston with a control pressure of 4 MPa without shear so that a quasi-1D system with a LJ liquid confined between fcc solid walls was achieved.After the equilibration, further relaxation run to achieve a steady shear flow was carried out for 10 ns by moving the particles in the outmost layers of both walls with opposite velocities of ±100 m/s in the x-direction, using the top wall as a piston with a control pressure of 4 MPa.Finally, steady shear flow simulation was carried out, keeping their z-position constant at the average position during the 2nd relaxation run, where the system pressure resulted in 3.61 MPa.
We tested the MoP expression (TABLE I) in the steady state, where the local density, velocity, advection term, and stress were obtained as the time-average of 200 ns on a grid with x-normal bin faces with a height ∆z = 0.150 nm and z-normal ones with a width ∆x = 0.145 nm.Assuming that the system is quasi-1D, the distribution in the x-direction was averaged for bins with identical z-positions.For a comparison, we also obtained the density and velocity distributions based on a standard volume average (VA), where the timeaverage in equally divided bin volumes parallel to the solid wall with a height of 0.15 nm were calculated.The density and velocity differences between MoP and VA are shown in Fig. 1 (c).The density difference was within 10 kg/m 3 , which is less than 1 % of the bulk density, and the diag. and off-diag.stress components, diag.stress component, ¦ zz (S z ), advection term §u x u x (S x ) and ¨xx (S x )−©u x u x (S x ) [MPa] Distributions of (a) the diagonal stress component τ zz (S z )(≡ τ int zz + τ kin zz ), advection term velocity difference was also within 0.5 m/s, showing that proposed MoP can extract the density and velocity distribution consistent with VA.
Figure 2 (a) shows the distributions of τ zz (≡ τ int zz +τ kin zz ), τ xx −ρu x u x [≡ τ int xx +(τ kin xx −ρu x u x )] and ρu x u x , where the first two were directly obtained with simple addition based on Eqs.(9)   and ( 21) as also listed in the top part of TABLE I, while ρu x u x was obtained from the density ρ and velocity u x .Note that τ zz − ρu z u z is shown as τ zz because u z is equal to zero in the whole area of the present system.Also note that the calculation of τ int kl in Eq. ( 9) was the same as in equilibrium systems without macroscopic flow.As clearly observed, τ xx − ρu x u x including the advection and τ zz are different away from the solid walls, indicating that the flow effect should be removed to properly evaluate the fluid stress.typically observed also in equilibrium systems, 10,22 because of the layered structure of the liquid as displayed in the density distribution in Fig. 1 (b).On the other hand, τ zz was constant except near the walls, where the solid-liquid (SL) interaction acts as the external force on the liquid.Regarding the off-diagonal components τ zx (= τ xz ), they were constant except just around the walls, where friction from the solid is included in the force balance even in the laminar flow.
In addition to the normal velocity component u k on the MoP plane S k , the calculation of u l (l = k) tangentially to S k is needed for the separation of τ kin kl (S k )−ρu l u k (S k ) in Eq. ( 24) to properly define the stress in general flows with u l = 0 and u k = 0. Including this tangential velocity, we compared the distributions of the density ρ, the mass flux ρu x and the velocity u x averaged on x-normal and z-normal bin faces as another numerical test in the present system in Fig. 1.More concretely, the density ρ(S x ) and ρ(S z ) averaged on x-normal and znormal bin faces S x and S z , respectively were calculated by Eq. ( 15) with setting k = x and k = z, whereas the macroscopic mass flux ρu x (S x ) and ρu x (S z ) were obtained by Eq. (17)   with l = x on S x and S z , respectively.With these definitions, u x (S x ) ≡ ρux (Sx) ρ(Sx) on S x as well as u x (S z ) ≡ ρux (Sz) ρ(Sz) on S z can be obtained as in Eq. (18).Note that in the present laminar flow system with u z = 0, the calculation of u x is practically not needed for the stress separation for τ kin zz , τ kin zx and τ kin xz in Eq. (24). Figure 3 shows the distributions of the (a) density ρ, (b) mass flux ρu x , and (c) velocity u x defined on x-normal and z-normal bin faces, in the system in Fig. 1, where the values averaged on each bin face of x-normal and z-normal are plotted with setting the z-position at the center of each bin face, respectively, i.e., they are staggered by ∆z/2.In the bulk, ρ, ρu x and resulting u x averaged on bin faces with different normal directions agreed well, indicating that the separation of the stress and advection terms in Eq. ( 24) is possible with the velocity values properly evaluated by the proposed method.The difference seen around the top and bottom is due to the layered structures around the two walls, i.e., the values on S z are the average on a surface parallel to the layered structure whereas those on S x are the average across the layers.

B. Quasi-2D shear flow with solid-liquid-vapor contact lines
The top panel of Fig. 4 shows the MD simulation system of a quasi-2D Couette-type flow, where the basic setups are the same as in the quasi-1D system.The periodic boundary condition was set in the x-and y-directions, and 20,000 LJ particles were confined between two parallel solid walls with a distance about 10.4 nm and a dimension of x × y = 3.92 × 39.2 nm 2 so that the LJ fluid may form two quasi-2D menisci with contact lines on the walls upon the preliminary equilibration at a control temperature T = 85 K without shear.
The static contact angle on both top and bottom walls were about 57 degrees.After the equilibration, further relaxation run to achieve a steady shear flow with asymmetric menisci were carried out for 10 ns by moving the particles in the outmost layers of both walls with opposite velocities of ±10 m/s in the x-direction.
After the relaxation run, the density, velocity and stress distributions were obtained by the present MoP expression in the steady state with the time-average of 500 ns on x-normal bin faces with a length of ∆z = 0.149 nm and z-normal ones with a length of ∆x = 0.150 nm.
The middle panel of Fig. 4 shows the distributions of density ρ calculated on the x-normal Indeed, these apparent flow features can be qualitatively visualized by another methods such as atomic stress 24 , but the present method provides the distributions of physical properties defined on a surface establishing a direct link with the conservation laws for arbitrary local volume as described in Introduction, and is generally applicable to a wide range of nanoscale systems with liquid flow.One of our future research targets is dynamic wetting [34][35][36] , for which we plan to examine the mechanical balance exerted on the fluid around a CV set around the moving contact line in Fig. 4. Through the comparison with the static case, 10 this would enable the analysis of advancing and receding contact angle from a mechanical point of view.

IV. CONCLUDING REMARKS
In this work, we showed a calculation method of local stress tensor applicable to non- indicating that the flow effect, i.e., the advection term, was removed to evaluate stress properly.Furthermore, we showed the density, velocity, and stress tensor distributions by the MoP even in a quasi-2D steady-state system with a moving contact line.In our method as opposed to VA, the density, mass and momentum fluxes are defined on a surface, which is essential to have consistency with the conservation laws in dynamic systems.

Figure 1 ( 5 FIG. 1 .
Figure 1 (b) shows the distributions of density ρ and macroscopic velocity in the xdirection u x calculated by the proposed MoP and standard VA as a reference.Note that these distributions by the MoP can be calculated both on x-normal bin faces and on z-normal and (b) diagonal and off-diagonal stress components τ xx (S x ), τ zz (S z ), τ zx (S z ) and τ xz (S x ).

Figure 2 (
b) displays the distributions of the stress component τ xx , τ zz , τ zx and τ xz , where τ zz and τ zx were calculated on z-normal bins whereas the others were obtained on x-normal bins.As explained above, the stress value τ xx was calculated by adding ρu x u x to τ xx − ρu x u x whereas the advection terms for the others can be neglected considering u z = 0.In the bulk region sufficiently away from the walls, τ xx = τ zz and τ zx = τ xz are satisfied as expected from the solution of a laminar Couette flow, and the former indicates that the stress τ xx is adequately calculated by the proposed MoP with the resulting value −τ xx (= −τ zz ) equal to the external pressure value of 3.61 MPa.The wall-tangential diagonal stress τ xx fluctuates near the walls as

FIG. 3 .
FIG. 3. Comparison of the time-averaged distributions of the (a) density ρ (b) mass flux ρu x and (c) velocity u x averaged on x-normal and z-normal bin faces S x and S z , respectively.
FIG. 4. Top: Quasi-2D Couette-type flow system of a Lennard-Jones liquid confined between two solid walls.Middle: Distributions of density ρ, velocity u, and off-diagonal stress component τ zx .Black arrow denotes the macroscopic velocity calculated by the proposed Method of Plane.Bottom: Distributions of diagonal stress components τ xx and τ zz .
MD systems based on the Method of Plane(MoP).From the relation between the macroscopic velocity distribution function and the microscopic molecular passage across a fixed control plane, we derived a basic equation to connect the macroscopic field variable and the microscopic molecular variable.Based on the connection, we derived a method to calculate the basic properties of the macroscopic momentum conservation law including the density, velocity and momentum flux as well as the interaction and kinetic terms of the stress tensor defined on a surface with a finite area.Any component of the streaming velocity can be obtained on a control surface, which enables the separation of the kinetic momentum flux into the advection and stress terms in the framework of the MoP.We verified the present method through the extraction of the density and velocity distributions by volume average (VA) and the MoP in a quasi-1D steady-state Couette flow system, seeing that the stress tensor distribution by the MoP satisfies the solution of a laminar Couette flow in the bulk,

TABLE I .
Microscopic expressions for the calculation of the corresponding macroscopic properties defined as the average on bin face S k in steady-state systems.The top four properties can be directly calculated from steady-state systems through the MoP procedure, whereas the others below are derived from the four.