# Numerical and Experimental Investigation of a Supersonic Flow Field around Solid Fuel on an Inclined Flat Plate

Authors: Uzu-Kuei Hsu
Journal: Modelling and Simulation in Engineering
 (2009)

---

## Abstract

This research adopts a shock tube 16 meters long and with a 9 cm bore to create a supersonic, high-temperature, and high-pressure flowfield to observe the gasification and ignition of HTPB solid fuel under different environments.  Also, full-scale 3D numerical simulation is executed to enhance the comprehension of this complex phenomenon.  The CFD (Computational Fluid Dynamics) code is based on the control volume method and the pre-conditioning method for solving the Navier-Stokes equations to simulate the compressible and incompressible coupling problem. In the tests, a HTPB slab is placed in the windowed-test section.  Various test conditions generate different supersonic Mach numbers and environmental temperatures. In addition, the incident angles of the HTPB slab were changed relative to the incoming shock wave.  Results show that as the Mach number around the slab section exceeded 1.25, the flowfield temperature achieved 1100 K, which is higher than the HTPB gasification temperature (930 K~ 1090 K).  Then, gasification occurred and a short-period ignition could be observed.  In particular, when the slab angle was 7∘, the phenomenon became more visible.  This is due to the flow field temperature increase when the slab angle was at 7∘.


---

## Body

## 1. Introduction

With  the development of the space shuttle and solar system exploration, hypersonic high technology in aviation will play an important role in the next-generation frontier [1]. However, if the engine of a hypersonic vehicle is a scramjet with the combustor inlet air supersonic, it will be much more useful and powerful than a ramjet [2].  In the 1960s, researchers indicated that the regression ratio is the key in mixer-rocket studies. For this reason, many test models have been developed in different combustion conditions. Marxman and Gilbert [3] considered that the optimal position of the flame should be at the top of the fuel surface, and the regression rate is the minimum in the turbulence layer at about 10–20%. Muzzy [4] also verified that the heat convection effect is very important in fuel consumption. Smoot and Price [5] indicated that the fuel regression ratio is proportional to the oxidizer flow rate and is in the 0.8 order in the lower oxidizer flow rate. During the 1990s, many basic studies were conducted. Greiner and Grederick  [6] expressed the fuel regression ratio as proportional to the oxidizer flow rate, and the pressure fluctuation decreased when raising the mixed region length. Chiaverini [7] indicated that the vapor temperature of the fuel surface was between 930 K and 1190 K, depending on the properties of several different HTPB composites.In order to investigate the ignition and combustion efficiency of a supersonic combustion ramjet and simplify the components, a 16-meter shock tube was established, as shown in Figure1. The device consists of a long tube divided into a high pressure and a low pressure section by an aluminum diaphragm. When the diaphragm in the high pressure driver section ruptured, a series of compression waves coalesced into a single shock front which compressed and heated to a high pressure the gas in the low pressure region, and created the supersonic gas flow condition. In the shock tube, there are many complex phenomena including the normal shock, contact surface, expansion wave, and so forth. We investigated the flowfield of the supersonic flow through the plate-like HTPB solid fuel was investigated under this unsteady condition. The difficult problem, however, was that it is not easy to create a supersonic condition for a long time. The best test period was about 10 milliseconds. Therefore, care had to be taken in experiments and operation. There have been several shock tube wind tunnels established for research since 1950, for example, Glass and Hall [8], Lukasiewicz [9], Nagamatsu [10], Bradley [11], Soloukhin [12], and so forth. In fact, the wind tunnel test is very important in classical aerodynamics. However, there is not a complete understanding of the full phenomena because of the limits of the experiment. Computational Fluid Dynamics (CFD) is a good tool to deal with the problems. In this study, both of these two methods are used. Using CFD simulates the supersonic flow through the HTPB slab, and its results are compared with the experimental data.Shock tube dimension (a)  high pressure sections (length: 290 cm, diameter: 28.5 cm), (b) convergent nozzle space (length: 10 cm, internal diameter in divergent section: 28.5 cm, internal diameter in convergent section: 9 cm), (c) test section (window:25cm×3cm), (d) dimension.
(a)
(b)
(c)
(d)


## 2. Experimental Apparatus

In Figure1, the length of the shock tube is 16 m. The high pressure region is higher than 147 psi, and the low pressure region is below 1 psi. An HTPB slab, 15cm (length)×3cm (spread)×0.5cm (thickness), is placed at 14.55 m from the starting point of the high pressure region. There are two pressure transducers placed at 14.55 m (no. 1) and 14.65 m (no. 2), separately.  The tube’s initial temperature is 300 K. In order to produce a snapshot, a 10 W pulse laser (t=50μs) and high speed CCD camera (t=200μs) are set up.


## 3. Numerical Model

### 3.1. Geometry and Grid System

Uniformly spaced grids are used to cover the flow field, and the stretching transformation clusters, using the Roberts generalized stretching transformation technique, are made near the boundary layer. The shock tube is symmetrical about the centre-plane and, therefore, only the right half of the shock tube and plate-like model needs to be modelled. The multiblock grid approach is used in the present study. The total number of cells was 2 094 750 with respect to half of the 3D shock tube, as shown in Figure2. The instantaneous solution was obtained by solving the time-dependent governing equations, and the residual was measured by the order of magnitude of the decay. The convergent solution was achieved when the residual had decayed by about 3 orders of magnitudes. Computation was performed on finer and coarser grids for a constant grid resolution, it, being found that the total grid sizes, especially in the y direction, depend on the turbulence model used. According to the present study, the average value of the first-layer cell, y+, closest to the boundary surface was 0.2 with the exact solution for turbulence models [13]. The DELL OPTIPLX GX270 workstation was used for the computation.Figure 2
Longitudinal cut view of the grid system.


