Numerical Analysis of Transonic Flow around Cones

Rozie Zangeneh

Department of Mechanical Engineering, Lawrence Technological University, Southfield, MI, USA.

DOI: 10.4236/ojfd.2020.104017                

Abstract

A new solver is presented for transonic flow around cone-cylinder, axisymmetric bodies. Ground experiments almost always suffer from uncertainty due to operating in the presence of high levels of facility noise. Besides, experimental measurements of these mechanisms are not available at high-speed flows. Direct Numerical Simulations have made it possible to compute details of the transonic mechanisms but still a significant challenge due to the cost. This study aims to present a new solver to model transonic flows. To assess the new solver, the surface Mach number and the drag coefficient are investigated as the freestream Mach number varies. The results are in excellent agreement with experimental data, indicating the new model is capable of accurately predicting the aerodynamics coefficients at transonic flow regimes.

Keywords

Axi-Symmetric Bodies, Transonic Flow, Surface Mach Number

1. Introduction

Transonic flow past certain two-dimensional bodies has been studied extensively, and therefore, the phenomena are well understood. Some of the earliest theoretical studies are done by Cole [1], Guderley [2], and Vincenti [3], which applied to two-dimensional wedge airfoils. The theory and experimental results conducted by Bryson [4] and Griffith [5] agreed well. Two-dimensional and axially symmetric bodies are of considerable theoretical and practical interest to study since these two cases are simplified cases of the problem around complex and arbitrary geometries.

The study of transonic flow around axisymmetrical bodies is not as complete as two-dimensional bodies. The similarity physicals of axially symmetric transonic flow studied by Von-Karman [6] and Oswatitsch [7] discuss general transonic flow past finite cones. The first theoretical results for supersonic flow past a cone were presented in 1929 by Busemann [8]. The Busemann’s solution predicts the smooth shock-free compression from supersonic to subsonic flow for specific combinations of cone angle and freestream Mach number. The conical solution also showed that for a given freestream Mach number and cone angle, the surface Mach number is less than the freestream Mach number right behind the conical shock wave; as the freestream Mach number decreases, the surface Mach number also decreases and finally changes from supersonic to subsonic values. It was concluded that the conical solution for the semi-infinite cone is valid for a finite cone only when the flow is supersonic everywhere. However, when the surface Mach number becomes subsonic, the perturbation due to the corner or the cone’s shoulder propagates forward through the subsonic field and therefore interrupts the conicity of the flow. Thereby, the conical solution applies only for large enough freestream Mach numbers so that the freestream Mach number is supersonic.

For the first time, Taylor and Maccoll [9] presented the numerical solution of the axisymmetric conical flow around semi-infinite cones and validated their experimental data results. The experiments and theoretical results showed notable discrepancies, especially in the shockwave form when the surface Mach number is subsonic. Yoshihara [10] computed the flow around a cone cylinder at the sonic Mach number, which was verified experimentally. The theoretical solution, however, has not been developed for transonic flow around finite cones. Solomon [11] reported the experiment results on several interesting characteristics of the transonic flow around finite cones. The experiment evaluated the deviation of Mach number from the predicted value of conical theory for the transonic freestream Mach number, which leads to an evaluation of drag coefficients.

In the design of aerodynamic vehicles such as missiles, rockets, space shuttles, etc., various shapes are used to reduce the aerodynamic drag to achieve the best performance. Investigating the different parameters on the flow/shock characteristics, such as the shape of the shock wave near the nose, shock detachment distance, and the local Mach number is of special importance in order to determine the parameters that provide minimum aerodynamic drag since drag reduction is essential for the better performance of the aerodynamic vehicles. The formation of the bow shock in the vicinity of the nose, shock detachment distance, shock layer, flow turning angle, etc. plays a substantial role in modifying the aerodynamic characteristic to achieve better cones’ performance.

