QDS Method in One-Dimension Flow

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 solver to 200. Simulation time is 0.1 s, and the walls are reflected. Figure 3-2 shows density results for the first to the third order QDS and the Riemann solver. We observed a great improvement between the first and the second orders. Although the second and the third order QDS results are closer to the benchmark than first order results, the improvement between the second and the third order is negligible.

3.2.2 Shock-blast wave interaction

The shock blast problem suggested as a benchmark in [Woodward et al., 1984]

included two blast waves contacted by strong shocks. By initial pressure jumps, this version of the problem will indicate a flow containing three distinct fluids.

The initial condition for simulation is expressed as:

(55) Each result was examined at time 0.038 s. The reflective boundary condition is applied to both ends. The comparable benchmark is an almost identical to shock tube

ρ,v,T

grids. The benchmark is used to measure errors by comparing the different methods considered here. Figure 3-3 shows a contrasting profile of the density computed by the QDS method with 400 cells using the MC limiter. As can be seen, the second order QDS method is much better than the first order, and the third order QDS method is slightly better than the second order. The difference between second and third order schemes can be seen easily at a range of 0.7 to 0.8 in the x direction. However, the pre-QDS method can be extended to calculate the case in which acoustic wave oscillation is quite large. Therefore, the next test case considers the pre-QDS method, which is an easy way to measure the oscillation accurately by ignoring the limiter.

3.2.3 Shock acoustic wave interaction

An ideal case for testing the general one-dimensional extension of the pre-QDS method is the interaction of a Mach 3 shock wave with an acoustic wave as proposed by Shuv [Shu, et al., 1988]. When a higher density shock wave contacts a smooth acoustic wave, an amplified wave with higher frequency results. The initial conditions are:

(56)

Results for various QDS configurations with 200 cells compared to a WENO [Huang, et al., 2009] benchmark with 2000 cells are shown in Figure 3-4. The result is obtained at time 1.8 and the limiter is not applied. The improvement from first to second order is the most significant; higher orders show only slight improvement. The general trend is in agreement with the WENO benchmark and third order QDS solutions that result in similar levels of dissipation to those of a fifth-order WENO solution.

3.3 QDS Method in Two-Dimension Flow ρ, ν, P

In this section, we present four test cases for comparing theN² and 2N methods.

The major discussions focus on the time cost and accuracy between the two methods in two-dimensional problems.

3.3.1Shock-Bubble Interaction

The strength of correct multi-dimension reconstruction is demonstrated in two dimensions with the shock/bubble interaction problem [Čada et al., 2009]. The initial conditions for this problem are shown in Figure 3-5. The simulation calculated a shock wave, moving from left to right with a velocity of Much number 2.85 in an ideal, inviscid gas and interacting with a pseudo bubble at x=0.3. The results are presented at t=0.2. The results of the numerical schlieren (gradients of density) are presented in Figure 3-6 for various QDS methods and the TVD schemes on a grid of 1700 × 500 cells. The application of correct multi-dimensional reconstructions results in a relatively high resolution of the circulation and reflected shock located at x = 0.6.

Figure 3-7 displays two QDS methods with different numbers of cells. We compare the N² method (Figure 3-7a) against the 2N method (Figures 3-7b, c and d) by using 300 × 100, 450 × 150, and 600 × 200 cells. For the sake of comparison, the limiter for each simulation is the monotonized central (MC) limiter. In Figures 3-7a and b the difference in resolution is clear despite the fact that both these methods employ the same number of cells. As the number of cells employed by the 2N method increases (shown up to 600 × 200 cells here), the results gradually approach those of the N² scheme with relatively few cells (1/4). Obviously, the multi-dimensional computation (N² method) achieves higher accuracy than the 2N scheme.

Further, we consider the computation time required by each method in this case.

The N² method is true-directional, in that each possible combination of discrete velocity

