By Terry Wood,2014-08-22 13:41
10 views 0

    Monte Carlo simulation of nonlinear Couette ;ow in a dilute gas

     )a Jos ?e Mar?ıa Montanero

    Departamento de Electr?onica e Ingenier?ıa Electromec?anica, Universidad de Extremadura, E-06071 Badajoz, Spain

     )b ?o)c Andr?es Santos and Vicente Garz

    Departamento de F?ısica, Universidad de Extremadura, 2000 E-06071 Badajoz, Spain

    (February 1, 2008) The Direct Simulation Monte Carlo method is applied to solve the Boltzmann equation in the Aug steady planar Couette ow for Maxwell molecules and hard spheres. Nonequilibrium boundary con-

    ditions based on the solution of the Bhatnagar-Gross-Krook (BGK) model for the Couette ow are 18 employed to diminish the inuence of finite-size efects. Non-Newtonian properties are characterized

    by five independent generalized transport coefcients: a viscosity function, a thermal conductivity

    function, two viscometric functions, and a cross coefcient measuring the heat ux orthogonal to the thermal gradient. These coefcients depend nonlinearly on the shear rate. The simulation results

    are compared with theoretical predictions given by the Grad method and the BGK and the ellip-

    soidal statistical (ES) models. It is found that the kinetic models present a good agreement with

    the simulation, especially in the case of the ES model, while the Grad method is only qualitatively

    reliable for the momentum transport. In addition, the velocity distribution function is also measured

    and compared with the BGK and ES distributions.

    PACS numbers: 47.50.+d, 51.10.+y, 05.20.Dd, 05.60.-k


    One of the most interesting states for analyzing transport processes far from equilibrium is the steady planar Couette [cond-mat.stat-mech] ow. The physical situation corresponds to a system enclosed between two infinite, parallel plates in relative motion and, in general, kept at diferent temperatures. These boundary conditions lead to combined heat and momentum transport. If x and y denote the coordinates parallel to the ow and orthogonal to the plates, respectively, then the corresponding steady hydrodynamic balance equations are

    @P @Pyy xy = (1) = 0; @y @y

    @u @qy x Pxy + (2) = 0; @y @y

    where u = ubx is the ow velocity, P is the pressure tensor, and q = q bx + q by is the heat ux. The presence of q in x x yyEq. (2) indicates that a thermal gradient @T=@y is induced by the velocity gradient, even if both plates are kept


    the same temperature. The balance equations (1) and (2) do not constitute a closed set unless the dependence of the pressure tensor and the heat ux on the hydrodynamic fields is known. If the gradients are small, the uxes P and q are described by the Navier-Stokes (NS) constitutive relations, which in this problem yield @u x ; Pxx = P = P; Pxy (3) =;). yyzz0 @y arXiv:cond-mat/0003364v3

    @T y ; qq = 0; (4) =;)1 x0 @y

    where . and 1 are the NS shear viscosity and thermal conductivity coefcients, respectively. As a consequence of 0 0 the absence of normal stress diferences in the NS description, the hydrostatic pressure p = (P + P + P =) 3 is xx yy zza constant, on account of the balance equation (1). Even in the linear regime described by the NS equations, one still needs to know the spatial dependence of the transport coefcients to obtain the exact solution of the hydrodynamic equations. The problem becomes tractable


