Numerical Solution of the Taylor-Quette Flow Problem: A Commodious Statistical Approach

Problem statement: The main problem in solving of the Boltzmann equat ion by statistical methods is its long computational time (CPU time). The problem of decreasing of CPU time usage in statistical solution methods of Boltzmann equation for rarefied vortex flows was studied. Approach: In this study the Boltzmann equation in a rarefied Tay lor-Quette flow regime was solved using the new Monte Carlo method which is officially named time R laxated Mont Carlo method that applied the equilibrium conditions in each time step. Results: The results obtained from time Relaxated Mont Carlo method for the problem at hand were compared with those from the usual Direct Simulation Monte Carlo method. This comparison showed good agr eement between the two sets of results. Conclusion: Whereas, the number of collisions and CPU usage ti m in the suggested method, as compared to the direct simulation Monte Carlo metho d, showed a significant decrease.


INTRODUCTION
The Taylor-Quette Flow is a flow between two coaxial cylinders, in which the outer cylinder is fixed and the inner one rotates. Taylor (1923) while studying this flow noticed that at certain angular speeds of the inner cylinder, vortex flows are created. Karman (1934) analyzed the Taylor-Quette flow in dense flow regimes and Kuhltau studied the Taylor-Quette flow in rarefied gas flow regimes experimentally (Kuhlthau, 1960). Since the Navier-Stockes equation can not describe rarefied gas flow regimes properly, it is essential to solve the Boltzmann equation which is valid for any flow regime (rarefied and dense). Within the scope that was developed by Kuhltau, Reichelmann and Nanbu numerically solved the Boltzmann equation by the Direct Simulation Monte Carlo method, known as the DSMC method and found that the Results from the DSMC were in good conformity with the experimental ones (Reichelmann and Nanbu, 1993). Stefanov and Cercignani (1993) used the DSMC method to study specially vortex flows between two coaxial cylinders. Following that, Bird, by considering more details, used the DSMC method and could obtain, in a more exact fashion, the stable solution for the Taylor-Quette flow (Bird, 1994).
The main problem in all previous solution methods of Taylor-Quette flow (Reichelmann and Nanbu, 1993;Stefanov and Cercignani, 1993;Bird, 1994) is their too long CPU time usage. So using of other appropriate method which leads to adequate result in less time for problem at hand is needed.
On the other hand, Carlen et al. (2000) solved the Boltzmann equation by using a method which was compatible with the application of the Maxwellian equilibrium conditions and Pareschi and Russo (2001) introduced the Time Relaxed Monte Carlo method that was proportional to the method introduced by Carlen et al. (2000). Pareschi's method is abbreviated as TRMC, in the references. Pareschi analyzed an outer flow problem with his new method. Since, Maxwellian equilibrium conditions applied in TRMC method decreases the calculation time, TRMC method can lead to appropriate results for the problem at hand in smaller time. In (Zeytoonian et al., 2009), for the first time, the Taylor-Quette flow between two coaxial cylinders was solved by the TRMC method and a comparison between some results with those obtained by the Bird's DSMC method done.
In this study, the authors have completed their before study and made a complete comparison between two methods. The complete comparison of the results from TRMC method with those obtained by the DSMC method, shows that obtained results are in good agreement with those of the DSMC and that the calculation time required to get stable results for the TRMC method, as compared to the same for the DSMC method, is significantly less.

