Discrete Boltzmann Method with Maxwell-Type Boundary Condition for Slip Flow?

Yu-Dong Zhang(張玉東)Ai-Guo Xu(許愛國(guó))? Guang-Cai Zhang(張廣財(cái)) and Zhi-Hua Chen(陳志華)

1Laboratory of Computational Physics,Institute of Applied Physics and Computational Mathematics,Beijing 100088,China

2Key Laboratory of Transient Physics,Nanjing University of Science and Technology,Nanjing 210094,China

3Center for Applied Physics and Technology,MOE Key Center for High Energy Density Physics Simulations,College of Engineering,Peking University,Beijing 100871,China

1 Introduction

In recent years,the development of natural science and engineering technology has moved towards miniaturization.One of the most typical examples is Micro-Electro-Mechanical System(MEMS).[1?4]It is extremely important to investigate the underlying physics of unconventional phenomena at the micro-scale.Those unconventional phenomena cannot be explained by the traditional macro-model and has become a key bottleneck limiting the further development of MEMS.Among these unconventional physical problems,the gas flow and heat transfer characteristics at the mico-scale are especially critical.

Due to the reduction of the geometric scale,the mean free path of gas molecules may be comparable to the length scale of the device.The Knudsen number(Kn),a dimensionless parameter used to measure the degree of rarefaction of the flow and de fined as the ratio of the mean free path of molecules to characteristic length of the device,may be larger than 0.001 and reaches slipflow regime(0.001

In practical,gas flow in a microchannel may encounter continuum,slip and transition regimes simultaneously.Traditionally macro-scale models cannot apply to such a broad range of Knudsen numbers by only one set of equations.Besides,the numerical solutions of higher order macro-equations are difficult to obtain because of the numerical stability problem.[3]It is commonly accepted that Direct Simulation Monte Carlo(DSMC)[4,8]is an accurate method for rare fied gas flow,which has been veri fied by experimental data.However,the computation cost in numerical simulations is too expensive for low speed gas flow.To reduce the huge ratio of the noise to the useful information,extremely large sample size is needed.Although,the information preservation(IP)method was presented to treat this problem,[9]the contradiction between the noise problem and sample size has not been well solved.

It has been know that rare fied gas dynamics are represented properly by the Boltzmann equation due to its kinetic nature.That is continuum,slip and transition regimes can be described by one equation.Unfortunately,the Boltzmann equation is a 6-dimensional problem and the computational cost is formidable to solve such an equation by numerical method directly.In order to alleviate the heavy computational burden of directly solving Boltzmann equation,a variety of Boltzmann equationbased methods,such as the uni fied gas-kinetic scheme(UGKS),[10?11]the discrete velocity method(DVM),[12]the discrete uni fied gas-kinetic scheme(DUGKS),[13?14]the lattice Boltzmann method(LBM),[2,15?26]have been presented and well developed.Recently,Discrete Boltzmann Method(DBM)has also been developed and widely used in various complex flow simulations,[27?29]such as multiphase flows,[30]flow instability,[31?32]combustion and detonation,[33?34]etc. From the viewpoint of numerical algorithm,similar to finite-different LBM,the velocity space is substituted by a limited number of particle velocities in DBM.However,the discrete distribution function in DBM satis fies more moment relations,which make it fully compatible with the macroscopic hydrodynamic equations including energy equation.The macroscopic quantities,including density,momentum,and energy are calculated from the same set of discrete distribution functions.From the viewpoint of physical modeling,beyond the traditional macroscopic description,the DBM presents two sets of physical quantities so that the nonequilibrium behaviors can have a more complete and precise description.One set includes the dynamical comparisons of nonconserved kinetic moments of distribution function and those of its corresponding equilibrium distribution function.The other includes the viscous stress and heat flux.The former describes the speci fic nonequilibrium flow state,the latter describes the in fluence of current state on system evolution.The study on the former helps understanding the latter and the nonlinear constitutive relations.[35]The new observations brought by DBM have been used to understand the mechanisms for formation and effects of shock wave,phase transition,energy transformation,and entropy increase in various complex flows,[30,34,36]and to study the in fluence of compressibility on Rayleigh–Taylor instability.[31?32]In a recent study,it is interesting to find that the maximum value point of the thermodynamic nonequilibrium strength can be used to divide the two stages,spinnodal decomposition and domain growth,of liquid-vapor separation.

Some of the new observations brought by DBM,for example,the nonequilibrium fine structures of shock waves,have been con firmed and supplemented by the results of molecular dynamics.[37?39]It should be pointed out that the molecular dynamics simulations can also give microscopic view of points to the origin of the slip near boundary,such as the non-isotropic strong molecular evaporation flux from the liquid,[40]which might help to develop more physically reasonable mesoscopic models for slip- flow regime.

In order to extend DBM to the micro- fluid,it is critical to develop a physically reasonable kinetic boundary condition.Many efforts have been made to devise mesoscopic boundary condition for LBM to capture the slip phenomenon.[2,22?26]However,the previous works are most suitable for two-dimensional(2D)models with a very small number of particle velocities and can not be directly applicable to the DBM.On the other hand,those boundary conditions fail to capture flow characteristics in the Knudsen layer so the effective viscosity or effective relaxation time approach needs to be adopted.[15,19]Besides,the results of LBM and DBM should be veri fied by the results of continuous Botlzmann equation.In 2009,Watari[16]gave a general diffuse re flection boundary for his thermal LB model[41]and investigated the velocity slip and temperature jump in the slip- flow regime.Then,in his sequent work,[42]he compared the relationship between accuracy and number of particle velocities in velocity slip.Two types of 2D models,octagon family and D2Q9 model,are used.It was found that D2Q9 model fails to represent a relaxation process in the Knudsen layer and the accuracy of the octagon family is improved with the increase in the number of particle velocities.However,all the boundary conditions were set as diffuse re flection wall and the tangential momentum accommodation coefficient(TMAC)was not taken into account.

Because of the dependence of the mean free path on microscopic details of molecular interaction,especially the collision frequency,the Knudsen number may have different values in various interaction models for the same macroscopic properties.

In this paper,we first clarify the de finitions of Knudsen number and the connection between the hard sphere model and BGK model for three-dimensional(3D)condition so that the results obtained from various models can be compared dynamically.Then a general Maxwell-type boundary condition for DBM is represented and accommodation coefficient is introduced to implement a gas-surface interaction model.Two kinds of gas flows,Couette flow and Poiseuille flow,in a microchannel are simulated.In the section of Couette flow,the relation between the analysis solutions based on hard sphere and BGK model are veri fied.The simulation results with various Knudsen numbers and accommodation coefficients are compared with analytical ones based on linear Boltzmann equation in the literature not only on the velocity slip but on the Knudsen pro files.While in the section of Poiseuille flow,the simulation results are compared with analytical solution based on Navier–Stokes equation and the second order slip boundary condition.

2 Models and Methods

2.1 De finition of Knudsen Number

The Knudsen number is de fined as the ratio of the free path of molecules(λ)to the characteristic length(L),

Throughout the paper,we consider the characteristic lengthLas unit,so the Knudsen numberKnis equal to the value ofλ.

For the hard sphere collision model,the molecules are considered as hard spheres with diameterd,the mean free path of moleculesλHScan be calculated by

wherenis the number density of molecules.[4]According to Chapman and Enskog,[43]the viscosity coefficientμof hard sphere molecules can be expressed by

wherekis the Boltzmann constant,mis the molecular mass,andTis the temperature.It should be noted that gas constantRcan be expressed byR=k/m.Combining the state equation of ideal gas(p=ρRT),we have the following relationship betweenλHSand macroscopic quantities:

For the BGK model,the mean free path of moleculesλBGKis de fined as

whereτis the reciprocal of collision frequency and called relaxation time,ˉcis the average thermal speed.[4]The definition ofKnin DBM is in accordance with the de finition here.

According to the kinetic theory of gas molecules,is expressed by

in 3D case.From the Chapman–Enskog expansion,we know thatτhas the following relation with macroscopic quantities:

The comparison of Eq.(4)and Eq.(8)yields the relationship between viscosity-based mean free pathλHSandλBGK,

2.2 Discrete Boltzmann Model

The 3D discrete Boltzmann model taking into account the effect of the external force was presented based on the thermal model represented by Watari.[41]The evolution of the discrete distribution functionfkifor the velocity particlevkiis given as

To recover the NS equations,the local equilibrium distribution function should retain up to the fourth order terms of flow velocity.The discrete local equilibrium distributioncontaining the fourth rank tensor is written as

The velocity particlesvkiconsist of a rest particle and 32 moving particles.Each moving particle has four speeds and can be obtained from the unit vectors in Table 1 multiplied by differenceck.The speedsckof moving particle is determined according to the method presented by Watari in Ref.[16].

Table 1 Discrete velocity model,where λ=

Table 1 Discrete velocity model,where λ=

i direaction Unit vector(vix,viy,viz)i=1–8 λ(±1,±1,±1)i=9–12 λ(0,±φ?1,±φ)i=13–16 λ(±φ,0,±φ?1)i=17–20 λ(±φ?1,±φ,0)i=21–24 ?(0,±φ,±1)i=25–28 ?(±1,0,±φ)i=29–32 ?(±φ,±1,0)

TheFkin Eq.(11)is the weighting coefficient for the particle velocityvkiand is determined byckusing the following equations:

2.3 Boundary Condition Models

To solve the evolution equation(10),finite-difference method is adopted.The spatial derivative is solved by the second-order upwind scheme and time derivative is solved by the first-order forward scheme.Then the evolution equation(10)can be rewritten as

Fig.1 Schematic of the space grid.

Forvkiα≥0,the evolution equation(17)with Eq.(18)is applied fromI=3 up to the right wall.However,at the nodeI=2,the second-order upwind scheme in Eq.(18)is not applicable.The first-order upwind scheme,

is applied there.Forvkiα<0,the evolution equation(17)with Eq.(18)is applied from the left wall to the nodeI=N?2.Likewise,at the nodeI=N?1,the firstorder upwind scheme,

is applied.

To solve the value of the distribution function on the left wall forvkiα>0 and on the right wall forvkiα<0,boundary condition models are required.

對(duì)于每一位學(xué)生而言都有實(shí)現(xiàn)以上目標(biāo)的可能性。成績(jī)回報(bào)、與外校同學(xué)建立人脈，是可以自我控制的。參加技能競(jìng)賽方面，一個(gè)院校往往至少派一個(gè)隊(duì)(四名同學(xué))參加，因此只要個(gè)人努力，參加競(jìng)賽的可能性還是很大的。

(i)Diffuse re flection boundary

The complete diffuse re flection model assumes that the molecules leaving the surface with a local equilibrium Maxwellian distribution irrespective of the shape of the distribution of the incident velocity.It can be expressed as

The equilibrium distribution functions,and,are determined from the wall conditions including the velocities and the surface temperatures.Using the zero-mass flow normal to the wall,[16]the densityandcan be respectively calculated by the following two equations:

As a result,the distribution function on the left wall(fki,1)forvkiα>0 and on the right wall(fki,N)forvkiα<0 are solved under the diffuse re flection boundary condition.

(ii)Specular re flection boundary

The specular re flection model assumes that the incident molecules re flect on the wall surface as the elastic spheres re flect on the entirely elastic surface.The component of the relative velocity normal to the surfaces reverses its direction while the components parallel to the surface remain unchanged.As an example,the direction normal to the wall surface parallels to thexaxis,then the molecules leave the surface with a distribution as Since the distribution functionforvkix<0 andcan be solved by Eq.(17)with Eq.(18),the distribution function on the right wall(fki,N)forvkix<0 and on the left wall(fki,1)forvkix>0 are easy calculated from Eqs.(25)and(26).

(iii)Maxwell-type boundary

In practice,the real re flection of molecules on the body surfaces cannot be described properly by complete diffuse re flection or pure specular re flection.So the Maxwell-type re flection model which is composed of the two re flection modes is needed.The TMAC,αis introduced to measure the proportion of complete diffuse re flection.[4]Theαportion of the incident molecules re flect diffusely and the other(1?α)portion re flect specularly.The value of TMAC is used to characterize the degree to which the re flected molecules has adjusted to the tangential momentum of the surface,

whereτiandτrare the tangential components of the momentum fluxes of the incident and re flected molecules,respectively.τwis the tangential momentum fluxes of the molecules in the wall.α=1 corresponds to the case of complete tangential momentum accommodation and the molecules re flect with the Maxwellian distribution under wall condtion,uwandTw.α=0 corresponds to the case when the incident molecules are entirely not adjusted to the conditions of the surface,τr=τi.

Under this boundary condition,the distribution function on the right wall,(fki,N),forvkix<0 and on the left wall,(fki,1),forvkix>0 are solved by the following equations,respectively,

3 Simulation Restlts

3.1 Couette Flow

Fig.2 Schemtic of typical velocity pro file for Couette flow with slip effect(right half).

Consider a gas flow between two parallel walls,one atx=?Land the other atx=L.The two plates are kept at uniform temperatureT0and move with velocity(0,?v,0)and velocity(0,v,0),respectively.Velocity slip becomes more signi ficant with the decrease of the distance between the two plates or with the increase of the mean free path of the molecules,more exactly,with the increase of Knudsen number.the NS flow area.The flow near the wall possesses pronounced non-equilibrium characteristics,and the corresponding flow layer is known as the Knudsen layer whose thickness is of the order of the mean free path.In Fig.2,the linear portionA-B,whose gradient is dv/dx,corresponds to the NS flow area and the portionB-Dcorresponds to the Knudsen layer.The lineB-Cis extended from the lineA-Band the pointCis the cross point of the extended line with the right wall.The slip velocityvslipis de fined as the difference between the velocity of the right wall(the value the pointW)and the velocity value of the pointC.Considering the complete diffuse re flection,Sone gave the relation betweenvslipand the mean free pathλin Ref.[44].

For the hard sphere model,

The typical velocity pro file between parallel plates in the slip- flow regime is depicted in Fig.2.Only right half of the pro file is shown because of its antisymmetry.The gas flow away from the wall can be described by NS equations,and the corresponding flow area is referred to as

The relationship betweenλHSandλBGKdeduced in Subsec.2.1 is veri fied by Eq.(30)and Eq.(31)since 1.2540λHS≈1.0162λBGK.

Knudsen pro file?vis de fined as the difference between the curves,B-DandB-C.Sone[44]gave also the relation between ?vandλ,by introducing the so-called Knudsen layer function,Y0(η),whereηis a coordinate transformed fromxthrough the following conversion:

The Knudsen layer function for hard sphere model,and for BGK model,,are both shown in Fig.3.The correction of the functionaccording to the relation in Eq.(9)is also plotted.It can be seen that,the pro file of the corrected function is in excellent agreement with the pro file.As a consequence,the Eq.(9)is revalidated.In addition,the results based on hard sphere model can be compared with those from BGK model under same macro conditions by using the relation of Eq.(9).

Considering the Maxwell-type boundary condition,Onishi[45]gave the expression of slip velocity and Knudsen layer functionY(α)0(η)under various TMACs as

Aiandksis the coefficient calculated by re fined moment methods.[45?46]According to Onishi,[45]the solutions ofN=7 are good approximations with high and sufficient accuracy to the exact ones.The coefficients for partial values ofαare listed in Table 2.It should be noted thatis in complete agreement withY0(η)shown in Fig.3 whenα=1.

Fig.3 Pro files of the function Y0(η).A local enlargement of the curves is shown in the inset.The solid line is for the Hard Sphere(HS)model.The dashed line is for the BGK,and the dash-dotted line is for the Hard Sphere model with correction(HS-corrected).

Table 2 Coefficients for several kinds of values of α.

Fig.4 DBM simulation results for different Knudsen numbers under complete diffusion re flection boundary conditions.(a)Velocity pro files for different Knudsen numbers.(b)Slip velocities from DBM and those from Sone’s formulas.

The DBM simulation results for the Couette flow with different Knudsen numbers under complete diffusion boundary condition are shown in Fig.4.Results for different Knudsen numbers are obtained by changing the relaxation-time constantτaccording to Eq.(5).

Figure 4(a)shows that the phenomena of velocity slip become more signi ficant with the increase ofKn.Comparison of the values of slip velocity normalized by dv/dxbetween the DBM results and Sone’s results is shown in Fig.4(b).The two kinds of results are in excellent agreement with each other.The DBM accurately capture the velocity slip.Besides,the Knudsen layer is also well de-scribed by DBM.As shown in Fig.5,comparison of normalized Knudsen pro files calculated from DBM are also in excellent agreement with Sone’s results.The Knudsen pro files?vin Fig.5 are normalized by

TakingtheTMAC (α)into consideration,the Maxwell-type boundary condition is adopted in the following simulation.The DBM simulation results with several different values ofαare shown in Fig.6.From Fig.6(a),we can see that the phenomena of velocity slip are more signi ficant with the decrease ofα.The values of velocity slip for different values ofαare compared with those given by Eq.(34).Figure 6(b)shows good agreement of DBM simulation results with those of Onishi.[45]The results offor various values ofαare compared in Fig.7.The DBM results also show good agreement with those of Eq.(35).

Fig.5 Comparison of normalized Knudsen pro files from DBM simulations and those from analysis.

Fig.6 DBM simulation results for different values of α under Maxwell-type re flection boundary conditions.(a)uy(x)pro files for different values of α.The various dashed lines are for DBM simulation results and the solid lines are for linear fitting results.A local enlargement is shown in the inset.(b)Comparison of normalized slip velocities from DBM simulations and those from analysis under different values of α.

3.2 Poiseuille Flow

Pressure driven gas flow,known as Poiseuille flow,in a microchannel is also very common in MEMS.In the slipflow regime,the Navier–Stokes equations with slip boundary condition are applicable.Accurate second-order slip coefficients are of signi ficantly since they directly determine the accuracy of the results given by Navier–Stokes equations.[47]The first second-order slip model was presented by Cercignani.Using the BGK approximation he obtained

It is clear that the pressure gradient and the viscosity coefficient vanish in the nondimensionalizedU(y).

Fig.7 Comparison of normalized Knudsen pro files from DBM simulations with those from analysis for different values of α.A local enlargement is shown in the inset.

Fig.8 Comparison of velocity pro files between DBM results and analysis for different Knudsen numbers.Symbols represent the DBM results while lines represent the analytic ones.

The simulation results by DBM for different Knudsen numbers are shown in Fig.8.It can be found that the velocity slip is signi ficant with the increase of Knudsen number.The nondimensional velocity has a higher maximum value for a smaller Knudsen number.The velocity pro file described by Eq.(41)withα=1 is also plotted for comparison.The simulation results show good agreement with the expression of Eq.(41).

Fig.9 Comparison of velocity pro files between DBM results and analysis for different values of α.Symbols represent the DBM results while lines represent the analytic results.

Considering the behaviours of velocity slip under various TMACs.Then the Maxwell-type boundary condition is used.It can be seen from Fig.9 that the velocity slip is more signi ficant with the decrease of TMAC and the nondimensional velocity has the highest maximum value when the complete diffuse re flection occurs.It is concluded that the effect of Knudsen and TMAC on velocity is in the opposite direction.The numerical results are in well agreement with Eq.(41)for different values ofαwhich verify the accuracy of the Maxwell-type boundary condition.

4 Conclusion

A discrete Boltzmann method with Maxwell-type boundary condition for slip flow is presented.The de finition of Knudsen number is clari fied for DBM.The relation between the Knudsen number based on hard sphere model and that based on BGK model is given.Two kinds of gas flows,including Couette flow and Poiseuille flow,are simulated to verify and validate the new model.The results show that the DBM with Maxwell-type can reasonably capture both the velocity slip and the flow characteristics in Kundsen layer under various Knudsen numbers and tangential momentum accommodation coefficients.

Acknowledgement

The authors would like to thank Drs.Wei Jiang,Hongwei Zhu,Wei Wang,and Ge Zhang for fruitful discussions on discrete Boltzmann modeling of slip flows.

[1]M.Gad-El-Hak,Mec.Indust.2(2001)313.

[2]Y.Zhang,R.Qin,Y.Sun,et al.,J.Stat.Phys.121(2005)257.

[3]F.Bao and J.Lin,Inter.J.Nonlinear Sci.Numer.Simul.6(2005)295.

[4]C.Shen,Rare fied Gas Dynamics:Fundamentals,Simulations and Micro Flows,Springer,Berlin Heidelberg(2005).

[5]D.Burnett,Proceedings of the London Mathematical Society 40(1935)382.

[6]W.Chen,W.Zhao,Z.Jiang,and H.Liu,Phys.Gases 1(2016)5.

[7]H.Grad,Commun.Pure Appl.Math.2(1949)331.

[8]G.Bird,Molecular Gas Dynamics and the Direct Simulation of Gas Flows,Clarendon Press,Oxford(2003).

[9]J.Fan and C.Shen,J.Comput.Phys.167(2001)393.

[10]K.Xu and Z.Li,J.Fluid Mech.513(2004)87.

[11]K.Xu,Commun.Comput.Phys.314(2016)305.

[12]L.Yang,C.Shu,J.Wu,and Y.Wang,J.Comput.Phys.306(2016)291.

[13]Z.Guo,K.Xu,and R.Wang,Phys.Rev.E:Stat.Nonlinear and Soft Matter Phys.88(2013)033305.

[14]Z.Guo,R.Wang,and K.Xu,Phys.Rev.E:Stat.Nonlinear and Soft Matter Phys.91(2015)033313.

[15]Z.L.Guo,B.C.Shi,and C.G.Zheng,EPL 80(2007)24001.

[16]M.Watari,Phys.Rev.E 79(2009)066706.

[17]X.Shan,X.Yuan,and H.Chen,J.Fluid Mech.550(2006)413.

[18]V.Sofonea and R.F.Sekerka,J.Comput.Phys.207(2005)639.

[19]G.H.Tang,Y.H.Zhang,X.J.Gu,and D.R.Emerson,EPL 83(2008)226.

[20]Q.Li,Y.L.He,G.H.Tang,and W.Q.Tao,Micro fluid.Nano fluid.10(2011)607.

[21]Q.Li,K.Luo,Q.Kang,et al.,Prog.Energy Combust.Sci.52(2016)62.

[22]S.Ansumali and I.V.Karlin,Phys.Rev.E:Stat.Nonlinear and Soft Matter Phys.66(2002)026311.

[23]S.Succi,Phys.Rev.Lett.89(2002)064502.

[24]X.D.Niu,C.Shu,and Y.T.Chew,EPL 67(2004)600.

[25]G.H.Tang,W.Q.Tao,and Y.L.He,Phys.Fluids 17(2005)329.

[26]Z.Guo,B.Shi,T.S.Zhao,and C.Zheng,Phys.Rev.E:Stat.,Nonlinear,Soft Matter Phys.76(2007)056704.

[27]A.Xu,G.Zhang,and Y.Ying,Acta Phys.Sin.Ch.Ed.64(2015)184701.

[28]A.Xu,G.Zhang,and Y.Gan,Mech.Eng.38(2016)361.

[29]A.G.Xu,G.C.Zhang,Y.J.Ying,and C.Wang,Sci.China:Phys.,Mech.Astron.59(2016)650501.

[30]Y.Gan,A.Xu,G.Zhang,and S.Succi,Soft Matter 11(2015)5336.

[31]H.Lai,A.Xu,G.Zhang,et al.,Phys.Rev.E 94(2016)023106.

[32]F.Chen,A.Xu,and G.Zhang,Front.Phys.11(2016)114703.

[33]C.Lin,A.Xu,G.Zhang,and Y.Li,Combust.Flame 164(2016)137.

[34]Y.Zhang,A.Xu,G.Zhang,et al.,Combust.Flame 173(2016)483.

[35]A.Xu,G.Zhang,Y.Li,and H.Li,Prog.Phys.34(2014)136.

[36]C.Lin,A.Xu,G.Zhang,et al.,Phys.Rev.E:Stat.Nonlinear&Soft Matter Phys.89(2014)013307.

[37]H.Liu,W.Kang,Q.Zhang,et al.,Front.Phys.11(2016)1.

[38]H.Liu,Y.Zhang,K.Wei,et al.,Phys.Rev.E 95(2017)023201.

[39]H.Liu,W.Kang,H.Duan,et al.,Sci.China:Phys.,Mech.Astron.47(2017)070003.

[40]W.Kang,U.Landman,and A.Glezer,Appl.Phys.Lett.93(2008)123116.

[41]M.Watari and M.Tsutahara,Physica A 364(2006)129.

[42]M.Watari,J.Fluids Eng.132(2010)101401.

[43]S.Chapman and S.T.Cowling,The Mathematical Theory of Non-Uniform Gases:An Account of the Kinetic Theory of Viscosity,Thermal Conduction and Diffusion in Gases,Cambridge University Press,Cambridge(1953).

[44]Y.Sone,Molecular Gas Dynamics:Theory,Techniques,and Applications,Springer Science&Business Media(2007)91.

[45]Y.Onishi,Osaka Prefecture University Bulletin 22(1974).

[46]Y.Sone and Y.Onishi,J.Phys.Soc.Jpn.35(1973)1773.

[47]N.G.Hadjiconstantinou,Phys.Fluids 15(2003)2352.