must be considered (nine instances with three discrete velocities per direction), while the 2N method employs approximate dimensional extension and only requires six discrete velocity computation in a two-dimensional computation. Moreover, for each particle, three space-averaging computations are required for each fraction falling into separate destination cells. Therefore, the N² method requires more computation time for the same number of cells. The computation time of the two solvers are summarized in Table 3-1. According to these data, the N² extension of QDS requires approximately three times the computation time as compared to that of the original 2N method for the same number of cells. However, for any given degree of accuracy, we find that the N² method provides an increase in computational efficiency of almost three times (300 × 100 vs. 600 × 200 for 2N vs. N²). Thus, the application of the N² method is justified over that of the 2N method for high-resolution solutions.

3.3.2 Euler-Four-Shocks problem

This test case was introduced in Serna[2006], which computed the numerical solution employing thepiecewise hyperbolic method-Marquina’s flux formula(PHM-MFF) and power PHM-MFF schemes. The test problem is initially divided into four quadrants sharing a common corner at 0.75 and 0.75 in the domain [0,1] × [0,1], as illustrated in Figure 3-8. These quadrants initially have the following different but uniform conditions:

(57)

Figure 3-9 shows four results of a comparison between the 2N and N² solvers at the time of 0.4. The Courant–Friedrichs–Lewy factor (CFL) is set as 0.5. We compare

ρ,u,υ, p

the results of the two QDS solvers using 100 × 100 to 300 × 300 cells. As can be seen, the result of the QDS N² solver obtained using a coarse grid (100 × 100 cells) is only approached by the 2N solver when employing considerably fine grids (300 × 300 cells).

Furthermore, the result we obtained using the QDS N² solver on a computational grid of 1000 × 1000 cells is similar to that obtained using the TVD-MUSCL scheme [Čada et al., 2009] (see Figure 3-10a).

An investigation of the computational expense of each scheme showed that the N² solver takes approximately four times longer to complete the simulation than the 2N solver. This comparison of computational expense is summarized in Table 3-2. The increase in computational time with the refinement of the computational grid is due to the constant “kinetic” CFL condition that we employ, which is defined as follows:

(58)

This basically ensures that particles in free flight do not travel further than the adjacent cells. Although the result take more time to compute using the N² solver, the accuracy is considerably better than that of the 2N solver, in fact, it is not clear that the 2N solver will ever approach the solution obtained using the N² solver, irrespective of the number of cells employed.

3.3.3 Euler-Four-Contacts problem

This test case involves the Euler-four-contacts interaction problem defined by Schulz-Rinne, Collins, and Glaz [Schulz-Rinne et al., 1993]. The same test case with a different higher-order method is also presented in [Salichs et al., 2005]. This Riemann problem briefly shows four constant states consisting of four quadrants and two shocks

CFL= Max

(

)

y)=(0.5, 0.5). A representation of the structure of the flow domain is shown in Figure 3-11. The initial flow condition is imposed by four difference shock waves and satisfies the following relation:

(59)

Figure 3-12 shows the numerical result of the second-order TVD-MUSCL method for a density contour profileon a 1000 × 1000 uniform grid, taken from [Čada et al., 2009].For the QDS method, the result obtained at the time of 0.8 on a 1000 × 1000 uniform grid can be seen in Figure 3-13. Two results are shown for the second-order method with the N² and 2N solvers. Both enforce a constant CFL number of 0.25. The contours of density are presented with levels of 0 to 2.4. In this case, a shock wave is generated and spirals from the contact point in an unsteady fashion. By comparing the two figures, we find that both the N² and the 2N solver results are symmetrical and that the result obtained using the N² solver is closer to the TVD-MUSCL result presented in Figure 3-12. As in the previous test cases, in the current test case, the accuracy of the N² method is superior to that of the 2N method. In this instance, however, the WENO results [Schulz-Rinne et al., 1993] are still superior to the N² results: this can be attributed to the small stencil employed for the estimation of the higher-order gradients, or to flux splitting employed and the inevitably associated numerical dissipation; this may require further investigation.

Further, we have compared the timings and the accuracy for this test problem with different N for both QDS-2N and QDS-N² with 1000 × 1000 cells since both schemes

ρ,u,υ, p

