Effect of intermolecular potential on compressible Couette flow in slip and transitional regimes

Effect of intermolecular potential on compressible Couette flow in slip and transitional regimes" (2014). Articles you may be interested in Molecular dynamics simulations of high speed rarefied gas flows AIP Conf. The sign change effect of the energy flux and other effects in the transitional regime for the Couette problem AIP Conf. Navier–Stokes modeling of a Gaede pump stage in the viscous and transitional flow regimes using slip-flow boundary conditions The effect of intermolecular potentials on compressible, planar flow in slip and transitional regimes is investigated using the direct simulation Monte Carlo method. Two intermolecular interaction models, the variable hard sphere (VHS) and the Lennard-Jones (LJ) models, are first compared for subsonic and supersonic Couette flows of argon at temperatures of 40, 273, and 1,000 K, and then for Couette flows in the transitional regime ranging from Knudsen numbers (Kn) of 0.0051 to 1. The binary scattering model for elastic scattering using the Lennard-Jones (LJ) intermolecular potential proposed recently [A. Venkattraman and A. Alexeenko, " Binary scattering model for Lennard-Jones potential: Transport coefficients and collision integrals for non-equilibrium gas flow simulations, " Phys. Fluids 24, 027101 (2012)] is shown to accurately reproduce both the theoretical collision frequency in an equilibrium gas as well as the theoretical viscosity variation with temperature. The use of a repulsive-attractive instead of a purely repulsive potential is found to be most important in the continuum and slip regimes as well as in flows with large temperature variations. Differences in shear stress of up to 28% between the VHS and LJ models is observed at Kn=0.0051 and is attributed to differences in collision frequencies, ultimately affecting velocity gradients at the wall. For Kn=1 where the Knudsen layer expands the entire domain, the effect of the larger collision frequency in the LJ model relative to VHS diminishes, and a 7% difference in shear stress is observed. C 2014 AIP Publishing LLC.


I. INTRODUCTION
The direct simulation Monte Carlo (DSMC) technique is a powerful numerical method for solving rarefied gas flow problems encountered in high-altitude hypersonic aerothermodynamics, 1, 2 in-space propulsion, 3 vacuum technology, 4,5 and microsystems. 6DSMC is a stochastic approach to solve the Boltzmann equation by simulating the motion and interaction between gas molecules based on a specified interaction potential.The first implementation of the DSMC method applied a simple hard sphere (HS) interaction 7 which belongs to the class of purely repulsive interactions.The HS model leads to a square-root dependence of viscosity coefficient on temperature which deviates from experimental observations for common gases. 8The variable hard sphere (VHS) model proposed by Bird 9 results in a more general power-law viscosity variation with temperature and has been widely used for DSMC simulations of single species gas flows due to its computational efficiency and ease of implementation.
Later, several variations were proposed to the VHS model including the variable soft sphere (VSS), 10 generalized hard sphere (GHS), 11 and generalized soft sphere (GSS) 12 which all belong to a class of purely repulsive interactions.The VSS model modifies the scattering law of the VHS model by using a parameter that allows reproduction of measured diffusion coefficients in addition to viscosity coefficient.The GHS model uses the same scattering law as the VHS model but applies a modified collision cross-section with parameters chosen to reproduce both viscosity and diffusion coefficients.The GSS model has a cross-section similar to GHS model but with a scattering law similar to the VSS model and has fewer parameters than the GHS model.Dimpfl et al. 13 used the Born-Mayer potential with a hard sphere scattering kernel to develop what is referred to as the Extended-VHS (EVHS) model to deal with hyperthermal gas flows.Other collision models such as the μ-DSMC method 14 can reproduce arbitrary viscosity variation with temperature by adjusting parameters of the hard sphere model in each cell based on the local time-averaged temperature.In summary, the parameters of all these models are determined in such a way that they match observed or theoretical bulk transport properties such as viscosity and diffusion coefficients.
Each of these models have limited validity due to the fact that they do not account for the attractive force between molecules at large distances.With the exception of the GSS model, most of these models are limited to a relatively narrow temperature range in which the viscosity variation is accurately reproduced.For problems involving a wide range of temperatures, this would be insufficient.Examples include flows with large temperature variations such as supersonic plume expansions in vacuum as encountered in space propulsion and in low-pressure deposition of thin film materials where the vapor temperature varies from the melting point of the material to very low temperatures due to rapid expansion into vacuum.In other applications, the detailed collision dynamics that includes the contribution of the attractive interactions between molecules becomes important.For example, this is the case when the orientation of a molecule incident on a solid surface should be accurately predicted along with the incident energy of the molecule for thin film growth modeling.
Implementation of realistic repulsive-attractive interaction potentials in DSMC simulations have been reported in the past 15,16 but has never been used widely.One of the popular attractive-repulsive interaction potential is the Lennard-Jones (LJ) potential which is known to accurately represent the interaction of molecules and was used by Koura and Matsumoto 16 in DSMC simulations to investigate the velocity distribution functions within an argon shock wave at free-stream temperatures much lower than the potential well depth of argon.Their implementation used numerical integration to obtain the elastic scattering angle for every collision.More recently, Valentini and Schwartzentruber 17 performed large scale molecular dynamics simulations to revisit the computation of velocity distribution functions within an argon shock wave.They reported significant differences between the DSMC computations with the VHS model and the molecular dynamics simulations with the LJ potential.Sharipov and Strapasson 18 demonstrated a data look-up table implementation in DSMC, applicable for arbitrary intermolecular potentials, and applied it to Couette and Fourier flows using the LJ potential.An approach presented recently 19 reduces the computational cost relative to direct implementation by representing the scattering angle as a polynomial expansion in non-dimensional collision parameters and is referred to as the LJ polynomial approximation (LJPA) model.In the aforementioned studies, [15][16][17][18][19] the focus has been either on demonstrating new numerical schemes for implementing the LJ potential in DSMC, or evaluating the accuracy of the VHS model in normal shock waves near continuum.The present paper aims to quantify the differences in transport properties between the widely used, repulsive VHS model and the attractive-repulsive LJ model for a variety of flow regimes.
The remainder of the paper is organized as follows.Section II provides the necessary background about existing DSMC collision models and the LJ potential.A verification of the DSMC implementation of the LJPA model is presented in Sec.III, while Sec.IV presents the DSMC simulations of Couette flow in the slip regime using the VHS and LJPA collision models.DSMC simulations of low-temperature Couette flows in the transitional regime using the VHS and LJPA collision models are presented in Sec.V. Finally, concluding remarks are summarized in Sec.VI.

