Transonic combustion: Steady and unsteady potential models

Abstract

Theoretical studies of combustion have mostly focused on low speed deflagrations as well as pure supersonic detonations; however, transonic conditions, have not been as widely studied. Transonic flow conditions naturally occur in a variety of applications such as: turbine engines, power generation systems, missiles, turbo-machinery, industrial reactors and many more. A small disturbance and irrotational model is proposed in [11]. The current work investigates relaxing the small disturbance assumption. In addition, the common practice is to treat the pressure equation as elliptic (pressure correction methods, pressure projection methods). Instead, a hyperbolic equation for the pressure is used in this study for unsteady problems. The current work investigates integration of the hyperbolic pressure equation as part of a reactive full potential flow model.

Keywords

Transonic flows potential flows combustion

1. Introduction

Increasingly computer simulations are being used in engineering, and simulations involving the reactive Euler equations are expensive and simplified models can reduce the calculations a design team needs to do. For this reason small disturbance, or potential models, which may be permissible under certain operating conditions are desirable. Rusak et al. demonstrated that steady state solutions of the compressible reacting flow problem with detonations behind the shock waves can be found. The model was used to study transonic combustion at various inlet Mach numbers, amounts of incoming reactant mass, reaction rates, and channel geometries [11]. The present work has two principal goals: (1) to extend the small disturbance combustion model proposed in [11] for transonic full potential flows, and (2) to develop a pressure-based unsteady Euler model that is capable of handling the effect of chemical reactions. Thus, in the present work we investigate a natural alternative model that fits between a reactive small disturbance model and a full reactive Euler model.

2. Combustion equations

The general governing equations, including viscous effects, are the two-dimensional compressible reactive flow Equations for a premixed chemical fuel species with Arrhenius type kinetics:

$\displaystyle\rho_{t}+(\rho u)_{x}+(\rho v)_{y}=0,$ $\displaystyle(\rho u)_{t}+(\rho u^{2})_{x}+(\rho vu)_{y}+p_{x}=\mu\left(u_{xx}% +u_{yy}+\frac{1}{3}\left(u_{x}+v_{y}\right)_{x}\right),$ $\displaystyle(\rho v)_{t}+(\rho uv)_{x}+(\rho v^{2})_{y}+p_{y}=\mu\left(v_{xx}% +v_{yy}+\frac{1}{3}\left(u_{x}+v_{y}\right)_{y}\right),$ $\displaystyle(\rho C_{p}T)_{t}+\left(\rho u\left[T+\frac{1}{2}\left(u^{2}+v^{2% }\right)\right]\right)_{x}+\left(\rho v\left[T+\frac{1}{2}\left(u^{2}+v^{2}% \right)\right]\right)_{y}$ $\displaystyle\qquad=\lambda\Delta T+\Phi_{\mu}+\Phi_{D}+Bw,$ $\displaystyle\left(\rho Y\right)_{t}+\left(\rho uY\right)_{x}+\left(\rho vY% \right)_{y}=D\Delta Y-w,$ $\displaystyle p=\rho R_{gas}T,$ $\displaystyle w=D_{a}\rho Ye^{-\frac{E_{a}}{R_{gas}T}}.$ (1)

Where $\mu$ is the viscosity, $\lambda$ is the thermal conductivity, $A$ is the pre-exponential reaction factor, and $D$ is the diffusivity coefficient, $\Phi_{\mu}$ is the viscous dissipation and $\Phi_{D}$ is the mass diffusion, and $B$ is the specific heat release. The choice of premixed chemistry with Arrhenius kinetics is to emphasize the effect of heat release with realistic reaction rates coupled to the flow. $Y$ is a reaction progress variable, $0\leqslant Y\leqslant 1$ , where $Y=0$ means the reactant is fully consumed, and $Y=1$ implies no reaction has yet taken place.

2.1 Nondimensionalization

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

$\displaystyle M_{\infty}=\frac{u^{2}_{\infty}}{\gamma R_{gas}T_{\infty}},\rho=% \frac{[\rho]}{\rho_{\infty}},p=\frac{[p]}{\rho_{\infty}u^{2}_{\infty}}=\frac{[% p]}{p_{\infty}\gamma M^{2}_{\infty}},$

