Transonic combustion: Steady and unsteady potential models

Consider a compressible flow of a reactive fluid in a two dimensional channel. Letting the ∞ -state denote the uniform flow conditions at the channel inlet, a non-dimensionalization is made by, where [ ⋅ ] denotes a dimensional variable, ∞ the inlet condition and L is half the width of the channel:

M ∞ = u ∞ 2 γ ⁢ R g ⁢ a ⁢ s ⁢ T ∞ , ρ = [ ρ ] ρ ∞ , p = [ p ] ρ ∞ ⁢ u ∞ 2 = [ p ] p ∞ ⁢ γ ⁢ M ∞ 2 ,

(2) T = [ T ] ( V ∞ 2 / C p ) , u = [ u ] u ∞ , v = [ v ] v ∞ , Y = [ Y ] Y ∞ , x = [ x ] L , y = [ y ] L .

Some authors non-dimensionalize T = [ T ] T ∞ , a conversion between this non-dimensionalization and ours is given by:

T = [ T ] ( V ∞ 2 / C p ) = [ T ] T ∞ ⁢ 1 γ - 1 ⁢ 1 M ∞ 2 .(3)

The dimensionless parameters are given by:

M ∞ 2 = u ∞ 2 γ ⁢ R g ⁢ a ⁢ s ⁢ T ∞ , A = D a u ∞ , θ = E ⁢ a C p ⁢ T ∞ ⁢ 1 γ - 1 ⁢ 1 M ∞ 2 , β = [ B ] u ∞ 2 = [ B ] C p ⁢ T ∞ ⁢ 1 γ - 1 ⁢ 1 M ∞ 2 .(4)

Now, the dimensionless governing equations, in the inviscid limit of the momentum equations, without viscous dissipation or mass diffusion become:

ρ t + ( ρ ⁢ u ) x + ( ρ ⁢ v ) y = 0 , f := ( ρ ⁢ u ) t + ( ρ ⁢ u 2 ) x + ( ρ ⁢ v ⁢ u ) y = - p x , g := ( ρ ⁢ v ) t + ( ρ ⁢ u ⁢ v ) x + ( ρ ⁢ v 2 ) y = - p y , ( ρ ⁢ T ) t + ( ρ ⁢ u ⁢ [ T + 1 2 ⁢ ( u 2 + v 2 ) ] ) x + ( ρ ⁢ v ⁢ [ T + 1 2 ⁢ ( u 2 + v 2 ) ] ) y = 1 P ⁢ r ⋅ R ⁢ e ⁢ Δ ⁢ T + β ⁢ w , ( ρ ⁢ Y ) t + ( ρ ⁢ u ⁢ Y ) x + ( ρ ⁢ v ⁢ Y ) y = 1 S ⁢ c ⋅ R ⁢ e ⁢ ( Y x ⁢ x + Y y ⁢ y ) - w , p = ρ ⁢ T ⁢ γ - 1 γ , w = A ⁢ ρ ⁢ Y ⁢ e - θ T .(5)

Where P ⁢ r = μ ⁢ C p k is the Prandtl Number, and S ⁢ c = μ ρ ⁢ D is the Schmidt Number. This is a model for the reactive Euler Equations for compressible flow with premixed Arrhenius combustion, see [1, 7] for more details.

3. Reactive transonic full potential model

Decomposition of the velocity field into a curl-free component and a divergence-free component gives:

u = ϕ x + ψ y ρ , v = ϕ y - ψ x ρ .(6)

This Helmholtz orthogonal decomposition of the velocity field gives us a potential and stream function formulation for the compressible Euler Eq. (2.1) without combustion ( A = 0):

( ρ ⁢ ϕ x ) x + ( ρ ⁢ ϕ y ) y + ρ t = 0 , ( ψ x ρ ) x + ( ψ y ρ ) y = - ω , ( ρ ⁢ u ⁢ ω ) x + ( ρ ⁢ v ⁢ ω ) y = F ⁢ ( ρ , ω ) , f x + g y = ρ t ⁢ t - ( p x ⁢ x + p y ⁢ y ) ,(7) ρ ⁢ T t + ρ ⁢ u ⁢ [ T + 1 2 ⁢ ( u 2 + v 2 ) ] x + ρ ⁢ v ⁢ [ T + 1 2 ⁢ ( u 2 + v 2 ) ] y = 1 P ⁢ r ⋅ R ⁢ e ⁢ Δ ⁢ T .