A. Purely repulsive interaction models
The various purely repulsive models proposed in the past are briefly summarized below.A typical collision for purely repulsive interactions such as the HS, VHS, VSS, etc. is shown in Figure 1.The single parameter that completely describes an elastic collision between two atoms or molecules by relating the pre-collisional and post-collisional velocities is the scattering angle, χ , which depends strongly on the nature of intermolecular forces between the two molecules.χ depends on the relative energies of the two colliding molecules, = m r c 2 r /2, and the impact parameter, b.In the expression for , m r is the reduced mass, m 1 m 2 /(m 1 + m 2 ), for collisions between molecules of mass m 1 and m 2 , and c r is the relative velocity between the colliding molecules.The purely repulsive interaction models result in a scattering angle that is non-negative and lies between 0 and 180 • .The scattering angle for the VHS model, 9 depends on the collision energy through the diameter, d V H S .This VHS diameter is varied with the relative velocity between the colliding molecules 9 via the form of Eq. (2), is similar in form to the VHS scattering angle, and in fact when α is unity the VHS scattering is obtained.The diameter dependence on relative velocity, d V SS , is the same as for d V H S .The GHS model also uses a scattering law similar to the VHS, but the total cross-section is instead varied with relative collision energy through the sum of two terms each depending on relative collision energy.The GHS total cross-section, 11 has a set of four constants: α 1 , α 2 , ω 1 , and ω 2 which are determined by fitting to experimental or theoretical viscosity and diffusion coefficients.Apart from the four fitting parameters, the two parameters from the Lennard-Jones (LJ) intermolecular potential (σ LJ and LJ described later) are also used, making it a total of six molecular model parameters.
The GSS model combines aspects of the GHS and the VSS models by using a total crosssection similar to the GHS and a scattering law similar to the VSS model.The GSS model simplifies the total cross-section calculation by reducing the number of model parameters to set.This is accomplished through curve fits to viscosity and diffusion cross-sections obtained using the LJ potential.All constants in the GSS model, four cross-section constants and one scattering parameter, are independent of the gas under consideration; thereby requiring only the two LJ potential parameters to implement the collision model.This results in a model that is easy to implement and is general enough in that it does not require the determination of additional parameters apart from the LJ potential parameters.The four cross-section constants have thus been computed for the GSS model, and the total cross-section is given by 12 with the scattering angle, 12 However, the curve fits for diffusion and viscosity cross-sections used to obtain the GSS parameters are valid only for 0.3 ≤ kT/ LJ ≤ 400 which may not be accurate for temperatures outside this range.Also, the rapid increase in the cross-section for low-energy collisions leads to a very large value of (σ T c r ) max when a low-energy collision is encountered, thereby making the method extremely inefficient without special treatment for the low-energy collisions.Previous researchers 20,21 who used the GSS model in DSMC typically opt for a constant or a linearly varying cross-section which does not have physical basis and could make the model inaccurate.Also, as pointed out by Kim et al., 22 the GSS model may not lead to accurate values for other collision integrals apart from the viscosity and diffusion cross-sections.The modified GSS model 22 which attempts to address this issue proposes a total cross-section of the form with constants C 1 , C 2 , ω 1 , and ω 2 depending on the gases under consideration.A table with constants for various combinations of common gases is provided. 22The use of gas-dependent constants offsets the advantage of the GSS model which uses only the LJ potential parameters to completely describe the scattering cross section.