scale differently with N. The results are essentially the same as those obtained for N = 3 when N increases to 5 or 9 for both the abovementioned methods, as shown in Figure 3-14. This is reasonable since the integration of a Gaussian function with a polynomial, having two or less degrees, becomes exact, if the number of Gaussian-Hermite integration points is 3 or more. Expectedly, the computation time increases roughly 3 times from N = 3 to N = 9 for both the methods. Further, Figure 3-15 shows the density contours when the grid resolution increases from 1000 × 1000 cells to 2000 × 2000 and 3000 × 3000 cells. In brief summary, for both the QDS-2N and the QDS-N² methods for solving the Euler equation, accuracy effectively increases with increasing grid resolution, while it is essentially the same when N ≥ 3.

3.3.4 Mach 3 Flow over a Forward Step

The fourth test problem for the second-order QDS method is the Mach 3 flow over a forward facing step in a high-speed wind tunnel. The main of this case is designed to solve complex oblique shock reflections. This problem was first introduced by Emery [Emery et al., 1968] and has since been used by Woodward and Colella to test a number of differencing schemes [Woodward et al., 1984]. The same problem using second order QDS scheme was used by Smith [2009] to demonstrate that QDS the effectively captures all features of the flow fluid and compare with Keats and Lien’s scheme [Keats et al., 2004]. Keats mentions that this case is well known problem to evaluate the robustness of numerical method. It is difficult to maintain positivity pressure and density using various numerical methods when the strong shock reflection at a step during first time step. Therefore, this test case has proven to be a well test for a long time by several of methods. For QDS method, it must be more careful to deal with particle actions facing the corner of the step. That will useful for us to expend a laugh of

blocks in a domain.

Here present the geometry and boundary conditions used in the problem (see Figure 3-16). As an initial condition, the flow is everywhere uniform at Mach 3 with a density of 1.4, a pressure of 1.0, gas constant R = 1.0 and a specific heat ratio c = 1.4.

The CFL number is set as 0.1. This is equivalent to an impulsively started flow and the simulations capture the unsteady development of the flow structure.

The QDS simulations utilized a second-order scheme and third-order with 4 QDS particles per cell on a uniform grid consisting of 100,800 square cells (the number of cells for total domain which included the step is 600 × 200). The number of contours is 30 from 0.2568 to 6.067. Figure 3-17 shows the density profile generated using two QDS solver – 2N solver and N² solver − at time = 4.0 s. The result is compared with that of Keats and Lien who employed a second order Godunov method on an adaptively refined mesh [Keats et al., 2004]. We observe the results in both figures are similar whenever the solver is 2N or N² solver.

3.3.5 Shock Wave Diffraction over a 90 degree sharp corner

This test case uses the same geometrical as previous test case which a block in the domain. The geometric configuration of this case is the forward-facing steps. This case is also an important case to observe shock wave diffraction which is designed to solve the Euler equation. The shock wave diffraction induces more phenomena in the perturbed region behind the shock [Skews, 1967]. Those included shock wave, vortex, terminator and an incident sound wave. The problem contain both computational and the mathematical studies. Secondly, the experimental results are available for variety of geometries and Mach number [Takayama et al., 1991]. More numerical studies can be found from Skews [1967], Schardin [Dyke, 1997], Bazhenova [Bazhenova et al., 1984].

Figure 3-18 shows structure of the perturbed region behind a diffracting shock wave presented form Skews [1967]. Skews performed experiments for a variety of Mach numbers and convex corner angles and has outlined the structure of the perturbed region.

The detail of the structure is described in appendix A.

The output format of this case shows in Figure 3-19[Takayama et al., 1991]. The incident shock Mach number is 1.5 which the normal shock moves to right. The isopycnics are to be displayed with each isopycnic corresponding to an increase of 4%

of the initial density. The geometry and boundary conditions of this test case can be of the sound of the gas in state 1. Figure 3-21illustrates the schematic of moving shock waves relating the W, up, state 1 and 2.The W is important for that is relates the wave velocity of the moving shock wave to the pressure ratio across the wave and the speed of sound of the gas into which the wave is propagating, shown as follows (derived in [Anderson, 1990]):

