Efficient Turbulent Compressible Convection in the Deep Stellar Atmosphere
Abstract
This paper reports an application of gaskinetic BGK scheme to the computation of turbulent compressible convection in the stellar interior. After incorporating the Subgrid Scale (SGS) turbulence model into the BGK scheme, we tested the effects of numerical parameters on the quantitative relationships among the thermodynamic variables, their fluctuations and correlations in a very deep, initially gravitystratified stellar atmosphere. Comparison indicates that the thermal properties and dynamic properties are dominated by different aspects of numerical models separately. An adjustable Deardorff constant in the SGS model and an amplitude of artificial viscosity in the gaskinetic BGK scheme are appropriate for current study. We also calculated the densityweighted auto and crosscorrelation functions in Xiong’s ([1977]) turbulent stellar convection theories based on which the gradient type of models of the nonlocal transport and the anisotropy of the turbulence are preliminarily studied. No universal relations or constant parameters were found for these models.
Vol.0 (200x) No.0, 000–000
1 Introduction
Turbulent convection has tight relations to the unsolved problems in the theory of stellar structure and evolution, especially, for the massive stars (Deng et al. [1996a]; [1996b] and references therein). These problems cannot be settled with pure analytical ways. Along with the development of computing science, numerical simulations become a powerful tool to investigate the hydrodynamic properties of astrophysical flows. It is widely used in the study of formation of cluster, accretion disk and evolution of galaxies. The convection in stellar interior have also been studied by many authors with numerical experiments. Due to the difficulties in this problem, the progress made in this field is limited. However, it is generally believed that numerical testing of some analytical models and local highresolution simulations can make a lot of sense to our understanding of stellar convection. So far, the numerical hydrodynamic scheme applied to the stellar convection are LaxWendroff scheme (Graham [1975], etc.), alternating direction implicit method on staggered mesh (ADISM) (Chan et al. [1982]; [1986], etc.), pseudospectrum scheme (Hossain & Mullan [1990]; [1991], etc.), piecewiseparabolic method (PPM) (Porter & Woodward [1994]; [2000], etc.), upwind scheme and so on. At the present time, the most suitable numerical scheme for the turbulent flow is the spectrum method, but it cannot handle the discontinuity in the motion of fluids.
Gaskinetic BGK scheme is a recently matured method for computational fluid dynamics and attached much attention in the practical problems (Xu [2001]). It is accurate and robust for computing supersonic unsteady flows. However its application to astrophysical flows is not popular attributed to its complex. Theoretical analysis shows that near the surface of stellar envelope, the motion of fluids becomes supersonic which cannot be selfconsistently treated by traditional mixinglength theory (MLT) (Deng & Xiong [2001]). Our original aim is to simulate the supersonic turbulent convection in the outer region of yellow giants by gaskinetic BGK scheme which would involve many efforts in different directions. We have already extend the BGK scheme to include the gravitational acceleration (Tian et al. [2007]). Before using it to compute the supersonic stellar convection, the turbulence model, radiation transfer model and realistic input physics must be correctly implemented.
Restricted by the capacity of digital computer, we cannot afford the very high resolution numerical experiments for the stellar type of turbulent convection. Large eddy simulation (LES) which calculates the large eddies explicitly while mimics the subgrid eddies by models may be the most feasible way in current stage. In current paper we implement the SGS turbulence model (Smagorinsky [1963]; Deardorff [1971]) into the BGK scheme and validate the threedimensional BGK code by calculating the turbulent compressible convection in a deep stellar atmosphere. For very high Reynolds number, the behaviors of the turbulent flows are greatly affected by the numerical and physical dissipation in the scheme. A investigation of these effects is very necessary before the code is applied to the practice. By varying the Deardorff number in the SGS model and artificial viscosity parameter introduced to capture the shock, we constructed three models which are similar to those studied by Chan & Sofia ([1989], hereafter CS89; [1996], hereafter CS96). The empirical relations derived by them were reexamined. A study of densityweighted auto and crosscorrelation function, anisotropy of turbulence and diffusive type of models of nonlocal turbulent transports in the turbulent stellar convection theory of Xiong ([1977]) was conducted too.
In the next section, we give a description of gaskinetic BGK scheme and mainly focus on the incorporation of SGS model. The computed physical models are formulated in Sect. 3. The numerical results are shown in Sect. 4 where the discussions are also presented. The conclusions are summarized in the last section.
2 Gaskinetic BGK Scheme
General numerical method for hydrodynamic problems is to directly discretize the NavierStokes equations,
(1)  
(2)  
(3) 
where is the density, is the velocity, is the pressure, is the gravitational acceleration and is the summation of internal energy and kinetic energy.
is the viscous stress tensor, where is the strain rate tensor, is the identity tensor, and are the dynamical and bulk viscosity coefficient, respectively. is the diffusive type of energy flux. Differently, the BGK scheme works on the BGK equation,
(4) 
which is an approximation of Boltzmann equation (Bhatnagar et al.[1954]). In above expression, is the gas distribution function in the phase space, is the particle velocity, is the collision time, and . The righthand side of equation (4) is the socalled relaxation model, which is a simplification of the complicated collision term in the Botlzmann equation (Vincent et al. [1965]). It physically means that the initially nonequilibrium distribution will approach the equilibrium state after the particles collide once. A larger corresponds to a further state from , i.e., the stronger nonequilibrium transport effects, such as viscosity and conduction. In our study, the equilibrium state, , in equation (4) is taken to be Maxwellian distribution. It can be proved mathematically that the solutions of NavierStokes equations (1–3) are automatically obtained through solving the BGK equation with the following definitions of dissipative coefficients:
(5) 
where is the thermal conductivity, is the internal degree of freedom for particles, is the molecule mass and is the Boltzmann constant.
The basic idea of BGK method is to find a local approximation to the nonlinear equation (4). Then evaluate the macro quantities (e.g., fluxes) using the micro distribution function . Finally, the cell average values are updated according to the conservation laws (finite volume method). In our first attempt, the BGK method (Xu [2001]) has been extended to include the external force. A detailed description of the three dimensional multidimensional gaskinetic BGK scheme for the NavierStokes equations under gravitational fields was given by Tian et al. ([2007]). Here, we only outline the new implements, i.e., the incorporation of turbulence model.
In the study of convection, the combined effects of viscosity, heat conduction and temperature difference on the instability of flows are measured by Rayleigh number: , where is the thermal expansion coefficient and is the kinematic viscosity. For the polytropic gas defined in Sect. 3, the Rayleigh number can be written as:
(6) 
where is the gas constant, is the depth of computational domain and is the ratio of specific heat, the meaning of other symbols can be found in Sect. 3. In the gaskinetic BGK scheme, the viscosity is controlled by collisions of particles. We can relate the Rayleigh number to the collision time by the following way. Suppose during each timestep , the particles collide times. Then we have
(7) 
where is Courant number, is the spatial resolution, is speed of fluids and is the sound speed. From equation (5), (6) and (7), we get
(8) 
where is the Mach number, is the vertical grids size and is Prandtl number. For current study, we have , , , , and . In efficient turbulent convection, the energy transfer by heat conduction is negligible which means the Prandtl number is very large. In all of our simulations, the is set to be . Therefore, approximately we have . Similarly, we have Reynolds number . While in the typical stellar convection, . Hence, is needed when we do a direct numerical simulation (DNS) of stellar convective flows where . These values can be reduced by increasing grids number which is very expensive for the present generation hardware. At the same time, a very small would introduce large computational error. An alternate way is the LES which simulate the large eddies directly and approximate the small eddies with models. There are a lot of approaches to perform the LES, the simplest way may be the SGS model (Smagorinsky [1963]; Deardorff [1971]).
In the BGK scheme, the viscosity is introduced through the collisions of particles. The natural way of implementing SGS model is to modify the collision time. Chen et al. ([2003]) included the renormalization group large eddy model into the BGK equation, where is the turbulent kinetic energy and is the turbulent dissipation. In the model, tow additional equations are needed to be solved. For sake of simplicity, we consider the Smagorinsky ([1963]) model. In current study, the collision time is defined as
(9) 
where the first term of righthand side represents the molecule viscosity and the second term is introduced to increase the numerical dissipation when there is a jump in pressure around the control volume boundary. is an adjustable constant, and are the reconstructed pressures at the left and right side of a cell interface (see sketch in Fig. 1). In the strong supersonic region, additional dissipation caused by this term is essential to stabilize the computing.
The last term in righthand side of equation (9) is implemented to account for the SGS viscosity,
(10) 
where is an adjustable constant, usually has the value: 0.1 0.2, the filter width is taken to be the local resolution, colon stands for the contract of tensor and . Above model is called Smagorinsky model or sometimes SmagoringskyLilly model. In our calculations, the eddy viscosity is computed in the control volume by staggered mesh strategy and then interpolated to the cell interface. In the BGK scheme, there is intrinsic diffusion caused by particle collisions. In current study, we are only interested in the turbulent properties. So, the molecule Prandtl number is set very large. And the following diffusive flux (CS96) is implemented explicitly into the BGK scheme,
(11) 
where is the specific entropy, is the specific heat at constant pressure and is the adiabatic gradient. In the stable layer, is set to make the diffusion carry out the input energy flux. In the convection zone, is very close to zero. represents the turbulent diffusion and is set to zero in the stale region, where is the effective Prandtl number of SGS turbulence and taken to be .
3 Physical Models
Our physical problem is very similar to those studied by CS89 and CS96. An ideal gas () in a rectangular box is considered with gravity in the vertical direction. The side boundaries are periodic and the top and bottom boundaries are impenetrable and stress free. In order to avoid boundary effects, a very thin stable layer is placed below the upper boundary and the diffusive flux is gradually enhanced near the lower boundary to make it carry out total flux at the lower boundary. A constant energy flux is fed at the bottom. At the top, the entropy is fixed. The system is initially static:
(12)  
(13)  
(14)  
(15) 
where is the normalized parameter, and is the polytropic gas index. The gravitational acceleration comes from the hydrostatic equilibrium . The subscripts and denote top and bottom values respectively. Above solutions to NavierStokes equations are not stable against small perturbations. In all of our calculations, the velocity field is slightly perturbed initially. After longtime thermodynamic relaxation, the system will reach a statistical steady state. We defined a series of runs to test the effects of numerical parameters. The numerical effects of a variety of parameters were studied by Chan & Sofia ([1986], hereafter CS86) and CS89 in details. The effects of changing turbulent Prandtl number was tested by Singh & Chan([1993]). Here, we just focus on the important parameter in the SGS model and new parameter appearing in the BGK scheme. The details are given in the second line of Table 1.
All the cases we computed use a mesh. The vertical grid decreases smoothly with height (about grids per PSH) and the horizontal grid is uniform. The aspect ratio (width/depth) of the box is .
4 Results
In this part, we show the results from the numerical simulations. All the runs were evolved numerical timesteps, corresponding to a dimensionless time around , before the statistical analysis is performed. The statistical steady state is indicated by the balance of the input energy flux from the bottom and the outgoing energy flux through the top. In our calculations, the spatial variation of averaged total energy flux from is within (see solid line in Fig. 2). Another criterion is the averaged vertical mass flux which is less then everywhere in all cases. Hence, the system will not undergo substantial adjustment any more. The statistical average covers numerical timesteps.
Except for the instant velocity fields, all the other quantities investigated here are the mean values. For an arbitrary quantity , represents its combined horizontal and temporal mean, denotes the deviation from , stands for the root mean square (rms) fluctuation from the .
In some figures, the integral pressure scale height (PSH) is shown by the vertical lines. For example, in Fig. 2, the vertical solid line denotes the location of stableunstable interface near the upper boundary. The second dashed line at the left side of the solid line is 2 PSHs away from the upper stableunstable interface. Although the numerical parameters are different for case A to case C, their relaxed thermal structure are nearly equivalent and the discrepancy between the integral PSH locations is very small. So, we can plot the results from all runs in the same figures.
4.1 Velocity Fields
In CS86, the threedimensional turbulent flow structure was well depicted by the pseudo stream lines. For sake of clarity, Fig. 3 shows the velocity fields projected in xz plane and xy plane. From the left panel of Fig. 3, it is evident to see that the high speed motions exist in the top region and are associated with the downward streams. In our calculations, the maximum Mach number occasionally exceeds one. In CS96, the SGS viscosity was enhanced by a factor to suppress the shocks occurring in the top region which would easily trigger the instability of the numerical computation, especially, during the early thermal relaxation. Based on fluxsplitting (see Fig. 1), such implement can be avoided in gaskinetic BGK scheme. If the nonlinear van Leer limiter is replaced by central interpolation, the supersonic motion cannot be handled correctly. Numerical tests also show that a less than will smooth the turbulence, the circulations become laminar. The solutions are almost unchanged when great than . In current study, we adopt . The networks of downward streams can be seen the right panel of Fig. 3.
4.2 Approximate Relations
Relation Parameter  A  B  C  CS89  Deviation  ID  Category 
0.2  0.25  0.2  0.2  (%)  
1  1  0      
0.69 ( 0.12 )  0.66 ( 0.11 )  0.62 ( 0.11 )  0.61 ( 0.05 )  8  R0  III  
0.57 ( 0.10 )  0.58 ( 0.10 )  0.64 ( 0.11 )  0.61 ( 0.05 )  2  R1  III  
0.83 ( 0.01 )  0.83 ( 0.02 )  0.83 ( 0.01 )  0.89 ( 0.04 )  7  R2  I  
0.57 ( 0.04 )  0.56 ( 0.04 )  0.57 ( 0.04 )  0.57 ( 0.07 )  1  R3  I  
0.90 ( 0.01 )  0.90 ( 0.01 )  0.90 ( 0.01 )  0.94 ( 0.03 )  4  R4  I  
0.48 ( 0.01 )  0.47 ( 0.01 )  0.48 ( 0.00 )  0.26 ( 0.01 )  83  R5  II  
0.88 ( 0.14 )  0.85 ( 0.13 )  0.87 ( 0.14 )  0.51 ( 0.03 )  70  R6  IV  
1.52 ( 0.15 )  1.51 ( 0.13 )  1.51 ( 0.15 )  0.90 ( 0.10 )  68  R7  IV  
0.98 ( 0.00 )  0.98 ( 0.02 )  0.98 ( 0.03 )  0.99 ( 0.01 )  1  R8  I  
0.93 ( 0.01 )  0.92 ( 0.01 )  0.92 ( 0.01 )  0.89 ( 0.06 )  4  R9  I  
0.83 ( 0.02 )  0.83 ( 0.02 )  0.82 ( 0.03 )  0.82 ( 0.05 )  1  R10  I  
0.59 ( 0.04 )  0.59 ( 0.04 )  0.59 ( 0.03 )  0.49 ( 0.05 )  20  R11  II  
0.78 ( 0.02 )  0.79 ( 0.02 )  0.77 ( 0.02 )  0.81 ( 0.03 )  4  R12  I  
0.76 ( 0.02 )  0.78 ( 0.03 )  0.75 ( 0.02 )  0.81 ( 0.03 )  6  R13  I  
0.64 ( 0.01 )  0.67 ( 0.01 )  0.63 ( 0.01 )  0.74 ( 0.03 )  13  R14  II  
1.00 ( 0.02 )  0.99 ( 0.02 )  1.00 ( 0.02 )  1.00 (    )  0  R15  I  
1.49 ( 0.06 )  1.44 ( 0.05 )  1.50 ( 0.06 )  1.24 ( 0.08 )  19  R16  II  
1.46 ( 0.06 )  1.41 ( 0.06 )  1.46 ( 0.06 )  1.26 ( 0.08 )  15  R17  II  
1.28 ( 0.05 )  1.24 ( 0.04 )  1.28 ( 0.05 )  1.20 ( 0.08 )  6  R18  I  
0.81 ( 0.08 )  0.85 ( 0.07 )  0.80 ( 0.09 )  0.58 ( 0.07 )  41  R19  IV  
1.49 ( 0.06 )  1.44 ( 0.05 )  1.50 ( 0.06 )  1.25 ( 0.08 )  18  R20  II  
1.20 ( 0.15 )  1.23 ( 0.13 )  1.20 ( 0.15 )  0.72 ( 0.06 )  68  R21  IV  
83  R22  II  
( )  ( )  ( )  ( )  
13  R23  II  
( )  ( )  ( )  ( )  
21  R24  II  
( )  ( )  ( )  ( ) 