B. Lennard-Jones interaction model
The LJ potential is one of the realistic intermolecular potentials which accounts for both the short-range repulsive and long-range attractive forces between neutral molecules through inverse power laws.Due to the inclusion of the attractive force, the LJ potential, has a minimum energy denoted by LJ at a distance r = σ LJ 2 1/6 .This minimum potential energy along with the intermolecular separation distance corresponding to zero potential energy, σ LJ , constitute the LJ model parameter set.The inclusion of the attractive component also results in negative scattering angles for larger impact parameters, and a typical collision is shown in Figure 2. It should be observed that the impact parameter is much greater than that in Figure 1.The viscosity variation with temperature, T, obtained from Chapman-Enskog theory 23,24 using the LJ potential is where m r is the reduced mass, k is the Boltzmann constant, and V and W (2) (2) are functions of kT/ LJ that are, for example, tabulated in the work of Hirschfelder et al. 23,24 Figure 3 compares the viscosity variation of argon obtained from experiments, [25][26][27] predictions using the LJ potential, 23 the most commonly used VHS model for argon 9 and a curve fit obtained by Maitland and Smith 8 using data from a number of experiments for argon.The Maitland-Smith fit is expected to lead to an error of less than 1% in the temperature range 80-600 K and 1.5% between 600 and 2200 K.As can be seen, the viscosity variation obtained using the LJ potential parameters agrees extremely well with the experiments and the Maitland-Smith curve fit; whereas the VHS model's performance is best near its reference temperature of 273 K.This is mainly due to the fact that the VHS model does not account for the attractive component of the force which becomes important at low temperatures.At higher temperatures, the VHS model predicts much smaller collision cross-sections and as a result over-predicts viscosity.It should be mentioned that use of a different VHS model in the low-or high-temperature regimes will lead to good agreement with the experimental data but a single VHS model may not be able to provide an accurate description of the viscosity behavior over a wide range of temperatures.Though the most fundamental method to obtain LJ potential parameters is to construct the potential energy surface for the interacting molecules, a first order estimate can be obtained using empirical relations based on the solid state molar volume and the melting temperature of the gas species under consideration.These relations are 28 where T m is the melting temperature in K, and V m,sol is the solid state molar volume in cm 3 /g mol.  9 LJ potential, 23,24 Maitland-Smith fit, 8 and experiments by Kestin et al. 25 and Vogel et al. 27 This article is copyrighted as indicated in the article.Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions.For the LJ potential, the scattering angle χ may be written as 16 6 , and u = r 0 /r.The distance of closest approach between the two molecules is denoted by r 0 ; while the reduced relative collision energy is a ratio of the relative collision energy to the LJ potential well depth, * = / LJ .A reduced form of the impact parameter is also used, where the impact parameter is non-dimensionalized by the intermolecular separation distance of zero energy, b * = b/σ LJ .Varying u from 0 to 1 varies r from ∞ to r 0 .z in the above integral depends on b * and * and is obtained by solving the implicit equation, For a given value of b * and * , Eq. ( 13) is solved using a bisection method, and the solution for z is then used to compute χ via numerical integration of Eq. ( 12).The numerical integration is performed using the Gauss-Chebyshev quadrature method, which approximates the integral of Eq. ( 12) as where w k are the weights for Gauss-Chebyshev quadrature, y k are the zeros of the Mth degree Chebyshev polynomial denoted by φ M , and I is the integrand of Eq. ( 12) given by The zeros, y k , of φ M are given by and the weights are all equal to For the results presented in this work, 800 quadrature points (M = 800) were used; with further increase in the number of quadrature points leading to negligible difference in the computed value of χ including regions near the singularity due to orbiting collisions.It should be mentioned that any other numerical integration scheme such as Gauss-Legendre could have been used and would have produced the same result, provided a sufficient number of points were used.