Where f and g stand for the inertia terms in the momentum equations. The above Poisson equation for the pressure is derived by taking the divergence of the momentum equations in primitive variables:

- ∇ ⋅ ∇ ⁡ p = ∂ ⁡ f ∂ ⁡ x + ∂ ⁡ g ∂ ⁡ y = ∇ ⋅ F → + ∂ ⁡ ∇ ⋅ ρ ⁢ u ∂ ⁡ t = ∇ ⋅ F → + ∂ 2 ⁡ ρ ∂ ⁡ t 2 .(8)

This model admits the incompressible flow solution as a special case when density is constant. Some simplifications of the above model are possible. For example, in the case of irrotational compressible flow, the stream function equation can be omitted, and in the case of irrotational incompressible flow, instead of the pressure equation, Bernoulli’s law can be used. We will assume the flow is irrotational (i.e. ω = ∇ × V → = 0 ), but we will account for heat release due to chemical reaction coupling as in the reactive Euler Equations. In particular, for the expanding duct flow under consideration in this paper, shock fronts are almost normal to the principal flow direction. Since Crocco relation relates vorticity to shock curvature [8], this suggests a framework where under such conditions a compressible potential flow model might be used.

For irrotational flows, a velocity potential function ϕ must exist such that q → = ∇ ⁡ ϕ . Under this assumption the velocity field is curl-free and reduces to: u = ϕ x , v = ϕ y . An irrotational unsteady compressible reactive full potential flow model is given by:

( ρ ⁢ ϕ x ) x + ( ρ ⁢ ϕ y ) y + ρ t = 0 , [ ( ρ ⁢ u 2 ) x + ( ρ ⁢ u ⁢ v ) y ] x + [ ( ρ ⁢ u ⁢ v ) x + ( ρ ⁢ v 2 ) y ] y - 1 P ⁢ r ⋅ R ⁢ e ⁢ Δ ⁢ ( ( ρ ⁢ u ) x + ( ρ ⁢ v ) y ) = ρ t ⁢ t - ( p x ⁢ x + p y ⁢ y ) , ρ ⁢ T t + ρ ⁢ u ⁢ [ T + 1 2 ⁢ ( u 2 + v 2 ) ] x + ρ ⁢ v ⁢ [ T + 1 2 ⁢ ( u 2 + v 2 ) ] y = 1 P ⁢ r ⋅ R ⁢ e ⁢ Δ ⁢ T + β ⁢ w , ρ ⁢ Y t + ρ ⁢ u ⁢ Y x + ρ ⁢ v ⁢ Y y = 1 S ⁢ c ⋅ R ⁢ e ⁢ Δ ⁢ Y - w , p = ρ ⁢ T ⁢ γ - 1 γ , w = A ⁢ ρ ⁢ Y ⁢ e - θ T .(9)

It should be noted that unlike the classical density-based approach, this is a pressure-based approach. For steady flow, the pressure equation can be treated as a Poisson’s equation. For unsteady flow, the following hyperbolic equation for pressure is introduced:

ρ t ⁢ t = ( p x ⁢ x + p y ⁢ y ) + f x + g y ⟶ γ γ - 1 ⁢ ( p T ) t ⁢ t = ( p x ⁢ x + p y ⁢ y ) + f x + g y .(10)

For the present work we simply choose:

1 S ⁢ c ⋅ R ⁢ e = 1 P ⁢ r ⋅ R ⁢ e = O ⁢ ( Δ ⁢ x ) .(11)

Using a small disturbance expansion around uniform inlet conditions gives a reactive transonic small disturbance model similar to that proposed in [11]:

[ 1 - M ∞ 2 - ( γ + 1 ) ⁢ M ∞ 2 ⁢ φ x ] ⁢ ∂ 2 ⁡ φ ∂ ⁡ x 2 + ∂ 2 ⁡ φ ∂ ⁡ y 2 = β ∼ ⁢ K δ ⁢ A ∼ ⁢ Y ∼ ⁢ e - θ ∼ / ( 1 + T ∼ ) ,

(12) T ∼ = - ( γ - 1 ) ⁢ φ x , Y ∼ x = A ∼ ⁢ Y ∼ ⁢ e - θ ∼ / ( 1 + T ∼ ) .