in the case of a low density gas, where the state of the system is completely specified by the velocity distribution 1function (r v;), which obeys the Boltzmann equation. A relevant dimensionless quantity in a dilute gas is the f; t

    Knudsen number Kn = 2=`, defined as the ratio of the mean free path 2 to the scale length of the hhydrodynamic

    gradients `. In many laboratory conditions, Kn;; 1 and so the Boltzmann equation can be solved by means of the h 2Chapman-Enskog method as an expansion of the distribution function in powers of the Knudsen number. The zeroth order approximation leads to the Euler hydrodynamic equations, while the first order approximation yields the NS equations with explicit expressions for the transport coefcients . and 1. The results show that the ratio . =1 is 0 0 0 0a constant. Consequently, it then follows from the NS hydrodynamic equations (1)-(4) that the ow velocity profile @u x is quasi-linear, . (5) = const; 0 @y

    and the temperature is quasi-parabolic,

    / 02 / 02 1@u@ 0 x 1 .0 (6) = const: T =;); 0 @y 0 @y .

    through the temperature. Anal- Note that the profile of u is not strictly linear, due to the space dependence of . x 0 ogously, the temperature profile is not strictly quadratic. In fact, the specific form of both profiles depends on the interaction potential under consideration. On the other hand, from Eqs. (5) and (6) it is easy to derive a nice

    result, namely that if the temperature T is seen as a function of u rather than as a function of the coordinate space y, then x

    one has 2.0 @T (7) =;) : @u2 x 0 1

    This is a sort of nonequilibrium "equation of state," according to which the temperature is a quadratic function of the ow velocity. Moreover, the "curvature" of the profile is practically universal, given the weak inuence of the

     2interaction potential on the Prandtl number Pr = 5k . =2m1 0;;?3 , where kB is the Boltzmann constant and m is B 0the mass of a particle.

    Since the mean free path is inversely proportional to the density, the Knudsen number, at a given value of the scale length `, increases as the gas becomes more rarefied. In general, when the Knudsen number is not small, the NS hrelations are not expected to hold and the transport must be described by nonlinear constitutive equations. In the )4 ,3 4special case of Maxwell molecules (particles interacting via an r potential), it has been shown that the Boltzmann equation admits a consistent solution in the nonlinear Couette ow characterized by a constant pressure p and profiles similar to those obtained in the NS regime, Eqs. (5)-(7), except that . and 1 are replaced by a generalized shear 00viscosity coefcient .(a) = . F (a) and a generalized thermal conductivity coefcient 1 (a) = 1F (a), respectively. 0 yy 0 ?Here, a = (. =p)@u=@y is a constant (dimensionless) shear rate and F and F are nonlinear functions of a. In 0 x? = q = In this problem, the hydrodynamic scale length can be identified as = andaddition, P;;;: P P; x xxyy zz p p )1 ;0. ;`h k T=m(@u k T=m(. , while the mean free path isThus, the reduced shear rate a ,;; B B ;;;,; 2represents the Knudsen number in this problem, i.e. a;;, Kn. Henceforth, we will use the reduced shear rate a to =p). =@y) 0xrefer to the Knudsen number Kn. The solution considered in Refs. 3,4 describes heat and momentum transport for arbitrary velocity and thermal gradients in the bulk domain, where boundary efects are negligible. On the other hand, the full nonlinear dependence of F (a) and F (a) is not explicitly known, since it involves the infinite hierarchy ? 4 2of moment equations. Their knowledge is limited to super-Burnett order and the result is F (a) = 1;) 3:111a and 2F (a) = 1;) 7:259a. ?Consequently, if one wants to get the transport properties for arbitrary values of the shear rate and the thermal gradient, one must resort to approximate schemes or to computer simulations. In the first alternative, explicit expressions for the nonlinear transport coefcients in the Couette ow have been obtained from exact solutions of 5, 697 the Bhatnagar-Gross-Krook (BGK) model and related models - for general interactions, as well as from the Grad 1010method. In the simulation side, Risso and Cordero have recently studied the shear-rate dependence of F and F ? by means of molecular dynamics simulations of a hard disk gas. Comparison between the diferent analytical results with those obtained from the simulation shows that the predictions given by kinetic models are in better agreement 9than those given by the Grad method, especially in the case of the thermal conductivity. Nevertheless, given the difculties associated with molecular dynamics simulations to achieve large shear rates in a dilute gas, the above comparison is restricted to a range of shear rates for which non-Newtonian efects are hardly significant. For instance, the shear viscosity has only decreased around 10% with respect to its Navier-Stokes value for the largest value of the 10 shear rate considered by Risso and Cordero. In order to overcome such difculties and extend the range of values


    11 of a, one may use the so-called Direct Simulation Monte Carlo (DSMC) method, which is known to qualify as an efcient and accurate method to numerically solve the Boltzmann equation.

    The aim of this paper is to solve the Boltzmann equation by means of the DSMC method for a gas subjected to

    the planar Couette ow. The motivation of this work is twofold. On the one hand, we want to test the reliability of

    the far from equilibrium results obtained from kinetic models and the Grad method by making a comparison with the Boltzmann solution in the case of Maxwell molecules, for which the form of the hydrodynamic profiles in the bulk region (far from the boundaries) is known. We will determine not only the hydrodynamic profiles but also the

    nonlinear transport coefcients and the velocity distribution function. On the other hand, we want to investigate whether the above results for a system of Maxwell molecules extend to other interaction potentials. This extension holds when the Boltzmann equation is replaced by kinetic models where, in terms of a conveniently scaled space variable, all the results are independent of the interaction law. Thus, we will also solve numerically the Boltzmann equation by the DSMC method for a hard-sphere gas.

    Since we are interested in describing transport properties in the bulk region, i.e., free from finite-size efects, we need to use appropriate boundary conditions in the simulations. In the conventional boundary conditions, 10,12 the gas is assumed to be enclosed between two baths at equilibrium in relative motion and, in general, at diferent temperatures. Under these conditions, a particle leaving the system is formally replaced by a particle coming from the bath, so the incoming velocity is sampled from a local equilibrium distribution. As a consequence, there exists a mismatch between the velocity distribution function of the reemitted particles and that of those particles of the gas located near

    the wall and moving along the same direction. In this case, in order to inhibit the inuence of boundary efects


    needs to take very large systems (normal distance between the plates much larger than the mean free path), what is not practical from a computational point of view. To overcome this difculty, a possibility is to assume that both baths are out of equilibrium in a state close to that of the actual gas. Since such a state is not known "a priori",

    in this paper we assume that the state of the baths is described by the BGK solution of the planar Couette 6 ow.

    Although the boundary efects are still unavoidable, we expect that the above mismatch between reemitted and gas particles will be much smaller. As a matter of fact, the use of these conditions has been shown to be appropriate to 13 analyze bulk transport properties in the special case of planar Fourier ow (both walls at rest). The plan of the paper is as follows. In Sec. II we give a brief description of the planar Couette ow and a summary of the main results obtained from the Boltzmann equation and kinetic models. The boundary conditions used in the simulations and the DSMC method are described in Sec. III. Section IV presents the main results of the paper, where special attention is paid to the nonlinear transport coefcients. A comparison with the analytical results derived from the BGK and the ellipsoidal statistical (ES) models and from the Grad method is also carried out. The comparison II. DESCRIPTION OF THE PROBLEM AND SUMMARY OF THEORETICAL RESULTS shows in general a good agreement of the kinetic models with computer simulations, even for large shear rates. In addition, the velocity distribution function obtained from the simulation in the bulk domain is compared with the Let us consider a dilute gas. In this case, a kinetic description is sufcient to characterize the state of the ones given by the kinetic models. It is shown again that the agreement is qualitatively good. We close the paper system in by means of the velocity distribution function f (r; v; t). This distribution function obeys the nonlinear Boltzmann Sec. V with some concluding remarks. 1 equation, which in the absence of external forces is given by @f (8) f = J[f; f ]; + v;);?@t

    where Z Z ;,,;;dbk gI(g; bk) [f (v)f (v );) f (v)f (v )] dv (9) J[f; f ] = 1 11

    is the collision operator. In this equation, I(g; bk) is the diferential cross section, g;;!;,v)v being the relative velocity, 1;,and (v; v are precollisional velocities yielding postcollisional velocities (v; v 1). From the distribution function, one 10

    may define the hydrodynamic quantities


     (10) n = dv f;

    Z 1 u = (11) dv v f; n


    Z m 3 2dv Vf; nk T (12) B2 2 =

    and the momentum and heat uxes


     (13) P = m dv VV f;

    Z m 2dv VV f: q = (14) 2

    Here, n is the number density, u is the ow velocity, T is the temperature, P is the pressure tensor, q is the heat ux, and V = v;) u is the peculiar velocity. In addition, the equation of state is that of an ideal gas, i.e., p = nk T . B14Most of the known solutions to Eq. (8) for spatially inhomogeneous states correspond to the special case of )4Maxwell molecules, namely, a repulsive potential of the form V (r);, r . For this potential, the collision rate gI(g; bk) is independent of the relative velocity and this allows the infinite hierarchy of velocity moments to be recursively solved in some specific situations. Furthermore, the NS transport coefcients . and 1 can be exactly obtained from 0 0 2the Chapman-Enskog method. They are given by

    p 15 k B 0 . . ;0 ; (15) 0 = 1= 4 4 m

    15 where 4 = /n, / being an eigenvalue of the linearized Boltzmann collision operator.

    In this paper we are interested in studying the planar Couette ow for a dilute gas. We consider a gas


    between two parallel plates in relative motion and maintained at diferent temperatures. Under these conditions, the system is driven out of equilibrium by the combined action of the velocity and thermal gradients along the


    normal to the plates. After a transient period, the gas is expected to reach a steady state and the corresponding @ Boltzmann equation reads v f(16) = J[f; f ]; @y y

    where we have chosen the axis y as the one normal to the plates. In general, this equation must be solved

    subjected 2to specific boundary conditions. Nevertheless, in the same spirit as in the Chapman-Enskog method, one may look for "normal" solutions in which all the space dependence of the distribution function occurs through a functional dependence on the hydrodynamic fields. The normal solution describes the state of the gas in the hydrodynamic regime, namely, for times much longer than the mean free time and for distances from the walls much larger than ,3 4the mean free path. As mentioned in the Introduction, it has been proved that, in the case of Maxwell


    Eq. (16) admits an exact solution corresponding to the planar Couette ow problem. This solution belongs to the normal class, so that no explicit boundary conditions appear. In other words, all the space dependence of f is given through the local density, the local velocity, the local temperature, and their gradients. The solution is characterized by hydrodynamic profiles that are a simple generalization of those predicted by the NS approximation, Eqs. (5)-(7),

    (17) p = const; namely

    @ux 1 (18) a = const; ; (y) @y 4

    1 22 @ 1 2m (19) a): (T =;)Pr (y) @y 4kB

    The dimensionless parameter (a) is a nonlinear function of the reduced shear rate a that, by construction, behaves 2as ;;: a=5 in the limit a;;; 0. Again, the temperature can be seen as a quadratic function of the ow velocity, 2but now the coefcient . =1 appearing in Eq. (7) is replaced by a shear-rate dependent coefcient (. 0=1 0)5(a)=a. 0 0Furthermore, in this solution the pressure tensor is independent of the thermal gradient and the heat ux vector is linear in the thermal gradient, but all these uxes are nonlinear functions of the shear rate. This nonlinear dependence