C. Comparison of VHS and LJ potential elastic scattering
Since there is no general relation between the parameters of the LJ potential and the VHS model, comparisons of the scattering angle for the two cases can be made only for a given gas.Figure 4 shows contours of the scattering angle obtained using the VHS model and the LJ potential computed using parameters for Argon.The relevant parameters were: σ LJ = 3.405 × 10 −10 m, LJ = 0.0103 eV , d ref = 4.17 × 10 −10 m at 273 K, and ω = 0.81.As can be observed, there are significant differences for a range of b * and * values.The differences are pronounced for low energy collisions, as can be observed from the contour lines close to the x-axis, and also for long-range collisions corresponding to values of b > d V H S .Figure 5 shows the variation of χ as a function of b * for four different values of * in order to compare the LJ potential, VHS model, and GSS model scattering angles more rigorously.Since the VHS and GSS models are purely repulsive models, negative values of scattering angle are not obtained while the LJ potential has a shallow well for χ even at relative energies of * = 5.0.For the lowest value of * shown, the GSS model leads to a very large collision

III. VERIFICATION OF LJPA MODEL
The LJPA model is verified through comparison of equilibrium collision frequencies computed at several temperatures to theory.Equilibrium collision frequency is computed as a product of number density and the mean of the product of total cross-section and relative velocity.The mean of the product, σ T c r , is averaged over the relative velocity distribution function, f c r , such that the equilibrium collision frequency has the form, 9 The total cross-section for the LJ potential would theoretically extend to infinity, and therefore a maximum impact parameter, 10,19 is set to a value beyond which scattering angles less than χ min are neglected.This variation of total cross-section with relative velocity also necessitates the use of numerical integration techniques, such as those described in Sec.II B. In this section, as well as in the remainder of the paper, a minimum scattering angle of 0.1 radians is used.
The DSMC implementation of the LJ potential itself follows the work of Koura and Matsumoto 16 very closely with the principal difference being the use of the polynomial expansion for the scattering angle instead of a concurrent numerical integration.To reiterate, the reduced impact parameter, b * , is uniformly distributed between 0 and b max once a collision pair is chosen, and then the elastic scattering angle for the chosen collision pair is obtained using the LJPA model described in detail in Ref. 19.The DSMC simulations were performed using the one-dimensional code of Bird, 9 DSMC1.FOR, after implementation of the LJPA model.
Collision frequencies are computed from DSMC using both the LJPA model and the direct LJ scattering from the integral, Eq. ( 14).A single cell of length, 1 mm, with 40 000 molecules and a number density corresponding to Kn = 10 is used.A time-step of τ 0 /10 is used, where τ 0 is the mean collision time as computed from theory (1/ν 0 ).DSMC sampling speeds of each LJ scattering angle implementation are recorded using a single processor on the Hansen compute cluster.The Hansen compute cluster has four 2.3 GHz 12-Core AMD Opteron 6176 processors per node, 10GB Ethernet interconnections, and 48.8 TeraFLOPS performance.Sampling speeds for the LJPA model ranged from 33.0 samples/s to 39.9 samples/s, while the direct LJ scattering implementation resulted in approximately half the sampling speed of 14.9-15.1 samples/s.Collision frequencies computed from DSMC and theory are reported in Table I along with their corresponding errors relative to theory.The statistical errors reported as a two-sided 95% confidence interval are much smaller than the errors in collision frequency relative to theory, thus providing confidence in the results.Both DSMC implementations of the LJ scattering angle have comparable collision frequencies at each of the simulated temperatures, and are in error by less than 0.02%.The good agreement verifies the implementation of the LJPA model in DSMC, and the errors are consistent with both those reported for the VHS model at 300 K 29 and for the VSS model in Table II.