This governing equation is of mixed type, elliptic and hyperbolic. Here T ∼ , Y ∼ denote small disturbance variables, and A ∼ , β ∼ , K δ , and θ ∼ are constant reaction control parameters, and M ∞ is the inlet Mach number. The velocity can be calculated from: u ∼ = φ x , v ∼ = φ y . Details of the derivation, scaling, and variations are available in [11]. The small disturbance approach is limited to small heat release, i.e. diluted premixtures, and small area variation. According to [11] the precise limitations are upstream flows with Mach number between 0.75 and 1.2, premixtures with reactant mass fractions up to 0.1, and channel cross-sectional area variation of 15 percent or less. The boundary conditions at the channel inlet ( - 1 ⩽ y ⩽ 1 ) are:

u ⁢ ( 0 , y ) = 1 + α ⁢ u 0 ⁢ ( y ) , v x ⁢ ( 0 , y ) = 0 , T ⁢ ( 0 , y ) = 1 + α ⁢ T 0 ⁢ ( y ) , Y ⁢ ( 0 , y ) = δ ⁢ Y 0 ⁢ ( y ) ,(13)

where . We will take α = 0 in this work. At the channel outlet ( y 𝑙𝑜𝑤𝑒𝑟 ⩽ y ⩽ y 𝑢𝑝𝑝𝑒𝑟 ):

v ⁢ ( x 𝑜𝑢𝑡 , y ) = 0 , ∂ ⁡ T ⁢ ( x 𝑜𝑢𝑡 , y ) ∂ ⁡ y = 0 , ∂ ⁡ Y ⁢ ( x 𝑜𝑢𝑡 , y ) ∂ ⁡ y = 0 .(14)

OPEN IN VIEWER

Figure 1.

Small Disturbance Model (TSD) for M a = 0.75, 0.80, 0.85, 0.90, 0.95 for A = 0, showing transition to shocked flow with Zeirup Singularity at 12% nozzle area variation parameter.

OPEN IN VIEWER

Figure 2.

Small Disturbance Model (TSD) T and Y profiles for M a = 0.75 for A = 0, 10000, 20000, 30000. Higher pre-exponential frequency factor A consumes more reactant (left), results in a lower thermal energy (middle) and higher Mach number (right) at the outlet.

OPEN IN VIEWER

Figure 3.

Small Disturbance Model (TSD) T and Y profiles for M a = 0.80, A = 20000 for range of channel divergence parameters ϵ 𝑤𝑎𝑙𝑙 = 3%, 6%, 9%, 12%, 15%. Higher channel divergence consumes more reactant (left) and results in a higher thermal energy (middle) and lower Mach number (right) at the outlet.

Along the channel walls:

v ≈ d ⁢ G 𝑏𝑜𝑑𝑦 d ⁢ x .(15)

As in [11] the nozzle geometry divergence is specified smoothly from a cosine continuation at x = 1 :

(16)

Typically ϵ 𝑤𝑎𝑙𝑙 = 0.12 = 12 % , unless otherwise specified, as in Fig. 3. We assume a symmetric flow, at the axis of symmetry, where y = 0 :

v ⁢ ( x , 0 ) = 0 , T y ⁢ ( x , 0 ) = 0 , Y y ⁢ ( x , 0 ) = 0 .(17)

In these calculations: β ∼ ⁢ K δ = 2 , θ ∼ = 12 , A ∼ = ϵ 2 / 3 ⁢ A . The mesh is 181 points in the x -direction from x = 0 to x = 3 , and 101 points in the y -direction from y = 0 to y = 1 . In this case, the 4-point operator of Murman-Cole is used to solve this mixed type elliptic hyperbolic equation with line successive under relaxation implicit in the y -direction [9]. The reaction ODE is solved with Simpson’s integration method as in [11]. The Mach number is calculated by:

M a = M ∞ ⁢ ( 1 + ( γ + 1 ) ⁢ u ∼ ) .(18)

Non-reactive and reactive validation test cases are shown below and are in good agreement with [11]. Example results showing the transition from fully subsonic to shocked flow with a supersonic bubble for the non-reactive flow are given in Fig. 1.

As reported in [11], it is found that higher pre-exponential frequency factor A consumes more reactant and results in a lower thermal energy and higher Mach number at the outlet, see Fig. 2. Higher peak Mach number as the duct is widening is also associated in this case with higher A and higher reactant consumption.

Similarly, it is found that variation of the channel curvature through the factor ϵ consumes slightly more reactant and results in a higher thermal energy and lower Mach number at the outlet, see Fig. 3.

