0 Introduction

Numerical solution of Euler equations in conservative form is discussed in Reference[2]using a generalization of Murman Scheme[1].In[3],Chattot constructed four operators for subsonic,supersonic,sonic and shock points using eigenvalues and eigenvectors of the Jacobian matrix and the compatibility relations to solve the conservation laws of a perfect gas in a convergent divergent nozzle.In this paper,a simplified scheme is introduced and the calculations for flows with shocks are presented.In the second part of the paper,compressible flow with friction and heat addition is considered followed by steady and unsteady detonation.

We use a predictor-corrector compact scheme that combines upwinding and central differencing.Consider 1D Non-Reactive Euler equations in a flux vector form,where =[ρ,ρu,ρE] gives:

点击浏览公式

The present numerical scheme is a two-step method.Central difference scheme is used in the subsonic regions and at the sonic points,backward difference is used in the supersonic regions,with a special shock point operator at shock points.The scheme is conservative almost everywhere and,unlike Chattot's calculations,no special sonic point operator is used in this work. The finite difference formulas and their stencils are shown in Figs. 1-3.

Fig. 1 Subsonic/sonic stencil and scheme.

图选项

Fig. 2 Supersonic stencil and scheme.

图选项

Fig. 3 Shock point stencil and scheme.

图选项

A first validation test case shows the sharpness of the shock,correct location (x=0.82) and flow through the sonic point,see Fig. 4.

Fig. 4 Compressible nozzle flow 1D validation test case.

图选项

Examples in 2D for transonic flow over a thin parabolic arc at incident M = 0.675 and M = 1.4 are in good agreement with Ni’s,who used second order accurate Lax-Wendroff scheme [4]. In Fig. 5 an Oswatitsch-Zeirep singularity at the root of the shock is observed. The present method shows better resolution of this feature and a sharper shock with higher peak than Lax-Wendroff scheme.

Fig. 5 (M∞=0.675) Present vs. 2-step Lax-Wendroff method over thin parabolic arc at different grids.

图选项

In Fig. 6,a grid refinement study is shown in the case of flow over a thin parabolic arc.

Fig. 6 Grid refinement study for present method showing representative grids of (121×40),(241×80) and (301×100).

图选项

Numerical cases for 2D transonic flow over a circular cylinder are also considered. Comparison with the results of 2-step Lax-Wendroff is shown in Fig. 7. The numerical results are in good agreement with the literature [5].

It should be remarked that in the two dimensional calculations,the present scheme is used only in the flow direction,while centered differences,augmented with artificial viscosity,are used for the lateral direction.

Fig. 7 (M∞=0.5) 2D flow over cylinder,comparison of Lax-Wendroff 2-step to present scheme.

图选项

For a hypersonic flow over a rectangular block is computed for Mach number M=10. This flow is comparable with the example problem and results published in [6]. In this calculation γ=1.17,the numerical mesh is 88×136 in the x and y directions respectively. Uniform far field boundary conditions are applied ,as well as symmetry conditions on the centerline axis. The computational domain is 0≤y≤34. The rectangular body starts at x=18 and the top of the body is at y=12,see Fig. 8.

Fig. 8 Hypersonic flow over a rectangular block is computed for M=10,present vs. 2-step Lax-Wendroff method.

图选项

The Quasi-1D Euler equations with friction and heat addition are given by the conservation laws for mass,momentum and energy respectively:

点击浏览公式

点击浏览公式

点击浏览公式

Together with the definitions:

A steady compressible perfect gas in a duct of constant cross sectional area is considered. There are two possibilities to get from the initially supersonic condition to the final sonic flow: by introduction of the heat source q or the friction source f. For compressible duct flows with friction and heat addition,several formulas are available in terms of Mach number,the friction factor and heat addition,see Reference [7]. Examples of choking due to friction and heat addition are calculated with the present scheme,see Fig. 9 and the results are in good agreement with those in Reference [7]. Although the calculation is very sensitive to the amount of heat and the shock moves to the right or the left with slight variation of the critical amount,the shock is very sharp,especially when compared with other artificial viscosity methods,for example Lax-Wendroff scheme.

Fig. 9 Compressible duct flow 1D with heat addition (up) and friction (down).

图选项

An ordinary differential equation for the species,together with algebraic relations derived from the shock jump conditions gives the steady Zeldovich-Von Nuemann-Doring (ZND) family of solutions to the Euler equations. Although detonations are funda-mentally unstable,the ZND solutions are suitable,in an averaged sense,for many cases. It is shown however by several authors,see for example [8, 9, 10, 11],that when solving the Euler equations,the solutions are predicted to be oscillatory and even chaotic for many examples. In general,as the parameter E_a,governing the activation energy for the chemical species,is increased,instability increases and the ZND solution becomes less predictive. However,for sufficiently small E_a,the true solution is nearly steady and the correspondence with the ZND solution agrees very well. Detonation in one dimension is governed by the reactive unsteady Euler equations:

点击浏览公式

点击浏览公式

点击浏览公式

点击浏览公式

Where m,f,q are source terms for mass,friction and heat respectively. The parameter k is a constant factor that relates the reaction and convection rates. The following non-dimensional variables are used:

点击浏览公式

Here u denotes axial velocity component and C_f,M is the coefficient of friction and Mach number respectively. In our calculations:

In characteristic coordinates,in the absence of source terms,we are solving the system:

点击浏览公式

点击浏览公式

点击浏览公式

点击浏览公式

where ξ=x-st. Initially the shock speed s is given from the ZND theory. Thereafter,the shock speed s=s(t),is corrected for unsteady flow in the characteristic coordinates ξ by requiring Y=0.5 remain at x=0. Letting L be a linear operator and N be the nonlinear reaction operator, w =[ρ,ρu,ρE,ρY],the reactive Euler Equations can be written compactly as:

点击浏览公式

In our computation an articial viscosity with ε= O(Δx) is added to stabilize the solution near the shock. We use an implicit third order upwinding scheme for conservation of mass,momentum and Energy

点击浏览公式

To solve this stiff system,an IMEX operator splitting method for the linear and nonlinear terms in the species equation is used,by updating Y first implicitly with third order upwinding,then explicitly updating the nonlinear terms.

点击浏览公式

Several validation computations are performed,for example the steady state ZND solution which is derived from the non-dimensionalized shock-jump conditions. At the shock front the following relations hold:

点击浏览公式

点击浏览公式

点击浏览公式

Downstream of the shock we then solve an ordinary differential equation for the chemical species Y,

点击浏览公式

The solution is achieved by combining the mass and momentum conditions to get ,and combining the momentum and energy conditions to get a quadratic equation for u:

点击浏览公式

The ZND model setup is shown in Fig. 10.

Fig. 10 Sketch of ZND solution.

图选项

Hugoniot curves for ZND Theory are shown in Fig. 11. As heat is released,supersonic flow decelerates continuously to the sonic condition,or we can go along the Rayleigh line to the Von-Neumann spike and then go back along the same Rayleigh line to the sonic condition.

Fig. 11 Sketch showing C-J detonation as well as strong and weak detonations.

图选项

Several verification computations are performed using the steady state ZND solution and are in agreement with results reported in [12]. An example showing several ZND temperature profiles for a range of E_a values is shown in Fig. 12.

Fig. 12 Several ZND temperature profiles for a range of Ea values.

图选项

The third order upwinding method is compared with several ZND solutions for verification. A representative comparison is shown in Fig. 13. The ZND solution is computed starting from the shock front while the Euler equations are integrated across the shock starting upstream with unreacted flow.

Fig. 13 ZND solution by Algebraic/ODE formulation vs. reactive Euler equations. Notice the scale difference between the upper and lower figures to emphasize the difference between the two solutions.

图选项

A Chapman-Jouquet detonation following [13] is computed in Fig. 14. We compute the shock speed to be s=5.441 9 which exactly matches the CJ-shock speed given in [13]. Additional comparisons with steady state computations are in good agreement with the results reported in [12, 13, 14].

Fig. 14 Comparison with C-J detonation.

图选项

The case of friction losses dominating heat transfer effects can sometimes occur when a detonation wave passes through a porous media [14]. In such a case,momentum loss due to friction can be much more influential than heat transfer effects. Following the work of [8] and [11],Fig. 15 shows that a source of friction in the momentum equation causes qualitatively the same effect as higher E_a.

Fig. 15 Ea = 25 (top),Ea= 28 (bottom),γ=1.2 showing a range of Cf coeffcients and transition to instability.

图选项

Our computations show that friction causes a momentum loss which decreases the peak temperature and local velocity,as shown in Fig. 16. Increase in the friction causing reduction in temperature was also reported in [8] and [11].

Fig. 16 Ea = 23,γ=1.2,comparing Cf = 0 with Cf=0.004,lower temperature,T,and velocity,u,are observed.

图选项

The effect of friction on detonation stability is also considered. Treating E_a as a bifurcation parameter,the peak pressure at the Von-Nuemann state is stored for several simulations as a time series. The attractor space is then constructed from the local maxima and minima of each time series. Representative examples are shown in Fig. 17. In our results,increasing friction causes period doubling for lower values of E_a,this trend is also reported in [8] and [11].

Fig. 17 Von-Neumann peak pressure for Cf = 0 (left),Cf=0.002 (center),Cf = 0.004 (right), γ= 1.2, showing translation and compression of period doubling regions for variable Ea.

图选项

In the present work we have demonstrated a simple scheme for computing transonic flows which can be considered as a generalization of Murman’s four point operator to solve the Euler equations. We have examined several transonic steady compressible 1D flows with friction and heat addition. Unsteady flows with weak dissipation and friction have been considered in the context of detonation theory for several reactive compressible cases. Computations with higher order IMEX schemes is also of interest [15].