(2) $\displaystyle T=\frac{[T]}{\left(V^{2}_{\infty}/C_{p}\right)},u=\frac{[u]}{u_{% \infty}},v=\frac{[v]}{v_{\infty}},Y=\frac{[Y]}{Y_{\infty}},x=\frac{[x]}{L},y=% \frac{[y]}{L}.$

Some authors non-dimensionalize $T=\frac{[T]}{T_{\infty}}$ , a conversion between this non-dimensionalization and ours is given by:

$\displaystyle T=\frac{[T]}{\left(V^{2}_{\infty}/C_{p}\right)}=\frac{[T]}{T_{% \infty}}\frac{1}{\gamma-1}\frac{1}{M^{2}_{\infty}}.$ (3)

The dimensionless parameters are given by:

$\displaystyle M^{2}_{\infty}=\frac{u^{2}_{\infty}}{\gamma R_{gas}T_{\infty}},A% =\frac{D_{a}}{u_{\infty}},\theta=\frac{Ea}{C_{p}T_{\infty}}\frac{1}{\gamma-1}% \frac{1}{M^{2}_{\infty}},\beta=\frac{[B]}{u^{2}_{\infty}}=\frac{[B]}{C_{p}T_{% \infty}}\frac{1}{\gamma-1}\frac{1}{M^{2}_{\infty}}.$ (4)

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

$\displaystyle\rho_{t}+(\rho u)_{x}+(\rho v)_{y}=0,\ \ \ f:=(\rho u)_{t}+(\rho u% ^{2})_{x}+(\rho vu)_{y}=-p_{x},$ $\displaystyle g:=(\rho v)_{t}+(\rho uv)_{x}+(\rho v^{2})_{y}=-p_{y},$ $\displaystyle\left(\rho T\right)_{t}+\left(\rho u\left[T+\frac{1}{2}\left(u^{2% }+v^{2}\right)\right]\right)_{x}+\left(\rho v\left[T+\frac{1}{2}\left(u^{2}+v^% {2}\right)\right]\right)_{y}=\frac{1}{Pr\cdot Re}\Delta T+\beta w,$ $\displaystyle\left(\rho Y\right)_{t}+\left(\rho uY\right)_{x}+\left(\rho vY% \right)_{y}=\frac{1}{Sc\cdot Re}\left(Y_{xx}+Y_{yy}\right)-w,$ $\displaystyle p=\rho T\frac{\gamma-1}{\gamma},$ $\displaystyle w=A\rho Ye^{-\frac{\theta}{T}}.$ (5)