MATERIALS AND METHODS
The appropriate mathematical model of the analysis of rarefied gas dynamics is the Boltzmann equation: In relation (1), f is the density distribution function (a non-negative function) and depends on the position x, the velocity v and the time t. In Eq. 1, ε is the Knudson number which is proportional to the mean free path between the collisions. Dual collision operator, which appears as Q (f, f) in Eq. 1 defines the collision between two molecules in a monoatomic gas as: In the relation above, ω is a vector on the unit sphere 2 3 S ⊂ ℜ which shows the direction of the molecule's collision. σ Kernel is a non-negative function which specifies the details of interaction between two molecules. Molecule's velocities after the collision ( ) The relative velocity of two molecules with respect to each other is defined as If f is the density in the phase space, then; the macroscopic density is defined as: Similarly, the macroscopic velocity is defined as: and the temperature and energy of the fluid can be obtained from: In the equilibrium state, where ( ) Q f , f 0 = , every local distribution function has a Maxwellian distribution form as: In relation (8), T , u ,ρ are the mean temperature, velocity and density of the gas, respectively. To solve the Boltzmann equation by Mont Carlo method, transient step (considering Q 0 = ) and collision step (considering x f 0 ∇ = ) will be solved separately. To begin the numerical solution of each steps simplest discritization, the first order time discretization of the Boltzmann equation, is discussed. After discretization, the main problem seems to appear in the way of solving the collision step. Therefore, it is discussed this issue in more details. So with considering only collision step: By using the simple Forward Euler method for discretization of relation (9) (which is used mostly for discretization) (Bird, 1994), we will have: It can be seen that the essential condition for the numeric scheme stability is that t ∆ should be in the same order of amplitude as the mean free path. So in small value of ε, this discretization method is not a proper one. In other hands, in the implicit discretization method, stability condition becomes less important, but other problems associated with solving the implicit collision operator (that is, insufficient data for the initial conditions required for the solving), makes this method, not a proper one, also. So, an appropriate solution for a nonlinear equation like (9), which, in addition to having unconditional stability, can still remain valid near the dense range, can be obtained as follows. Considering a variable change such as ( ) ( ) Eq. 9 becomes (Carlen et al., 2000): Initial condition in (11) is similar to (1), 0 µ ≠ is a constant parameter and P is a linear two-directional operator. By changing the variables and defining new variables (Carlen et al., 2000), we get: In which τ is the Relaxed Time. Also Eq. 11 changes to the simple form as follow: Considering Eq. 13 as a Qushi equation, it can be solved by the following power series: In which, k f functions can be calculated from the relation (15): By substituting the main variables (the ones before the variable change), we get the following equation which is the solution of the Qushi Eq. 9 as: The (16) where,

( )
M v is the local Maxwellian distribution function. Considering this case, the following solution, which is based on the Maxwellian equilibrium for m 1 ≥ , can be obtained (Pareschi and Russo, 2001): In (18), ( ) n f f n t = ∆ and t ∆ is a small time step.
Considering 4 µ = π σ ρ , the following result will be derived (Pareschi and Russo, 2000) (to have a non negative result, such assumption is essential): Here, σ is the maximum collision section in the flow field and is calculated according to (20). In relation (20), σ is the collision section between two specified molecules: Considering the result from (18), It should be noted that, this method can be derived through different weight functions which make the best estimate for high order coefficient in Wild extension (Pareschi and Russo, 2000); so, in the overall form: f k is calculated from (15) One function which satisfies the above conditions and also is based on relation (18) It should be noted that another weight function could be selected, but it must be checked to make sure that it is an optimal choice (Pareschi and Russo, 2000). Now with considering the TRMC and DSMC solution methods (both of them are based on the Monte Carlo method), for solving the previous equations, At first, it is presented a Direct Simulation of the Monte Carlo method which corresponds to Eq. 10 and is compatible with the Nanbu-Babovsky algorithm. It should be pointed out that the molecule here means a model molecule which, by itself, represents many more molecules. Considering f as a probability distribution of the particles, we have: Having used the variable change, the Boltzmann's new form was the relation (11). With discritization and using the Forward Euler method and using ( ) It should be noted that if Thus, the probabilistic Interpretation of Eq. 27 is that, the position of a molecule at n 1 f + , could be determined through the sum total of partial probabilities ( ) 1 t / − µ∆ ε of n f plus t / µ∆ ε of ( ) n n P f , f / µ .
It should be mentioned that the interpretation of Eq. 27, when the value of t / ∆ ε is too large, is not true (because the f n coefficient on the right hand side of the equation may come out negative). So, in a semi-dense flow, due to the mean free path of the flow being small, the time step of the solution should be reduced to satisfy the restriction condition. Consequently, this method of solution takes much time and practically is not usable. Velocities of the molecules after collision are obtained as: Components of ω vector on the unit sphere are according to: 1 ξ and 2 ξ are random numbers which will be generated in a random manner on [0,1] interval. Further information on that is available in (Bird, 1994). On the other hand, the TRMC method can also be used to determine the f distribution function. The TRMC method is based on considering three terms of the Eq. 21 right hand according to the following: Since relation (30) satisfies the condition of (23), the probabilistic expression of this relation is that, the position of a molecule at n 1 f + could be determined through the summation of partial probabilities 0 A of n f plus 1 A of ( ) 1 f P f , f / = µ plus A 2 of the Maxwellian function M. for calculate of 0 1 2 A , A , A amounts, it is essential to evaluate τ from Eq. 12 and also use of Eq. 25. Relation 30 is valid on whole range of step is based on Maxwellian distribution sample taking. In other words it can be said that n 1 f + convert to distribution function on equilibrium condition (Pareschi and Russo, 2000). Figure 1 shows the problem's geometries. The conditions presumed for solving the problem are similar to the conditions used by Stefanov and Cercignani (1993) and Bird (1994) all of whom solved the Taylor-Quette problem by the DSMC method. The radius of the outer cylinder, r 2 , is twice the radius of the inner cylinder and the gas between the two cylinders is a monoatomic gas with hard spheres model (in the simulation, the molecular properties of Argon has been considered). The gas, at the initial conditions, is static and uniform. Density is so selected that the mean free path becomes ( ) 2 1 r r / 50 − . The selected length along the rotation axis is 4 times the distance between the two cylinders. The two end plates along the rotation axis are assumed with the axial symmetry boundary condition and the cylinder surfaces are assumed with the reflection boundary condition in equilibrium. At t = 0, the inner cylinder rotates with a speed that is three times the most probable speed. Hence, in these conditions, The Knudson number is 0.02 and the velocity ratio (sec) is 3.0.

