Specific Objectives of this Thesis

Based on the preceding discussion of studies related to the QDS-2N method, it is clear that further numerical study is needed to improve the accuracy of the QDS-2N method and, consequently, lead to more effective applications.

In this thesis, we extend the second-order QDS algorithm (QDS-2N) [Smith et al., 2009] to flux reconstruction through true-direction polynomial multi-dimensional reconstruction of conserved properties across each cell width; this method is called QDS-N².

The specific objectives and organization of this thesis can be summarized as follows:

1. To improve the QDS-2N method by first studying its advantages and disadvantages. Flux reconstruction calculated using QDS-2N neglects neighbourhood cells when calculating diagonally cells. (Chapter 2)

2. To develop a QDS-N²method for solving Euler’s equation for inviscid fluid flow.

The fluxes of conserved properties are calculated by a sum of weighted moments over the polynomial spatial reconstruction of mass, momentum, and energy across the cell width. The particle properties are updated, considering the average value of the conserved quantity between the region bounds, which are required in translational directions and the application of splitting. (Chapter 2) 3. To develop the QDS-N² method in three-dimensions. (Chapter 2)

4. To apply the QDS-N² method to a numerical problem and an experimental case.

Various numerical methods, including Riemann solvers and total variation diminishing (TVD) methods, were compared as a benchmark. (Chapter 3)

5. To apply the QDS-N² method in the three-dimensional flow. (Chapter 3)

6. To develop the parallel QDS-N²method for large-scale applications such that the computational time problem is reduced efficiently. The scheme includes weak scaling and strong scaling. (Chapter 3)

7. To analyzedifferences associated with the 2Nmethod and compare QDS-2Nto QDS-N².This involves measurements to reveal differences between the two methods and provides information to determine which problem will be considered first. (Chapter 4)

8. The concluded by summarizing the major findings found in this thesis and outlining the recommendations for the future work. (Chapter 5)

Chapter 2

Numerical Methods

2.1 Overview of Euler Equation Solver

The Euler equations describe how the density, velocity, and pressure of a moving fluid are related. The Euler equations directly represent conservation of mass, momentum, and energy, and correspond to the Navier-Stokes equations without viscosity and heat conduction terms. Eq. (1) shows a two-dimensional formulation of the energy equation has been included in the Euler equations in fluid dynamics literature [Anderson, 1995].

2.1.1 Computational Fluid Dynamics

Engineers made further approximations and simplifications to the equation set until they had a group of equations that they could solve. Recently, high speed computers have been used to solve approximations to the equations using a variety of techniques,

such as finite difference, finite volume, finite element, and spectral methods. This area perfect gas, internal energy is only dependent on temperature. Pressure is presented by p

= (γ - 1)ρe with γ = cp/cυ the ratio of specific heats (γ = 7/5 for air).

Most numerical schemes for solving Euler equations are built around the Riemann solver. For example, Godunov technique provided two and three dimension applications in new finite volume numerical schemes and total variation diminishing (TVD) properties [Toro, 2009; LeVeque, 2004], and achieved second-order accuracy. Although, such numerical methods were generally accurate, they incurred high computational costs. However, approximate numerical schemes that reduce computational costs are not only less accurate and less robust but are also based on solutions of a Riemann problem.

More extensive introductions to numerical method for the Euler equations ware given by Godlewski and Raviart [1996], Kroner [1997], Laney [1998], Majda [1984], Toro [1997], Smoller [1983], and Hirsch [1990].

2.1.2 Kinetic Method

The classical kinetic theory of gas emerged from a combination of mechanics and statistics. The motions of molecules are described by probability rather than their

pressure, temperature, and a generalized equation of state for gases, and 2) transport properties (velocity, thermal conductivity, diffusion coefficients) based on first principles.

The kinetic method can be introduced by the Boltzmann equation, which is used in the study of a collection of particles in non-equilibrium statistical mechanics. The Boltzmann equation was devised by Ludwig Boltzmann in 1872 [Lerner et al., 1991].

