Finite-difference Discretization of the Two-Dimensional MBE

We consider a class of model Boltzmann equations of the form

wheref( vrv ,v, t)is the velocity distribution function which depends on space, rv , molecular velocity, vv , and time, t; ν is the collision frequency and f is an ^N appropriate distribution function depending on the model selected. The number density, macroscopic flow velocity, and temperature of the gas are the first three

moments of the distribution function

∫

where k is the Boltzmann constant and δ_ij is the Kronecker delta. The heat flux

vector q is

The elastic collision frequency is of the form µ

ν =nkT (2.23)

where µ is the viscosity and is assumed to have a temperature dependence µ χ

µ _∞ = (T T_∞) (2.24)

Here χ is a constant for a given gas. If we assume the dependence of the viscosity on the temperature as for the Chapmann-Enskog gas of inverse ζ power law, we have 2

(

-1

)

3 ζ

χ= ζ ⁺ . For Maxwell molecules, ζ =5 then χ =1.; thus the collision

frequency is independent of temperature. The viscosity coefficient µ is related to _∞ the freestream mean free path λ by the relation _∞

In this study we consider two kinetic models for f ; one is the BGK model and the ^N other is the Shakov model. For the BGK model, we have f equal to the Maxwellian ^N distribution f : ^M

For the Shakov model, we have

5)/(5pRT)]

Here, Pr is the prandtl number and is equal to 3

2 for a monatomic gas.

We note that the derivation of the continuum Navier-Stokes equations from the BGK model or the Shakov model can be obtained using a Chapmann-Enskog procedure.

To illustrate the numerical approach, we describe in detail the relevant equations for two-dimensional problem, for the purpose of reduction in computer storage requirements, the following reduced distribution functions (Chu 1965) are introduced:

z Where L is the reference length, T_∞ is the reference temperature. Then, the non-dimensional variables can be defined as follows:

L

After the process of non-dimensionalization and integrating out the v_z dependence in Eq. (2.16) using Eq. (2.28), the single model Boltzmann equation in three space dimensions reduces to the following two simultaneous equations in two

space dimensions and then we neglect the signal "∧ . "

(

G -g

)

2pT)]

Without causing any confusion we shall drop the hat in the equations in the following.

The macroscopic moments are found as follows:

y

strong conservation law from with stiff source terms as follows:

To treat general geometry we consider the conservation equation of the two-dimensional rarefied gasdynamics in general coordinates (ξ,η)

F S The metric Jacobian and the metric terms are given by

ξ

variables (for two-dimensional case). To remove the functional dependency on the velocity space of the equations, the discrete ordinate method [10] is applied. This method, which consists of replacing the integration over velocity space of the distribution function by an appropriate quadrature, requires the values of the distribution function only at certain discrete velocities. The choice of the discrete values of velocity point is dictated by the consideration that our final interest is not in the distribution functions themselves but in the moments. Hence, the macroscopic moments given by integrals over molecular velocity space can be evaluated by the same quadrature. The discrete ordinate method is then applied the set of Eq. (2.36) for the )(v_x,v_y velocity space. That is the value of g(ξ,η,t,v_x,v_y) become

hyperbolic partial differential equations with source terms in the physical space

, , ,

discrete velocity point (v_σ,v_δ) respectively, where σ =- N₁,...,-1,1, ..., N₁ , and

2 2,...,-1,1,...,N N

=-

δ . Also to apply the discrete ordinate method, the integrals

appeared in Eq. (2.33) are expressed as finite sums according to the quadrature define as

W_σ are the corresponding weights of the Gauss-Hermite quadrature. Both full-range and half-range Gauss-Hermite quadrature are needed. It can be shown that above quadrature formula is equivalent to approximate the Maxwellian distribution by the

discrete distribution

Where δ is the Dirac delta function. This can be considered as the optimum discrete approximation in the sense that the first 2N moments of the Maxwellian can be

exactly duplicated Where Γ represents the usual Gamma function. The discrete velocity points and the

Giddens (1968) and by Shizgal (1981).

Once the discrete distribution functions g_σ,δ and h_σ,δ are solved, one can obtain all the moment integrals as

)