can be characterized through five generalized transport coefcients. First, the shear stress P defines a generalized xy

    shear viscosity .(a) as

    @u @ux x : (20) P =;).(a) !;). F (a) ;xy 0 @y @y

    Analogously, the component of the heat ux parallel to the thermal gradient defines a generalized thermal conductivity coefcient 1 (a): yy

    @T @T : (21) (a) q =;)1 !;)1F (a) ;yyy 0 ?@y @y

    The dimensionless functions F (a) and F (a) are the most relevant quantities of the problem. They are related to ? 2the function (a) by (a) = aF (a)=5F (a). Normal stress diferences are diferent from zero and are measured by ?a) defined by (two viscometric functions ( ,21

    P ) P; 2yy xx (a)a; = ((22) 1p ;

     Pzz 2) P;= ( (a)a: yy (23) 2 p ;

    . To NS 3(a) =;)P=P Another interesting quantity related to the pressure tensor is the friction coefcientxyyy we simply have 3(a) = a. In the non-Newtonian regime, we generalize this coefcient as 3(a) = aF (a), where the order, µfriction function F (a) is µ

    F (a) : F(a) = (24) µ1;) [((a);) ((a)]a=3 2 1 2

    Finally, there exists a nonzero component of the heat ux orthogonal to the thermal gradient given by

    @ @ T: (25) (a) q =;)1 ;;!;)1~(a)a Txxy 0@y @y

    a) and ~(a) are generalizations of Burnett coefcients. In fact, ((0) =;)14=5, ((0) = 4=5, (The three functions ( 1 2 ,21 2and ~(0) =;)7=2 for Maxwell molecules and for hard spheres in the first Sonine approximation. The determination of , and ~ would imply the solution of an infinite hierarchy that the nonlinear transport coefcients F , F , ( 2, ?13 cannot This hierarchy can only be solved step by step when one performs a perturbation be solved in a recursive way. 4expansion in powers of the shear rate. In fact, Tij and Santos determined the solution up to super-Burnett order:

    4 2F (a) = 1;) 3:111a +;?(a ); (26)

    4 2(a) = 1;) 7:259a +;?(a ): F (27) ? 2 2 2 43 (a) = (a=5)[1 + 4:148a +;?(a)]. In addition, F µ(a) = 1;) 1:911a +;?(a ). From Eqs. (26) and (27) it follows that

    Although the above analyses are valuable, they have two main limitations. On the one hand, they are restricted to the special case of Maxwell molecules. For other interaction potentials (e.g., hard spheres), the collisional moments involve all the moments of the distribution function and, as a consequence, the hydrodynamic profiles (17)-(19) are not strictly true. On the other hand, even for Maxwell molecules, the above perturbative solution is not useful for finite shear rates. These two limitations can be overcome, in the context of analytical methods, by introducing 10 ,975 additional approximations, such as the Grad method, or by describing the system by means of kinetic models. - In these approaches, one looks for a solution having the same hydrodynamic profiles as in the case of the Boltzmann

    equation, cf. Eqs. (17)-(19). As before, this solution describes the properties of the system in the bulk region, which is insensitive to the details of the boundary conditions. From the Grad method and from the BGK and ES kinetic models, one explicitly gets the full nonlinear shear-rate dependence of the transport coefcients in a non-perturbative way. Moreover, the results are universal , namely the functions F (a), F (a), . . . are independent of the interaction ?potential, provided that the reduced shear rate is defined as in Eq. (18) with 4 = p=. . 0The thirteen-moment Grad method consists of replacing the actual distribution by


    7 1/ 0 28 mV m 1 ; 1 + 2 P (28) (f;;; f ) 1 ;V;) q + ) pffi)V V;L ij ij ij 2 n (k 5k T ) 2 Bwhere T B/ 0 / 0 /23 mV m f = n exp 2 (29) );;L26k T 2k B

    is the local equilibrium distribution function. When the approximation (28) is inserted into the Boltzmann equation T B(8) and velocity moments are taken, one gets a closed set of equations for n, u, P, and q. According to the geometry of the planar Couette ow, there are eight independent moments instead of the original thirteen moments appearing in Eq. (28). Risso and Cordero 10,16 found that the set of independent moment equations, neglecting nonlinear terms in the uxes, admits a solution consistent with the profiles (17)-(19). In addition, they obtained explicit expressions for the transport coefcients introduced above as nonlinear functions of the reduced shear rate a. Those expressions are displayed in the Appendix.

    Now we consider the kinetic model approach. The basic idea is to replace the detailed Boltzmann collision operator by a much simpler term that otherwise retains the main physical properties of the original operator J[f; f ]. The most 1 familiar choice is the BGK model,

    (30) J[f; f ];;;)4(f;;) f); L

    where 4 is an efective collision frequency that depends on the temperature, according to the interaction potential. The NS transport coefcients of the BGK model are . = p=4 and 1 = 5pk =2m4. Thus the BGK model has 0 0 Bthe

    drawback that it predicts an incorrect value for the Prandtl number, Pr = 1. This is a consequence of the fact that 4 is the only adjustable parameter in the model. This deficiency is corrected by the so-called ellipsoidal statistical (ES) 1model, in which case (31) J[f; f ];;;)-(f;;) f); R