Despite several studies that have been done on flow past different cones numbers, detailed numerical investigation of transonic and supersonic flow with attached shock wave past a circular cone at zero incidences is of fundamental importance both from a theoretical point of view and for practical applications. The present work deploys a new numerical solver for the compressible Navier-Stokes equations to model the transonic flow past cone cylinders. The model is based on a central differencing scheme developed by Kennedy and Gruber [12] for complex geometries. Unlike the existing methods, it does not involve Riemann solvers or characteristic decomposition; therefore, suitable for complex geometries. This unique is particularly essential to achieving an approach to incorporate in standard industrial solvers.

2. Computational Methodology

2.1. A Hybrid-Energetic Numerical Model

To address the conflict between turbulence modeling and shock capturing schemes, a hybrid algorithm is a natural solution, where a dissipative scheme can capture discontinuities at shocks, and the non-dissipative scheme resolves the small scales at the turbulent region. Here an energy conserving scheme is achieved in a finite volume frame-work through hybridization of convective terms, which is evaluated as the inviscid flux, within the Navier-Stokes equations as follows:

${\varnothing }_{\text{Hybrid-Energetic}}=\left(1-\in \right){\varnothing }_{\text{Central-nondissipative}}+\in {\varnothing }_{\text{Shock-capturing}}$ (1)

The non-dissipative component is based on an improved skew-symmetric formulation developed by Kennedy and Gruber [12] for convection terms. In this respect, for the quadratically nonlinear terms, the following single generalized expression is used:

$\nabla \left(u\vartheta \right)=\alpha \nabla \left(u\vartheta \right)+\left(1-\alpha \right)\left[\vartheta \nabla u+u\nabla \vartheta \right]$ (2)

where u is the velocity vector and $\vartheta$ denotes a generic transport variable. The skew-symmetric form is obtained by setting $\alpha =1/2$ to minimize aliasing errors during numerical simulation within the Navier-Stokes equations. In order to reduce to a minimum unphysical energy transfer and growth at the high wavenumbers, convective formulations with the $\alpha$ set to values near $\alpha =1/2$ should be employed. This form is proved to minimize the aliasing errors that are associated with low-order non-dissipative schemes. Here, a modified form is used to compute the convection terms as:

${\varnothing }_{\text{Central-nondissipative}}=\frac{u_P+u_N}{2}\frac{{\vartheta }_P+{\vartheta }_N}{2}$ (3)

where P and N are cell centroids (Figure 1). The cubically non-linear terms can be discretized in the same manner. It is noted that using second-order methods are common within applications in complex geometries, even for LES applications [13].

A key role in shock-capturing schemes is played by “shock sensors,” that must be defined to confine numerical dissipation in shocked regions effectively and, at the same time, not to effect smooth parts of the flow field [14]. The pertinence of

the choice of a sensor based on the pressure gradient, ${\delta }_i=\frac{\left|P_{i+1}-2P_i+P_{i-1}\right|}{\left|P_{i+1}+2P_i+P_{i-1}\right|}$, to

capture shock discontinuities usually found in aerodynamics is discussed by Swanson and Turkle [15]. Ducros et al. [13] developed a new correction to the

Figure 1. Computational domain with structured mesh. The gridding is coarsened for visual clarity.

sensor by multiplying the standard sensor by the local function $\kappa$, which is defined as the following:

${\kappa }_i=\frac{{\left(\nabla \cdot u_i\right)}^2}{{\left(\nabla \cdot u_i\right)}^2+{\left(\nabla \times u_i\right)}^2+\epsilon }$ (4)

where $\epsilon =1\text{e}-23$ is a small positive real number chosen to prevent numerical divergence in regions where $\nabla \cdot u_i^2+\nabla \times u_i^2$ is zero. The sensor is able to evaluate the smoothness if the numerical solution by $\kappa =0$ in the smooth zones, and $\kappa =1$ in the presence of shocks [13]. This function alters between 0 for weakly compressible regions to about 1 in shock regions. Here the following modified shock sensor is used as follows:

$\in =\max \left({\kappa }_i\times {\delta }_i,{\kappa }_{i+1}\times {\delta }_{i+1}\right)$ (5)

The modified sensor shows smooth correction, proportional to the degree of local compressibility and is proved to predict the right decay of turbulence kinetic energy in turbulent regions out of the shock [13].

2.2. Central-Upwind Scheme for Compressible Flows

The ${\varnothing }_{\text{Shock-capturing}}$ component in the hybrid algorithm in Equation (1) is based on a model introduced by Kurganov et al. [16] [17]. The main advantage of these central schemes is the high resolution, due to the minimal numerical dissipation, and the simple application. Because there are no Riemann solvers and characteristic decomposition involved, and this makes it an appropriate model for various types of applications. The shock capturing component is based on the use of accurate information about the local speeds of propagation. Figure 1 shows the grids of polyhedral cells with an arbitrary number of faces and vertices.

An owner cell and a neighbor cell were assigned for each face. $S_f$ is a vector normal to the face surface pointing out of the owner cell. The magnitude of the vector is the area of the surface. In such a collocated system, all properties and dependent variables are located at each cell centroid, p. The discretization of a general dependent tensor field $\Psi$ is described by interpolation of values ${\Psi }_p$ at cell centroid to values ${\Psi }_f$ at cell faces. The Integral of convective term over a control volume is linearized as [18]:

${\displaystyle {\int }_V\nabla \cdot \left[u\Psi \right]\text{d}V}={\displaystyle {\int }_S\text{d}S\cdot \left[u\Psi \right]}\approx {\displaystyle {\sum }_f{S}_f\cdot u_f{\Psi }_f}={\displaystyle {\sum }_f{\varphi }_f{\Psi }_f}$ (6)

where ${\varphi }_f$ is volumetric fluxes. For compressible flows, however, fluid properties are transported by the propagation of waves in addition to transport by the flow. Therefore, the flux interpolation should be stabilized based on transports in any direction [19]. Since the interpolation is done to a given face only from neighboring cell values, the original form of Kurganov and Tadmor (KT) and Kurganov, Noelle, and Petrova (KNP) methods are used [19]. The interpolation procedure is based on splitting into two directions corresponding to flow outward and inward of the face owner cell. The discretization is as follows:

${\displaystyle {\sum }_f{\varphi }_f{\Psi }_f}={\displaystyle {\sum }_f\alpha {\varphi }_{f+}{\Psi }_{f+}}+\left(1-\alpha \right){\varphi }_{f-}{\Psi }_{f-}+{\omega }_f\left({\Psi }_{f-}-{\Psi }_{f+}\right)$ (7)

The last term in the above equation is volumetric fluxes associated with the local speed of propagation. It is an extra diffusion term related to the maximum speed of propagation of any discontinuity that may occur at a face. Therefore, the value is interpolated in the $f-$ and $f+$ directions. The diffusive volumetric flux is calculated according to:

${\omega }_f=\{\begin{array}{l}\alpha \max \left({\phi }_{f-},{\phi }_{f+}\right)\text{forKTmethod}\\ \alpha \left(1-\alpha \right)\left({\phi }_{f-}+{\phi }_{f+}\right)\text{forKNPmethod}\end{array}$ (8)

where ${\phi }_{f-}$ and ${\phi }_{f+}$ are defined as:

$\{\begin{array}{l}{\phi }_{f+}=\max \left(C_{f+}\left|S_{f+}\right|+{\varphi }_{f+},C_{f-}\left|S_{f-}\right|+{\varphi }_{f-},0\right)\\ {\phi }_{f-}=\max \left(C_{f+}\left|S_{f+}\right|-{\varphi }_{f+},C_{f-}\left|S_{f-}\right|-{\varphi }_{f-},0\right)\end{array}$ (9)

Here $C_f=\sqrt{\gamma R{T}_f}$ is the speed of sound of the gas at the face, with respect to outward and inward of the owner cell. A flux limiter function, $\beta \left(r\right)$, is used in the interpolation procedure to switch between low and high-order schemes where r represents the ratio of successive gradients of the interpolated variable. It can be described according to:

$r=2\frac{d\cdot {\left(\nabla \Psi \right)}_p}{{\left({\nabla }_d\Psi \right)}_f}-1$ (10)

where ${\left(\nabla \Psi \right)}_p$ is the full gradient calculated at the owner cell, P, and ${\left({\nabla }_d\Psi \right)}_f={\Psi }_N-{\Psi }_P$. Then, the $f+$ interpolation of $\Psi$, for example, is evaluated as:

${\phi }_{f+}=\left(1-g_{f+}\right){\Psi }_P+g_{f+}{\Psi }_N$ (11)

where $g_{f+}=\beta \left(1-{\omega }_f\right)$, and the symmetric TVD limiter is $\beta \left(r\right)=\frac{r+\left|r\right|}{1+r}$, from

van Leer et al. [20]. The resolution of this semi-discrete central-upwind scheme can be further improved, especially of the contact waves, by adding a correction term, $q_f$, which is anti-diffusion to numerical fluxes in Equation (7) as:

${\displaystyle {\sum }_f{\varphi }_f{\Psi }_f}={\displaystyle {\sum }_f\alpha {\varphi }_{f+}{\Psi }_{f+}}+\left(1-\alpha \right){\varphi }_{f-}{\Psi }_{f-}+{\omega }_f\left({\Psi }_{f-}-{\Psi }_{f+}-q_f\right)$ (12)

The “correction” term, $q_f$ is, in fact, a built-in anti-diffusion term, can be computed as [21]:

$q_f=\gamma \min \mathrm{mod}\left({\Psi }_{f+}-w_f,w_f-{\Psi }_{f-}\right),\gamma \in \left(0,1\right)$ (13)

where:

$w_f=\alpha {\Psi }_{f+}+\left(1-\alpha \right){\Psi }_{f-}-\left\{{\varphi }_{f+}{\Psi }_{f+}-{\varphi }_{f-}{\Psi }_{f-}\right\}$

These terms help to reduce the numerical dissipation present in the original form of the semi-discrete central-upwind scheme. Finally, the Laplacian term on the polyhedral cell is discretized by splitting into orthogonal and non-orthogonal components as following:

${\displaystyle {\int }_S\text{d}S\cdot \left[\text{г}\nabla \Psi \right]}\approx {\displaystyle {\sum }_f\text{г}{S}_f\cdot {\left(\nabla \Psi \right)}_f}={\displaystyle {\sum }_f\text{г}\left\{A\left({\Psi }_N-{\Psi }_P\right)+a\cdot {\left(\nabla \Psi \right)}_f\right\}}$ (14)

where $\text{г}$ is the diffusion coefficient which is interpolated linearly from the cell center values, $A={\left|S_f\right|}^2/\left(S_f\cdot d\right)$ and $a=S_f-Ad$.

2.3. Implementation in OpenFOAM

The new framework is developed using C++ code and linked to the existing density-based library (rhoCentralFoam) for compressible flows in the open-source toolbox, OpenFOAM [22], as Hybrid-Energetic algorithm. For time advancement, the original Euler implicit scheme is preserved. However, a higher-order Runge-Kutta scheme can be employed in the future for more accurate temporal discretization.

rhoCentalFoam comprises a dissipative model based on original central-upwind KT schemes, which initially is developed for laminar compressible flow problems and hence ineffectual for resolving turbulence. Therefore, the new hybrid model can potentially be a practical implementation for modeling interactions of shocks and turbulence for OpenFOAM users. The model is already compared with the original solver [23] for transonic flows and validated on numerous cases [24] [25].

2.4. Computational Domain

The flow around the axisymmetrical cone-cylinder is investigated with semi-angles of 20 and 25 deg. The base diameter is 20 mm. And the total length of the model is 50 mm. For the discretization of the computational domain, structural grids are employed. The grid size uniformly expands in the wall-normal direction in order to cluster grid points near the wall; therefore, the grid size changes from $y^{+}=1.3$ near the wall to $y^{+}=5$ at the top of the domain. The transonic Mach number, $M_{\infty }=1.4$, changes from 0.6 to values close to 1.4, similar to the experimental setup [11]. The Reynolds number varies from $Re=50000$ to $Re=500000$. The flow is assumed to be Newtonian and the specific heat ratio and Prandtl number are $\gamma =1.4$ and $Pr=0.72$, respectively. The no-slip and adiabatic conditions are enforced on the surface of the cones. The origin of the coordinate system is the leading edge of the cone. As shown in Figure 1. Only the axisymmetry model about z-axis was simulated due to symmetry.

3. Results

Figure 2 shows supersonic flow at infinity, $M_{\infty }=1.4$, past the cylindrical cone for the semi-angle cone of 25 deg with an attached curved shockwave and subsonic flow between the initial portion of the shock wave and the cone.

3.1. Local Mach Number Contours

Figure 3 (Left) shows the example of the sonic-line location when the flow at the infinity is supersonic with a detached shockwave and subsonic flow between the cone and the shockwave. It is apparent that the sonic-line location originates slightly upstream the cone shoulder which is due to the effect of boundary-layer rounding the cone shoulder. Figure 3 (Middle) illustrates the case with nearly attached curved shockwave with a small supersonic to subsonic compression region on the front section of the sonic-line. The sonic-line initiates at the cone shoulder before it terminates on the shockwave. The freestream Mach number, $M_{\infty }=1.273$, is slightly less than the detachment Mach number predicted by the exact conical theory [22]. Figure 3 (Right) depicts the sonic-line location for the attached shockwave when the flow at the infinity is supersonic, $M_{\infty }=1.41$, a mixed supersonic and subsonic flow regime is present between the cone and the shockwave. The sonic-line again originates at the corner and now terminates at the cone tip and not on the shock wave, unlike the two previews cases. In this case, a shock-free supersonic to subsonic compression happens at the front section of the sonic-line.

Figure 2. The 25 deg semi-angle cone at $M_{\infty }=1.41$.

Figure 3. The location of the sonic-line for 25 deg semi-angle cone cylinder at $M_{\infty }=1.229$ (Left), the location of the sonic-line for 25 deg semi-angle cone cylinder at $M_{\infty }=1.273$ (middle), the location of the sonic-line for 25 deg semi-angle cone cylinder at $M_{\infty }=1.41$ (right).

Figure 4 shows the sonic-line for a 20 deg semi-angle of the cone. The problem of the smooth shock-free supersonic to subsonic compression is subjected to discussion for many decays. The numerical results prove that such a smooth flow transition occurs on the transonic cone. The condition is a sonic-line bounds the zone of subsonic flow enclosed by supersonic flow. It is noted that the transonic flow on an airfoil is an example of the non-shock free supersonic to subsonic compression. A smooth compression through sonic condition does not occur since the shockwave terminates the local supersonic flow on the airfoil.

3.2. Surface Mach Number Distribution

The distribution of the surface Mach number, on a 20 deg cone versus various values of the freestream Mach number at different stations, x/c, where c is the length of the model, is shown in Figure 5 and for the 25 deg cone in Figure 6. The values of Surface Mach number close to the tip agrees well with the experimental values (Solomon, 1955 [11] ). At the corner of the cone, the Mach number should approach sonic, although the boundary layer can affect that. As the freestream Mach number reaches the sonic condition, the surface Mach number approaches a constant value. The behavior is more apparent for the cone angle of 25 deg. This concept is proven from conical theory and is established for axisymmetric flows. The stationary concept can be generalized to three-dimensional bodies for some selective ranges of near sonic velocity.

Figure 4. Sonic-line location for 20 deg semi-angle cone cylinder at $M_{\infty }=1.247$ (Left), Sonic-line location for 20 deg semi-angle cone cylinder at $M_{\infty }=1.293$ (right).

Figure 5. Surface Mach number on 20 deg semi-angle cone.

Figure 6. Surface Mach number on 25 deg semi-angle cone.

Figure 7. The drag coefficient for cone cylinders.

3.3. Drag Coefficient

Figure 7 shows the variation of the drag force coefficient of the cone cylinder with the freestream Mach number. The slope of the drag coefficient agrees with the experimental value (Solomon, 1955 [11] ). The drag coefficient’s numerical values in the transonic region are well predicted, especially for the semi-angle cone of 25 deg.

4. Conclusion

Numerical results are presented for transonic flow over axisymmetric cone cylinder. The smooth transition from the supersonic to subsonic compression is determined by the numerical results in agreement with the experimental observation. The surface Mach number and the drag coefficient are evaluated as the freestream Mach number is varied, and the results are compared with the experimental data. The variation of the drag force with freestream Mach number is calculated from the numerical values. Results show as the freestream Mach number reaches the sonic, the surface Mach number approaches a constant value. The stationary concept can be generalized to three-dimensional bodies for some selective ranges of near sonic. The results are in excellent agreement with experimental data, indicating the new model is capable of accurately predicting the aerodynamics coefficients at transonic flow regimes. The model can be used to investigate different parameters on the flow/shock characteristics, such as the shape of the shock wave near the nose, shock detachment distance, and the local Mach number to determine the parameters that provide minimum aerodynamic drag since drag reduction is essential for the better performance of the aerodynamic vehicles.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Cole, J.D. (1951) Drag of a Finite Wedge at High Subsonic Speeds. Journal of Mathematics and Physics, 30, 79-93. https://doi.org/10.1002/sapm195130179
[2]	Yoshihara, H. and Guderley, G. (1950) The Flow over a Wedge Prole at Mach Number 1. Journal of the Aeronautical Sciences, 17, 723-735. https://doi.org/10.2514/8.1782
[3]	Wagoner, C.B. and Vincenti, W.G. (1952) Transonic Flow past a Wedge Cone with Detached Bow Wave. NACA-TR1095.
[4]	Bryson, A. (1952) An Experimental Investigation of Transonic Flow past Two-Dimensional Wedge and Circular-Arc Sections Using a Mach-Zehnder Interferometer. NACA-TR1094.
[5]	Griffith, W. (1952) Shock-Tube Studies of Transonic Flow over Wedge Cones. Journal of the Aeronautical Sciences, 19, 249-257. https://doi.org/10.2514/8.2242
[6]	Von Karman, T. (1947) The Similarity Law of Transonic Flow. Journal of Mathematics and Physics, 26, 182-190. https://doi.org/10.1002/sapm1947261182
[7]	Berndt, S. and Oswatitsch, K. (1950) Aerodynamic Similarity at Axisymmetric Transonic Flow around Slender Bodies. In: Oswatitsch, K., Ed., Contributions to the Development of Gasdynamics, Vieweg+Teubner Verlag, German, 68-75.
[8]	Busemann, A. (1929) Drücke auf kegelformige Spitzen bei Bewegung mit überschallgeschwindigkeit. Journal of Applied Mathematics and Mechanics, 9, 496-498. https://doi.org/10.1002/zamm.19290090619
[9]	Maccoll, J.W. and Taylor, G.I. (1933) The Air Pressure over a Cone Moving at High Speeds. Proceedings of the Royal Society, 139, 278-311. https://doi.org/10.1098/rspa.1933.0017
[10]	Yoshihara, H. (1952) On the Flow over a Cone-Cylinder Body at Mach Number One. WADC TR q2-295.
[11]	Solomon, G.E. (1953) Transonic Flow Past Cone-Cylinders. REPORT 1242, California Institute of Technology, California.
[12]	Kennedy, C. and Gruber, A. (2008) Reduced Aliasing Formulations of the Convective Terms within the Navier-Stokes Equations for a Compressible Fluid. Journal of Computational Physics, 227, 1676-1700. https://doi.org/10.1016/j.jcp.2007.09.020
[13]	Ducros, F., Ferrand, V., Nicoud, F., Weber, C., Darracq, D., Gacherieu, C. and Poinsot, T. (1999) Large-Eddy Simulation of the Shock/Turbulence Interaction. Journal of Computational Physics, 152, 517-549. https://doi.org/10.1006/jcph.1999.6238
[14]	Pirozzoli, S. and Turkel, E. (2011) Numerical Methods for High-Speed Flows. Annual Review of Fluid Mechanics, 43, 163-194. https://doi.org/10.1146/annurev-fluid-122109-160718
[15]	Swanson, R.C. and Turkel, E. (1992) On Central-Difference and Upwind Schemes. Journal of Computational Physics, 101, 292-306. https://doi.org/10.1016/0021-9991(92)90007-L
[16]	Kurganov, A. and Tadmor, E. (2000) New High-Resolution Central Schemes for Nonlinear Conservation Laws and Convection-Diffusion Equations. Journal of Computational Physics, 180, 241-282. https://doi.org/10.1006/jcph.2000.6459
[17]	Kurganov, A., Noelle, S. and Petrova, G. (2001) Semi-Discrete Central-Upwind Scheme for Hyperbolic Conservation Laws and Hamilton—Jacobi Equations. SIAM Journal on Scientific Computing, 23, 707-740. https://doi.org/10.1137/S1064827500373413
[18]	Greenshields, C.J., Weller, H.G., Gasparini, L. and Reese, J.M. (2009) Implementation of Semi-Discrete, Non-Staggered Central Schemes in a Colocated, Polyhedral, Finite Volume Framework, for High-Speed Viscous Flows. International Journal for Numerical Methods in Fluids, 63, 1-21. https://doi.org/10.1002/fld.2069
[19]	Kurganov, A. and Petrova, G. (2005) Central-Upwind Schemes on Triangular Grids for Hyperbolic Systems of Conservation Laws. Numerical Methods for Partial Differential Equations, 21, 536-552. https://doi.org/10.1002/num.20049
[20]	Van Leer, B. (1974) Towards the Ultimate Conservative Difference Scheme, II: Monotonicity and Conservation Combined in a Second Order Scheme. Journal of Computational Physics, 14, 361-370. https://doi.org/10.1016/0021-9991(74)90019-9
[21]	Kurganov, A. and Lin, C.T. (2007) On the Reduction of Numerical Dissipation in Central-Upwind Schemes. Computer Physics Communications, 2, 141-163.
[22]	Wu, J.M. (1973) A Theory For Subsonic And Transonic Flow Over A Cone-With Or Without Small Yaw Angle. AD-776 374, Tennessee University Space Institute, Tullahoma.
[23]	Zangeneh, R. (2020) Development of a New Algorithm for Modeling Viscous Transonic Flow on Unstructured Grids at High Reynolds Numbers. Journal of Fluids Engineering, 143, Article ID: 024504. https://doi.org/10.1115/1.4048611
[24]	Zangeneh, R. (2020) Evaluation of Reattaching Shear-layer in Compressible Turbulent Flows; A Large-Eddy Simulation Approach. Proceedings of the ASME 2020 Fluids Engineering Division, Virtual, Online, 13-15 July 2020, No. V003T05A035.
[25]	Zangeneh, R. (2020) A New Framework for Modeling Shock-Turbulence Interactions. SAE International Journal of Aerospace. https://doi.org/10.4271/2020-01-5092