RESULT AND DISCUSSION
In both methods, 400 calculation cells are considered along the axial direction and 100 cells along the radial direction with cells having similar dimensions to each other, the fluid field is divided into 40,000 cells. The time step value is selected as 7 t 1 0 − ∆ = for the two solution models and each molecule of the model stands for 10 14 real molecules.
Because of the existing high temperature differences and high density gradients created by the supersonic rotational speed of the cylinder, the local and effective values of many no dimensional parameters, such as the Knudsen number, are different from their nominal values at stable flow conditions (beginning and end conditions). Therefore, at the beginning of the problem's solution, many small vortices are generated which together, gradually, form larger vortices and after about 30 rotations of the inner cylinder, three large vortices form between the two coaxial cylinders. It can be seen from Fig. 2 that the system's stable flow answer obtained from the TRMC method, takes the form of three vortices whose size, geometry and physical properties are similar to each other and the only difference is their direction of rotation which changes for every other vortex. The result, thus obtained, totally agrees with the results from (Bird, 1994). Also, in Fig. 2, density contours are presented. Since the three vortices in Fig. 2 have similar flow properties, for a more precise analysis, only the middle vortex is considered. In Fig. 3 and 4, density results obtained from the TRMC and DSMC methods are shown. As it was predicted, the density increases in the radial direction, but, in the axial direction, because of existing vortex flows, it decreases and increases alternately.
Due to the high speed of the inner cylinder, gas temperature increases significantly. The field temperature increase continues until the heat transferred to the surface equals the study done by the surface of the rotating inner cylinder. As it can be observed in Fig. 5 and 6 which represent the results obtained from the TRMC and DSMC methods, the maximum temperature occurs near the inner cylinder. As the flow moves away from the inner cylinder's surface, its temperature decreases. Also, it can be seen that the temperature gradient inside the vortex flow, in the axial direction is greater than the temperature gradient on the outside of the vortex. In Fig. 7 and 8, no dimensional rotational speed contours ( m w / v ) along the radial and axial directions of the cylinders and in the distance between the two cylinders are shown. As can be expected, the maximum rotational speed of the flow occurs near the rotating inner cylinder. In Fig. 9, the numbers of collisions for ten succeeding time steps are shown. CPU time required for the DSMC is 63 h and for the TRMC is 43.1 h. Calculations are performed on a Pentium 4, 2.8 GHz processor. CPU time for the calculations based on the TRMC method is 63% of the CPU time required for the DSMC. This shows that the TRMC is faster than the DSMC. The main reason for the increase in calculation speed is that the number of collisions in each time step is reduced. Fig. 3, 5 and 7 with Fig. 4, 6 and 8, respectively and in the same order, it is observed that the results obtained from the TRMC and the DSMC methods have many similarities. In general, due to the application of the Maxwellian equilibrium conditions in each time step, the results obtained from the TRMC are more uniform and show less fluctuations. Therefore, the TRMC method, by eliminating these fluctuations, provides more fitting results, when compared to the DSMC scheme. However, the main advantage of the TRMC method over the DSMC, is in the reduction of the time required to process the program algorithm. The main reason for this time reduction in the TRMC method is that the number of intermolecular collisions (which engage the bulk of the computations), especially in more dense areas of the flow, decreases that is because of exerting of Maxwellian equilibrium condition in each time step.

ACKNOWLEDGMENT
The Computing Centre of Mechanical Department of Islamic Azad University of Marivan is acknowledged for allowing us to use their advance computers. The authors also gratefully acknowledge the software support of Mechanical Department of Islamic Azad University of Marivan.