where - is again an efective collision frequency and the reference function f is R_ )1 ;)3/2 ;)1/2 ;R (32) = n6 (det A) : VV ; fexp;;;)A

     )1 ;)1 ;)1 ;ij = (2k T=m)Pr where A ffi)2(Pr )1)P =mn. The NS coefcients are now . = p=(-Pr ;) and 1 = 5pk =2m-. B ij ij 0 0 B)1 If, as in Eq. (15), we define a collision frequency 4 = p=. , then 4 = -Pr in the ES model. Note that the ES 0model reduces to the BGK model if we formally make Pr = 1. Therefore, the ES model can be seen as an extension of the

    BGK model to account for the correct Prandtl number, which can be particularly relevant in the Couette ow where heat transport and momentum transport are coupled. The results derived from the BGK and the ES models for the Couette ow problem are also given in the Appendix. Apart from obtaining the transport properties, the velocity 6 distribution function can be explicitly written. In particular, the BGK distribution is given by / 0Z t /23 3/2) 1 ff(1 + ff 2m )2 5/dt 2t;) (1;) ff)f (r; v) = n 2t 26k , ;;,; t0 y5 " #) ( / 2 01 + ff T 1;) t 2aff 1;) t 2ff B2 2 : x + + 5 + (33) );;);;5 exp ; y2 1 + ff ,5 2t;) (1;) ff)t y 1 + ff , 5z

     /Here, (t 0; t 1) = (0; 1) if 5 y > 0 and (t 0; t 1) = [1; 2=(1;) ff)] if 5 y < 0. Besides, ?;;! (m=2k BT ) V,21

    , 1/2; (34) ff = ( + 2,

    8) and

    / 0 /21 2k 1 @ (35) , = T lnT 4 @y B

    m is a (local) reduced thermal gradient. Equation (33) shows that the distribution function is a highly nonlinear function of the reduced gradients a and ,.


    In Ref. 9 a comparison between the analytical results obtained from the Grad method and the BGK and ES models 10with those obtained from molecular dynamics simulations for hard disks was carried out. It was found that the three theories reproduced quite well the simulation data for F , but the Grad method failed for F and ~. The ?latter quantity was reproduced better by the ES model than by the BGK model. Notwithstanding this, more definite conclusions could not be reached because the simulation data were restricted to rather small shear rates, namely

     a < :2. In this range of shear rates, non-Newtonian efects are not especially significant. In addition, although the ; 0

    molecular dynamics simulations correspond to a very dilute gas (area fraction _;: 1%), the collisional contributions to the transport coefcients (absent in a Boltzmann description) are not strictly zero. Finally, as a minor point, conclusions drawn in the context of two-dimensional systems should not be extrapolated without caution to the more realistic case of three dimensions. As said in the Introduction, the aim of this paper is to solve numerically the Boltzmann equation by means of the DSMC method for the planar Couette ow and compare the results with the theoretical predictions. We will consider three-dimensional systems of Maxwell molecules and hard spheres subjected to shear rates as large as a;? 1:2. In addition to the nonlinear transport coefcients, the comparison will be extended to the level of the velocity distribution function itself.


    A. Boundary conditions

    The goal now is to solve numerically the Boltzmann equation corresponding to the planar Couette ow by using 11the successful DSMC method. The gas is enclosed between two parallel plates located at y = 0 and y = L, which are moving along the x-direction with velocities U = Ubx and U = Ubx, respectively. In addition, they are kept 0 0L Lat temperatures T and T, respectively. In order to solve Eq. (16), we need to impose the corresponding boundary 0 L;,conditions. They can be expressed in terms of the kernels K 0,L (v; v) defined as follows. When a particle with velocity ;,;,v hits the wall at y = L, the probability of being reemitted with a velocity v within the range dv is K(v; v)dv; L;,17 the kernel K(v; v) represents the same but at y = 0. The boundary conditions are then 0Z ,;;,;;,;;?()v )v;;,f (y =;,0; L?; v) = ?()v ) dv;;,v y y y 0,L (v; v ) ;K y,;;?()v f( y =;,0;L?; v ;t): ;(36) y0 ;,(v; v) = In this paper we consider boundary conditions of complete accommodation with the walls, so that K ;0,L ;,(v) does not depend on the incoming velocity v and can be written K0,L as Z )1 ;K (37) (v) = A (v); A0,L = (v); ?()v )v dv ?()v )v ,L 0,L 0,L 00,L ; y y

    _ _ yywhere _ (v) represents the probability distribution of a fictitious gas in contact with the system at y =;;,0; L?. 0,L Equation (37) can then be interpreted as meaning that when a particle hits a wall, it is absorbed and then replaced (v) must be consistent with the imposed wall by a particle leaving the fictitious bath. Of course, any choice of _ 0,L velocities and temperatures, i.e., Z

    dv v _ (38) U0,L = (v); 0,L x

    Z 1 ) k m 0,L = (39) (v): dv (v;) U 0,L 0,L 23 _ T B

    The simplest and most common choice is that of a Maxwell-Boltzmann (MB) distribution:

    / 0 1 2 2 /23 ) m m (v;) U 0,L MB _ : (v) = (40) exp;;;); 2k 0,L 26k 0,L 0,L T BT Under these conditions, the system is understood to be enclosed between two independent bathsB at equilibrium. While 1 18 ,the conditions (40) are adequate for analyzing boundary efects, they are not very efcient when one is interested in obtaining the transport properties in the bulk region. In order to inhibit the inuence of boundary efects, it is more convenient to imagine that the two fictitious baths are in nonequilibrium states resembling the state of the actual


    gas near the walls. More specifically, we can assume that the fictitious gases are described by the BGK equation, whose exact solution for the steady planar Couette ow is given by Eq. (33). In this case, the probability distributions _0,L (v) are

    /2 Z 31 t 2_/2 m ff0,L 0,L) (1 + ff )5BGK )3/2 ;_ (v) = 6 t ) 2;) (1;)dtt ff 0,L 0,,L v k T 0,L 0,L t0 ( B / 0 ,; /21 1 + ff ym k 2ff 21;) t 0,L 0,L );;);;0,L exp ;T )0,L , Bt;) (1;) ff 22k 0,L 0,L 0,Lv 1 + ff 2t y #) "/ m 2 0T ;,2aff B1;) t 0,L 22 + v + ; 0,L + (41) v ;) U y;; x1 + ff 0,L ,0,L v z ;,1/2 =[. Here, a)] (+ 8 )] if v < 0, and = , where (t ; t ) = (0; 1) if v > 0 and (t ; t ) = [1; 2=(1;) ff y 0 0 1 0,L 0,L 0,L 0,L BGK y 2, 1 ;, a is the estimated value of the reduced shear rate, as predicted by the BGK model for specific values of the boundary ff ;,a) is obtained from Eq. (A6). Given the values of the four independent boundary (, and parameters U and T 0,L 0,L BGK ;,(as well as the distance L), the shear rate a and the local thermal gradients , are fixed parameters U and T 0,L 0,L 0,L by the conditions (17)-(19). Therefore,

    U ) U;L 0 a ; (42) ?;;=

    0 11 2 / /2,;;TL mBGK 1 2k T ) T;0,L 2 0 (a )Pr ;, : 0,L = (43) );;B ?; ?; T k 0,L ;m 0,L ;T BIn these equations,;;? is related to the actual separation L between the plates through the nonlinear equation

    Z ;;;?;ds ; (44) L = (s) 40

    where s is a variable in terms of which the temperature is a quadratic function, namely T (s) = T ,0 [1 + a )Prs =k T], and the s-dependence of the collision frequency 4(s) appears only through (0 B 01/2;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;, 2 (m=2k T) s;) m B 0BGK a , , , T , and L. Nevertheless, from a practical point of view, it is more convenient , and;;? for any choice of0,L 0,L 0,L the temperature (taking into account that p = const). The solution of the nonlinear set of equations (42)-(44) gives ;, U , a, and , to fix U, T as independent parameters. Without loss of generality we take U = 0 and T = 1. 0,L 0 00L In addition, we will choose , = 0. This implies that, if boundary efects were absent, the temperature would have a 0 maximum at the lower plate y = 0. Thus,

     1 2 /21 ? ;,= a?; U ; (45) ? = ;L a) (2Pr ;,BGK

    ? : (46) , =;)2 L?;; =. Finally, the actual distance L is given T;) 1 and we have taken m = 1 and k In the above equations ?;! 0B 21 by Eq. (44). However, from the simulation point of view, it is more useful to express L in terms of the collision frequency 4 corresponding to the temperature T and the average density n. In other words, instead of Eq. (44) we L use 1 Z 1 Z;;;?; L n (s) ds (47) n = dy n(y) = L L (s) 40 0

     ?to determine L. For the sake of concreteness, let us consider repulsive potentials, for which 4 = 4(n=n)(T=T) , L(21 hard spheres). In that case, where ! ranges from 0 (Maxwell molecules) to

    1 Z ;;;?;? ;): (48) L = ds [T (s)] ;4 0 ?; ?;;)1 For Maxwell molecules, this simply reduces to L =;?=4, while for hard spheres one has L = (?=4) tan ?. ( ?);, a and ?, the separation L between the plates and the velocity of the upper plate In summary, given the values of= ;U are uniquely determined. In addition, the thermal gradient at the upper plate , is also determined, while the LL value at the lower plate is fixed = 0. The knowledge of these boundary parameters allows one to obtain the BGK as , distributions _ (v), according to Eq. (41). 0 0,L


    B. The DSMC method