An important consequence of the heat release due to combustion is the capability to cause shocks as shown in Fig. 4. In this figure, for M a = 0.85, without chemical reaction ( A = 0) the flow is completely subsonic and smooth, and with chemical reaction ( A = 20000) the flow contains a supersonic bubble which abruptly becomes subsonic through a normal shock. The outlet Mach number is higher and the outlet coefficient of pressure ( C p ) is lower for the reactive flow case where heat is released.

OPEN IN VIEWER

Figure 4.

Small Disturbance Model (TSD), M a = 0.85, A = 0 and A = 20000, shows supersonic bubble due to chemical reaction and small Zierup Singularity.

5. Numerical methods for the full potential model

To solve the full potential equation we use the artificial compressibility methods [4, 13]. We consider the density biasing and flux biasing formulations.

( ρ ⁢ ϕ x ) x + ( ρ ⁢ ϕ y ) y = 0 ⟶ ( ρ ¯ ⁢ ϕ x ) x + ( ρ ¯ ⁢ ϕ y ) y = 0 .(19)

A modified artificial density method that is based on flux biasing is given by:

ρ ¯ ( i + 1 / 2 , j + 1 / 2 ) = ρ ( i + 1 / 2 , j + 1 / 2 ) - [ u q 2 ⁢ Δ ⁢ x ⁢ ∂ ⁡ F ∂ ⁡ x + v q 2 ⁢ Δ ⁢ y ⁢ ∂ ⁡ F ∂ ⁡ y ] ( i + 1 / 2 , j + 1 / 2 ) .(20)

After calculating the local Mach number at the cell centers:

M ( i + 1 / 2 , j + 1 / 2 ) 2 = q 2 a 2 = u 2 + v 2 1 M ∞ 2 + 1 2 ⁢ ( γ - 1 ) ⁢ [ 1 - u 2 - v 2 ] .(21)

The flux F is calculated at the cell-centers from:

(22)

Where the following definitions are used:

a * = ( 2 ⁢ γ - 1 γ + 1 ⁢ [ 1 γ - 1 ⁢ 1 M ∞ 2 + 1 2 ] ) 1 / 2(23) ρ * = ( 1 - M ∞ 2 ⁢ ( γ - 1 ) 2 ⁢ ( ( a * ) 2 - 1 ) ) 1 / ( γ - 1 ) .(24)

OPEN IN VIEWER

OPEN IN VIEWER

OPEN IN VIEWER

OPEN IN VIEWER

OPEN IN VIEWER

OPEN IN VIEWER

OPEN IN VIEWER

OPEN IN VIEWER

OPEN IN VIEWER

OPEN IN VIEWER

In the density biasing scheme, instead of biasing the flux F := ρ ⁢ q , now the density alone is biased:

ρ i + 1 / 2 , j + 1 / 2 = ρ i + 1 / 2 , j + 1 / 2 - μ ⁢ ( ρ i + 1 / 2 , j + 1 / 2 - ρ i - 1 / 2 , j + 1 / 2 ) .(25)

Where the upwind shifting parameter μ is defined in terms of the local Mach number [5]:

μ = τ ⋅ max ⁡ ( 0 , 1 - M c 2 M i + 1 / 2 , j + 1 / 2 2 , 1 - M c 2 M i - 1 / 2 , j + 1 / 2 2 ) .(26)

Where τ , M c 2 are user specified parameters. The energy and species equations are solved using a conservative mixed upwind/centered scheme with backward difference in the x -direction. A staggered grid is used as show in Fig. 5.

For the case of steady state computations, all the time dependent terms are removed from the governing equations. The pressure Poisson equation is solved for p :

- [ ( ρ ⁢ u 2 ) x + ( ρ ⁢ u ⁢ v ) y ] x - [ ( ρ ⁢ u ⁢ v ) x + ( ρ ⁢ v 2 ) y ] y = ( p x ⁢ x + p y ⁢ y ) ,(27)

and the density is updated via the perfect gas law by line successive under relaxation until convergence, where n denotes the iteration count:

ρ n + 1 = ρ n + ω 𝑆𝑂𝑅 ⁢ ( p T ⁢ γ γ - 1 - ρ n ) .(28)

Since the the pressure, temperature and species equations are elliptic, they are iteratively solved until convergence before updating the potential function Φ again, and then recalculating the velocity field.