the velocity Up behind the wave in state 2 is defined as:

1 velocity in the state 2 is set to u . Table 3-3 summarizes the temperature ratio, density

ratio and up initial conditions for the moving shock Mach number Ms, determined form [Anderson, 1990]. The simulation time t is depend on the location of the incident shock, which is approximately W × t. For example in this test case, the incident shock is 1.3L as shown in Figure 3-19 , the initial normal shock wave is set at the 0.3L where near from inlet boundary on the edge of the step, the simulation time is 1.464.

Figure 3-22 shows three density results that obtain using the second-order TVD extension of Godunov method [Takayama et al., 1991], N² and 2N method with 400 × 400 cells. The vortex obtained using both QDS method below the slipstream is close to the benchmark. The position of the vortex, incident shock and slipstream are perfectly in the correct please as the second-order TVD extension of Godunov method [Takayama et al., 1991].It is able to gauge the ability of the QDS method to detect shock, contact and expansion regions. To compare with experimental result made by Ritzerfeld et al. [Takayama et al., 1991], the schlieren result is easily to identify that the vortex obtained using N² method is clearer to observe as we can see from Figure 3-23.

According the results in this test case, the accuracy obtain using N² method is considerably better than that of the 2N method.

3.3.6 Advection of Vortical Disturbance

The final test case consists of an inviscid unsteady flow in which a vortex is located at the canter of a uniform domain (xc, yc). The mean flow for this case uses Mach number M^∞=0.1. The case tests the capabilities of the QDS method compared to the exact solution taken from [Visbal et al., 1999] in order to accurately advent vertical disturbance. This problem also appears in Tutkun and Bdis [Tutkun et al., 2010]. The initial conditions are shown as follows:

(63)

where u, υ and Rc determine the Cartesian velocity components and the vortex radius. C is the vortex strength parameter, defined as follows:

(64)

The density is assumed to be constant and the vortex radius Rc is taken to be 1.0 in this case.

Figure 3-24 shows the vorticity contours of the N² and the 2N methods with 800 × 800cells using the second-order method. The limiter in this case is the monotonized central (MC) method. A constant CFL number (0.1) is enforced such that the non-dimensional time step size is . The result of the N² solver is essentially the same as the exact solution and shows a perfect circular shape of the vorticity distribution while that of the 2N solver does not. The result of the 2N solver shows more significant dissipation and anisotropy errors as compared to that of the N² solver. Figure 3-25 shows the vorticity distributions of various simulations along a horizontal line (at y = 8.0) passing through the vortex center in Figure 3-24. We have compared the results obtained by using the two solvers (2N and N²) on a uniform grid containing three different levels of resolution (160 × 160, 800 × 800, and 1600 × 1600 cells). The result obtained using the N² solver in the case of 800 × 800 cells is in excellent agreement with the exact solution and radial symmetry, while the results obtained using the 2N solver are far from the correct solution even in the case of 1600 × 1600 cells. Thus, the influence of multi-dimensional reconstruction is significant for the

u= U_∞−C y

(

)

QDS, particularly on the numerical accuracy of the solution for a gas flow field as in the current problem.

The investigation of the computational expense again reveals a trade-off between computational time and accuracy. The computational time of the N² solver in terms of calculation time is approximately 3~4 times less than that of the 2N solver for any given computational grid although the accuracy of the former is considerably better than that of the latter. This leads to a question whether the use of the N² method is worthwhile or not. Thus, we compare the results obtained using the N²method using160 × 160 cells with those obtained using the 2N method using1600 × 1600 cells, as shown in Figure 3-25. The results show that they are essentially the same for the same level of accuracy;

thus, the proposed the N² solver is approximately 25 times faster than the 2N solver in this case. Once again, we are unsure whether the 2N result will ever converge to the analytical solution, thereby justifying the application of the N² solver and its proposed multi-dimensional reconstruction of QDS particles.