IV. COUETTE FLOW SIMULATIONS IN THE SLIP REGIME
The LJPA model is evaluated by initially using it in DSMC simulations of compressible subsonic and supersonic Couette flow problems similar to those used by Bird 9 to verify the VHS model and Macrossan 14 to verify the μ-DSMC technique.The three cases each have all the same initial conditions and numerical parameters, and only differ by the specified wall temperature.Use of the same number density and wall velocity results in differing Knudsen numbers, Kn, and Mach numbers, M, between the cases, respectively.The flow conditions for each of the three cases as well as the numerical parameters used are summarized in Table III.
The Knudsen numbers reported in Table III are defined based on the mean free path, obtained using the mean free path theory. 30The viscosity, μ, is the viscosity corresponding to the wall temperature, ρ is the initial density, and c is the mean thermal velocity based on the wall temperature.Apart from the theoretical viscosity, the viscosity in DSMC simulations may be computed given shear stresses and velocity gradients.
Viscosity is computed from DSMC simulations through its relationship to shear stress and velocity gradient.Since both shear stresses, τ DSMC =< ρc x c y > , (24)   and velocity gradients, ∂v/∂x, are directly calculated from microscopic properties, viscosity may then be determined through the relationship The brackets, . . ., in the expression for shear stress (Eq.( 24)) denote an average value.In Secs.IV A-IV C, reported average viscosity and shear stress values are averaged over the central 60% of the domain.This averaging procedure is used in order to exclude the Knudsen layer, which extends several mean free paths from the walls.The sampled viscosities and shear stresses have statistical errors inherent to the DSMC method, but are estimated to be below 0.2% and 0.001% for viscosity and shear stress, respectively.

A. Case 1: Moderate-temperature, subsonic, slip flow
For the compressible, subsonic Couette flow problem, a moderate wall temperature of 273 K was chosen.This temperature corresponds to the reference temperature for the VHS model, and as such the VHS model is expected to perform well in this case.The Kn for this wall temperature is 0.012 and is therefore in the slip regime.
Figure 6 compares the variation of normalized temperature obtained using the VHS model and the LJPA model.The agreement between the two models is excellent with the maximum difference  is expected to be small.The average values of shear stress, including statistical error estimates, obtained using the VHS model and LJPA model are 6.28 × 10 −3 Pa and 6.60 × 10 −3 Pa, which corresponds to a difference of about 5%.Viscosities obtained from DSMC simulations are compared to the theoretical viscosities obtained using the LJ intermolecular potential and the VHS model at the local temperature, T DSMC in Figure 7.The LJPA model agrees with the theoretical value obtained using the LJ intermolecular potential with an average error of 1.0% in regions outside the Knudsen layer where the DSMC viscosity computed using Eq. ( 25) is not expected to match the theoretical value.The DSMC viscosity obtained using the VHS model agrees with the theoretical value with an average error of 1.15%.Statistical errors are estimated to be less than 0.1% for this case and therefore have negligible effects on the DSMC viscosity comparisons to theory.
Comparisons of theoretical viscosity to the experimental fit of Maitland and Smith 8 are also shown Figure 7, where the fit corresponds to the local temperatures, T DSMC at each of the spatial locations.The VHS model is in better agreement with the experimental fit than LJ.While VHS is in error by 4.9%, the LJ model has a 9.5% error.This is expected since the VHS model parameters are based on viscosity measurements made at 273 K.
The difference in shear stress of about 5% for the VHS model and LJPA model reported earlier is a combination of the models deviating to small degrees from the theoretical value and also the LJ theoretical viscosity being about 4% higher than the VHS viscosity at 273 K.The other key aspect that can be observed in Figure 7 is the Knudsen layer where the DSMC viscosity deviates from theoretical viscosity.It can also be seen that the Knudsen layer is smaller for the LJPA model due to smaller mean free path.
For Case 1, the time taken for 1 000 000 sampling time-steps using the VHS model was 2288 s whereas for the LJPA model, the time taken was 3623 s which is about one and a half times larger than that of the VHS model.However, it should be mentioned that most of this increase is contributed by the higher number of collision events due to the long-range nature of the LJ potential.A collision pair is selected at a rate of approximately 6 for every 10 collision attempts, per the acceptancerejection method, 9 for both VSS and LJPA models.The GSS model which is the closest to the LJPA model in terms of fidelity was also implemented in DSMC1.FOR for comparison with LJPA and VHS models.The ratio of collision events to selections remains about the same for all cases, but the collision frequency significantly changes.Thus, the low-energy collisions in Case 2 result in less collisions occurring each sampling time-step-reducing the computational time required.Applying the same reasoning, Case 3 requires more computational time as a result of the increased collision frequency.