inviscid Euler equations. In this work, we not only need to solve the discrete distribution functions (not in equilibrium) but also to use them to evaluate the

macroscopic moment by numerical quadratures. The selection of the discrete velocity point and the range of velocity space in the Newton-Cotes formulas are somewhat

2.5 Numerical Algorithm in solving the MBE

In this section we describe the numerical algorithm for solving the set of Eq.

(2.41). Both the time-accurate explicit method using operator splitting for unsteady flow problems and implicit method using lower-upper (LU) factorization for steady-state calculations are considered. We follow and extend our previous high resolution non-oscillatory scheme for hyperbolic system of conservation laws to include a source term. Some general remarks can be given here. When explicit methods are used to integrate the equations for g_σ,δ and h_σ,δ, one can decouple the equations and solve them separately. When implicit methods are employed, the equations in general are coupled through the jacobian of the source terms since the source terms are functionals of g_σ,δ and h_σ,δ. In the following we still keep the equations in vector-matrix form and with the understanding that they can be coupled into scalar form and solved in scalar manner.

Define a uniform computational mesh system (ξ_j,η_k) with mesh sizes characteristic variables in the local ξ-direction and η-direction respectively as

k

where J (J_j,k J_j1,k) / 2

The time-splitting method described above closely patterns the procedure first proposed by Bird and used in particles schemes, in which free molecular motion and the intermolecular collisions are two independent stages of the algorithm that update the particle position and velocity.

In terms of operator form we have the time integration schemes as

The time integration of the governing equations is carried out on each pair of discrete velocity point(v_σ,v_δ)with finite difference approximations. Without causing any ambiguity, we omit the subscripts (σ ,δ) in the time integration operators

η

method

using simple averages. The components of ^N

k

]

The class of schemes covered by Eq. (2.54) includes the total variation diminishing (TVD) and essentially nonoscillatory (ENO) scheme. For ω =0 and

=0

ϑ , one has a third-order ENO scheme, denoted as ENO3. a first-order upwind

scheme, denoted as UW1, can be deduced from Eq. (2.54) by setting all the elements

l k

ej, and d^l_j,k equal to zero. The accuracy and Fourier stability of schemes defined by Eq. (2.54) can be analyzed by looking at different possible combinations of the

arguments in the m and m limiter functions.

Equation () can be approximately factored in several different ways. Here we adopt the lower-upper method and Eq. () become

n

In Eq. (2.66) denote the backward and forward difference operators, respectively.

The split jacobian matrices are Λ^± =diag{λ^l±}, where λ^l± =( λ^l ± λ^l )/2 . The

± are defined analogously by Eq. (2.53). for steady-state calculation, the use of Eq. (2.58) and (2.59) can lead to the undesirable results that the steady state depends on the time step ∆t and causes slow convergence. We use the following approximation which still maintains the spatial accuracy:

2 (z) (z) 1ψ

σ =

⎪⎪

and it is implemented in the sequence:

n

which can be show to produce the least amount of error among several possible factorizations, particularly when the norms of the source terms are large. The collision source term, S, of the model equation in general is a functional of the reduced distribution functions g_σ,δ and h_σ,δ. The excat evaluation of the jacobian matrix of the source term, C, is difficult. In this work, we approximate the jacobian of the source term by

With this simplified approximation the equations become diagonal and completely decoupled and the solution procedure becomes rather simple and can be solved scalarly. The numerical experience indicates that such an approximation works well.

Boundary Conditions

To specify the interaction of the molecules with the solid surface, it is assumed that molecules which strike the surface are subsequently emitted with a Maxwellian velocity distribution characterized by the surface temperature T and zero net _w tanfential velocity. The two-stream concept is also applied here by defining the half-range distribution functions wall, the wall distribution function is given by

0 priori and may be found by applying the condition of zero mass flux normal to the surface at the wall. One has

y

The farfield boundary condition at the freestream is given by the Maxwellian distribution

The inflow and outflow boundary conditions are treated using characteristics-

