Chapter 4 Difference Analysis of QDS-2NMethod as Compared to QDS-
4.3 Results and Discussion
4.3.2 Example with Large Difference
4.3.2.4 Shock Wave Diffraction over a 90-degree Sharp Corner
Different from the preceding geometric configuration, the oblique shock does not occur in this test case instead of secondary physical phenomena such as a second shock, contact surface, and vortex. The last case discussed in this chapter involves the use of the same geometrical and boundary conditions as those considered in the case discussed in Section 3.3.5. The initial conditions of the flow for Ms are listed in Table 3-3.
According to the explanation of the structure of the perturbed region by Skews [1967], the location of the vortex is well defined for Ms < 1.5. Therefore, in order to observe more secondary physical phenomena in this test case, we discuss the case of shock wave diffraction by using the initial flow velocity of Ms = 2.4. The experimental result obtained by Schardun [Dyke, 1997] is selected as the benchmark. The flow structure around a perturbed region is outlined in Appendix A. Because the comparison has an experimental result as a benchmark, the resolutions of the simulation are calculated using a considerably fine cell: 1000 × 1000 uniform grids. The CFL number is set to 0.2.
The simulation time is calculated using the same principle as that discussed in Section 3.3.4 and is t = 0.775.
Figure 3-13 shows three schlieren results using the second-order 2N method, N2 method, and the experimental result separately. These results can be used for gauging the ability of these methods to detect the expansion region by juxtaposing them with the three results shown in Figure 3-18. The shape of the primary shock wave shown in Figure 3-13b and c matches the experimental result. The secondary shock waves obtained using the two methods are accurately located at the correct position behind the
wall and between the slipstream and the contact surface corresponding to the experimental result. The accuracy of the second shock wave in both the results is clear.
However, the phenomena in the vortex and the contact surface of both the results are contrary to the second shock. It is obvious that the vortex obtained using the N2 method is presented in considerable detail than that obtained using the 2N method. The contact wave is considerably diffused as compared to that in the case of the N2 method and the experimental result. The density, temperature, and velocity contours in this case are shown in Figures 4-14 to 4-16. These results show that the vortex belongs to the region with a low density, low temperature, and high velocity. The contact surface in the results is shown in the region of low density, low temperature, and low velocity. This is reasonable for supporting the theory in Section 4.3.1 that a more significant difference between the two methods shows the same trend as the large discrepancy is in the region of low density, low temperature, and low velocity.
4.4 Brief Summary
The major findings of the study of the difference analysis of the QDS-2N method presented in this chapter can be summarized as follows:
1. There are two regions of flow properties where a large discrepancy of conservative fluxes occurs between the QDS-2N and the QDS-N2 methods. The first is in the region of low density (down to 1.0), low temperature (down to 2.5), and high velocity (up to 2.0). The second is in the region of low density and low velocity (down to 1.2).
2. The conservative fluxes of the QDS-2N method that move along the diagonal direction exhibit a considerably large difference as compared to the QDS-N2
method. In contrast, the conservative fluxes that move along the horizontal or the vertical direction exhibit a significantly smaller discrepancy.
3. The normal shock and oblique shock waves in the resolution obtained by using the two methods are approximately located in the same region, and most of the shock waves, as predicted by the QDS-2N method, move in the diagonal direction.
4. Because the properties of an after-expansion wave are low density, low temperature, and high velocity, we have found a large discrepancy between the QDS-2N method and the QDS-N2 method.
Chapter 5
Conclusion and Recommendations of Future Work
5.1 Summary
In this thesis, a true-direction multi-dimensional higher-order extension of the QDS method, referred to as the QDS-N2 method, for solving the inviscid Euler equation is investigated numerically and theoretically. The major findings of this thesis are summarized in the following two sections in turn.
5.1.1 Numerical Investigation of QDS-N2 Method
1. The results of the one-dimensional shock and acoustic wave interaction problem demonstrate an improvement for higher orders of accuracy (up to third-order) of the QDS-N2 method.
2. The QDS-N2 method improves the solution in the flows unaligned with the computational grid as compared to the QDS-2Nmethod.
3. The QDS-N2 method significantly reduced the amount of numerical dissipation within the solution as compared to the QDS-2Nmethod.
4. Despite the additional computational expense associated with the QDS-N2 method for the same computational grid, for any given degree of accuracy, the proposed solver was found to be several times (up to 25 times in the case of the advection of vortical disturbances) faster than the original QDS-2N method.
5. Of particular interest is the test case of the advection of vortical disturbances, where the QDS-N2 method improves the radial symmetry of the result
approaching the analytical solution, while the QDS-2N method failed to converge to the analytical solution even when a very fine grid is used.
6. The results are essentially the same when N≥ 3 because the integration of the Gauss function with a polynomial (degree ≤ 2) using the Gauss-Hermite integration technique becomes exact.
7. Parallel performance studies, including strong and weak scaling, show that the parallel efficiency of shock bubble interaction for a large-scale problem (0.5, 2, and 12.5 million cells) can reach up to 75%, 68.5%, and 65.5% respectively using 256 processors at the APLS cluster of National Center for High-Performance Computing, Taiwan.
8. Parallel performance of weak scaling shows that the average efficiency of shock bubble interaction using 20,000 cells per processor is about 1.2, which the ideal efficiency is 1.0.
5.1.2 Theoretical Analyses of Conservation Fluxes of QDS-2N Method and QDS-N2 Method
1. There are two regions of flow properties where a large discrepancy of conservative fluxes occurs between the QDS-2N and the QDS-N2 methods. The first is in the region of low density (down to 1.0), low temperature (down to 2.5), and high velocity (up to 2.0). The second is in the region of low density, and low velocity (down to 1.2).
2. The conservative fluxe of the QDS-2N method that moves along the diagonal direction has considerably large difference as compared to the QDS-N2 method.
On the contrary, the conservative fluxes move along horizontal or vertical direction has a much smaller discrepancy.
3. The normal shock and oblique shock wave in the flow obtained using two methods are approximately located in the same region, even the direction of the shock wave, predicted by the QDS-2N method, moves to be more in diagonal direction.
4. Because the properties of an after expansion wave are in the region of low density, low temperature, and high velocities, we have found that a large discrepancy occurs between the QDS-2N method and the QDS-N2 method.
5.2 Recommendations of Future Work
In this thesis, we have demonstrated that the QDS-N2 method is a very fast numerical method without being subject to convergence problem like conventional CFD methods. However, there are several areas need to be done in pushing the method forward, The areas outlined below should be examined in the future:
• To further reduce the computational time for large-scale multi-dimensional problems, the method should be implemented on multiple extension graphics processing units (GPU).
• To further reduce the numerical diffusion, a high-order stencil in calculating the conservative fluxes may be considered.
• To implement an adaptive mesh function in the region where large gradient of flow properties occurs.
• To extend the QDS-N2 method for modelling the Navier-Stokes equation by employing the Chapmann-Engskog expansion theory to account for the viscous effect.
References
[1] Albright B.J. Don S. Lemons, Michael E. Jones, and Dan Winske., 2002, “Quiet direct simulation of Eulerian fluids,” Phys. Rev. E, Vol. 65, pp.1-4.
[2] Albright B.J, W. Daughton, D.S. Lemons, D. Winske, and M.E. Jones, 2002,
“Quiet direct simulation of plasmas,” Phys. Plasma, Vol. 9, Issue 5,pp. 1898.
[3] Anderson J. D.,1990, Modern Compressible Flow. McGraw-Hill, New York.
[4] Bazhenova T. V., Gvozdeva L. G., and Nettleton M. A., 1984, “Unsteady interactions of shock waves.” Prog. Aerospace Sci., Vol. 21, pp. 249–331.
[5] Bird G.A., 1994, Molecular Gas Dynamics and the Direct Simulation of Gas Flows, Clarendon Press, Oxford.
[6] Bird G.A., 1963,“Approach to translational equilibrium in a rigid sphere gas,”
Phys. Fluids, Vol. 6, pp. 1518-1519.
[7] Boyd I.D., Y. Jafry and J.V. Beukel, 1994, “Particle simulations of helium micro thruster flows,” J. Spacecraft Rockets, Vol. 31, Issue 2, pp. 271-277.
[8] Cave H.M., C.-W. Lim, M.C. Jermy, S.P. Krumdieck, M.R. Smith Y.-J. Lin, and J.-S. Wu, 2010, “Multi-Species Fluxes for the Parallel Quiet Direct Simulation (QDS) Method,” in: Proceedings of27th International Symposium on Rarefied Gas Dynamics.
[9] Čada M., Torrilhon M., 2009,“Compact third-order limiter functions for finite volume methods,” J. Compt. Phys.,Vol.228, pp. 4118-4145.
[10] Coronell D.G. and K.F. Jensen, 1992, “Analysis of Transition Regime Flows in Low Pressure Chemical Vapor Deposition Reactors Using the Direct Simulation Monte Carlo Method,” J. Electro chemical. Soc., Vol. 8, pp. 2264-2273.
[11] Chen W.F., 2003, “Study of hybrid DSMC/EPSM method,” Chin. J. Comput.
[12] Cook G., 1998, “High accuracy capture of curved shock fronts using the method of space-time conservation element and solution element,” in:
Proceedings of 37th AIAA Aerospace Sciences Meeting and Exhibit.
[13] Culbert B. Laney, 1998, “Computational Gasdynamics,” Cambridge University Press.
[14] Dyke M. van, 1997, An Album of Fluid Motion. Parabolic Press, Stanford.
[15] Gardiner C.W., 1985, “Handbook of Stochastic methods,” 2nd ed., Springer, New York, pp.106-107.
[16] Gombosi T.I., 1994, Gaskinetic Theory, Cambridge University Press, Cambridge.
[17] Godlewski, Edwige, 1996, “Numerical approximation of hyperbolic systems of conservation laws,” Springer, New York.
[18] Huang J.-C., Lin Herng, Yang J.-Y., 2009, “Implicit preconditioned WENO scheme for steady viscous flow computation,” J. Comp. Phys., Vol. 228, pp.
420-438.
[19] Hirsch C., 1990, “Numerical Computation of Internal and External Flows, Vol.
2,” John Wiley & Sons, New York.
[20] Jameson A., W. Schmidt and E. Turkel, 1981, “Numerical solution of the Euler equations by finite volume methods using Runge-Kutta time-stepping schemes,”
the 14th Aerospace Sciences Meeting, AIAA-81-12S9.
[21] Jameson A., T. J. Baker and N. P. Weatherill, 1986, “Calculation of inviscid transonic flow over a complete aircraft,” AIAA24th Aerospace Sciences Meeting, AIAA -86-0103.
[22] John D. Anderson Jr., 1995, Computational Fluid Dynamics.
[23] Jiang G., Shu C.-W., 1996, “Efficient implementation of weighted ENO schemes,” J. Comp. Phys., Vol. 126,pp. 202–228.
[24] Keats W.A., Lien F.S., 2004, “Two-dimensional anisotropic Cartesian mesh adaptation for the Euler Equations,” Int. J. Numer. Meth. Eng., Vol. 46, pp. 1099.
[25] Kröner D, 1997, “Numerical Schemes for Conservation Laws.” Wiley Teubner.
[26] Le Beau G. J. and F. E. Lumpkin III, 2001, “Application highlights of the DSMC Analysis Code (DAC) software for simulating rarefied flows,” Comp.
Meth. Appl. Mech. Eng. Vol. 191, Issue 6, pp. 595-609.
[27] Lim Chin Wai, Hadley M. Cave, Mark C. Jermy, Susan P. Krumdieck and Jong-Shinn Wu, 2010, “An Approximate Method for Solving Unsteady Transitional and Rarefied Flow Regimes in Pulsed Pressure Chemical Vapor Deposition Process using the Quiet Direct Simulation Method,” in: Proceedings of 27thInternational Symposium on Rarefied Gas Dynamics, California.
[28] LeVeque R. J., 2004, “The dynamics of pressureless dust clouds and delta waves,” J. Hyperbolic Differential Equations, Vol 1, pp. 315-327.
[29] Karniadakis G.E. and A. Beskok, 2002, “Micro flows: fundamentals and simulation” Springer, New York.
[30] Kwon M.-K. and Y.-K. Hwang, 2006, “A Study on the Pumping Performance of a MolecularDrag Pump in the Rarefied Gas Flow,” in: Proceedings of 25th International Symposium on Rarefied Gas Dynamics, St. Petersburg, Russia.
[31] MacCormack R. W. and B. S. Baldwin, 1975, “A numerical method for solving the Navier-Stokes equations with application to shock boundary layer interaction”, the13th Aerospace Sciences Meeting, AIAA-7s-1.
[32] Macrossan M.N., 1989, “The Equilibrium Flux Method for the calculation of flows withnon-equilibrium chemical reactions,” J. Comp. Phys. Vol. 80, Issue 1, pp. 204-231.
[33] Macrossan M.N., 2001, “A particle only hybrid method for near continuum flows,” in: Proceedings of 22nd International Symposium on Rarefied Gas Dynamics, Vol. 585, pp. 426.
[34] Majda A., 1984, “Compressible fluid flow and systems of conservation laws in several space variables,” Springer-Verlag, New York.
[35] Pullin D.I., 1980, “Direct simulation methods for compressible inviscid ideal gas flow,” J. Comp. Phys. Vol. 34, pp. 231-244.
[36] Salichs S. Serna, 2005, High Order Accurate Shock Capturing Schemes for Hyperbolic Conservation Laws based on a new Class of Limiters, Ph.D. Thesis, Doctoral Theses in Applied Mathematics, Universitat de Valencia.
[37] Schulz-Rinne C.W., Collins J.P., and H.M. Glaz, 1993, ‘Numerical Solution of the Riemann problem for two-dimensional gas dynamics,” SIAM J. Sci. Compt., Vol. 14, pp. 1394–1414.
[38] Serna S., 2006, “A class of extended limiters applied to piecewise hyperbolic methods,“ SIAM J. Sci. Compt., Vol. 28, pp. 123–140.
[39] Shu C.-W., Osher S., 1988, “Efficient implementation of essentially non-oscillatory shock-capturing schemes,” J. Comp. Phys., Vol. 77, pp. 439–471.
[40] Skews B.W., 1967, “The perturbed region behind a diffracting shock wave.” J.
Fluid Mech., Vol. 29, Issue 4, pp.705–719.
[41] Smith M.R., 2003, “Hybrid methods in near-continuum flows,” M.Phil. Thesis, University of Queensland, Australia.
[42] Smith M.R., M. N. Macrossan, M. M. Abdel-Jawad, 2008, “Effects of direction decoupling in flux calculation in finite volume solvers,” J. Comp. Phys. Vol. 227, pp. 4142-4161.
[43] Smith M.R., 2008, “The true direction equilibrium flux method and its application”, PhD Thesis, University of Queensland, Australia.
[44] Smith M.R., H.M. Cave , J.-S. Wu , M.C. Jermy , Y.-S. Chen, 2009, “An improved Quiet Direct Simulation method for Eulerian fluids using a second-order scheme,” J. Comp. Phys. Vol. 228, pp. 2213–2224.
[45] Smoller, Joel, 1983, “Shock waves and reaction-diffusion equations,” Springer-Verlag, New York.
[46] Sun M. and K. Takayama., 1997, “The formation of a secondary shock wave behind a shock wave diffracting at a convex corner.” Shock Waves, Vol. 7, pp.
287–295.
[47] Takayama K. and O. Inoue., 1991, “Shock wave diffraction over a 90 degree sharp corner.” Shock Waves, Vol. 1, pp. 301–312.
[48] Teshima K. and M. Usami, 2001, “DSMC simulation of axisymmetric supersonic free jets,” Comp. Fluid Dyne. J., Vol. 10, pp. 525-530.
[49] Toro E.F., 1999, Riemann Solvers and Numerical Methods for Fluid Dynamics, second ed., Springer, Berlin, Germany.
[50] Toro E. F., 1991, “A linearized Riemann solver for the time-dependent Euler equations of gas dynamics,” Proc. Roy. Soc. London Ser. A, 434, pp. 683-693.
[51] Toro E. F., 2009, “Riemann solvers and numerical methods for fluid dynamics:
a practical introduction.” Springer-Verlag, New York.
[52] Tutkun, B., Edis, F., O., 2010, “A GPU Application for High-Order Compact Finite Difference Scheme,” 22nd Int. Conference on Parallel CFD 2010, Kaohsiung, Taiwan.
[53] Van Leer B., 1977, "Towards the ultimate conservative difference scheme III.
Upstream-centered finite-difference schemes for ideal compressible flow", J.
Comp. Phys., Vol. 23, Issue 3, pp. 263–275.
[54] Visbal Miguel R., Gaitonde Datta V., 1999, “High-Order Accurate Methods for Complex Unsteady Subsonic Flows,” AIAA Journal, Vol. 37, pp. 10.
[55] Wagner W., 1992, “A convergence proof for Bird's direct simulation Monte Carlo method for the Boltzmann equation,” J. Stat. Phys., Vol. 66, pp. 1011-1044.
[56] Woodward P., Colella P., 1984, ”The numerical simulation of two-dimensional fluid flow with strong shock,” J. Comp. Phys., Vol. 54, pp.115.
[57] Wu M.G., 2003, “Study on hybrid EPSM/EPSM method for chemical reacting flow,” Chin. J. Comp. Phys., Vol. 20, Issue 3, pp. 564.
Appendix A
Flow Structures of the Shock Wave Diffraction
Keats [Keats et al., 2004] describes the theory of the shock wave diffraction that including experimental and computational result from Skews [1967], and summarized the secondary physical phenomena in the perturbed region behind the shock wave.
Skews performed experiments for a variety of Mach numbers and convex corner angles, and has outline the structure of the perturbed region; the structure is shown in Figure 3-5.1. Skews determined experimentally and tabulated the following correlations:
• The slipstream angle variation with the shock Mach number Ms.
• The terminator angle variation with Ms.
• The relationship between Ms and the velocity of the secondary shock.
• The contact surface velocity variation with Ms.
• The variation of the vortex angle and velocity with Ms. The flow structures can be described as follows:
• Incident shock: Diffracts in a similar way to a sound wave: its radius of curvature is approximately Wt.
• Reflected sound wave: Propagates upstream and marks the start of the curvature of the incident shock.
• Slipstream: Due to separation, it separates high-velocity gas on the upper side from almost stationary gas on the lower side. It represents the outermost characteristic of the Prandtl-Meyer expansion fan.
• Terminator: The first characteristic of the Prandtl-Meyer expansion; the single separating the terminator from the horizontal increases with rising Ms.
• Second shock: The region between the slipstream and the terminator is a uniform flow region parallel to the slipstream, and the second shock is a normal shock caused when the flow in this region exceed Mach 1.0 [Sun et al., 1997].
• Vortex: Located just below the slipstream, its location is well defined for Ms≤ 1.5. The angle between the vortex and the slipstream decreases as Ms increases.
• Contact surface: Originates at the intersection point of the reflected sound wave and the incident shock, but is highly diffuse in this region. It becomes better-defined as it nears the region containing the rest of the flow structures.
Tables
Table 2-1 The value of weight and abscissas for the Gaussian quadrature.
Number of QDS particles Weight (wJ) Abscissas (qJ)
2 1
2 2
± 1
2 π
3 0 2
3 π 1 6
±2 1
6 π
4 1 3 6
2 2
± − 4 3
(
−π 6)
1 3 6
2 2
± + 4 3
(
+π 6)
Table 3-1 Comparison of computational expenses for QDS schemes using 2N and N2 dimensional reconstruction.
Number of cells
QDS solvers
2N N2
300x100 8.41 min 23.15 min
450x150 28.9 min 78.3 min
600x200 68.56 min 183.6 min
1000x500 478.4 min 1282.6 min
Table 3-2 QDS scheme time cost in Euler-4-shocks interaction case.
Number of cells QDS solvers
2N N2
1000x1000 13.29 hours 55.6 hours
100x100 45 (s) 189 (s)
200x200 375(s) 1520 (s)
300x300 1307 (s) 5199 (s)
Table 3-3 pre- and post-shock fluid initial conditions.
Mach number
1.5 2.4
ρ2/ρ
1 1.862 3.212
T2/T1 1.32 2.04
Up 0.8215 1.956
Table 3-4 Parallel Performance for a 2D shock-bubble problem with 2.4 millioncomputational cells.
Number of processors Computation time (sec.) Number of cells
1 361.613 20000
4 348.995 80,000
9 428.812 180,000
16 410.294 320,000
25 425.229 500,000
36 444.135 720,000
49 472.599 980,000
Table 3-5 Parallel computation times for shock-bubble problem with 500,000 cells at 2000 time steps in simulation time 0.2.
Number of processors Computation time (s) Number of Cells for one processor
1 8220.19 500,000
2 4299.89 250,000
4 2761.14 125,000
8 1210.08 62,500
16 762.245 31,250
25 425.229 20,000
32 347.833 15,625
64 158.845 7,813
80 132.77 6,250
100 106.164 5,000
256 42.9935 1,953
Table 3-6 Parallel computation times for shock-bubble problem with 2 million cells in 2000 time steps.
Number of processors Computation time (sec.)
Number of cells for one processor
1 29740.59 2,000,000
4 8744.45 500,000
8 4978.56 250,000
16 2890.74 125,000
25 1827.92 80,000
32 1348.79 62,500
64 679.745 31,250
128 340.572 15,625
256 172.599 7,813
Table 3-7 Parallel computation times for shock-bubble problem with 12.5 million cells in 2000 time steps.
Number of processors Computation time (sec.)
Number of cells for one processor
1 187130.1 12,500,000
8 27532.8 1,562,500
16 19929.3 781,250
32 8613.3 390,625
64 4463.82 195,313
100 2806.29 125,000
128 2225.91 97,656
256 1115.02 48,828
Figures
Figure 2-1. Schematic showing the way fluxes of conserved quantities between source and destination cells are calculated using the “overlap” function in QDS [Smith et al.,
2009].
Destination cell
Source cell (Area, AS)
vyk∆t A
vxj∆t
Overlap area, A= vxjvyk∆t2 Mass transported from source to destination cell is proportional to A/AS
Figure 2-2.Flowchart describing QDS-2N particle computation with gradient inclusion.
Figure 2-3. The special reconstruction convention for current amount of conserved
quantity Q in one cell.
Figure 2-4. QDS flux procedure within a general (arbitrary) spatial reconstruction of conserved quantity Q.
Figure 2-5. Flowchart describing QDS particle computation with gradient inclusion.
Figure 2-6. Two-dimensional motion of a single QDS particle showing “sub-particle”
contributions.
Figure 2-7. Three-dimensional motion of a single QDS particle showing “sub-particle”
contributions.The green parallelogram is presented the concept in two-dimensional.
DESTINATION A
DESTINATION B
DESTINATION C
DESTINATION D
DESTINATION E
DESTINATION F
DESTINATION G
SOURCE
Figure 3-1. Waves generated in shock tube following rupture of diaphragm [Anderson, 1990].
Contact surface moving at the velocity of the gas behind the shock when diaphragm ruptured.
Normal shock wave propagating
Expansion wave propagating
outflow outflow
Figure 3-2. The shock tube problem as computed by pre-QDS method with a uniform grid of 200 zones. The results were discussed thedifference to the QDS 1st ~3rdmethod
and Riemann solver using MINMOD limiter at time 0.1.
Figure 3-3. The interaction of two blast wave computed by the QDS method with 400 grids at t = 0.0038. The solid black line is WENO (fifth order) scheme with 10,000
grids.
Figure 3-4. Density profile of the shock-acoustic-wave case at t = 1.8. The solid black line is WENO-3 (fifth order) with 2000 grids compared with QDS method which
without limiter form 1st order to 3rdorder.
0.5
Figure 3-5. The structure of shock bubble interaction.
(a) (b)
(c)
Figure 3-6. Zoom of shock-bubble Schlieren image with 1000×500 cells at time of 0.2.
QDS 2nd order (a) 2N method with van Leer’s limiter, (b) N2 method, and (c) 2nd order TVD result presented in [Čada, 2009] using the same resolution.
Figure 3-7. Zoom of Schlieren image of shock bubble problem at time of 0.2; (a) QDS-N2 method with 300×100 cells; (b) QDS-2N method with 300x100 cells; (c) QDS-2N
method with 450×150 cells; (d) QDS-2N scheme with 600×200 cells.
Figure 3-8. The initial conditions for the first problem of Euler-4 shocks interaction.
(a) (b)
(c) (d)
Figure 3-9. Zoom of density contour line of Euler-four-shocks problem. Comparing the second-order QDS-N2 method (a) using 100×100 grids with MC limiter and 2N method
Figure 3-9. Zoom of density contour line of Euler-four-shocks problem. Comparing the second-order QDS-N2 method (a) using 100×100 grids with MC limiter and 2N method