New relation defined in current paper.

.

These are standard deviations of the leastsquares fits.
Our BGK code has been extensively tested for laminar flows (Tian et al. [2007]). In order to validate the incorporation of SGS model, we reestimate quantitatively some approximate relations among the thermodynamic variables and their fluctuations. The results are given in Table 1 where the results from CS89 are also listed for comparison. Our computational models are partially similar to CS89 and partially to CS96. Models of CS89 undergo substantially adjustment near the boundaries and we cannot afford the highresolution desired in CS96. During the data analysis, the correlation function of quantity and is defined as . The standard deviations of these approximations () are given in the brackets. The deviation of from CS89 is calculated by . We only concentrate on the middle region in the convection zone, namely, 1 PSH from the bottom and 2 PSHs from the upper stableunstable interface. The investigated layer expands about 3 PSHs.
The approximate relations and correlations in Table 1 are classified roughly into four categories which are indicated by Roman numerals in the last column according to the goodness of fit and discrepancy from CS89. These categories are: Category I: and %; Category II: and %; Category III: and %; Category IV: and %.
Category I contains the best relations where the thermal variables are mainly involved, i.e. , , and their fluctuations. These relations should be weekly affected by numerical scheme and details of the computational models. In the work of Kim et al. ([1995]) where a realistic equation of state (EOS) was used, these relations are different from current studies. Therefore, Category I may be dominantly determined by EOS.
Most relations of Category II are functions of pressure , velocity and their fluctuations. The amplitude of the fluctuations of the components of velocity in our computation is lower than those in CS89 (Fig. 1 therein) and CS96 (Fig. 2 therein) while the amplitude of the relative fluctuations of thermal quantities are around two times larger than those in CS89 (see Fig. 2 therein) and in Singh’s work ([1993]). Hence when the fluctuations of thermal variables are expressed in terms of , , etc., the coefficients are nearly doubled, e.g., and . This can also be seen in Fig. 3a in CS96.
Ratios of and do not deviate from CS89 very much, but the fitting is not accurate. This may be caused by the effects of the boundary and transition layers. It is obvious from Fig. 4 that the height distribution of is not flat as CS89. There is a small hump below the upper transition layer (one PHS from the unstablestable interface) which can also be found in CS96. We believe that the unstablestable transition is responsible for this. The large difference in Category IV should be caused by the high order powers of velocity fluctuations.
are commonly approximated in MLT. In current study, these relations are fitted by linear approximation very well. However, the slopes are different from the CS86 and CS96. The situation of can be explained as above. For , the reasons may lie in the fact that in the upper convective region, the amplitude of from CS96 is larger than ours and therefor we need a smaller . CS96 got a smaller (0.78) than CS89. In current study, it is even smaller which implies a smaller superadiabatic gradient .
The cause of most of the discrepancies between current study and CS89 may be the larger and smaller . In the traditional theory of stellar convection, the enthalpy flux is proportional to the fluctuations of temperature and velocity, i.e., and the kinetic flux is totally neglected. When the kinetic flux is comparable to the total flux, this kind of proportional relation is not exactly held any more. Based on the following facts:

, , (Singh [1993], Fig. 1 and Fig. 5);

, , (CS89, Fig. 1 and Fig. 2);

, , (CS96, Fig. 2);

, , (current, Fig. 4);

, , (Kim [1995], Fig. 7 and Fig. 8),
we can conclude that the and are proportional to in a nonlinear way. Note that current study use a different numerical scheme and Kim adopted a realistic EOS. It is beyond the scope of current study to perform a quantitative analysis of such relations.
Above comparison shows that the dynamical properties and thermal properties for stellar type of convection may be affected by different aspects of the numerical models separately. Thermal structure is mainly determined by physical parameters while dynamic motions can easily affected by numerical parameters.
4.3 Anisotropic Turbulence
Xiong’s nonlocal timedependent stellar convection theory is based on Reynolds stress method. It is dynamic theory of auto and crosscorrelation functions of turbulent velocity and temperature fluctuations. These fluctuations are defined as the derivation from the density weighted average, i.e.,
(16) 
The start point of Xiong’s theory is a set of partial differential equations of
(17) 
where , and the summation convention is adopted. The numerical results of , and are given in Fig. 5. In the closure models of Xiong’s theory, three adjustable parameters, i.e., , and are introduced. is used to describe the anisotropic turbulent motions. In Deng’s ([2006]) work, was related to turbulent velocity by in the fully unstable zone. In the upper overshooting region, they proposed that and is independent of . In the lower overshooting zone, and decrease as decrease. There is no enough room for overshooting in our models. The ratio of from numerical simulations is given in Fig. 6 from which we can see this ratio is slightly dependent on and . In the upper efficientinefficient convection interface (about 1 PSH from unstablestable interface), this value approximately equals to . It takes its maximum at the location where the turbulent convection starts to become inefficient near the bottom. Its maximum is affected evidently by and . Current study cannot give a definite solution for anisotropic turbulence. Here, we present preliminary suggestions. Suppose that is held in current models, we can conclude that is infinitely large at the upper efficientinefficient interface and decreases as the distance from the boundaries increases. The minimum of is about .
4.4 Nonlocal Transport Models
Generally, a Reynolds stress method suffers the socalled closure problem. CS96 numerically studied the popular closures and found they were poor. Some third order moments representing the nonlocal transport effects were approximated by gradient models in Xiong’s work ([1989]), i.e.,
(18)  
(19)  
(20) 
with , where is the Lagrangian integral length scale of turbulence. The nonlocal transports and coefficients from numerical simulations for all cases are shown in Fig. 7 . From panel (a) of Fig. 7, we can see that the is about one order larger than and is three orders less than . Hence, the nonlocal transports are dominated by turbulent kinetic energy. Panel (b)(d) of Fig. 7 show that the , and gradually increase with the distance from the top boundary and change this trend near the interface where the turbulent convection become inefficient. The variation of is slow in the lower half convection zone around the value of 1.8. But it is affected by and obviously. is about ten times larger than . Both of and vary rapidly which suggests that there may not exist the universal constant for these closure models.
4.5 Effects of Numerical Parameters: and
In current LES, the local effective grid Reynolds number is enlarged by SGS model whose amplitude is controlled by Deardorff constant . An inadequately small would cause the buildingup of the kinetic energy at the twogrid level and make the computation crashed. If is too large, the turbulent motions will be damped down. The proper value of is resolution and methoddependent since the numerical dissipation also plays an important role in the behavior of turbulence. Deardorff ([1971]) suggested a value of for which was used in the CS89. In CS96, this value was increased to .
In the fluxsplitting method, the discontinuity is introduced at the cell interface (see Fig. 1) by limiters. Additional dissipation (the second term in the left hand side of equation (9)) is employed to handle the strong shock waves near the sharp jumps of pressure. The typical value of is 1. In the smooth region this kind of discontinuity should be very small. However, our simulations have lower resolution and the computing zone extends about seven PSHs. So it is necessary to check if these values are adequate for studying the stellar type of convection.
Relations and in Table 1 show that the effects of changing these parameters on the thermal structure are really slight. They both affect the eddy properties with very small amplitude which can be seen from Fig. 8 where the autocorrelations of the vertical velocity are shown. Their profiles are nearly symmetrical except in the upper stable and transition zone. The half width at half maximum (HWHM) of these profiles should be sensitive on the viscosity (CS86). and in Table 1 are the and , respectively. They should be nearly equal to each other for isotropic turbulent flows. In our study, for case A (, ), the discrepancy between them is clear which becomes slight for case B (, ) and further slight for case C (, ). This kind of anisotropy comes from the initially perturbation and can also bee found in Fig.7 and Fig.11 in CS86. So it seems that is more suitable for current study. In the deep stellar convection, the shock waves are mild. Hence, the fluxsplitting may enough to handle it as in our tests. We suggest that it is better to make as small as possible because it would enhance the anisotropy which can be diminished by larger .
5 Conclusions
In this paper, we present a preliminary application of gaskinetic BGK scheme to the simulation of turbulent convection in stellar atmosphere. The approximate relations among thermodynamic variables, their fluctuations and correlations were examined. The anisotropy and diffusive models of nonlocal transport were investigated too. The effects of varying numerical parameters were also tested. The main conclusions are summarized as follows:

The behavior of the thermal variables and dynamic variables are affected by different aspects of the models and numerical scheme. For example, the fluctuations of density and pressure are dominantly determined by physical models while the fluctuations of velocity are sensitively dependent on numerical parameters, e.g., and .

There is no constant ratio of for anisotropic turbulence in current models. We suggest that take an infinite value at the boundary and approach its minimum () in the deep convective region.

The diffusive models for nonlocal transport are not applicable since the coefficients for different quantities are dramatically different. The best situation is the turbulent transport of turbulent kinetic energy where a roughly flat region exists for .

For current resolutions, is better than . And should be set as small as possible in any cases. A fluxsplitting technique is needed to stabilized the shock waves near the top. A Rayleigh number less than may smear the turbulent motions.
Our simulations may suffer from the lower resolutions, aspect ratio and the number of testing cases. But it is enough for the purpose of validation and get some preliminary results. A further study will be performed in the future.
Acknowledgements.
We wish to acknowledge the contributions of K. Xu to the developing of the hydrodynamic code. We also thank the Department of Astronomy at Peking University for providing computer time on their SGI Altix 330 system. This work was partially funded by the Chinese National Natural Science Foundation (CNNSF) through 10573022, 10773029 and by the national 973 program through 2007CB815406. KLC thanks Hong Kong RGC for support.References
 Bhatnagar P.L., Gross E.P., Krook M., 1954, Phys.Rev., 94, 511
 Chan K.L., Wolff C.L., 1982, J.Comput. Phys., 47, 109
 Chan K.L., Sofia S., 1986, \apj, 307, 222 (CS86)
 Chan K.L., Sofia S., 1989, \apj, 336, 1022 (CS89)
 Chan K.L., Sofia S., 1996, \apj, 466, 372 (CS96)
 Chen H., Kandasamy S., Orszag, S., Shock, R., Succi,S., Yakhot, V., 2003, Science, 301, 633
 Deardorff J.W., 1971, J.Comput. Phys., 7, 120
 Deng L., Bressan A., Chiosi C., 1996a, \aap, 313, 145
 Deng L., Bressan A., Chiosi C., 1996b, \aap, 313, 159
 Deng L., Xiong D.R., 2001, ChJAA, 1, 50
 Deng L., Xiong D.R., Chan K.L., 2006, \apj, 643, 426
 Graham E., 1975, J. Fluid. Mech., 70, 689
 Hossain M., Mullan D.J., 1990, \apj, 354, L33
 Hossain M., Mullan D.J., 1991, \apj, 380, 631
 Kim Y.C., Fox P.A., Sofia S., Demarque P., 1995, \apj, 442, 422
 Porter D.H., Woodward P.R., 1994, \apjs, 93,309
 Porter D.H., Woodward P.R., 2000, \apjs, 127,159
 Singh H. P., Chan K. L., 1993, \aap, 279, 107
 Smagorinsky J.S., 1963, Mon. Weather Rev., 91, 99
 Tian C.L., Xu K., Chan K.L., Deng,L.C., 2007, J.Comput. Phys, 226, 2003
 Xiong D.R., 1977, Acta Astron. Sinica, 18, 86
 Xiong D.R., 1985, \aap, 150, 133
 Xiong D.R., 1989, \aap, 209, 126
 Xu K., 2001, J. Comput. Phys., 171, 289
 Vincent, W.G. and Kruger Jr., C.H., 1965, Introduction to Physical Gas Dynamics, New York, Wiley