For the case of unsteady computations a stability term of the form: ν ⁢ ρ t , with , is added to the pressure equation to help with convergence:

[ ( ρ ⁢ u 2 ) x + ( ρ ⁢ u ⁢ v ) y ] x + [ ( ρ ⁢ u ⁢ v ) x + ( ρ ⁢ v 2 ) y ] y = ρ t ⁢ t + ν ⁢ ρ t - ( p x ⁢ x + p y ⁢ y ) .(29)

This term is a natural choice because when small viscosity is present in the momentum equations:

ν ⁢ Δ ⁢ [ ( ρ ⁢ u ) x + ( ρ ⁢ v ) y ] ≈ - ν ⁢ ρ t .(30)

Where 0 ⩽ ν ⩽ 1 . One way to proceed is that the pressure equation is solved for ρ and the gas law is used to update p . The other, which has the advantage of being able to treat operators implicitly, is to solve the pressure equation for p and update the density from the gas law. In this case, one must first update the temperature to the current ( n + 1 ) time step, T n + 1 , then, using the gas law, p n + 1 is solved according to the nonhomogeneous wave equation:

γ γ - 1 ⁢ [ ( P T ) n + 1 - 2 ⁢ ( P T ) n + ( P T ) n - 1 Δ ⁢ t 2 + ν ⁢ ( P T ) n + 1 - ( P T ) n - 1 2 ⁢ Δ ⁢ t ] - ( p x ⁢ x + p y ⁢ y ) n + 1 = f x n + g y n .(31)

Boundary conditions, in terms of the standard momentum equations, must be used to eliminate spurious solutions.

A validation test case is performed in 1D and compared against the quasi-1D Euler equations following the example given in [2]. We write the energy equation in terms of temperature using the relation:

E = T γ + 1 2 ⁢ u 2 .(32)

Then the energy equation becomes:

( A ⁢ ρ ⁢ E ) t + ( A ⁢ ρ ⁢ u ⁢ H ) x = 0 , ( A ⁢ ρ ⁢ ( T γ + 1 2 ⁢ u 2 ) ) t + ( A ⁢ ρ ⁢ u ⁢ ( T + 1 2 ⁢ u 2 ) ) x = ϵ ⁢ ( A ⁢ ρ ⁢ E ) x ⁢ x ,(33) ( A ⁢ ρ ⁢ T ) t + γ 2 ⁢ ( ( A ⁢ ρ ⁢ u ) 2 A ⁢ ρ ) t + γ ⁢ ( A ⁢ ρ ⁢ u ⁢ ( T + 1 2 ⁢ u 2 ) ) x = γ ⁢ ϵ ⁢ ( A ⁢ ρ ⁢ E ) x ⁢ x .

Differentiation of conservation of mass with respect to t and conservation of momentum with respect to x , inclusive of artificial viscosity terms, gives:

- ( A ⁢ ρ ) t ⁢ t = ( A ⁢ ρ ⁢ u ) x ⁢ t - ϵ ⁢ ( A ⁢ ρ ) 𝑥𝑥𝑡 ,(34) ( A ⁢ p ) x ⁢ x = - ( A ⁢ ρ ⁢ u ) t ⁢ x - ( A ⁢ ρ ⁢ u 2 ) x ⁢ x + ( A x ⁢ p ) x + ϵ ⁢ ( A ⁢ ρ ⁢ u ) 𝑥𝑥𝑥 .

Adding these equations together implies:

( A ⁢ p ) x ⁢ x - ( A ⁢ ρ ) t ⁢ t = - ( A ⁢ ρ ⁢ u 2 ) x ⁢ x + ( A x ⁢ p ) x - ϵ ⁢ ( A ⁢ ρ ) 𝑥𝑥𝑡 + ϵ ⁢ ( A ⁢ ρ ⁢ u ) 𝑥𝑥𝑥 .

Therefore, the present pressure based model for the quasi-1D Euler Equations and small artificial viscosity (with strength ϵ ) is given by:

( A ⁢ ρ ) t ⁢ t = ( A ⁢ p ) x ⁢ x + ( A ⁢ ρ ⁢ u 2 ) x ⁢ x - ( A x ⁢ p ρ ) x + ϵ ⁢ ( A ⁢ ρ ) 𝑥𝑥𝑡 - ϵ ⁢ ( A ⁢ ρ ⁢ u ) 𝑥𝑥𝑥 , ( A ⁢ ρ ⁢ u ) x = - ( A ⁢ ρ ) t + ϵ ⁢ ( ( A ⁢ ρ ⁢ u ) x ⁢ x + ( A ⁢ ρ ) x ⁢ x ) , ( A ⁢ p ) t = - γ - 1 2 ⁢ ( A ⁢ ρ ⁢ u 2 ) t - ( γ - 1 ) ⁢ ( A ⁢ ρ ⁢ u ⁢ ( T + 1 2 ⁢ u 2 ) ) x ⁢ ( γ - 1 ) ⁢ ϵ ⁢ ( A ⁢ p ) x ⁢ x , T = γ γ - 1 ⁢ p ρ .(35)

This system is hyperbolic and converges with Δ ⁢ t = O ⁢ ( Δ ⁢ x ) . Alternatively, one can solve the continuity equation for a potential function ϕ , where u = ϕ x , without loss of generality as is proposed earlier.

( A ⁢ ρ ⁢ ϕ x ) x + ( A ⁢ ρ ) t = 0 .(36)

A comparison showing the quasi-1D Euler equations (left) and the quasi-1D full potential model (right) is shown in Fig. 6. The correct shock location ( x = 0.82) and flow through the sonic point is acheived in both models, see [2].

The small disturbance model is essentially mixed type in the x -direction only. In contrast, the full potential can handle mixed type flows in multiple principal directions. As in the small disturbance model, a validation case showing the transition from fully subsonic to shocked flow with a supersonic bubble for a non-reactive flow case is shown first in Fig. 7.

When the full scheme is used, the results are slightly different from the small disturbance model, although the shock location prediction as well as peak Mach numbers are similar. The pre-shock solution and shock location are very similar, the post shock conditions are seen to vary between the TSD and full potential model, as shown in Fig. 8.

The profiles for Y and T are markedly different from the TSD model. In particular, where the TSD model predicts that temperature trends like C p , in the full model, the temperature increases with heat release. The Y profiles are similar to that of the TSD model, but in our calculations there is a more noticeable effect of the geometry creating variation in reactant consumption rate.

In Fig. 9 the case of M ∞ = 0.8 for a range of pre-exponential frequency factors A = 0, 10000, 20000, 30000 is computed. This result is in good agreement with the small disturbance model.

In Fig. 10 the case of M ∞ = 0.85 for a range of pre-exponential frequency factors A = 0, 10000, 20000, 30000 is computed. In this case, the results show the combustion induces a supersonic flow bubble. The trend with respect to Mach number and pressure coefficient are in agreement with the small disturbance model. As in the case with the TSD model, without chemical reaction ( A = 0) the flow is completely subsonic and smooth, and with chemical reaction a supersonic bubble is shown to develop, the outlet Mach number is higher and the outlet coefficient of pressure ( C p ) is lower for the reactive flow cases where heat is released.

In Fig. 11 the case of incident ρ -biased full potential model for M ∞ = 0.80, 0.85, 0.90, 0.95 with A = 20000 is computed.

An unsteady example has been computed as a validation test for M ∞ = 0.80 with A = 0. The results are in good agreement with the steady state model and are shown in Fig. 12.

Also the case where combustion causes a shock in the flow for M ∞ = 0.85 with A = 20000 has been computed. The results are shown in Fig. 13 and indicates the unsteady model produces a slightly sharper shock that is in nearly the same location as the steady model.

The time dependent transient history for M ∞ = 0.85 with A = 20000, starting from uniform quiescent initial data show the effect of combustion on the peak Mach number and fuel consumption in Fig. 14. This result shows that the steady state exists and that startup can cause shocked flow, but the long time behavior of the solution is not dependent on this transient behavior.

8. Conclusions

It is clear that the small disturbance approach is limited to small heat release and small area variation. Although the small disturbance model is limited in its range of applicability, within that range it is simpler to apply and can produce quite reasonable results, see [11]. The full potential model is more general and should be able to handle a larger heat release and area variation as indicated in the text. In addition, it should be straightforward to include a more detailed reaction mechanism without significant modifications to the model.

In the present work we have developed a full potential model with chemical reactions and numerical methods for its solution, including cases with strong sharp shocks. Moreoever, we introduced a new scheme based on pressure correction using a hyperbolic equation and the results have been verified for an unsteady one dimensional problem. The merits of the new formulation compared to the standard elliptic equation for the pressure are discussed.