B. Case 2: Low-temperature, supersonic, slip flow
The comparison performed for the subsonic compressible Couette flow was repeated for a much lower wall temperature of 40 K, at which differences are expected between the VHS model and the LJ potential.Sound travels much slower at low temperatures, and therefore the Mach number is higher under these conditions than in the other two cases.The Mach number is 2.55 and hence a significant region of the flow in the gap is supersonic as opposed to the other two cases in which the flow never becomes supersonic in the entire gap.The Knudsen number of 0.0051 is approximately half the value for Case 1.All flow conditions and numerical parameters for Case 2 are summarized in Table III.
Figure 6 shows a comparison of the normalized temperature and velocity variation in the gap obtained using the VHS and LJPA models.Temperature variation across the gap is more significant for this case due to increased viscous dissipation.While a maximum difference of only 0.4% is observed in the normalized temperature profile between VHS and LJ models, the velocity profiles for both VHS and LJPA models significantly differ from that of incompressible Couette flow.A linear velocity profile is the analytic solution to an incompressible Couette flow, but this case has compressible, supersonic flow.Thus, the deviation from the analytic solution is greatest in this case.The velocity slip at the wall is the least in this case as a result of the smaller mean free path, and the VHS and LJPA models are in worse agreement than in Case 1. Figure 7 compares the shear stress variation in the gap obtained using the two models and the higher VHS viscosity at temperatures around 40 K leads to a significantly higher shear stress in the gap.A slight increase in shear stress across the gap is observed in both models as a result of the nonlinear velocity profile typical of compressible Couette flow.The average shear stress obtained using the VHS model and the LJPA are 1.22 × 10 −3 Pa and 1.56 × 10 −3 Pa, respectively.This corresponds to a 28% higher shear stress predicted by the VHS model when compared to the LJPA model.In order to verify that the viscosity variation predicted theoretically is reproduced by the DSMC simulations, Figure 7 compares μ theory with μ DSMC for both VHS and LJPA models for the supersonic Couette flow.The agreement with theory is good again with the average error being 1.64% for the VHS model and 1.26% for the LJPA model.A viscosity measurement by Kestin et al. 31 at 50 K is also illustrated in Figure 7 (shown at a location corresponding to a local temperature of 50 K), and the LJPA model is observed to be in better agreement than the VHS model.
One of the important aspects to be considered in order to evaluate the practical applicability of this method is the computational overhead associated with the polynomial computations.The use of pre-computed scattering angles makes the implementation very efficient with insignificant overhead as compared to the VHS model.The time taken for 1 000 000 sampling time-steps is 2786 s for the LJPA model, while for the VHS model only 1961 s is required.Again, the additional time required for the LJPA model is a result of the higher number of collision events due to the long-range nature of the LJ potential.

C. Case 3: High-temperature, subsonic, slip flow
The behavior of the LJPA and VHS models are studied at higher temperatures in this case.A wall temperature of 1000 K is set with the same wall velocity as in previous cases.The Knudsen number is 0.017, which indicates the flow is in the slip regime similar to Case 1.The Mach number for this case is 0.51, and is the least of the three cases.All flow conditions and numerical parameters are summarized in Table III.
The variation in normalized temperature for both LJPA and VHS models are observed in Figure 6 to be in excellent agreement with each other, and have the least variation across the gap relative to the previous two cases.Figure 6 shows a more significant difference in velocity profiles between the two models at the walls.The velocity slip for this case is the largest as a result of the larger mean free path, with the VHS model predicting a larger velocity slip than the LJPA model.A mostly linear velocity profile is also observed in Figure 6 for most of the domain due to the lower Mach number.
Figure 7 shows the shear stress and viscosity variations across the gap, respectively, using the LJPA and VHS models.The shear stresses remain constant across the gap, with a mean value of 1.56 × 10 −2 Pa for LJPA and 1.74 × 10 −2 Pa for VHS models.The larger Knudsen layer observed for the VHS model in Figure 7 is due to the larger mean free path, and in turn causes larger shear stresses and viscosities.Agreement between the VHS model viscosities computed from DSMC and theory is better than that of the LJPA model.The average error in VHS model viscosity relative to theory is 1.09% while for the LJPA model the error is 3.45%.The Maitland and Smith 8 fit of viscosity measurements shown in Figure 7 indicate the LJPA model viscosities are in better agreement than the VHS with errors of 2.06% and 9.51% for LJPA and VHS models, respectively.The reason for the larger deviation in viscosity computed from the LJPA model in DSMC relative to theory is as yet unclear.
Several additional cases have been analyzed in order to determine the source of this error including: varied number of simulated particles per cell, varied Knudsen number, varied minimum cut-off angle, and varied time step and cell size.The DSMC computed viscosities and viscosity errors relative to theory are presented in Table IV for each of these diagnostic cases.Most of the cases show little effect on the viscosity errors between DSMC and theory, but a sufficiently small time step and cell size appears to be an important factor.Decreasing the cell size to 1.43 μm, or 7% of the mean free path, and the time step to 0.90 ns resulted in a -6.64% error between DSMC and theoretical viscosity.This is about a percent lower in magnitude, and is therefore significant.Furthermore, decreasing both the cell size and time step to just 5% of the mean free path and mean  collision time, respectively, reduced the viscosity error to -2.37%.Although the ratios of cell size to mean free path and time step to mean collision time were kept constant between 900 K and 1500 K cases, the no-time counter 9 (NTC) method requires these smaller cell sizes and time steps.Case 3 requires the most computational time of the three cases due to the higher number of collision events.One million sampling time-steps take 4016 s for the LJPA model and nearly half that, 2370 s, for the VHS model.