The equation is a phase space of system that contains seven variables: three coordinates for position coordinates x, y, z, where each coordinate is parameterized by time t and three for each momentum component ͘px, ͘py, ͘pz and each coordinate is parameterized by time t. The volume element for position r and momenta ͘p can be expressed as follows:

d³rd³͘p = dx dy dz d͘px d͘py d͘pz (3) For one chemical species, the Boltzmann equation can be written as follows:

i i particles in the d³rd^3͘p phase space volume element around the phase space point (r, ͘p).

∂F/∂t is the total time derivative of the phase space distribution function, and(δF^/δt)coll

= rate of change of the phase-space distribution function due to collisions.

The Euler equation can be described using by kinetic theory. Here the conservation equation of mass, momentum and energy are derived as follows:

( ) ( )

where mo is the mass of the molecular object and no is the number of molecule per unit configuration-space volume. Because the Euler equations are non-linear hyperbolic equations, the shock waves are generally described by these equations. Several highly successful algorithms have been developed to solve such problems [Jameson, 1986;

MacCormack et al., 1975; Jameson et al., 1981].

In addition to the Boltzmann equation, Maxwell-Boltzmann distribution also contributed to the kinetic theory of gasses. This distribution was first carried out in 1859 and was named after James Clerk Maxwell and Ludwig Boltzmann. The distribution function can be expressed as follows:

(

^x

^{, ,}

^{y z}

)

³²

⁽

^{vx2 vy2 vz2}

⁾

distribution is close to thermodynamic equilibrium. An understanding of the Maxwell-Boltzmann distribution function is essential when studying the QDS method. The detail will be discussed in the next section (subchapter 2.2).

2.2 Quiet Direct Simulation (QDS) Method [Albright et al., 2002]

The normal random variable N(0,1) is defined by the probability density:

(7) by using a Gaussian quadrature approximation, the integral of Eq.(7) over its limits can be approximated by:

p x

( )

^=e

−x² 2

2π

( ) ( )

where wJ and qJ are the weights and abscissas of the Gaussian quadrature (also known as the Gauss-Hermite parameters) shown in Table 1, and N is the number of terms. The abscissas are the roots of the Hermite polynomials, which can be defined by the

The moment of the form are represented as

( )

²²

The particle simulation of fluid behaviour involved random variables which governed by stochastic differential equations of motion. For example, the one-dimensional Ornstein-Uhlenbeck (OU) equations describe the random dynamics of a particle of mass relaxing at a rate γ to the local fluid velocity u and temperature

can be solved [Gardiner, 1985] from following equation:

( )

^t

(

0

)

1 ^t

( )

0,1

x t u e ^γ x u _υ e ^γ N

υ Δ = + ^{− Δ} υ − +σ − ^{− Δ} (14)

In the thermalization γ∆t >>1, eq. (13) can be described as ^u+συ^N

( )

^0,1 drawn H_n+1(q)= 2qH_n − 2nH_n−1

As same as DSMC calculations to split particle transport and particle thermalization into two distinct operations, we presented operations by differential OU process which denoted with subscript tr and th respectively as follows:

( )

²

The transport differential operator describes particle free streaming, while thermalization operator drives particle velocities toward ^u+συ^N

( )

^0,1 without changing their positions. In the QDS algorithm, the part of tr preforms particle properties i.e. masses, momenta, special internal energies in each mash; The part of th represents each particle which is advanced to a new position. Those parts are established local thermodynamic equilibrium throughout the fluid.

The net fluxes of mass, momentum and energy of a cell are given by the sum of individual flux contributions from all the particles flowing in and out as follows:

1 1 outflow particles respectively into the cell under consideration. Each of the individual contribution (with first order spatial accuracy) can be described by the expressions, e.g., in one-dimensional case:

where ρ is the density, u is the bulk (or mean) flow velocity, and σ = (RT)^1/2 in a given source cell. Note R is the universal gas constant, and T is the gas temperature. The total number of degrees of freedom ξ is defined as and Ω is the number of simulated degrees of freedom (e.g., Ω = 1 for one dimensional flow). In the existing QDS-2N [Smith et al., 2009], the values of ρ, u, and σ employed in QDS particle initialization are taken from reconstructions based on linear variations. Despite fluxes being true direction in nature, the reconstructions performed by previous implementations are direction decoupled – i.e. a flux is computed through the product of (separate) fluxes previously computed (for 2D flow) in the x and y directions. For the 2D case, the particle mass and velocities in Eq. (9) become:

J K variables are the same as those in 1D case. The internal energy remains identical to the 1-D case, allowing for a corresponding increase in Ω to account for the extra simulated dimension. The fluxes from sources cell to any arbitrary destination cell can be calculated by the particle position distributions. The fluxes of mass, momentum and energy, which are based on the proportion of the overlapped area to the area of the original cell, are given by:

MASS JK

2.3 QDS-2N Method

In the QDS-2N algorithm used in the present study, the concept of QDS “particles”

whose properties are interpolated onto a grid (as used by Albright et al. [2002] in their original development of the technique) is replaced by the concept of fluxes of a large number of particles uniformly distributed across the cell, as described by Smith et al.

[Smith et al., 2008]. In this finite volume approach, quantities such as mass, momentum and energy are exactly conserved by tracking fluxes from source volumes to destination volumes. If the particle position distributions (i.e. gradients of density in the flow) are known, the flux from the source region to any arbitrary destination volume can be calculated.

In the present implementation of QDS-2N, the flux scheme employed by Smith et al. [2009] is used for the efficient calculation of two dimensional, true direction fluxes.

Here the Nx fluxes in each coordinate direction are computed separately requiring the calculation of 2N fluxes for the two-dimensional case.

In the second order scheme the gradients of cell velocity in the x-coordinate direction (du/dx), can be used to update the flux velocity:

2 2

J L v J

v u du x q

dx σ

= + Δ + (21)

where ΔxL represents the location in the cell from where the flow properties are taken.

Fluxes moving to the right are assumed to take their quantities from the reconstructed state ΔxL = 0.5(Δx – vxj∆t) to the right of the cell canter, where ∆x is the cell size and ∆t is the simulation time step. This corresponds to the displacement of the centre of mass of the flux which moves into the destination cell. Left moving fluxes have properties constructed in a similar fashion for which ΔxL = 0.5(–Δx – vxj∆t). The flux then moves in free flight, justifying the use of a linear interpolation routine.

The total mass and energy associated with the particles in the particular “bucket”

for the second-order case in the x-direction can be determined from the cell’s density (ρ), energy (E) and their respective gradients by:

L J

where wj are weights of the Gauss-Hermite quadrature, and:

( )

²

where Ω=2 for two-dimensional simulations. Any unused translational and other non-translational degrees of freedom are thus treated as internal structural degrees of freedom.

The amount of mass which fluxes to the new cell can be determined by multiplying Eq. (22) by vJ∆t /∆x. The gradients used in Eq. (21) to (23) are determined using the MINMOD (Minimum Modulus) and the MC (Monotonized Central Difference) scheme [Van Leer, 1977]. Using density in the x-direction as an example, the gradient using the MC slope limiter is:

(24)

where the MINMOD scheme is:

(25)

It should be noted that when non-uniform grids are employed (for example, when adaptive mesh refinement and coarsening is employed) the fluxes must be calculated