Where $Pr=\frac{\mu C_{p}}{k}$ is the Prandtl Number, and $Sc=\frac{\mu}{\rho 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:

$\displaystyle u=\phi_{x}+\frac{\psi_{y}}{\rho},v=\phi_{y}-\frac{\psi_{x}}{\rho}.$ (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):

$\displaystyle(\rho\phi_{x})_{x}+(\rho\phi_{y})_{y}+\rho_{t}=0,\left(\frac{\psi% _{x}}{\rho}\right)_{x}+\left(\frac{\psi_{y}}{\rho}\right)_{y}=-\omega,$ $\displaystyle(\rho u\omega)_{x}+(\rho v\omega)_{y}=F(\rho,\omega),f_{x}+g_{y}=% \rho_{tt}-(p_{xx}+p_{yy}),$ (7) $\displaystyle\rho T_{t}+\rho u\left[{T+\frac{1}{2}\left(u^{2}+v^{2}\right)}% \right]_{x}+\rho v\left[T+\frac{1}{2}\left(u^{2}+v^{2}\right)\right]_{y}=\frac% {1}{Pr\cdot Re}\Delta 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:

$\displaystyle-\nabla\cdot\nabla p=\frac{\partial f}{\partial x}+\frac{\partial g% }{\partial y}=\nabla\cdot\vec{F}+\frac{\partial\nabla\cdot\rho u}{\partial t}=% \nabla\cdot\vec{F}+\frac{\partial^{2}\rho}{\partial 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. $\omega=\nabla\times\vec{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 $\phi$ must exist such that $\vec{q}=\nabla\phi$ . Under this assumption the velocity field is curl-free and reduces to: $u=\phi_{x}$ , $v=\phi_{y}$ . An irrotational unsteady compressible reactive full potential flow model is given by:

$\displaystyle(\rho\phi_{x})_{x}+(\rho\phi_{y})_{y}+\rho_{t}=0,$ $\displaystyle\left[(\rho u^{2})_{x}+(\rho uv)_{y}\right]_{x}+\left[(\rho uv)_{% x}+(\rho v^{2})_{y}\right]_{y}-\frac{1}{Pr\cdot Re}\Delta\left((\rho u)_{x}+(% \rho v)_{y}\right)=\rho_{tt}-(p_{xx}+p_{yy}),$ $\displaystyle\rho T_{t}+\rho u\left[T+\frac{1}{2}\left(u^{2}+v^{2}\right)% \right]_{x}+\rho v\left[T+\frac{1}{2}\left(u^{2}+v^{2}\right)\right]_{y}=\frac% {1}{Pr\cdot Re}\Delta T+\beta w,$ $\displaystyle\rho Y_{t}+\rho uY_{x}+\rho vY_{y}=\frac{1}{Sc\cdot Re}\Delta Y-w,$ $\displaystyle p=\rho T\frac{\gamma-1}{\gamma},$ $\displaystyle w=A\rho Ye^{-\frac{\theta}{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:

$\displaystyle\rho_{tt}=(p_{xx}+p_{yy})+f_{x}+g_{y}\longrightarrow\frac{\gamma}% {\gamma-1}\left(\frac{p}{T}\right)_{tt}=(p_{xx}+p_{yy})+f_{x}+g_{y}.$ (10)

For the present work we simply choose:

$\displaystyle\frac{1}{Sc\cdot Re}=\frac{1}{Pr\cdot Re}=O(\Delta x).$ (11)

4. Reactive transonic small disturbance model

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

$\displaystyle\left[1-M^{2}_{\infty}-(\gamma+1)M^{2}_{\infty}\varphi_{x}\right]% \frac{\partial^{2}\varphi}{\partial x^{2}}+\frac{\partial^{2}\varphi}{\partial y% ^{2}}=\tilde{\beta}K_{\delta}\tilde{A}\tilde{Y}e^{-\tilde{\theta}/(1+\tilde{T}% )},$

(12) $\displaystyle\tilde{T}=-(\gamma-1)\varphi_{{x}},\tilde{Y}_{x}=\tilde{A}\tilde{% Y}e^{-\tilde{\theta}/(1+\tilde{T})}.$

This governing equation is of mixed type, elliptic and hyperbolic. Here $\tilde{T},\tilde{Y}$ denote small disturbance variables, and $\tilde{A},\ \tilde{\beta},\ K_{\delta},$ and $\tilde{\theta}$ are constant reaction control parameters, and $M_{\infty}$ is the inlet Mach number. The velocity can be calculated from: $\tilde{u}=\varphi_{x},\tilde{v}=\varphi_{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\leqslant y\leqslant 1$ ) are:

$\displaystyle u(0,y)=1+\alpha u_{0}(y),v_{x}(0,y)=0,T(0,y)=1+\alpha T_{0}(y),Y% (0,y)=\delta Y_{0}(y),$ (13)

where $0\leqslant\alpha<<1$ . We will take $\alpha=0$ in this work. At the channel outlet ( $y_{\textit{lower}}\leqslant y\leqslant y_{\textit{upper}}$ ):

$\displaystyle v(x_{\textit{out}},y)=0,\frac{\partial T(x_{\textit{out}},y)}{% \partial y}=0,\frac{\partial Y(x_{\textit{out}},y)}{\partial y}=0.$ (14)

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.

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.

Figure 3.
Small Disturbance Model (TSD) T and Y profiles for $M_{a}=$ 0.80, $A=$ 20000 for range of channel divergence parameters $\epsilon_{\textit{wall}}=$ 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:

$\displaystyle v\approx\frac{dG_{\textit{body}}}{dx}.$ (15)

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

$\displaystyle G_{\textit{body}}(x)=\left\{{{\begin{array}[]{ll}1,&\text{for }x% <1\\ \frac{1}{2}\left[2+\epsilon_{\textit{wall}}\left(1+\cos\left(\pi x\right)% \right)\right],&\text{for }1\leqslant x\leqslant 2\\ 1+\epsilon_{\textit{wall}},&\text{for }2<x\leqslant 3.\end{array}}}\right.$ (16)

Typically $\epsilon_{\textit{wall}}=0.12=12\%$ , unless otherwise specified, as in Fig. 3. We assume a symmetric flow, at the axis of symmetry, where $y=0$ :

$\displaystyle v(x,0)=0,T_{y}(x,0)=0,Y_{y}(x,0)=0.$ (17)

In these calculations: $\tilde{\beta}K_{\delta}=2$ , $\tilde{\theta}=12$ , $\tilde{A}=\epsilon^{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:

$\displaystyle M_{a}=M_{\infty}\sqrt{\left(1+(\gamma+1)\tilde{u}\right)}.$ (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 $\epsilon$ 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.

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.

$\displaystyle(\rho\phi_{x})_{x}+(\rho\phi_{y})_{y}=0\longrightarrow(\bar{\rho}% \phi_{x})_{x}+(\bar{\rho}\phi_{y})_{y}=0.$ (19)

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

$\displaystyle\bar{\rho}_{(i+1/2,j+1/2)}=\rho_{(i+1/2,j+1/2)}-\left[\frac{u}{q^% {2}}\Delta x\frac{\partial F}{\partial x}+\frac{v}{q^{2}}\Delta y\frac{% \partial F}{\partial y}\right]_{(i+1/2,j+1/2)}.$ (20)

After calculating the local Mach number at the cell centers:

$\displaystyle M^{2}_{(i+1/2,j+1/2)}=\frac{q^{2}}{a^{2}}=\frac{u^{2}+v^{2}}{% \frac{1}{M^{2}_{\infty}}+\frac{1}{2}(\gamma-1)\left[1-u^{2}-v^{2}\right]}.$ (21)

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

$\displaystyle F=\left\{{\begin{array}[]{ll}\rho q,&\text{for }M>1\\ \rho^{*}q^{*}=\rho^{*}a^{*},&\text{for }M\leqslant 1\end{array}}\right.$ (22)

Where the following definitions are used:

$\displaystyle a^{*}=\left(2\frac{\gamma-1}{\gamma+1}\left[\frac{1}{\gamma-1}% \frac{1}{M^{2}_{\infty}}+\frac{1}{2}\right]\right)^{1/2}$ (23) $\displaystyle\rho^{*}=\left(1-\frac{M^{2}_{\infty}\left(\gamma-1\right)}{2}((a% ^{*})^{2}-1)\right)^{1/(\gamma-1)}.$ (24)

Figure 5.

Staggered grid showing the location of the duct, centerline and model variables.

Figure 6.

Compressible nozzle flow 1D validation test case. Quasi-1D Euler equations (left) and the quasi-1D full potential model (right).

Figure 7.

$\rho$ -biased full potential model for $M_{a}=$ 0.75, 0.80, 0.85, 0.90, 0.95 for $A=$ 0, showing transition to shocked flow at 12% nozzle area variation parameter.

Figure 8.

Flux Biasing and $\rho$ -Biasing vs. TSD with $A=$ 0 for $M_{a}=$ 0.8 (left) and for $M_{a}=$ 0.95 (right).

Figure 9.

$\rho$ -biased full potential model with $M_{\infty}=$ 0.8 for $A=$ 0, 10000, 20000, 30000.

Figure 10.

$\rho$ -biased full potential model with $M_{\infty}=$ 0.85 for $A=$ 0, 10000, 20000, 30000.

Figure 11.

$\rho$ -biased full potential model for $M_{\infty}=$ 0.80, 0.85, 0.90, 0.95 with $A=$ 20000.

Figure 12.

$\rho$ -biased unsteady full potential model for $M_{\infty}=$ 0.80 with $A=$ 0 validating the unsteady model converges to the same steady state solution.

Figure 13.

$\rho$ -biased unsteady full potential model for $M_{\infty}=$ 0.85 with $A=$ 20000 showing the steady state and unsteady solution overlayed.

Figure 14.

Unsteady peak Mach number (right) and minimum Y value (left) along the wall for $M_{\infty}=$ 0.85 with $A=$ 0, 20000. It is shown that the unsteady code converges from quiescent initial conditions to the same solution as the steady code.

In the density biasing scheme, instead of biasing the flux $F:=\rho q$ , now the density alone is biased:

$\displaystyle\rho_{i+1/2,j+1/2}=\rho_{i+1/2,j+1/2}-\mu\left(\rho_{i+1/2,j+1/2}% -\rho_{i-1/2,j+1/2}\right).$ (25)

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

$\displaystyle\mu=\tau\cdot\max\left(0,1-\frac{M^{2}_{c}}{M^{2}_{i+1/2,j+1/2}},% 1-\frac{M^{2}_{c}}{M^{2}_{i-1/2,j+1/2}}\right).$ (26)

Where $\tau$ , $M^{2}_{c}$ 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$ :

$\displaystyle-\left[(\rho u^{2})_{x}+(\rho uv)_{y}\right]_{x}-\left[(\rho uv)_% {x}+(\rho v^{2})_{y}\right]_{y}=(p_{xx}+p_{yy}),$ (27)

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

$\displaystyle\rho^{n+1}=\rho^{n}+\omega_{\textit{SOR}}\left(\frac{p}{T}\frac{% \gamma}{\gamma-1}-\rho^{n}\right).$ (28)

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

For the case of unsteady computations a stability term of the form: $\nu\rho_{t}$ , with $0<\nu<<1$ , is added to the pressure equation to help with convergence:

$\displaystyle\left[(\rho u^{2})_{x}+(\rho uv)_{y}\right]_{x}+\left[(\rho uv)_{% x}+(\rho v^{2})_{y}\right]_{y}=\rho_{tt}+\nu\rho_{t}-(p_{xx}+p_{yy}).$ (29)

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

$\displaystyle\nu\Delta\left[(\rho u)_{x}+(\rho v)_{y}\right]\approx-\nu\rho_{t}.$ (30)

Where $0\leqslant\nu\leqslant 1$ . One way to proceed is that the pressure equation is solved for $\rho$ 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:

$\displaystyle\frac{\gamma}{\gamma-1}\left[\frac{\left(\frac{P}{T}\right)^{n+1}% -2\left(\frac{P}{T}\right)^{n}+\left(\frac{P}{T}\right)^{n-1}}{\Delta t^{2}}+% \nu\frac{\left(\frac{P}{T}\right)^{n+1}-\left(\frac{P}{T}\right)^{n-1}}{2% \Delta t}\right]-(p_{xx}+p_{yy})^{n+1}$ $\displaystyle\quad=f^{n}_{x}+g^{n}_{y}.$ (31)

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

6. 1D Euler vs. 1D unsteady full potential

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:

$\displaystyle E=\frac{T}{\gamma}+\frac{1}{2}u^{2}.$ (32)

Then the energy equation becomes:

$\displaystyle(A\rho E)_{t}+(A\rho uH)_{x}=0,$ $\displaystyle\left(A\rho\left(\frac{T}{\gamma}+\frac{1}{2}u^{2}\right)\right)_% {t}+\left(A\rho u\left(T+\frac{1}{2}u^{2}\right)\right)_{x}=\epsilon(A\rho E)_% {xx},$ (33) $\displaystyle(A\rho T)_{t}+\frac{\gamma}{2}\left(\frac{(A\rho u)^{2}}{A\rho}% \right)_{t}+\gamma\left(A\rho u\left(T+\frac{1}{2}u^{2}\right)\right)_{x}=% \gamma\epsilon(A\rho E)_{xx}.$

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

$\displaystyle-(A\rho)_{tt}=(A\rho u)_{xt}-\epsilon(A\rho)_{\textit{xxt}},$ (34) $\displaystyle(Ap)_{xx}=-(A\rho u)_{tx}-(A\rho u^{2})_{xx}+(A_{x}p)_{x}+% \epsilon(A\rho u)_{\textit{xxx}}.$

Adding these equations together implies:

$\displaystyle(Ap)_{xx}-(A\rho)_{tt}=-(A\rho u^{2})_{xx}+(A_{x}p)_{x}-\epsilon(% A\rho)_{\textit{xxt}}+\epsilon(A\rho u)_{\textit{xxx}}.$

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

$\displaystyle(A\rho)_{tt}=(Ap)_{xx}+\left(A\rho u^{2}\right)_{xx}-\left(\frac{% A_{x}p}{\rho}\right)_{x}+\epsilon(A\rho)_{\textit{xxt}}-\epsilon(A\rho u)_{% \textit{xxx}},$ $\displaystyle(A\rho u)_{x}=-(A\rho)_{t}+\epsilon((A\rho u)_{xx}+(A\rho)_{xx}),$ $\displaystyle(Ap)_{t}=-\frac{\gamma-1}{2}\left(A\rho u^{2}\right)_{t}-(\gamma-% 1)\left(A\rho u\left(T+\frac{1}{2}u^{2}\right)\right)_{x}(\gamma-1)\epsilon(Ap% )_{xx},$ $\displaystyle T=\frac{\gamma}{\gamma-1}\frac{p}{\rho}.$ (35)

This system is hyperbolic and converges with $\Delta t=O(\Delta x)$ . Alternatively, one can solve the continuity equation for a potential function $\phi$ , where $u=\phi_{x}$ , without loss of generality as is proposed earlier.

$\displaystyle(A\rho\phi_{x})_{x}+(A\rho)_{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].

7. Results of the present two-dimensional full potential model

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_{\infty}=$ 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_{\infty}=$ 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 $\rho$ -biased full potential model for $M_{\infty}=$ 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_{\infty}=$ 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_{\infty}=$ 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_{\infty}=$ 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.

References

Buckmaster

and Ludford

, Theory of Laminar Flames, Cambridge University, 1982.

Chattot

J.J.

, Computational Aerodynamics and Fluid Dynamics, Springer-Verlag Berlin Heidelberg New York, 2002.

Cole

J.D.

and Cook

L.P.

, Transonic Aerodynamics, Elsevier Science Publishers B.V., 1986.

Hafez

M.M.

South

and Murman

E.M.

, Artificial Compressibility Methods For Numerical Solution of Transonic Full Potential Equation, AIAA Journal 17(8) (1979), 838–844.

Hafez

M.M.

Habashi

W.G.

and Kotiuga

P.L.

, Conservative calculations of non-isentropic transonic flows, International Journal for Numerical Methods in Fluids 5 (1985), 1047–1057.

Holst

T.L.

, Numerical Computation of Transonic FLow Governed by the Full-Potential Equation, (N83-19708) NASA Technical Memorandum 84310, 1983.

Kuo

K.K.

, Principles of Combustion, John Wiley & Sons, Inc, 2005.

Liepmann

H.W.

and Roshko

, Elements of Gasdynamics, Dover, 2012.

Murman

E.M.

and Cole

J.D.

, Calculation of Plane Steady Transonic flow, AIAA S. 9(1) (1971), 114–121.

10.

Osher

Hafez

M.M.

and Whitlow

, Jr., Entropy Condition Satisfying Approximations for the Full Potential Equation of Transonic Flow, Mathematics of Computation 44(169) (1985), 1–29.

11.

Rusak

Lee

J.C.

and Choi

J.J.

, A small-disturbance model of transonic combustion, Combustion Theory and Modeling 12(1) (2007), 93–113.

12.

Rusak

and Lee

J.C.

, Parametric Investigation of Combustion in Compressible Flows, AIAA-2005-4805: 4th AIAA Theoretical Fluid Mechanics Meeting, 2005.

13.

Volpe

and Jameson

, Transonic Potential Flow Calculations by Two Artificial Density Methods, AIAA Journal 26(4) (1988).