V. COUETTE SIMULATIONS IN THE TRANSITIONAL REGIME
The low-temperature, supersonic, Couette flow simulations of Sec.IV B are repeated here for Knudsen numbers ranging from slip to transitional regimes in order to study the importance of a realistic potential as Knudsen number increases.The wall temperatures are kept at a constant 40 K for each of the cases considered in this section with the Knudsen number varied through the initial number density.A domain length of 1 m and 20 particles per cell are again used, and the time-step and cell width are set to be approximately 15% of the mean collision time and 15% of the mean free path, respectively.These case conditions are summarized in Table V.
In the free-molecular limit, there are no collisions and hence the choice of intermolecular potential is inconsequential.The differences in shear stress between the LJ and VHS models are shown in Table VI to be tending towards zero as Knudsen number is increased-in agreement with the previous statement.
Another measure for quantifying the effect of the intermolecular potential on Couette flow is the difference in velocity distribution function relative to the equilibrium distribution function, The original 3D velocity distribution function has been integrated over the velocity in the z-direction, w, to produce the 2D distribution function, f (u, v).The shear stress is anyways in the xy-plane, and the w velocities are not affected.The equilibrium distribution function in Eq. ( 26) may be written A 1D velocity distribution, f (v), may similarly be obtained by further integration over the xvelocities, u.Doing so allows for the computation of statistical moments such as mean v-velocities and skewness of the distribution function, which also indicate the degree of non-equilibrium.The skewness, is a third-order moment where the expectations are integrations over the distribution function, The distribution function was computed for both VHS and LJ models in Case A and Case C using more than 700× 10 6 samples in 150 × 500 × 150(u × v × w) velocity bins at the spatial location: X = 10 −3 m.Several observations can be made from Figure 8. First, both L-J and VHS models have similar errors relative to the u-velocity at the lower Knudsen number of 0.0051.Also, the VHS model errors are shifted to higher v-velocities relative to the L-J errors.This is an indication of higher skewness in the distribution function for v-velocities, and is confirmed when the skewness is computed from Eq. ( 28).The skewness in the f (v) distribution function for the VHS model is 0.577nearly four times greater than that of the L-J model.The equilibrium distribution function has zero skewness, and therefore larger values of skewness indicate larger deviation from equilibrium.This result is in agreement with previous discussions about the relative sizes of the Knudsen layer-the VHS model is predicting a larger degree of nonequilibrium and therefore has a larger Knudsen layer.More can be deduced from the errors in the 2D velocity distribution functions at a Knudsen number of 0.051.One obvious difference in the errors at the higher Knudsen number is the distortion in the contours for both L-J and VHS models.This indicates that under these conditions both U-and Vvelocities are affected such that the distribution function is skewed in both directions.The skewness is larger for this case indicating a larger degree of nonequilibrium; with the skewness for the L-J and VHS models being 0.605 and 0.723, respectively.However, both models are similarly predicting high degrees of nonequilibrium, and the contours of distribution function errors in Figure 8 show close agreement between the two models.
A relationship between the degree of nonequilibrium and intermolecular potentials is now clear.The differences in degree of nonequilibrium between the L-J and VHS models is greatest at lower Knudsen numbers and thus results in the largest differences in shear stress and viscosity.As the Knudsen number is increased towards the free-molecular limit both models will indicate the same degree of nonequilibrium, and the choice of intermolecular potential on shear stress and viscosity will be negligible.