simulations, the amount of flux from one cell to another can be calculated trivially by determining the overlap area A=vJvk∆t² (where there are k = 1,…,N fluxes in the y-direction and vk is calculated in the same manner as Eq. (21) divided by the source cell area AS = ∆x∆y, as shown in Figure 2-1. The mass m and the energy ɛ are thus given energy Eflux and momentum in each coordinate direction px,flux and py,flux which must be added to the destination cell and subtracted from the source cell are given by:

flux

Friedrichs–Levy (CFL) number in the domain below a desired value (usually ≤1). It is important to note that this CFL restriction is to maintain physical realism and is not related to the numerical stability of the scheme. For a two-dimensional or axisymmetric simulation, the CFL number is given by:

( )

where qj(max) is the maximum value of the particle abscissas (i.e. the value which gives the maximum particle thermal velocity).

In the current implementation boundary conditions are handled using ghost cells.

These cells can be used to represent walls, stream boundaries, inflow boundaries and zero-gradient outflow boundaries. The interaction of a gas with a wall is identical to the interaction of that flow with an adjacent cell having the same flow properties but a reversed flow direction normal to the wall. The basic description of the simulation processes for QDS-2N method is available in Figure 2-2.

2.4 QDS-N² Method

2.4.1 Spatial Reconstruction and Flux Calculation

In the current study, referring to Figure 2-3, the general extension to higher order in QDS in one-dimensional case is performed using a spatial reconstruction of the form:

(33)

where Q(x) is the value of a conserved property (mass, momentum, or energy) at a distance x from the left hand side of the cell, and integer n indicates the order of the reconstruction. Note Qc is the value of Q(x) at the cell centre. This value is calculated from Q(x) integrating over the cell width divided by the cell width equalling to the existing average value of the source cell Q_s, presented below:

Q x

( )

^{= Q}c+ dQ

( )

By using our revised reconstruction, the bounds of integration are XL=0 and XR=∆x.

Then, Eq. (34) leads to:

Alternatively, Qc can be expressed as follows:

2 2 4 4 1 1

where n is assumed to be an odd number. This shows that the correction is only required when the scheme is third order (n = 3) accurate or higher, otherwise Q_c =Q_s (e.g., n=2).

Thus, the complete correct form of the higher order reconstruction of Q(x) using Q_s contains additional terms on every even derivative:

( ) ( ) ( )

Specifically, as n=2, the above is reduced to the following form because ofQ_c =Q_s^, as shown below:

The above reduces to exactly the same form as in QDS-2N [Smith et al., 2009].

Next, referring to Figure 2-4, the outgoing flux value of average conserved quantity successfully moving from the source cell into the destination cell (denoted by Q ) is: _tr

( ) ( ) ( ) ( )

Now, the transition of mean values Q can be used to calculate particle properties. _tr Assigning the flux out of average conserved properties Q_1tr, Q_2tr and Q_3tr as the mass, momentum, and energy, respectively, the resulting QDS particle properties for particle J are: reconstruction, it is important how the flux limiting is coupled. According to the value of conserved property Q(x) (see Eq. (33)), the gradient of Qc is defined in flux limiting during the reconstruction process. In each cell, we employ the monotonized central difference (MC) limiter to the effective gradients of conserved properties, as described below:

where is the equivalent flux limiter and F is the gradient calculated using forward differences. The theta is the ratio of the first order gradient calculated using forward and backward differences:

(43)

Therefore, an alternate representation of the variation of Q(x) over space must be:

(44)

2.4.2 One-Dimensional Flux Calculation and Implementation

In QDS simulations, we require flux from a volume to another volume. Since fluxes are split, the qualities of the flux depend entirely on the region from which they originated. The flux calculation is described as a flowchart in Figure 2-5, and summarized briefly as follows:

1. The gradients of conserved properties Q are first calculated using standard finite difference approximations in each cell i. For example, for a 5^th order accurate reconstruction, one might use the stencils like:

(45)

a. Calculate the approximate particle velocity based on the current cell Q , _s

particles to successfully move into the destination region.

d. Calculate the particle properties based on the average valuesQ . _tr

e. Calculate the fluxes of conserved properties to neighbouring cells following the standard QDS algorithm in subsection 2.2.

2.4.3 Two-Dimensional Flux Calculation and Implementation

Multi-dimensional extension is performed using the same principles applied for a one-dimensional reconstruction. The variation of conserved quantity Q(x, y) over two-dimensional space is assumed to be:

( )

2 2 4 4 1 1

Following this, the average value of conserved quantity in the region bound by [XL, YB] – [XR, YT] is formulated as:

Since the average requires bounding regions in both translational directions, application of splitting (as applied to TDEFM to improve computational efficiency) is impossible, and the full N² number of particles (i.e. nine when three particles are used per direction, sixteen for four, etc.) are required for a complete flux computation.

Previous extensions required only the 2N particles. Unlike the one-dimensional reconstruction, each particle carries three separate fluxes (for three different destination cells) and so any single QDS particle possesses three “sub-particles” based on different integral bounds. This concept is demonstrated in Figure 2-6, showing each unique sub-region (A – C). The area of the sub-sub-region A is u × v × dt as described earlier in Section 2.2 for the QDS-2N method.

2.4.4 Three-Dimensional Flux Calculation and Implementation

As previous section, the three-dimensional flux calculation is followed the same way for a one-dimensional reconstruction. The variation of conserved quantity Q(x, y) over three-dimensional space is assumed as follows:

( ) The subsequent cell centred value of Qc is:

2 2 4 4 1 1

Following this, the average value of conserved quantity in the region bound by [XL, YB, ZI] – [XR, YT, ZO] is formulated as:

( )( )( ) ( )

Therefore, each particle must be completed calculation for flux reconstruction in a time step. Especially the calculation for diagonal cell, the particle has 7 sub-regions (A – G) to be considered which are based on the difference bound. The concept is same as section 2.4.5 demonstrated in Figure 2-7.

2.5 Brief Summary

The characterization of the QDS-N² method developed in this chapter can be summarized as follows:

1. The QDS method replaces the random sampling method used in the DSMC method. Particles are permitted to move in physically realistic directions;

therefore flux exchange is not limited to cells sharing an adjacent interface as in conventional direction decoupled finite volume solvers.

2. In QDS method, the particles are recreated deterministically from the properties stored on the mesh using Gauss–Hermite quadrature weights and abscissas.

3. The QDS-N² and QDS-2N methods use the same procedure to calculate flux reconstruction in one-dimension.

4. In the QDS-N² method, article properties are updated, considering the average value of the conserved quantity between the region bounds, which are required in translational directions and the application of splitting.

5. With parallel implementation, an extension to three-dimensional QDS method is also demonstrated.

Chapter 3

Complex Gas Flow Simulations using the QDS-N

²

Method

3.1 Introduction

To determine the effectiveness of the QDS algorithm, several test cases, including three cases of one-dimensional domain, six of two-dimensional domain, and one of three-dimensional domain were used. The test cases were chosen for the following reasons:

• The cases are suitable for solving the Euler equation.

• The benchmark is well known and many numerical methods or experimental results are available for comparison.

• All qualitative, quantitative, computational, and experimental results are readily available.

The QDS-2N and QDS-N² methods are compared for both accuracy and computational time. We also parallelize the QDS-N² method to increase calculation speed. As well, the literature related to the test cases is surveyed to determine useful benchmarks to which the present code can be compared.

3.2 QDS Method in One-Dimension Flow 3.2.1 Shock Tube

The shock tube is an important application in unsteady wave motion, the study of high-temperature gases in physics and chemistry, and the testing of supersonic bodies

and hypersonic entry vehicles. Figure 3-1 shows the features of a shock tube after the diaphragm has been broken. The region to the left of the diaphragm is the driver section and the region to the right is the driven section.

We consider the shock tube problem to validate the accuracy of the QDS code, especially a Riemann solver [Toro, 1999], which represents the majority of solution methods. The initial condition for the simulation consists of two constants:

(54)

To compare the results easily, we set the number of cells for QDS and the Riemann

To compare the results easily, we set the number of cells for QDS and the Riemann

在文檔中 應用真實流向概念重建通量法則於靜直接模擬法速解尤拉方程式的研究 (頁 29-0)