### 3.2. Governing Equations

The Navier-Stokes equations in integral form for an arbitrary control volumeΩ with differential surface area dA⃗ can be expressed as(1)
∂
∂
t
∭
W
d
Ω
+
∬
[
F
⃗
-
G
⃗
]
·
d
A
⃗
=
0
,
where W is the dependent vector of conservative variables and F⃗, G⃗ are the inviscid and viscous flux vectors in standard conservation form(2)
W
=
[
ρ
ρ
u
ρ
v
ρ
w
ρ
E
]
,
F
⃗
=
[
ρ
V
⃗
ρ
V
⃗
u
+
p
i
⃗
ρ
V
⃗
v
+
p
j
⃗
ρ
V
⃗
w
+
p
k
⃗
ρ
V
⃗
E
+
p
V
⃗
]
,
G
⃗
=
[
0
τ
x
i
τ
y
i
τ
z
i
τ
i
j
v
j
+
q
]
,
where ρ,V⃗,E, and p are the density, velocity vector, total energy per unit mass, and pressure, respectively. The term τ is the viscous stress tensor, q is the heat flux vector. The total energy is related to the total enthalpy H by E=H-p/ρ, where H=h+V2/2, h=CpT, and T is the temperature. The system is closed by the equation of state, which typically has a form of  ρ=ρ(p,T).Now two of the newer low-diffusion flux-splitting methods, AUSM+ and AUSMDV, are presented. For both methods, the inviscid interface fluxFi+1/2 in the x direction is split into a convective contribution F1/2c plus a pressure contribution F1/2p. The convective term is defined as(3)
F
1
/
2
c
=
(
ρ
u
)
1
/
2
[
1
u
v
w
H
]
i
/
i
+
1
,
where state i is chosen for the column vector (1,u,v,w,H)T if the interface mass flux (ρu)1/2 is nonnegative and state i+1 is chosen if (ρu)1/2 is negative. The pressure term is defined as(4)
F
1
/
2
p
=
(
ρ
u
)
1
/
2
[
0
p
1
/
2
0
0
0
]
.The interface quantities(ρu)1/2  and p1/2 are expressed in terms of sets of polynomials in Mach number, defined as(5)
M
i
=
u
i
a
1
/
2
,
where a1/2 is an interface speed of sound. These sets of polynomials are required:(6)
M
(
1
)
±
=
1
2
(
M
±
|
M
|
)
,
M
(
4
)
±
=
{
±
1
4
(
M
±
1
)
2
±
1
8
(
M
2
-
1
)
2
,
|
M
|
<
1
,
1
2
(
M
±
|
M
|
)
,
otherwise,
P
(
5
)
±
=
{
1
4
(
M
±
1
)
2
(
2
∓
M
)
±
3
16
M
(
M
2
-
1
)
2
,
|
M
|
<
1
,
1
2
M
(
M
+
|
M
|
)
,
otherwise.The numerals in the subscripts ofM and P indicate the degree of the polynomials. With these, the interface quantities are defined as(7)
(
ρ
u
)
1
/
2
=
a
1
/
2
(
ρ
i
m
1
/
2
+
+
ρ
i
+
1
m
1
/
2
-
)
,
p
1
/
2
=
P
(
5
)
+
(
M
i
)
ρ
i
+
P
(
5
)
-
(
M
i
+
1
)
ρ
i
+
1
.
The differences between AUSM+ and AUSMDV are due to the definitions of m1/2±.For AUSM+,(8)
m
1
/
2
=
O
(
4
)
′
+
(
M
i
)
+
O
(
4
)
′
-
(
M
i
+
1
)
,
m
1
/
2
-
=
1
2
(
m
1
/
2
±
|
m
1
/
2
|
)
.For AUSMDV,(9)
m
1
/
2
+
=
o
(
1
)
′
+
(
M
i
)
+
ω
+
o
(
4
)
′
+
(
M
i
)
-
ω
+
o
(
1
)
′
+
(
M
i
)
,
m
1
/
2
-
=
o
(
1
)
′
-
(
M
i
+
1
)
+
ω
-
o
(
4
)
′
-
(
M
i
+
1
)
-
ω
-
o
(
1
)
′
-
(
M
i
+
1
)
,
where(10)
ω
+
=
(
p
i
ρ
i
)
/
[
1
2
(
p
i
ρ
i
+
p
i
+
1
ρ
i
+
1
)
]
,
ω
-
=
(
p
i
+
1
ρ
i
+
1
)
/
[
1
2
(
p
i
ρ
i
+
p
i
+
1
ρ
i
+
1
)
]
.As shown later, the inclusion ofω± in the AUSMDV flux splitting scheme includes a pressure-diffusion term in the interface mass flux. The presence of this term ensures monotonicity in the capturing of nongrid-aligned shock waves but can, under certain conditions, give rise to the carbuncle instability for bluff-body flows. The mass flux representation for AUSM+ does not contain a pressure-diffusion term and is not susceptible to the instability.  To save the monotonicity in the capturing of nongrid-aligned shock waves only the AUSMDV scheme was considered.


### 3.3. Preconditioning System

The original system of equations from the conservative variablesW can be transformed into the primitive variables Q=(p,u,v,w,T)T as follows:(11)
Γ
∂
∂
t
∭
Q
d
V
+
∬
[
F
⃗
-
G
⃗
]
·
d
A
⃗
=
0
.
The precoinditioning matrix Γ is based on that proposed by Choi and Merkle [14] and extended further by Weiss and Smith [15]. Given that(12)
U
ref
2
=
min
⌊
a
2
,
max
(
|
V
|
2
,
K
|
V
∞
|
2
)
⌋
,
Θ
=
1
U
ref
2
-
1
a
2
,
H
=
C
p
T
+
1
2
(
u
2
+
v
2
+
w
2
)
,
where |V| is the local velocity magnitude, |V∞| is a fixed reference velocity, a is the sound speed, and  K is a constant; the preconditioning matrix Γ takes the following form:(13)
Γ
=
[
Θ
+
1
R
T
0
0
0
-
ρ
T
u
(
Θ
+
1
R
T
)
ρ
0
0
-
ρ
u
T
v
(
Θ
+
1
R
T
)
0
ρ
0
-
ρ
v
T
w
(
Θ
+
1
R
T
)
0
0
ρ
-
ρ
w
T
H
(
Θ
+
1
R
T
)
-
1
ρ
u
ρ
v
ρ
w
ρ
[
C
p
-
(
H
T
)
]
]
.As indicated, the Weiss-Smith preconditioner is formed by the addition of the vectorΘ[1,u,v,w,T]T to the Jacobian matrix ∂W/∂Q, where W is the vector of conservative variables. The choice of the reference velocity Uref is crucial for success of the method. If the magnitude of the reference velocity is equal to the local speed of sound, the system (11) becomes identical to the nonpreconditioned system (1). It is our goal to make all eigenvalues of the same order of magnitude, so the reference velocity must be of the same order as a local velocity. The constant value of the reference velocity through entire flowfield can be used. However it is attractive to compute the reference velocity locally, because it is more flexible for the flows with significantly different speeds. In the present study, K is fixed at 0.25. The eigenvalues of Γ-1A(A=∂F/∂W) are  u, u′±a′, where u is the velocity component in the x direction and(14)
u
′
±
a
′
=
1
2
[
(
1
+
M
ref
2
)
u
±
a
(
1
-
M
ref
2
)
2
M
2
+
4
M
ref
2
]
,
M
ref
2
=
U
ref
/
a
2
2
,
M
=
u
a
.


### 3.4. High-Resolution Scheme

The numerical scheme, using the preconditioning finite volume method, was introduced to solve the governing flow equations. A second-order scheme was initially applied, so the left and right states were chosen to be the cell average values on the left and right of the cell faces. In a high-resolution scheme, in order to raise the order of accuracy of upwind differencing, all that is needed is to raise the order of accuracy of the initial-value interpolation that yields the zone-boundary data. Such schemes are labelled as high-resolution schemes as opposed to Total Variation Diminishing (TVD) schemes, which completely eliminate any of those spurious oscillations when applied to one dimensional nonlinear hyperbolic conservation laws and linear hyperbolic systems. The van Leer kappa-scheme, in which the kappa number is one-third, was selected to obtain the high-resolution upwind differencing [16–18].An optimal multistage scheme was used for the time integration, and the multistage coefficients were modified by Tai et al. [19], and redefined using the Courant number for multidimensional use. Also, a residual smoothing method was imposed to accelerate convergence and to improve numerical stability.


## 3.1. Geometry and Grid System

Uniformly spaced grids are used to cover the flow field, and the stretching transformation clusters, using the Roberts generalized stretching transformation technique, are made near the boundary layer. The shock tube is symmetrical about the centre-plane and, therefore, only the right half of the shock tube and plate-like model needs to be modelled. The multiblock grid approach is used in the present study. The total number of cells was 2 094 750 with respect to half of the 3D shock tube, as shown in Figure2. The instantaneous solution was obtained by solving the time-dependent governing equations, and the residual was measured by the order of magnitude of the decay. The convergent solution was achieved when the residual had decayed by about 3 orders of magnitudes. Computation was performed on finer and coarser grids for a constant grid resolution, it, being found that the total grid sizes, especially in the y direction, depend on the turbulence model used. According to the present study, the average value of the first-layer cell, y+, closest to the boundary surface was 0.2 with the exact solution for turbulence models [13]. The DELL OPTIPLX GX270 workstation was used for the computation.Figure 2
Longitudinal cut view of the grid system.


## 3.2. Governing Equations

The Navier-Stokes equations in integral form for an arbitrary control volumeΩ with differential surface area dA⃗ can be expressed as(1)
∂
∂
t
∭
W
d
Ω
+
∬
[
F
⃗
-
G
⃗
]
·
d
A
⃗
=
0
,
where W is the dependent vector of conservative variables and F⃗, G⃗ are the inviscid and viscous flux vectors in standard conservation form(2)
W
=
[
ρ
ρ
u
ρ
v
ρ
w
ρ
E
]
,
F
⃗
=
[
ρ
V
⃗
ρ
V
⃗
u
+
p
i
⃗
ρ
V
⃗
v
+
p
j
⃗
ρ
V
⃗
w
+
p
k
⃗
ρ
V
⃗
E
+
p
V
⃗
]
,
G
⃗
=
[
0
τ
x
i
τ
y
i
τ
z
i
τ
i
j
v
j
+
q
]
,
where ρ,V⃗,E, and p are the density, velocity vector, total energy per unit mass, and pressure, respectively. The term τ is the viscous stress tensor, q is the heat flux vector. The total energy is related to the total enthalpy H by E=H-p/ρ, where H=h+V2/2, h=CpT, and T is the temperature. The system is closed by the equation of state, which typically has a form of  ρ=ρ(p,T).Now two of the newer low-diffusion flux-splitting methods, AUSM+ and AUSMDV, are presented. For both methods, the inviscid interface fluxFi+1/2 in the x direction is split into a convective contribution F1/2c plus a pressure contribution F1/2p. The convective term is defined as(3)
F
1
/
2
c
=
(
ρ
u
)
1
/
2
[
1
u
v
w
H
]
i
/
i
+
1
,
where state i is chosen for the column vector (1,u,v,w,H)T if the interface mass flux (ρu)1/2 is nonnegative and state i+1 is chosen if (ρu)1/2 is negative. The pressure term is defined as(4)
F
1
/
2
p
=
(
ρ
u
)
1
/
2
[
0
p
1
/
2
0
0
0
]
.The interface quantities(ρu)1/2  and p1/2 are expressed in terms of sets of polynomials in Mach number, defined as(5)
M
i
=
u
i
a
1
/
2
,
where a1/2 is an interface speed of sound. These sets of polynomials are required:(6)
M
(
1
)
±
=
1
2
(
M
±
|
M
|
)
,
M
(
4
)
±
=
{
±
1
4
(
M
±
1
)
2
±
1
8
(
M
2
-
1
)
2
,
|
M
|
<
1
,
1
2
(
M
±
|
M
|
)
,
otherwise,
P
(
5
)
±
=
{
1
4
(
M
±
1
)
2
(
2
∓
M
)
±
3
16
M
(
M
2
-
1
)
2
,
|
M
|
<
1
,
1
2
M
(
M
+
|
M
|
)
,
otherwise.The numerals in the subscripts ofM and P indicate the degree of the polynomials. With these, the interface quantities are defined as(7)
(
ρ
u
)
1
/
2
=
a
1
/
2
(
ρ
i
m
1
/
2
+
+
ρ
i
+
1
m
1
/
2
-
)
,
p
1
/
2
=
P
(
5
)
+
(
M
i
)
ρ
i
+
P
(
5
)
-
(
M
i
+
1
)
ρ
i
+
1
.
The differences between AUSM+ and AUSMDV are due to the definitions of m1/2±.For AUSM+,(8)
m
1
/
2
=
O
(
4
)
′
+
(
M
i
)
+
O
(
4
)
′
-
(
M
i
+
1
)
,
m
1
/
2
-
=
1
2
(
m
1
/
2
±
|
m
1
/
2
|
)
.For AUSMDV,(9)
m
1
/
2
+
=
o
(
1
)
′
+
(
M
i
)
+
ω
+
o
(
4
)
′
+
(
M
i
)
-
ω
+
o
(
1
)
′
+
(
M
i
)
,
m
1
/
2
-
=
o
(
1
)
′
-
(
M
i
+
1
)
+
ω
-
o
(
4
)
′
-
(
M
i
+
1
)
-
ω
-
o
(
1
)
′
-
(
M
i
+
1
)
,
where(10)
ω
+
=
(
p
i
ρ
i
)
/
[
1
2
(
p
i
ρ
i
+
p
i
+
1
ρ
i
+
1
)
]
,
ω
-
=
(
p
i
+
1
ρ
i
+
1
)
/
[
1
2
(
p
i
ρ
i
+
p
i
+
1
ρ
i
+
1
)
]
.As shown later, the inclusion ofω± in the AUSMDV flux splitting scheme includes a pressure-diffusion term in the interface mass flux. The presence of this term ensures monotonicity in the capturing of nongrid-aligned shock waves but can, under certain conditions, give rise to the carbuncle instability for bluff-body flows. The mass flux representation for AUSM+ does not contain a pressure-diffusion term and is not susceptible to the instability.  To save the monotonicity in the capturing of nongrid-aligned shock waves only the AUSMDV scheme was considered.


## 3.3. Preconditioning System

The original system of equations from the conservative variablesW can be transformed into the primitive variables Q=(p,u,v,w,T)T as follows:(11)
Γ
∂
∂
t
∭
Q
d
V
+
∬
[
F
⃗
-
G
⃗
]
·
d
A
⃗
=
0
.
The precoinditioning matrix Γ is based on that proposed by Choi and Merkle [14] and extended further by Weiss and Smith [15]. Given that(12)
U
ref
2
=
min
⌊
a
2
,
max
(
|
V
|
2
,
K
|
V
∞
|
2
)
⌋
,
Θ
=
1
U
ref
2
-
1
a
2
,
H
=
C
p
T
+
1
2
(
u
2
+
v
2
+
w
2
)
,
where |V| is the local velocity magnitude, |V∞| is a fixed reference velocity, a is the sound speed, and  K is a constant; the preconditioning matrix Γ takes the following form:(13)
Γ
=
[
Θ
+
1
R
T
0
0
0
-
ρ
T
u
(
Θ
+
1
R
T
)
ρ
0
0
-
ρ
u
T
v
(
Θ
+
1
R
T
)
0
ρ
0
-
ρ
v
T
w
(
Θ
+
1
R
T
)
0
0
ρ
-
ρ
w
T
H
(
Θ
+
1
R
T
)
-
1
ρ
u
ρ
v
ρ
w
ρ
[
C
p
-
(
H
T
)
]
]
.As indicated, the Weiss-Smith preconditioner is formed by the addition of the vectorΘ[1,u,v,w,T]T to the Jacobian matrix ∂W/∂Q, where W is the vector of conservative variables. The choice of the reference velocity Uref is crucial for success of the method. If the magnitude of the reference velocity is equal to the local speed of sound, the system (11) becomes identical to the nonpreconditioned system (1). It is our goal to make all eigenvalues of the same order of magnitude, so the reference velocity must be of the same order as a local velocity. The constant value of the reference velocity through entire flowfield can be used. However it is attractive to compute the reference velocity locally, because it is more flexible for the flows with significantly different speeds. In the present study, K is fixed at 0.25. The eigenvalues of Γ-1A(A=∂F/∂W) are  u, u′±a′, where u is the velocity component in the x direction and(14)
u
′
±
a
′
=
1
2
[
(
1
+
M
ref
2
)
u
±
a
(
1
-
M
ref
2
)
2
M
2
+
4
M
ref
2
]
,
M
ref
2
=
U
ref
/
a
2
2
,
M
=
u
a
.


## 3.4. High-Resolution Scheme

The numerical scheme, using the preconditioning finite volume method, was introduced to solve the governing flow equations. A second-order scheme was initially applied, so the left and right states were chosen to be the cell average values on the left and right of the cell faces. In a high-resolution scheme, in order to raise the order of accuracy of upwind differencing, all that is needed is to raise the order of accuracy of the initial-value interpolation that yields the zone-boundary data. Such schemes are labelled as high-resolution schemes as opposed to Total Variation Diminishing (TVD) schemes, which completely eliminate any of those spurious oscillations when applied to one dimensional nonlinear hyperbolic conservation laws and linear hyperbolic systems. The van Leer kappa-scheme, in which the kappa number is one-third, was selected to obtain the high-resolution upwind differencing [16–18].An optimal multistage scheme was used for the time integration, and the multistage coefficients were modified by Tai et al. [19], and redefined using the Courant number for multidimensional use. Also, a residual smoothing method was imposed to accelerate convergence and to improve numerical stability.


## 4. Results and Discussion

### 4.1. Initial Conditions

In order to understand the velocity and temperature statute in the shock tube, the 1D shock tube theorem is applied to determine the shock speed, temperature, and action time as follows:(15)
P
4
P
1
=
147
(
psi
)
1
(
psi
)
=
147
,
∵
P
1
P
4
=
P
1
P
2
[
1
-
γ
-
1
2
γ
(
(
P
2
/
P
1
)
-
1
)
1
+
(
(
γ
+
1
)
/
2
γ
)
(
(
P
2
/
P
1
)
-
1
)
]
2
γ
/
γ
-
1
,
where γ is  the ratio of specific heats, so P2/P1=7.0 from a normal shock table, it shows Ms=2.48, and hence M2=0.5149,T2=635.4K.The history of the action time is shown in Figure3. The main shock wave was reflected to the test model before the contact surface arrived. Therefore, the period of the test time is only 3.2  milliseconds.Figure 3
1D theorem determined result.


### 4.2. Full-Scale Shock Tube Simulation

The time step in an unsteady simulation from the above determined 1D theorem can be set.  In Figure4, P1 and P2 are in the position of the pressure transducer (no. 1 and no. 2) individually.  The result shows the shock speed is 787.4 milliseconds in the test section. The main shock arrived at point no. 1 at t=11.08 milliseconds, and was 2.32 milliseconds faster than the theorem determined.  The numerical model simulates the real case when the calibre is reduced to 1/3 from the high to the low pressure region.  For this reason, the shock speed is rapidly increased.  The calculation of the theory in terms of the main straight tube neglects the heat effect and boundary layer effect because their response time is relatively slow compared with the test event interval.  Figures 5 and 6 show the pressure and temperature change at different time steps: (a) shows the initial state after the diaphragm was broken, and (b) and (c) show the incident shock wave propagation phenomena.  Figure 6(b) shows the temperature distribution, and the contact surface can clearly be found.  Figures 6(d) and 6(e) show the details of the reflected shock wave.  Because the wind tunnel is a closed type, the pressure drops naturally in the high pressure section after the main shock is released, as shown in Figure 5.Pressure and temperature profile at point no. 1 and no. 2.
(a)
(b)Figure 5
Snapshots of the pressure distribution in the tube.Figure 6
Snapshots of the temperature distribution in the tube.


### 4.3. Different Fuel Slab Angle Analysis

After understanding the phenomena of a full-scale shock tube, a fuel slab was placed in the test section(X=14.55~14.70M) to investigate the physical phenomena of the fuel slab surface at different angles of attack (AOA=0°,7°, and 10°).  Figures 7 to 9 show the slab surface pressure, temperature, and Mach number.  Due to the symmetric shape across the upper and lower area of the slab at AOA=0°, the distribution curves merged as one line. But  at AOA=7° and 10°, the temperature and pressure curve across the upper and lower surface appears different, as in Figures 7 to 9.  In Figure 10, at AOA=0°, the blunt shape of the leading edge causes a supersonic bow shock with another shock wave following to make the temperature rise further near the leading edge.  The channel flowed through a convergent nozzle between the upper and lower passage because of the boundary layer effect.  Also because of the passage flow wall effect, the shock wave was reflected, which made the speed rapidly reduce. Therefore, the Mach number decreased from 1.25 to M=0.5 when the air flow attached to a reflection shock at X/C=0.1.  Following this, the flow speed was raised to X/C=0.15~0.4 because of the convergent passage that the boundary layer affects.  However, the energy after the reflecting shock wave was reduced too fast, so the Mach number tended towards stability after X/C=0.5.  At 7 degrees of angle of attack, the impact effect on the upper surface increased so the temperature was higher than the 0 degrees case.  The supersonic flow reduced the speed and raised pressure because of the geometric convergent passage.  The shock wave which occurred in the leading edge of the upper surface reattached at the trailing edge through the reflected shock.  The main shock wave has the deflection characteristic due to the decline in the slab, so the air flow suppresses the boundary layer.  For this reason, the temperature increased on the upper surface.  The supersonic flow speed rose and pressure reduced when the air flowed over the lower surface which was like in a divergent nozzle. Therefore, at AOA=7° and when the HTPB fuel was mounted on the upper surface of the slab, the high temperature helped the gasification and ignition in HTPB plate.  In Figure 9, it was found that the air flow speed to be still above  M=1, either on the upper or lower surface at AOA=7°.  Therefore, this situation can satisfy the ignition in supersonic flow.  As AOA=10°, the average temperature was not as high as in AOA=7°, but the speed was still above M=1.  Because of the reflection of the main shock wave at the leading edge, the pressure was reduced, the temperature was raised, and the Mach number was reduced after the flow through the reflected shock (Figures 7 to 9, X/C=0.3).Figure 7
Pressure profile of HTPB surface at different AOA (C: length of the HTPB slab).


### 4.4. Experiment Visualization

Slab atAOA=0
Figure11 shows the time history snapshot of the shock wave using a CCD camera.  The shock wave arrived at the slab, AOA=0 deg, at t=11.1 milliseconds.  The bow shock occurred because of the blunt shape configuration.  The secondary oblique shock wave was induced at t=11.7 milliseconds when the flow went over the HTPB.  Compared with the CFD results, the phenomenon was matched, as shown in Figure 12.  The flow speed behind the shock rapidly reduced to subsonic, and the speed of the stagnation point on the leading edge reduced to zero.  Although the speed was reduced behind the bow shock, the secondary shock was induced at 1.17 cm from the leading edge because of the boundary effect, as shown in Figure 13.
The reflected shock was induced in the trailing edge att=14.3 ms as shown in Figure 12.  The period was about 3 ms between the shock’s arrival and the reflected shock attached to the test section.  The result fit the CFD calculations.  Figure 14 shows the time history of the pressure transducer record at 14.55 M (#1) and 14.65 M (#2).  Because the HTPB fuel was not burned, the change of the pressure transducer record was not clear until the impact of the main shock and reflected shock.  It was found that the pressure fluctuated after a rapid peak, as shown in Figure 14, because of the interaction of shocks in the tube.  The total pressure was not reduced until the pressure stabilized after the peak value.Slab atAOA=7
Figure15 shows the time history snapshot of the shock wave motion when the slab was at AOA=7 deg.  The shock wave arrived at the slab at t=11.1 milliseconds, and the secondary oblique shock wave was induced at t=11.7 milliseconds when it flowed over the HTPB.  Compared with the CFD results, as shown in Figure 16, the flow speed behind the shock remained transonic on the upper side, and the speed of the lower side increased to M=2.  Although the speed was reduced behind the oblique shock on the upper side, the temperature was also proportional to the length, as shown in Figure 8.
The reflected shock was also induced at the trailing edge att=14.1 milliseconds, as shown in Figure 11(d).  Because of the interaction of shocks in the tube, it was also found that the pressure fluctuated after a rapid peak, as shown in Figure 17.
It was observed that gasification exits as the shock wave moves across the HTPB.  From the recovered fuel, melting was observed on the leading edge, as shown in Figure18.  This region is a stagnation zone for a bow shock, where the flow speed was slow and the temperature was high.  The temperature reached 1100 K, as shown in Figure 8.  The temperature was already higher than the surface gasification temperature (930 K 1190 K) of the HTPB, and there was enough time for burning.  The flow speed in other places was too fast, so the temperature did not reach the gasification criterion, and therefore was unable to burn.Figure 8
Temperature profile of HTPB surface at different AOA  (C: length of the HTPB slab).Figure 9
Mach no. profile of HTPB surface at different AOA (C: length of the HTPB slab).Iso-density contour (a)AOA=0∘, (b) AOA=7∘, (c) AOA=10∘.
(a)
(b)
(c)Snapshots of shock wave (HTPB slab @AOA=0°)  (a) before event (t=0 second), (b) shock wave arrived (dir.: →,t=11.1 milliseconds),  (c) shock wave arrived (dir.:→, t=11.7 milliseconds),  (d) ref. Shock (dir.:←t=14.3 milliseconds).
(a)
(b)
(c)
(d)Figure 12
Comparison Exp. and CFD result (HTPB slab @AOA=0°).Figure 13
Mach No. distribution (HTPB slab @AOA=0°).Figure 14
Pressure measurement (HTPB slab @AOA=0°).Snapshots of shock wave (HTPB slab @AOA=7°) (a) before event (t=0 second), (b) shock wave  arrived (dir.:→, t=11.1 milliseconds), (c) shock wave through (dir.:→, t=11.7 milliseconds), (d) ref. Shock (dir.:←, t=14.1 milliseconds).
(a)
(b)
(c)
(d)Figure 16
Mach No. distribution (HTPB slab @AOA=7°).Figure 17
Pressure measurement (point #1, HTPB slab @AOA=7°).Comparison before and after the action (a) without action, (b) with action.
(a)
(b)


## 4.1. Initial Conditions

In order to understand the velocity and temperature statute in the shock tube, the 1D shock tube theorem is applied to determine the shock speed, temperature, and action time as follows:(15)
P
4
P
1
=
147
(
psi
)
1
(
psi
)
=
147
,
∵
P
1
P
4
=
P
1
P
2
[
1
-
γ
-
1
2
γ
(
(
P
2
/
P
1
)
-
1
)
1
+
(
(
γ
+
1
)
/
2
γ
)
(
(
P
2
/
P
1
)
-
1
)
]
2
γ
/
γ
-
1
,
where γ is  the ratio of specific heats, so P2/P1=7.0 from a normal shock table, it shows Ms=2.48, and hence M2=0.5149,T2=635.4K.The history of the action time is shown in Figure3. The main shock wave was reflected to the test model before the contact surface arrived. Therefore, the period of the test time is only 3.2  milliseconds.Figure 3
1D theorem determined result.


## 4.2. Full-Scale Shock Tube Simulation

The time step in an unsteady simulation from the above determined 1D theorem can be set.  In Figure4, P1 and P2 are in the position of the pressure transducer (no. 1 and no. 2) individually.  The result shows the shock speed is 787.4 milliseconds in the test section. The main shock arrived at point no. 1 at t=11.08 milliseconds, and was 2.32 milliseconds faster than the theorem determined.  The numerical model simulates the real case when the calibre is reduced to 1/3 from the high to the low pressure region.  For this reason, the shock speed is rapidly increased.  The calculation of the theory in terms of the main straight tube neglects the heat effect and boundary layer effect because their response time is relatively slow compared with the test event interval.  Figures 5 and 6 show the pressure and temperature change at different time steps: (a) shows the initial state after the diaphragm was broken, and (b) and (c) show the incident shock wave propagation phenomena.  Figure 6(b) shows the temperature distribution, and the contact surface can clearly be found.  Figures 6(d) and 6(e) show the details of the reflected shock wave.  Because the wind tunnel is a closed type, the pressure drops naturally in the high pressure section after the main shock is released, as shown in Figure 5.Pressure and temperature profile at point no. 1 and no. 2.
(a)
(b)Figure 5
Snapshots of the pressure distribution in the tube.Figure 6
Snapshots of the temperature distribution in the tube.


## 4.3. Different Fuel Slab Angle Analysis

After understanding the phenomena of a full-scale shock tube, a fuel slab was placed in the test section(X=14.55~14.70M) to investigate the physical phenomena of the fuel slab surface at different angles of attack (AOA=0°,7°, and 10°).  Figures 7 to 9 show the slab surface pressure, temperature, and Mach number.  Due to the symmetric shape across the upper and lower area of the slab at AOA=0°, the distribution curves merged as one line. But  at AOA=7° and 10°, the temperature and pressure curve across the upper and lower surface appears different, as in Figures 7 to 9.  In Figure 10, at AOA=0°, the blunt shape of the leading edge causes a supersonic bow shock with another shock wave following to make the temperature rise further near the leading edge.  The channel flowed through a convergent nozzle between the upper and lower passage because of the boundary layer effect.  Also because of the passage flow wall effect, the shock wave was reflected, which made the speed rapidly reduce. Therefore, the Mach number decreased from 1.25 to M=0.5 when the air flow attached to a reflection shock at X/C=0.1.  Following this, the flow speed was raised to X/C=0.15~0.4 because of the convergent passage that the boundary layer affects.  However, the energy after the reflecting shock wave was reduced too fast, so the Mach number tended towards stability after X/C=0.5.  At 7 degrees of angle of attack, the impact effect on the upper surface increased so the temperature was higher than the 0 degrees case.  The supersonic flow reduced the speed and raised pressure because of the geometric convergent passage.  The shock wave which occurred in the leading edge of the upper surface reattached at the trailing edge through the reflected shock.  The main shock wave has the deflection characteristic due to the decline in the slab, so the air flow suppresses the boundary layer.  For this reason, the temperature increased on the upper surface.  The supersonic flow speed rose and pressure reduced when the air flowed over the lower surface which was like in a divergent nozzle. Therefore, at AOA=7° and when the HTPB fuel was mounted on the upper surface of the slab, the high temperature helped the gasification and ignition in HTPB plate.  In Figure 9, it was found that the air flow speed to be still above  M=1, either on the upper or lower surface at AOA=7°.  Therefore, this situation can satisfy the ignition in supersonic flow.  As AOA=10°, the average temperature was not as high as in AOA=7°, but the speed was still above M=1.  Because of the reflection of the main shock wave at the leading edge, the pressure was reduced, the temperature was raised, and the Mach number was reduced after the flow through the reflected shock (Figures 7 to 9, X/C=0.3).Figure 7
Pressure profile of HTPB surface at different AOA (C: length of the HTPB slab).


## 4.4. Experiment Visualization

Slab atAOA=0
Figure11 shows the time history snapshot of the shock wave using a CCD camera.  The shock wave arrived at the slab, AOA=0 deg, at t=11.1 milliseconds.  The bow shock occurred because of the blunt shape configuration.  The secondary oblique shock wave was induced at t=11.7 milliseconds when the flow went over the HTPB.  Compared with the CFD results, the phenomenon was matched, as shown in Figure 12.  The flow speed behind the shock rapidly reduced to subsonic, and the speed of the stagnation point on the leading edge reduced to zero.  Although the speed was reduced behind the bow shock, the secondary shock was induced at 1.17 cm from the leading edge because of the boundary effect, as shown in Figure 13.
The reflected shock was induced in the trailing edge att=14.3 ms as shown in Figure 12.  The period was about 3 ms between the shock’s arrival and the reflected shock attached to the test section.  The result fit the CFD calculations.  Figure 14 shows the time history of the pressure transducer record at 14.55 M (#1) and 14.65 M (#2).  Because the HTPB fuel was not burned, the change of the pressure transducer record was not clear until the impact of the main shock and reflected shock.  It was found that the pressure fluctuated after a rapid peak, as shown in Figure 14, because of the interaction of shocks in the tube.  The total pressure was not reduced until the pressure stabilized after the peak value.Slab atAOA=7
Figure15 shows the time history snapshot of the shock wave motion when the slab was at AOA=7 deg.  The shock wave arrived at the slab at t=11.1 milliseconds, and the secondary oblique shock wave was induced at t=11.7 milliseconds when it flowed over the HTPB.  Compared with the CFD results, as shown in Figure 16, the flow speed behind the shock remained transonic on the upper side, and the speed of the lower side increased to M=2.  Although the speed was reduced behind the oblique shock on the upper side, the temperature was also proportional to the length, as shown in Figure 8.
The reflected shock was also induced at the trailing edge att=14.1 milliseconds, as shown in Figure 11(d).  Because of the interaction of shocks in the tube, it was also found that the pressure fluctuated after a rapid peak, as shown in Figure 17.
It was observed that gasification exits as the shock wave moves across the HTPB.  From the recovered fuel, melting was observed on the leading edge, as shown in Figure18.  This region is a stagnation zone for a bow shock, where the flow speed was slow and the temperature was high.  The temperature reached 1100 K, as shown in Figure 8.  The temperature was already higher than the surface gasification temperature (930 K 1190 K) of the HTPB, and there was enough time for burning.  The flow speed in other places was too fast, so the temperature did not reach the gasification criterion, and therefore was unable to burn.Figure 8
Temperature profile of HTPB surface at different AOA  (C: length of the HTPB slab).Figure 9
Mach no. profile of HTPB surface at different AOA (C: length of the HTPB slab).Iso-density contour (a)AOA=0∘, (b) AOA=7∘, (c) AOA=10∘.
(a)
(b)
(c)Snapshots of shock wave (HTPB slab @AOA=0°)  (a) before event (t=0 second), (b) shock wave arrived (dir.: →,t=11.1 milliseconds),  (c) shock wave arrived (dir.:→, t=11.7 milliseconds),  (d) ref. Shock (dir.:←t=14.3 milliseconds).
(a)
(b)
(c)
(d)Figure 12
Comparison Exp. and CFD result (HTPB slab @AOA=0°).Figure 13
Mach No. distribution (HTPB slab @AOA=0°).Figure 14
Pressure measurement (HTPB slab @AOA=0°).Snapshots of shock wave (HTPB slab @AOA=7°) (a) before event (t=0 second), (b) shock wave  arrived (dir.:→, t=11.1 milliseconds), (c) shock wave through (dir.:→, t=11.7 milliseconds), (d) ref. Shock (dir.:←, t=14.1 milliseconds).
(a)
(b)
(c)
(d)Figure 16
Mach No. distribution (HTPB slab @AOA=7°).Figure 17
Pressure measurement (point #1, HTPB slab @AOA=7°).Comparison before and after the action (a) without action, (b) with action.
(a)
(b)


## 5. Conclusion

From the numerical simulation results, the HTPB slab at an angle of attack of 7 degrees, had higher temperature and pressure in the upper surface than at 0 degrees and 10 degrees, and the flow speed of the upper and lower surfaces was kept at supersonic and contributed to gasification ignition.  According to the experimental data and the numerical results, the test period was both about 3 milliseconds, and resulted in the melting of the leading edge of the tested HTPB slab.  In other areas, although the gasification criterion was reached, it was unable to burn because the flow speed was too great and the test time was limited to only 3 ms in the shock tube facility.



---

*Source: 823874-2009-06-23.xml*