Now, we briey describe the numerical algorithm we have employed to solve the Boltzmann equation by means

    of 11? and positions;,y? of a sufciently large number of particles N . Given the geometry represented by the velocities;,v the so-called Direct Simulation Monte Carlo (DSMC) method. In this method, the velocity distribution function is i i

    of the problem, the physical system is split into layers of width ffiy, sufciently smaller than the mean free path.


    velocities and coordinates are updated from time t to time t+ffit, where the time step ffit is much smaller than the mean free time, by applying a streaming step followed by a collision step. In the streaming step, the particles are moved from the corresponding probability distribution K ballistically, y;;; y + v ffit. In addition, those particles crossing the boundaries are reentered with velocities sampled i i iy0,L (v). Suppose that a particle i crosses the lower plate between according to the following rules. First, a velocity e (with times t and t + ffit, i.e., y + v ffit < 0. Then, regardless of the incoming velocity v , a new velocity e is iv i iy iiv(v). If one is considering the MB boundary conditions, this velocity is accepted directly. Otherwise, the above vassigned iye > 0) is sampled with a probability proportional to MB v _ y0acts as a "filter" to optimize the acceptance-rejection procedure and the velocity e is accepted with a probability ivBGKMBproportional to the ratio _ (v)=_ (v). If the velocity is rejected, a new velocity e is sampled and the 00ivprocess

    is repeated until acceptance. The new position is assigned as v (ffit + y=v ). The process is analogous in the case iye i iyof the upper plate.

    The collision step proceeds as follows For each layer ff, a pair of potential collision partners i and j are chosen at random with equiprobability. The collision between those particles is then accepted with a probability equal to the corresponding collision rate times ffit. For hard spheres, the collision rate is proportional to the relative velocity v ;;) v , while it is independent of the relative velocity for Maxwell molecules (an angle cut-of is needed in i jthe

    latter case). If the collision is accepted, the scattering direction is randomly chosen according to the interaction law and post-collisional velocities are assigned to both particles, according to the conservation of momentum and energy. N N nL After the collision is processed or if the pair is rejected, the routine moves again to the choice of a new pair ? X n = n = y); (49) (?? ? iuntil the (N=L)ffiNffi=1 irequired number of candidate pairs has been taken. y y In the course of the simulations, the following "coarse-grained" local quantities are computed. The number density where ?(y) is the characteristic function of layer ff, i.e., ?(y) = 1 if y belongs to layer ff and is zero ? ?in layer ff is otherwise.

    Similarly, the ow velocity, the temperature, the pressure tensor, and the heat ux of layer ff are N 1 X ? = y)v ; (50) (u ?? i iN ? =1 i

    N p m ? 2 X k T = (51) ? = ?? (y)(v ;) u ) ; i i ?B 3N n? ? =1 i

    N L X m ? = (52) P ?? (y)(v ;) u )(v i;) u i i ?Nffi=1 i); ?y

    N m X L 2 ? = (53) ?? q (y)(v ;) u ) (v i;) u i i ?Nffiy i=1 ): ?2

    From these quantities one can get local values of the gradients and of the transport coefcients. For instance, the reduced shear rate is

    / 0? n TL u?+1,x ) u;?,x a ? = (54) ffi? 4n T? y and the viscosity function is


    P?,xy : (55) F =;); ,?a p ? ? 21 . It remains to fix the As said before, in the simulations we take units such that m = 1, T = 1, and k L B unit or, equivalently, the length unit. The standard definition of mean free path in the case of hard spheres time 2 is

    (56) 1 n68 22 2 =;;? ;

    where 8 is the diameter of the spheres. The Navier-Stokes shear viscosity is (in the first Sonine approximation) 2 1/2 .0 = 5(mk and the mean free path 2 are =168 . Consequently, the efective collision frequency 4 = p=.;? 0 =2 1/2) related as 4 = (8=5 6)(2k BT=m T=6) . As usual, we choose the mean free path corresponding to the average density B?; n as the length unit. This in turn implies that 4 = 8=5 6;;? 0:903. For convenience, we take the latter value

    for 5Maxwell molecules as well. The typical values of the simulation parameters are N = 2;) 10 particles, a layer width

     = 0:02, and a time step ffit = 0:003. ffiy

    The procedure to measure the relevant quantities of the problem is as follows. First the values of the imposed shear ;,rate a and the temperature diference ? are chosen. This choice fixes the system size L, as well as the upper velocity

     and the upper thermal gradient ,, according to Eqs. (45), (46), and (48). Starting from an equilibrium initial UL Lstate with T = T, the system evolves driven by the boundary conditions described before. After a transient period 0(typically up to t = 25), the system reaches a steady state in which the quantities uctuate around constant values. In this state, the balance equations predict that the quantities u, P , P , and uP + q are spatially uniform and yxyyy xxy ythis is used in the simulations as a test to make sure that the steady state has been achieved. Once the steady state is reached, the local quantities (49)-(55) are averaged over typically 100 snapshots equally spaced between t = 25 and

     = 55, which corresponds to about 60 collisions per particle in the case of hard spheres. t

    C. Test of the numerical algorithm

    Before closing this Section, it is worthwhile carrying out a test of the reliability of the numerical method. To that

    end, we have simulated the BGK equation by a DSMC-like method similar to the one described in Ref. 19. If

    the ;,consider the hard-sphere situation with a = 1 and ? = 5. This corresponds to = 0 :248, L = 1:81, U = :17, 3boundary conditions are implemented correctly and the simulation parameters are well chosen, then the simulation BGKL

    =;)3:15, where in this case Pr = 1. Figure 1 shows the marginal velocity distributions of particles reemitted results should agree with the theoretical BGK predictions. We have checked that this indeed the case. As an example, and , L by the walls,

    / 0 Z1 Z;;;(; /2(;;2k T 0,L dv dvx (57) (v); ?0,L ( ) ,L 0B z y5)(;;)(;;K m = ;/2 1vy 0,L) . The case 5 > 0 (5 < 0) corresponds to particles that are reemitted from the as functions of 5 = (m=2k T y B y ylower (upper) plate. The agreement with the imposed distribution is excellent. Note the strong asymmetry between the distributions corresponding to both plates, in contrast to the symmetry of the MB distributions obtained from (40). The temperature and velocity profiles are shown in Fig. 2. The simulation values overlap, within statistical uctuations, with the theoretical predictions. Apart from the profiles, we have verified that the generalized transport coefcients obtained from the uxes also agree with the theory. A more stringent test is provided in Fig. 3, where the marginal distribution function

     Z1 / 0 / 0Z;;;(; /2(;; /21 m 1 k T 2vdvx 5y y ' (5 ) f (v); (58) dv B y zn 2k )(;;)(;;= m = T is plotted at the point y=L = 0:5, which corresponds to a local thermal gradient , =;)0:60. It is apparent again that Ban excellent agreement exists between simulation and theory.


    This Section is devoted to a comparison between the simulation results for the Boltzmann equation obtained by the simulation method described in the previous Section and the theoretical predictions provided by the Grad method and the BGK and ES kinetic models. The comparison will be carried out at the levels of the transport coefcients and the velocity distribution, both for Maxwell molecules and hard spheres. Before that, the hydrodynamic profiles obtained from the simulations by using the two types of boundary conditions considered in Sec. III are presented.


Report this document

For any questions or suggestions please email