VI. CONCLUSIONS
Effects of intermolecular interaction models are evaluated in DSMC through the analysis of compressible Couette flow in the slip and transitional regimes.Flows in the continuum and slip regimes involving wide temperature ranges are expected to be most affected by the choice of molecular models.This includes flows with moderate to strong shock waves, in-space propulsion and vacuum technology.The theory of polynomial approximations was used to represent the scattering angle for the LJ potential using a two-variable polynomial of degree 7. The formulation for the scattering angle was first verified via equilibrium collision frequency calculations for argon.The agreement for viscosity between the two models is within 5% at a temperature of 273 K. On the other hand, at a temperature of 40 K the higher viscosity of the VHS model leads to a shear stress that is about 28% higher than the LJPA shear stress due to the importance of the attractive component of the force between molecules.At 1000 K, the VHS model predicts a 12% higher shear stress than the LJPA model, and again over-predicts the viscosity relative to experimental measurements.The DSMC simulations performed using the LJPA model were shown to reproduce the theoretical LJ potential viscosity of all three cases to within 3.5%.
The large difference in shear stress between the VHS and LJ models observed for the lowtemperature Couette flow diminishes as the Knudsen number is increased.Only a 7% difference is observed for a Knudsen number of 1, and this trend is in agreement with what is expected as the flow approaches the free-molecular limit.Comparisons of 2D velocity distribution functions relative to the equilibrium distribution function for the VHS and LJ models also indicate more significant differences at lower Knudsen numbers.A larger skewness in the velocity distribution function for the VHS model provides more evidence that the VHS model experiences a higher degree of nonequilibrium at the lower Knudsen numbers.In the free-molecular limit, there are no molecular collisions and hence the intermolecular potential plays no role in the velocity distribution functions.
The LJPA model has a computational cost that is comparable to the VHS and GSS models with an increase largely attributed to the higher number of collisions due to long-range nature of the LJ potential.The polynomial expansion approach enhances the fidelity of DSMC simulations and may be applied to other complex intermolecular potentials.
FIG.1.Schematic of a typical collision for purely repulsive models such as the HS, VHS models.
) with the purpose of improving agreement to viscosity measurements.In the previous expression, k is the Boltzmann constant and denotes the Gamma function.The three VHS parameters are: d ref the reference diameter, T ref the reference temperature, and ω the viscosity temperature exponent.A value of 0.5 for ω corresponds to the HS model.The VSS model uses a total of four molecular model parameters which includes the three VHS parameters d ref , T ref , ω, and an additional scattering parameter, α.The scattering angle for the VSS model, 9 FIG.2.Schematic of a typical collision at large impact parameter for models such as the Lennard-Jones potential which include the long-range attractive component of force between neutral molecules.

FIG. 3 .
FIG.3.Comparison of argon viscosity obtained using the VHS model,9 LJ potential,23,24 Maitland-Smith fit,8 and experiments by Kestin et al.25 and Vogel et al.27 FIG. 4. Comparison of scattering angle contours for the (a) LJ potential and (b) VHS model.

FIG. 6 .
FIG. 6.Comparison of (a) normalized temperatures for Cases 1-3, (b) velocity variation for Case 1, (c) velocity variation for Case 2, and (d) velocity variation for Case 3, in the gap obtained using VHS and LJPA models.

FIG. 7 .
FIG. 7. Comparison of (a) shear stress for Case 1, (b) viscosity variation for Case 1, (c) shear stress for Case 2, (d) viscosity variation for Case 2, (e) shear stress for Case 3, and (f) viscosity variation for Case 3, in the gap obtained using VHS and LJPA models.

Figure 7
compares the shear stress variation in the gap, which has to be constant across the gap for the Couette flow, obtained using the VHS model and the LJPA model.Since the viscosity obtained using VHS model and LJ potential agree very well at 273 K the difference in shear stress This article is copyrighted as indicated in the article.Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.210.206.145On: Wed, 23 Dec 2015 21:08:35 107102-12 Weaver, Venkattraman, and Alexeenko Phys.Fluids 26, 107102 (2014)

TABLE I .
LJ Collision frequencies computed from DSMC and theory.

TABLE II .
VSS collision frequencies computed from DSMC and theory.

TABLE III .
Summary of flow conditions and numerical parameters used for the subsonic and supersonic Couette flow cases.

TABLE IV .
Errors in viscosities computed from L-J model in DSMC relative to theory at 1500 K and τ eq = 8.49ns, λ = 20.1μm.

TABLE V .
Summary of flow conditions and numerical parameters used for the Couette flow in the transitional regime.

TABLE VI .
Shear stresses and viscosities computed from LJ and VHS models in the slip and transitional regimes.