based boundary conditions which are in accord with the upwind nature of the interior point scheme. For problems with symmetry, only half plane is computed and the symmetry condition is assigned to the distribution function for

) v - , v t, , , g(

) v , v t, , - ,

g(ξ η _{x y} = ξ η _{x y} (2.75)

Chapter 3. Preliminary Results and Discussion

3.1 Comparison of MBE Results with DSMC Numerical Method

The computed results using MBE are found to compare well with those using the

DSMC(direct simulation Monte Carlo) in the Near-continuum regime flows .

We using DSMC to simulated for M=2, and Kn=0.0033. In Fig. 3.114, we show the simulated results include, number density, temperature, Mach number, u-velocity, v-velocity, and velocity streamlines.

To compare with MBE simulated results. In Fig. 3.114 (a), an ultra high-density region appears at the very right-hand upper corner due to the moving plate at the top of the cavity, but smaller than MBE simulated results. In Fig. 3.114 (b), there is a temperature increased region in the cavity, the right-hand upper corner which temperature increased obvious as a result of density increased, but smaller than MBE simulated results. In Fig. 3.114 (c), (d), (e), and (f), The MBE simulated results about u-direction and v-direction negative speed also are all bigger slightly than the DSMC simulated results. Therefore, we may discover that smaller secondary eddies have been created at the two bottom corners and the center of the top vortex is moved slightly to right-side and toward the moving plane.

3.2 Driven Cavity Flows

3.2.1 Problem Description and Test Conditions

Fig. 3.1 sketch of the 2D square (L/H=1) driven cavity flow with moving top plate. Initially, we use argon gas at rest inside the cavity and at the same uniform temperature 300K. At time t=0 the upper plate begins moving instantaneously at speed Ma=0.5-2 and Kn=10-0.0033 based on the mean free path of wall temperature and size of the cavity. Table Ⅰ-Ⅲ shows the all cases which we simulates with different parameters .

3.2.2 Grid Convergence Tests

In order to must test the different grid to the result influence. We have used four kind of different grids (Fig. 3.2-5 to show) in the case (M=0.9, Kn=0.01). In Fig. 3.6, we found that the convergence is better while grid size becomes smaller towards the diffuse wall. However, when grid size nearby the diffuse wall is the same, the large number of mesh doesn't result in better convergence. In Fig. 3.7, we compare BGK and the Shakov these two kind of different models. Although in each kind of parameter result comparison all almost, but the Shakov MODEL convergence speed is faster.

3.2.3 Effects of Knudsen Number

In this section, we were observed effects of Mach number in different Knudsen numbers (Kn=10-0.0033). First, we were showed general simulation results include density, temperature, Mach number, u-direction, v-direction and streamline. Second, we were showed property distributions across cavity geometric center for x =0.5m, y=

0 to -1m and y=-0.5m, x=0 to 1m. Third, we showed property distributions near the solid walls. Finally, we were observed the recirculation center position in different cases.

3.2.3.1 General Simulation Results

3.2.3.1.1 Subsonic Moving Plate (M=0.5, 0.9)

Fig. 3.7 shows that number density contour for Ma=0.5 and Knudsen number 10, 1, 0.1, 0.01, and 0.0033 respectively. Driven plate takes particles to the right-hand upper corner. An ultra high-density region appears at the very right-hand upper corner due to the moving plate at the top of the cavity. Therefore the particles are larger than initial value. In addition, there are low densities at the left-hand upper corner. In the other series, fixed M=0.9 Fig. 3.12.

Fig. 3.8 show that temperature contour for Ma=0.5, and Knudsen number 10, 1, 0.1, 0.01, and 0.0033 respectively. We normalize temperature to divide the initial

temperature 300K. When the Knudsen number is decrease, there are two ultra temperatures region in the cavity. One of them is the right-hand upper corner which temperature increased as a result of density increased; the other one is left-hand upper corner which temperature increased due to high vertical speed. Therefore, the temperature increase more seriously as Knudsen number decreased. In the other series, fixed M=0.9 is showed in Fig. 3.13.

Fig. 3.9 show that Mach contour for Ma=0.5, and Knudsen number 10, 1, 0.1, 0.01, and 0.0033 respectively. In the other series, fixed M=0.9 is showed in Fig. 3.14.

Fig. 3.10 show that u-velocity contour for Ma=0.5, and Knudsen number 10, 1, 0.1, 0.01, and 0.0033 respectively. The maximum u-velocity values are 0.35, 0.4, 0.6, 0.9, and 0.9 with Knudsen number 10, 1, 0.1, 0.01, and 0.0033, respectively. Because of rarefaction effect caused slip phenomenon and the slip velocity along the solid walls increase with Knudsen number at the same Mach number. We normalize u-velocity to divide the upper plate velocity. The velocity is more and more decrease when Knudsen number increase. In the other series, fixed M=0.9 is showed in Fig.

3.15.

Fig. 3.11 show that v-velocity contour for Ma=0.5, and Knudsen number 10, 1, 0.1, 0.01, and 0.0033 respectively. We normalize v-velocity to divide the upper plate velocity. An ultra high-speed region appears at the left-hand and right-hand upper

region. The velocity of right-hand upper region is increase when Knudsen number decrease. In the other series, fixed M=0.9 is showed in Fig. 3.16.

As mentioned above, we can be briefly summarized as follows:

1. The slip velocity is more and more decrease when Knudsen number increase.

2. An ultra high-density region appears at the very right-hand upper corner due to the moving plate at the upper of the cavity

3. There are two ultra temperatures region in the right and left upper corner in cavity.

3.2.3.1.2 Supersonic Moving Plate (M=1.1, 2)

Fig. 3.17 shows that number density contour for Ma=1.1 and Knudsen number 10, 1, 0.1, 0.01, and 0.0033 respectively. Driven plate takes particles to the right-hand upper corner. An ultra high-density region appears at the very right-hand upper corner due to the high-speed moving plate at the top of the cavity. Therefore the particles are larger than initial value. In addition, there are low densities at the left-hand upper corner. In the other series, fixed M=2 Fig. 3.22.

Fig. 3.18 show that temperature contour for Ma=1.1, and Knudsen number 10, 1, 0.1, 0.01, and 0.0033 respectively. We normalize temperature to divide the initial temperature 300K. When the Knudsen number is decrease, there are two ultra

temperatures region in the cavity. One of them is the right-hand upper corner which temperature increased as a result of density increased; the other one is left-hand upper corner which temperature increased due to high vertical speed. Therefore, the temperature increase more seriously as Knudsen number decreased. In the other series, fixed M=2 is showed in Fig. 3.23.

Fig. 3.19 show that Mach contour for Ma=1.1, and Knudsen number 10, 1, 0.1, 0.01, and 0.0033 respectively. In the other series, fixed M=2 is showed in Fig. 3.24.

Fig. 3.20 show that u-velocity contour for Ma=1.1, and Knudsen number 10, 1, 0.1, 0.01, and 0.0033 respectively. The maximum u-velocity values are 0.35, 0.4, 0.6, 0.9, and 0.9 with Knudsen number 10, 1, 0.1, 0.01, and 0.0033, respectively. Because of rarefaction effect caused slip phenomenon and the slip velocity along the solid walls increase with Knudsen number at the same Mach number. We normalize u-velocity to divide the upper plate velocity. The velocity is more and more decrease when Knudsen number increase. In the other series, fixed M=2 is showed in Fig. 3.25.

Fig. 3.21 show that v-velocity contour for Ma=1.1, and Knudsen number 10, 1, 0.1, 0.01, and 0.0033 respectively. We normalize v-velocity to divide the upper plate velocity. An ultra high-speed region appears at the left-hand and right-hand upper region. The velocity of right-hand upper region is increase when Knudsen number decrease. In the other series, fixed M=2 is showed in Fig. 3.26.

As mentioned above, we can be briefly summarized as follows:

1. The slip velocity is more and more decrease when Knudsen number increase.

2. An ultra high-density region appears at the very right-hand upper corner due to the high-speed moving plate at the upper of the cavity

3. There are two ultra temperatures region in the right and left upper corner in cavity.

3.2.3.2 Property Distributions Across Cavity Centroid

Fig. 3.27 present the profiles of the number density along vertical line through geometry center(x/L=0.5) for Ma=0.5-2, and Knudsen number 10, 1, 0.1, 0.01, and 0.0033 respectively. In the other series, the number density have replaced by temperature, u-velocity and velocity is showed in Fig. 3.28-30.

Fig. 3.31 present the profiles of the number density along vertical line through geometry center(y/L=-0.5) for Ma=0.5-2, and Knudsen number 10, 1, 0.1, 0.01, and 0.0033 respectively. In the other series, the number density have replaced by temperature, u-velocity and velocity is showed in Fig. 3.32-34.

3.2.3.3 Property Distributions Near Solid Walls

Fig. 3.35 present the profiles of the number density along vertical line through

geometry center(x/L=0) for Ma=0.5-2, and Knudsen number 10, 1, 0.1, 0.01, and 0.0033 respectively. In the other series, the number density have replaced by temperature, u-velocity and velocity is showed in Fig. 3.36-38.

Fig. 3.39 present the profiles of the number density along vertical line through geometry center(x/L=1) for Ma=0.5-2, and Knudsen number 10, 1, 0.1, 0.01, and 0.0033 respectively. In the other series, the number density have replaced by temperature, u-velocity and velocity is showed in Fig. 3.40-43.

Fig. 3.44 present the profiles of the number density along vertical line through geometry center(y/L=-1) for Ma=0.5-2, and Knudsen number 10, 1, 0.1, 0.01, and 0.0033 respectively. In the other series, the number density have replaced by temperature, u-velocity and velocity is showed in Fig. 3.45-Fig. 3.47.

Fig. 3.48 present the profiles of the number density along vertical line through geometry center(y/L=-0) for Ma=0.5-2, and Knudsen number 10, 1, 0.1, 0.01, and 0.0033, respectively. In the other series, the number density have replaced by temperature, u-velocity and velocity is showed in Fig. 3.49-Fig. 3.51.

3.2.3.4 Recirculation Center Position

Fig. 3. 51-54 show plots of the velocity streamlines for M=0.5, 0.9, 1.1, and 2, respectively, and in each figure results are provided for various values of Knudsen

number. The corresponding velocity streamlines for M=0.5, 0.9, and 1.1, are given in Fig. 3.51-53. It is seen that while for Kn=10, 1, and 0.1 there is only one vortex, for Kn=0.01, and 0.0033, under the main vortex secondary eddies have been created at the two bottom corners. In Fig. 3.54 , it is seen that while for Kn=1, and 0.1 there is only one vortex, for Kn=10 one additional vortex, for Kn=0.01, and 0.0033 two additional vortices, under the first one, have been developed. As the Knudsen number is increased further, these secondary eddies grow under the first one.

In Fig. 3.55, we show the relative vertical distance (x/L) of the center of the top vortex of the cavity in term of M=0.9, 1.1, and 2 of Kn=10, 1, 0.1, and 0.01. It is seen that in these cases, as Kn is decreased, the center of the top vortex is moved slightly toward the right wall, when Kn=0.0033 is opposite. However, M=0.5 of Kn=10, 1, 0.1, 0.01, and 0.0033. It is seen that in these cases, as Kn is decreased, the center of the top vortex is moved slightly toward the right wall.

In Fig. 3.56, we show the relative vertical distance (y/L) of the center of the top vortex of the cavity in term of M=0.5, 0.9, and 1.1 of Kn=10, 1, 0.1, and 0.01. It is seen that in these cases, as Kn is increased, the center of the top vortex is moved

In Fig. 3.56, we show the relative vertical distance (y/L) of the center of the top vortex of the cavity in term of M=0.5, 0.9, and 1.1 of Kn=10, 1, 0.1, and 0.01. It is seen that in these cases, as Kn is increased, the center of the top vortex is moved

在文檔中 使用波茲曼模型方程式模擬自由分子流到近連續流的正方形空穴流場 (頁 34-0)