# Impact of the Surface Morphology on the Combustion of Simulated Solid Rocket Motor

Authors: Abdelkarim M. Hegab; Hani Hussain Sait; Ahmad Hussain
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

An advanced and intensive computational solution development is integrated with an asymptotic technique, to examine the impact of the combustion surface morphology on the generated rotational flow field in a solid rocket chamber with wide ranges of forcing frequencies. The simulated rectangular chamber is closed at one end and is open at the aft end. The upper and lower walls are permeable to allow steady and unsteady injected air to generate internal flow mimicking the flow field of the combustion gases in real rocket chamber.  The frequencies of the unsteady injected flow are chosen to be very close or away from the resonance frequencies of the adapted chamber. The current study accounts for a wide range of wave numbers that reflect the complexity of real burning processes. Detailed derivation for Navier-Stokes equations at the four boundaries of the chamber is introduced in the current study. Qualitative comparison is performed with recent experimental work carried out on a two-inch hybrid rocket motor using a mixture of polyethylene and aluminum powder. The higher the percentage of aluminum powder in the mixture, the more the corrugations of the combustion surface. This trend is almost similar to the computational and analytical results of a simulated solid rocket chamber.


---

## Body

## 1. Introduction

Nowadays, the designers of solid and hybrid rocket systems have a big challenge to understand the dynamics of the complex rotational fluid flow inside the rocket’s chamber, nozzle, and the exhaust plume. In order to have a stable rocket system with desirable performance, individual study for the combustion model that incorporates the burning of microscale propellant and the nonreactive model that accounts the complex fluid flow inside the chamber and the nozzles must be examined first. Some recent studies that consider the burning of microscale solid rocket propellant [1–6] concluded that the injection conditions of the cold models in [7–14] are found to be different than that considered with the reactive models [1–6]. As a result, the current study is conducted to examine the unsteady flow field inside the solid rocket motor chamber and the complex generated wave and vorticity patterns and their impact on the burning of real solid propellant surfaces. We try in this study to assign unsteady velocity injection with axial distributions to simulate the real variation of composite propellants.Recent studies of the rotational flow field generating from transient side wall mass injection have been surveyed in detail in [12–14]. These studies revealed that the coexistence of the acoustic wave patterns and the injection mass addition are the reason behind the generating of the vorticity dynamics in the combustion chamber. However, many researches have been directed to analyze the acoustic wave patterns including the acoustic velocity, temperature, and pressure time history without any consideration for the existence of rotational flow (vorticity) transients. A competitive linear [15–17] and nonlinear theory [11, 13, 14] are formulated to show that the precise descriptions of the fluid flow in the simulated rocket motor chamber require the presence of the acoustic and the rotational flow fields.Numerical approaches for the presence of vorticity in an injected internal flow are presented in [19, 20] by solving the Navier-Stokes equation. Moreover, the studies in [9, 10] presented detailed comparison of asymptotic and numerical approaches without any consideration to the thermal effect. Intensive numerical procedures for the effect of exothermic reactions very close to the sidewall have been carried out in [21, 22].Finally, the root mean square value for the axial distribution for the coexistence of rotational and irrotational flow field is calculated in [8, 23].In the present study, asymptotic and computational results are obtained to consider the unsteady variation of the vorticity and temperature dynamic in the prescribed chamber. Results are given for steady state solution generated from spatially uniform injection velocity. The convergent solution is oscillated with similar amplitude unsteady mass additions with wide ranges of injection wave numbers. The magnitude of the superimposed injection velocity is chosen to be large enough to make sure that the fluid is always injected from the side wall and reduce the probability of sucking outside the chamber.


## 2. The Physical and Mathematical Formulation

The present numerical study is performed to study the effect of the complicated surface morphology of the combustion surface on the generated vorticity in a chamber of lengthL
′ and half-height H
′. The aspect ratio in this study δ
=
L
′
/
H
′ is taken to be twenty. Pressure node is assigned at the chamber exit, while closed end is considered at the head end. Air is injected from the side wall with a characteristic velocity V
y
o
′ to simulate the combustion products from real burning of propellant. This modeling generates an axial flow characterized by U
z
o
′
=
δ
V
y
o
′ as in Figure 1.Figure 1
The physical domain and the boundary conditions for the simulated solid rocket motor chamber.Thesteady and the unsteady injection velocity are denoted by v
r
w
s
′
x, v
^
o
′
wall
(
x
′
,
t
′
), while the exit pressure and the temperature of the fluid injected from the side wall are denoted by P
o
′ and T
o
′.The parabolized Navier Stokes equations for unsteady, two dimensional, compressible fluid flow are solved numerically in the prescribed chamber and may be written as follows:(1)
∂
Q
∂
t
+
∂
E
∂
x
+
∂
F
∂
y
=
0
,
where
(2)
Q
=
ρ
,
ρ
u
,
ρ
v
,
E
T
T
,
E
=
M
ρ
u
,
M
ρ
u
2
+
1
γ
M
2
p
,
M
ρ
u
v
,
M
E
T
+
γ
-
1
p
u
T
,
F
=
M
ρ
v
,
M
ρ
u
v
-
δ
2
M
R
e
u
y
,
M
ρ
v
2
+
δ
2
γ
M
p
,
M
E
T
+
γ
-
1
p
v
-
γ
δ
2
M
R
e
P
r
T
y
T
.
E
T is the total energy (ρ
[
T
+
γ
(
γ
-
1
)
M
2
(
u
2
+
(
v
/
δ
)
2
]) and the equation of state for a perfect gas is
(3)
p
=
ρ
T
.
In this study, the specific heats, the viscosity, and conductivity coefficients are considered to be constant, since the change in the temperature in the cold model is relatively small. Nondimensional variables are given by
(4)
p
=
p
′
p
o
′
,
ρ
=
ρ
′
ρ
o
′
,
T
=
T
′
T
o
′
,
u
=
u
′
U
z
o
′
,
v
=
v
′
U
z
o
′
/
δ
,
x
=
x
′
L
′
,
y
=
y
′
H
′
,
t
=
t
′
t
a
′
.The axial flow Mach number, Prandtl number, flow and acoustic Reynolds numbers are written as:(5)
M
=
U
z
o
′
C
o
′
,
P
r
=
μ
o
′
C
p
o
′
k
o
′
,
R
e
=
ρ
o
′
U
z
o
′
L
′
μ
o
′
,
R
e
A
=
R
e
M
.
The characteristic Mach number is chosen to be O(10−1) while the characteristic Reynolds number is O(105-106). These orders of magnitudes approximate the characteristics of the fluid flow in similar real rocket engine. The nondimensionalized time t
a
′ is defined as the axial acoustic time t
a
′
=
L
′
/
C
o
′.


### 2.1. The Initial and Boundary Conditions

Fort
≤
0, time-independent side wall transverse velocity with V
wall
=
-
V
r
w
s
(
x
) is used to generate steady flow solution, while, for t
>
0, V
wall
=
V
r
w
s
(
x
,
t
)
=
-
V
r
w
s
(
x
)
-
ε
φ
wall
(
x
)
(
1
-
cos
⁡
ω
t
), as the unsteady longitudinal variation, for the injection transverse at resonance or away from the resonance case. The function φ
wall represents the longitudinal variation of the velocity at the injection surface. Moreover, the specified conditions at the four boundaries can be written as follows:
(6)
x
=
0
;
u
=
0
,
v
=
0
,
x
=
1
;
p
=
1
,
“pressure node”
y
=
0
;
v
=
0
,
∂
u
∂
y
=
∂
T
∂
y
=
∂
p
∂
y
=
∂
ρ
∂
y
=
0
,
y
=
1
;
u
=
0
,
T
=
1
,
v
=
-
v
r
w
s
x
,
t
≤
0
,
-
v
r
w
s
x
-
ε
φ
r
w
1
-
cos
⁡
ω
t
,
t
>
0
,
and the initial condition is
(7)
t
=
0
;
u
=
v
=
0
,
p
=
T
=
1
.
The initial conditions for the unsteady flow computations are the convergence solutions for the steady state conditions.


### 2.2. The Analytical Approach

#### 2.2.1. Steady and Unsteady Flow

Hegab [7] has concluded that the total axial velocity component includes two main parts. The first part represents the steady state component and the other part represents unsteady state component and can be written as follows:
(8)
u
,
v
,
P
,
ρ
,
T
=
u
s
,
v
s
,
P
s
,
ρ
s
,
T
s
+
u
^
,
v
^
,
P
^
,
ρ
^
,
T
^
.
Using the expended variables used in Hegab [7] and inserting them in (1), the side wall steady vorticity generation (Ω
o
s
wall) and the accompanying heat transfer (∂
T
o
s
/
∂
y
wall) are derived and written as follows:
(9)
Ω
os
wall
=
-
1
γ
v
r
w
s
∂
∂
x
γ
π
2
4
∫
x
1
v
r
w
s
(
x
^
)
∫
0
x
^
v
r
w
s
(
τ
)
d
τ
d
x
^
+
∂
v
r
w
s
x
∂
x
,
∂
T
os
∂
y
wall
=
0
.Equations (9) show that nonzero steady vorticity and zero heat transfer are found with the use of no-slip and constant injection boundary conditions.The sequence of the numerical solution consists of two steps. The first one is the initial steady state computation which arises from constant injection with isothermal condition. The second one is the unsteady flow computation based on initial conditions as computed from the first step. The perturbation theory introduced by Staab et al. [11] and applied by Hegab [7] leads to the following weakly rotational equations:
(10)
t
=
0
,
P
^
o
=
0
,
∂
P
^
o
∂
t
=
0
,
x
=
0
,
∂
P
^
o
∂
x
=
0
,
x
=
1
,
P
^
o
=
0
,
y
=
0
;
v
^
0
=
0
,
y
=
1
;
v
^
o
=
-
ε
φ
r
w
(
1
-
cos
⁡
ω
t
)
t
>
0
.
The final nonresonance solution for the flow in the simulated chamber in Figure 1 is
(11)
P
^
o
x
,
t
=
∑
n
=
0
∞
2
ε
γ
ω
∫
0
1
φ
r
w
(
x
)
cos
⁡
β
n
x
ω
2
-
β
n
2
×
ω
β
n
sin
β
n
t
-
sin
ω
t
cos
⁡
β
n
x
,
u
^
o
x
,
t
=
∑
n
=
0
∞
2
ε
γ
ω
∫
0
1
φ
r
w
(
x
)
cos
⁡
β
n
x
ω
2
-
β
n
2
×
ω
β
n
cos
⁡
β
n
t
-
1
-
β
n
ω
cos
⁡
ω
t
-
1
×
sin
β
n
x
,
v
^
o
x
,
y
,
t
=
-
y
φ
r
w
x
1
-
cos
⁡
ω
t
.
Equations (11) show the time history of the perturbed values variables P
^
o
,
u
^
o
,
v
^
o at any axial locations.Following the same policies, the resonance solutions are derived and formulated as(12)
P
^
o
x
,
t
=
γ
C
~
n
*
β
n
*
sin
β
n
*
t
-
β
n
*
t
cos
⁡
β
n
*
t
cos
⁡
β
n
*
x
,
u
^
o
x
,
t
=
C
~
n
*
-
2
β
n
*
cos
⁡
β
n
*
t
-
1
-
t
sin
β
n
*
t
×
cos
⁡
β
n
*
x
,
where C
^
n
*
=
2
∫
0
1
φ
r
w
(
x
)
cos
⁡
β
n
*
x.


### 2.3. The Computational Approach

In the present study, higher order accuracy difference equations are employed to solve the parabolized Navier-Stokes equations in the interior of the prescribed chamber and at the boundaries. The numerical multisteps scheme in [24] is used to solve the 2D, unsteady, and compressible Navier-Stokes equations for the points located near the boundaries, while the two-four explicit, predictor-corrector scheme in [25] is employed to solve the adapted equations at the interior points. The latter method is a fourth-order accuracy variant of the fully explicit MacCormack scheme. Uniform grids are used. The accuracy criterion employed by MacCormack [24] is used to determine the size of the time step.Here, the NSCBC technique applied by [7] is presented in detail. The full derivation for the parabolized Navier-Stokes at the upper and lower boundaries and at the head and exit ends is formulated to specify the nonphysical conditions at these boundaries. Moreover, this method is applied to reduce the generated reflected waves arising from the numerical techniques.Based on the idea of characteristics for Euler equations, the first derivatives normal to the boundaries appearing in the convective terms of the Navier-Stokes equations may be expressed in terms of the amplitude of incoming, outcoming, and entropy waves{
℘
i
’s
}. The fluid dynamics equations can be written in nondimensional wave modes form as follows:
(13)
∂
ρ
∂
t
+
1
C
2
℘
2
+
1
2
℘
4
+
℘
1
+
M
∂
ρ
v
∂
y
=
0
,
(14)
∂
u
∂
t
+
1
2
ρ
C
℘
4
-
℘
1
+
M
v
∂
u
∂
y
=
χ
x
ρ
,
(15)
∂
v
∂
t
+
℘
3
+
M
v
∂
v
∂
y
+
A
2
γ
M
1
ρ
∂
p
∂
y
=
χ
y
ρ
,
(16)
∂
p
∂
t
+
1
2
℘
4
+
℘
1
+
M
v
∂
p
∂
y
+
γ
M
p
∂
v
∂
y
=
χ
E
,
where χ
x and χ
y represent the viscous terms and χ
E refers to the combination of the viscous and conduction terms in the axial and transverse directions and C
=
C
′
/
C
o
′. The amplitudes of the left-running waves (℘
1), the right-running waves (℘
4), and the entropy waves (℘
2
,
℘
3) in terms of the characteristic velocities λ
i’s are
(17)
℘
1
=
λ
1
1
γ
∂
p
∂
x
-
M
ρ
C
∂
u
∂
x
;
λ
1
=
u
M
-
C
,
(18)
℘
2
=
λ
2
C
2
∂
ρ
∂
x
-
1
γ
∂
p
∂
x
;
λ
2
=
M
u
,
(19)
℘
3
=
λ
3
∂
v
∂
x
;
λ
3
=
M
u
,
(20)
℘
4
=
λ
4
1
γ
∂
p
∂
x
+
M
ρ
C
∂
u
∂
x
;
λ
4
=
u
M
+
C
.Two approaches are used to infer the values of the incoming amplitude-variations from the boundary conditions. One is the local one-dimensional inviscid relations (LODI) technique and the second is the Navier-Stokes characteristics boundary conditions (NSCBC). For the latter one, the complete description of the numerical boundary conditions at all boundaries is discussed below.


### 2.4. Exit Boundary Conditions

In this case, a pressure node is assumed at the exit boundary. The NSCBC relation (16) suggests that the amplitude of the reflected wave should be
(21)
℘
1
=
2
χ
E
-
℘
4
,
and the remaining conditions are
(22)
℘
2
=
M
u
C
2
∂
ρ
∂
x
-
1
γ
∂
p
∂
x
,
℘
3
=
M
u
∂
v
∂
x
,
℘
4
=
u
M
+
c
1
γ
∂
p
∂
x
+
M
ρ
C
∂
u
∂
x
.
The nondimensional equations can be written as follows:
(23)
∂
ρ
∂
t
+
M
u
∂
ρ
∂
x
-
M
u
γ
C
2
∂
p
∂
x
+
M
v
∂
ρ
∂
y
+
A
2
M
γ
-
1
R
e
P
r
1
C
2
∂
2
T
∂
y
2
=
0
,
∂
u
∂
t
+
1
ρ
C
u
+
C
M
1
γ
∂
p
∂
x
+
M
ρ
C
∂
u
∂
x
+
M
v
∂
u
∂
y
-
A
2
γ
-
1
R
e
P
r
1
ρ
C
∂
2
T
∂
y
2
-
M
A
2
R
e
1
ρ
∂
2
u
∂
y
2
=
0
,
∂
v
∂
t
+
M
u
∂
v
∂
x
+
v
∂
v
∂
y
=
0
,
d
p
d
t
=
0
.Here again, the nondimensional variables (ρ, u, and v) can be computed using one-sided finite-difference approximations.


### 2.5. Head End Boundary Condition

The velocity component normal to the boundary(
u
(
0
,
y
,
t
)
) is zero since the head end is impermeable. In contrast the transverse velocity component cannot be specified because the axial viscous terms in the transverse momentum equation have been suppressed. As the normal velocity drops to zero at x
=
0, the amplitudes ℘
2 and ℘
3, in (18) and (19), are zero and the characteristic velocities are λ
1
=
-
C and λ
4
=
+
C. The reflected wave ℘
4 enters the domain, while the wave ℘
1 propagates out through the boundary. The NSCBC relation (19) indicates that the amplitude of the incoming wave ℘
4 should be
(24)
℘
4
=
℘
1
,
and the remaining conditions are
(25)
℘
1
=
-
c
1
γ
∂
p
∂
x
+
M
ρ
C
∂
u
∂
x
,
℘
2
=
0
,
℘
3
=
0
.
The equations can be written in nondimensional form as follows:
(26)
∂
ρ
∂
t
-
M
ρ
∂
u
∂
x
+
1
γ
C
∂
p
∂
x
+
M
∂
ρ
v
∂
y
=
0
,
∂
u
∂
t
=
0
,
∂
v
∂
t
+
M
v
∂
v
∂
y
=
0
,
∂
p
∂
t
-
c
∂
p
∂
x
+
M
ρ
C
∂
u
∂
x
+
γ
M
p
∂
v
∂
y
-
γ
A
2
M
γ
-
1
P
r
R
e
∂
2
T
∂
y
2
=
0
.
More details about the NSCBC techniques are found in [23, 26].


### 2.6. Code Validations

The computational codes are designed, constructed, and validated for steady and unsteady flow as shown in Figures2 and 3. For the latter one, the constructed codes run for many wave cycles to detect the whole unsteady rotational flow dynamic across the chamber.Figure 2
Normalized axial velocity profiles at the midlength of the chamber for the steady state flow.Figure 3
Comparison between the computational and the asymptotic solutions for the perturbed temperature time history withω
=
1 (nonresonance), ε
=
0.4, M
=
0.02, and R
e
=
3
*
10
5 at the half-length of the chamber.In Figure2, comparison between the current numerical and asymptotic results dealing with the steady state flow solutions is carried out. The steady state flow is generated by adding constant side wall mass injection (
V
r
w
=
-
1
) from the injection surface of the prescribed chamber. Similar experiment in [18, 27] was carried out and compared with the current mathematical models in Figure 2. The asymptotic analysis solution represents the solid line, while the dotted line represents the converged steady state numerical solution for viscous flow with wall Reynolds number R
e
w
=
5
*
10
3. The nonfilled points refer to the experimental measured data. The experimental results are recorded for flow in a channel with fully transpired wall. Three different porous surfaces using three different sizes of honeycomb cells are used, [7]. In this figure, a comparison between the analytical and the numerical solutions at the midlength of the chamber is illustrated. This comparison shows a reasonable agreement between the three approaches. The comparison between the theoretical and the computational approaches shows little bit deviation. This deviation is logically expected since the used honey combs have three different sizes. The smaller size of honey combs, the smaller deviations between the two approaches.The second set of the validations show the solution for the unsteady flow computations at nonresonance forced mode. In this set of the results, the computational results are validated by comparing with the asymptotic analysis. These results represent the asymptotic solution when the sidewall injection velocity is represented byv
r
w
(
x
,
1
,
t
)
=
-
1
+
ε
φ
r
w
(
1
-
cos
⁡
ω
t
), φ
r
w
=
cos
⁡
β
n
x. In the current work and the work by Hegab [7], the disturbance forcing frequency is taken to be ω
=
1 (nonresonance) with the amplitude of ε
=
0.4 and the eigenvalues (β
n
=
n
π
/
2, n
=
1,3
,
5,7
,
9
,
…
,
etc.). The temperature time history (T
^) at the medlength location is presented in Figure 3 for R
e
=
3
*
10
5, M
=
0.02, δ
=
20, and n
=
1.0. For the nonresonance case, the solution shows stable conditions but it is not purely harmonic. This interesting behavior shows the ability of the constructed code with the good treatment to the boundary conditions to detect the non-harmonic trend as detected by the analytical approach. A comparison of the perturbed temperature (T
^) (27) with the numerical results in Figure 3 shows good agreement between the two approaches especially for t
<
20, while the computational solution shows little less deviation than the analytical results for t
>
20 due to the nonlinearity effects:
(27)
T
^
o
x
,
t
=
∑
n
=
0
∞
2
ε
ω
γ
-
1
∫
0
1
φ
r
w
(
x
)
cos
⁡
β
n
x
ω
2
-
β
n
2
×
ω
β
n
sin
β
n
t
-
sin
ω
t
cos
⁡
β
n
x
.Generally, the coexistence of the forced and the natural modes of the system is detected in the adapted computational and analytical approaches as seen in the earlier investigation by Hegab [7]. In contrast, the computational techniques in [10] showonly pure harmonic oscillations. The source of these deficiencies may be arising from the treatment of the boundary conditions especially at the head and aft ends of the simulated chamber.


## 2.1. The Initial and Boundary Conditions

Fort
≤
0, time-independent side wall transverse velocity with V
wall
=
-
V
r
w
s
(
x
) is used to generate steady flow solution, while, for t
>
0, V
wall
=
V
r
w
s
(
x
,
t
)
=
-
V
r
w
s
(
x
)
-
ε
φ
wall
(
x
)
(
1
-
cos
⁡
ω
t
), as the unsteady longitudinal variation, for the injection transverse at resonance or away from the resonance case. The function φ
wall represents the longitudinal variation of the velocity at the injection surface. Moreover, the specified conditions at the four boundaries can be written as follows:
(6)
x
=
0
;
u
=
0
,
v
=
0
,
x
=
1
;
p
=
1
,
“pressure node”
y
=
0
;
v
=
0
,
∂
u
∂
y
=
∂
T
∂
y
=
∂
p
∂
y
=
∂
ρ
∂
y
=
0
,
y
=
1
;
u
=
0
,
T
=
1
,
v
=
-
v
r
w
s
x
,
t
≤
0
,
-
v
r
w
s
x
-
ε
φ
r
w
1
-
cos
⁡
ω
t
,
t
>
0
,
and the initial condition is
(7)
t
=
0
;
u
=
v
=
0
,
p
=
T
=
1
.
The initial conditions for the unsteady flow computations are the convergence solutions for the steady state conditions.


## 2.2. The Analytical Approach

### 2.2.1. Steady and Unsteady Flow

Hegab [7] has concluded that the total axial velocity component includes two main parts. The first part represents the steady state component and the other part represents unsteady state component and can be written as follows:
(8)
u
,
v
,
P
,
ρ
,
T
=
u
s
,
v
s
,
P
s
,
ρ
s
,
T
s
+
u
^
,
v
^
,
P
^
,
ρ
^
,
T
^
.
Using the expended variables used in Hegab [7] and inserting them in (1), the side wall steady vorticity generation (Ω
o
s
wall) and the accompanying heat transfer (∂
T
o
s
/
∂
y
wall) are derived and written as follows:
(9)
Ω
os
wall
=
-
1
γ
v
r
w
s
∂
∂
x
γ
π
2
4
∫
x
1
v
r
w
s
(
x
^
)
∫
0
x
^
v
r
w
s
(
τ
)
d
τ
d
x
^
+
∂
v
r
w
s
x
∂
x
,
∂
T
os
∂
y
wall
=
0
.Equations (9) show that nonzero steady vorticity and zero heat transfer are found with the use of no-slip and constant injection boundary conditions.The sequence of the numerical solution consists of two steps. The first one is the initial steady state computation which arises from constant injection with isothermal condition. The second one is the unsteady flow computation based on initial conditions as computed from the first step. The perturbation theory introduced by Staab et al. [11] and applied by Hegab [7] leads to the following weakly rotational equations:
(10)
t
=
0
,
P
^
o
=
0
,
∂
P
^
o
∂
t
=
0
,
x
=
0
,
∂
P
^
o
∂
x
=
0
,
x
=
1
,
P
^
o
=
0
,
y
=
0
;
v
^
0
=
0
,
y
=
1
;
v
^
o
=
-
ε
φ
r
w
(
1
-
cos
⁡
ω
t
)
t
>
0
.
The final nonresonance solution for the flow in the simulated chamber in Figure 1 is
(11)
P
^
o
x
,
t
=
∑
n
=
0
∞
2
ε
γ
ω
∫
0
1
φ
r
w
(
x
)
cos
⁡
β
n
x
ω
2
-
β
n
2
×
ω
β
n
sin
β
n
t
-
sin
ω
t
cos
⁡
β
n
x
,
u
^
o
x
,
t
=
∑
n
=
0
∞
2
ε
γ
ω
∫
0
1
φ
r
w
(
x
)
cos
⁡
β
n
x
ω
2
-
β
n
2
×
ω
β
n
cos
⁡
β
n
t
-
1
-
β
n
ω
cos
⁡
ω
t
-
1
×
sin
β
n
x
,
v
^
o
x
,
y
,
t
=
-
y
φ
r
w
x
1
-
cos
⁡
ω
t
.
Equations (11) show the time history of the perturbed values variables P
^
o
,
u
^
o
,
v
^
o at any axial locations.Following the same policies, the resonance solutions are derived and formulated as(12)
P
^
o
x
,
t
=
γ
C
~
n
*
β
n
*
sin
β
n
*
t
-
β
n
*
t
cos
⁡
β
n
*
t
cos
⁡
β
n
*
x
,
u
^
o
x
,
t
=
C
~
n
*
-
2
β
n
*
cos
⁡
β
n
*
t
-
1
-
t
sin
β
n
*
t
×
cos
⁡
β
n
*
x
,
where C
^
n
*
=
2
∫
0
1
φ
r
w
(
x
)
cos
⁡
β
n
*
x.


## 2.2.1. Steady and Unsteady Flow

Hegab [7] has concluded that the total axial velocity component includes two main parts. The first part represents the steady state component and the other part represents unsteady state component and can be written as follows:
(8)
u
,
v
,
P
,
ρ
,
T
=
u
s
,
v
s
,
P
s
,
ρ
s
,
T
s
+
u
^
,
v
^
,
P
^
,
ρ
^
,
T
^
.
Using the expended variables used in Hegab [7] and inserting them in (1), the side wall steady vorticity generation (Ω
o
s
wall) and the accompanying heat transfer (∂
T
o
s
/
∂
y
wall) are derived and written as follows:
(9)
Ω
os
wall
=
-
1
γ
v
r
w
s
∂
∂
x
γ
π
2
4
∫
x
1
v
r
w
s
(
x
^
)
∫
0
x
^
v
r
w
s
(
τ
)
d
τ
d
x
^
+
∂
v
r
w
s
x
∂
x
,
∂
T
os
∂
y
wall
=
0
.Equations (9) show that nonzero steady vorticity and zero heat transfer are found with the use of no-slip and constant injection boundary conditions.The sequence of the numerical solution consists of two steps. The first one is the initial steady state computation which arises from constant injection with isothermal condition. The second one is the unsteady flow computation based on initial conditions as computed from the first step. The perturbation theory introduced by Staab et al. [11] and applied by Hegab [7] leads to the following weakly rotational equations:
(10)
t
=
0
,
P
^
o
=
0
,
∂
P
^
o
∂
t
=
0
,
x
=
0
,
∂
P
^
o
∂
x
=
0
,
x
=
1
,
P
^
o
=
0
,
y
=
0
;
v
^
0
=
0
,
y
=
1
;
v
^
o
=
-
ε
φ
r
w
(
1
-
cos
⁡
ω
t
)
t
>
0
.
The final nonresonance solution for the flow in the simulated chamber in Figure 1 is
(11)
P
^
o
x
,
t
=
∑
n
=
0
∞
2
ε
γ
ω
∫
0
1
φ
r
w
(
x
)
cos
⁡
β
n
x
ω
2
-
β
n
2
×
ω
β
n
sin
β
n
t
-
sin
ω
t
cos
⁡
β
n
x
,
u
^
o
x
,
t
=
∑
n
=
0
∞
2
ε
γ
ω
∫
0
1
φ
r
w
(
x
)
cos
⁡
β
n
x
ω
2
-
β
n
2
×
ω
β
n
cos
⁡
β
n
t
-
1
-
β
n
ω
cos
⁡
ω
t
-
1
×
sin
β
n
x
,
v
^
o
x
,
y
,
t
=
-
y
φ
r
w
x
1
-
cos
⁡
ω
t
.
Equations (11) show the time history of the perturbed values variables P
^
o
,
u
^
o
,
v
^
o at any axial locations.Following the same policies, the resonance solutions are derived and formulated as(12)
P
^
o
x
,
t
=
γ
C
~
n
*
β
n
*
sin
β
n
*
t
-
β
n
*
t
cos
⁡
β
n
*
t
cos
⁡
β
n
*
x
,
u
^
o
x
,
t
=
C
~
n
*
-
2
β
n
*
cos
⁡
β
n
*
t
-
1
-
t
sin
β
n
*
t
×
cos
⁡
β
n
*
x
,
where C
^
n
*
=
2
∫
0
1
φ
r
w
(
x
)
cos
⁡
β
n
*
x.


## 2.3. The Computational Approach

In the present study, higher order accuracy difference equations are employed to solve the parabolized Navier-Stokes equations in the interior of the prescribed chamber and at the boundaries. The numerical multisteps scheme in [24] is used to solve the 2D, unsteady, and compressible Navier-Stokes equations for the points located near the boundaries, while the two-four explicit, predictor-corrector scheme in [25] is employed to solve the adapted equations at the interior points. The latter method is a fourth-order accuracy variant of the fully explicit MacCormack scheme. Uniform grids are used. The accuracy criterion employed by MacCormack [24] is used to determine the size of the time step.Here, the NSCBC technique applied by [7] is presented in detail. The full derivation for the parabolized Navier-Stokes at the upper and lower boundaries and at the head and exit ends is formulated to specify the nonphysical conditions at these boundaries. Moreover, this method is applied to reduce the generated reflected waves arising from the numerical techniques.Based on the idea of characteristics for Euler equations, the first derivatives normal to the boundaries appearing in the convective terms of the Navier-Stokes equations may be expressed in terms of the amplitude of incoming, outcoming, and entropy waves{
℘
i
’s
}. The fluid dynamics equations can be written in nondimensional wave modes form as follows:
(13)
∂
ρ
∂
t
+
1
C
2
℘
2
+
1
2
℘
4
+
℘
1
+
M
∂
ρ
v
∂
y
=
0
,
(14)
∂
u
∂
t
+
1
2
ρ
C
℘
4
-
℘
1
+
M
v
∂
u
∂
y
=
χ
x
ρ
,
(15)
∂
v
∂
t
+
℘
3
+
M
v
∂
v
∂
y
+
A
2
γ
M
1
ρ
∂
p
∂
y
=
χ
y
ρ
,
(16)
∂
p
∂
t
+
1
2
℘
4
+
℘
1
+
M
v
∂
p
∂
y
+
γ
M
p
∂
v
∂
y
=
χ
E
,
where χ
x and χ
y represent the viscous terms and χ
E refers to the combination of the viscous and conduction terms in the axial and transverse directions and C
=
C
′
/
C
o
′. The amplitudes of the left-running waves (℘
1), the right-running waves (℘
4), and the entropy waves (℘
2
,
℘
3) in terms of the characteristic velocities λ
i’s are
(17)
℘
1
=
λ
1
1
γ
∂
p
∂
x
-
M
ρ
C
∂
u
∂
x
;
λ
1
=
u
M
-
C
,
(18)
℘
2
=
λ
2
C
2
∂
ρ
∂
x
-
1
γ
∂
p
∂
x
;
λ
2
=
M
u
,
(19)
℘
3
=
λ
3
∂
v
∂
x
;
λ
3
=
M
u
,
(20)
℘
4
=
λ
4
1
γ
∂
p
∂
x
+
M
ρ
C
∂
u
∂
x
;
λ
4
=
u
M
+
C
.Two approaches are used to infer the values of the incoming amplitude-variations from the boundary conditions. One is the local one-dimensional inviscid relations (LODI) technique and the second is the Navier-Stokes characteristics boundary conditions (NSCBC). For the latter one, the complete description of the numerical boundary conditions at all boundaries is discussed below.


## 2.4. Exit Boundary Conditions

In this case, a pressure node is assumed at the exit boundary. The NSCBC relation (16) suggests that the amplitude of the reflected wave should be
(21)
℘
1
=
2
χ
E
-
℘
4
,
and the remaining conditions are
(22)
℘
2
=
M
u
C
2
∂
ρ
∂
x
-
1
γ
∂
p
∂
x
,
℘
3
=
M
u
∂
v
∂
x
,
℘
4
=
u
M
+
c
1
γ
∂
p
∂
x
+
M
ρ
C
∂
u
∂
x
.
The nondimensional equations can be written as follows:
(23)
∂
ρ
∂
t
+
M
u
∂
ρ
∂
x
-
M
u
γ
C
2
∂
p
∂
x
+
M
v
∂
ρ
∂
y
+
A
2
M
γ
-
1
R
e
P
r
1
C
2
∂
2
T
∂
y
2
=
0
,
∂
u
∂
t
+
1
ρ
C
u
+
C
M
1
γ
∂
p
∂
x
+
M
ρ
C
∂
u
∂
x
+
M
v
∂
u
∂
y
-
A
2
γ
-
1
R
e
P
r
1
ρ
C
∂
2
T
∂
y
2
-
M
A
2
R
e
1
ρ
∂
2
u
∂
y
2
=
0
,
∂
v
∂
t
+
M
u
∂
v
∂
x
+
v
∂
v
∂
y
=
0
,
d
p
d
t
=
0
.Here again, the nondimensional variables (ρ, u, and v) can be computed using one-sided finite-difference approximations.


## 2.5. Head End Boundary Condition

The velocity component normal to the boundary(
u
(
0
,
y
,
t
)
) is zero since the head end is impermeable. In contrast the transverse velocity component cannot be specified because the axial viscous terms in the transverse momentum equation have been suppressed. As the normal velocity drops to zero at x
=
0, the amplitudes ℘
2 and ℘
3, in (18) and (19), are zero and the characteristic velocities are λ
1
=
-
C and λ
4
=
+
C. The reflected wave ℘
4 enters the domain, while the wave ℘
1 propagates out through the boundary. The NSCBC relation (19) indicates that the amplitude of the incoming wave ℘
4 should be
(24)
℘
4
=
℘
1
,
and the remaining conditions are
(25)
℘
1
=
-
c
1
γ
∂
p
∂
x
+
M
ρ
C
∂
u
∂
x
,
℘
2
=
0
,
℘
3
=
0
.
The equations can be written in nondimensional form as follows:
(26)
∂
ρ
∂
t
-
M
ρ
∂
u
∂
x
+
1
γ
C
∂
p
∂
x
+
M
∂
ρ
v
∂
y
=
0
,
∂
u
∂
t
=
0
,
∂
v
∂
t
+
M
v
∂
v
∂
y
=
0
,
∂
p
∂
t
-
c
∂
p
∂
x
+
M
ρ
C
∂
u
∂
x
+
γ
M
p
∂
v
∂
y
-
γ
A
2
M
γ
-
1
P
r
R
e
∂
2
T
∂
y
2
=
0
.
More details about the NSCBC techniques are found in [23, 26].


## 2.6. Code Validations

The computational codes are designed, constructed, and validated for steady and unsteady flow as shown in Figures2 and 3. For the latter one, the constructed codes run for many wave cycles to detect the whole unsteady rotational flow dynamic across the chamber.Figure 2
Normalized axial velocity profiles at the midlength of the chamber for the steady state flow.Figure 3
Comparison between the computational and the asymptotic solutions for the perturbed temperature time history withω
=
1 (nonresonance), ε
=
0.4, M
=
0.02, and R
e
=
3
*
10
5 at the half-length of the chamber.In Figure2, comparison between the current numerical and asymptotic results dealing with the steady state flow solutions is carried out. The steady state flow is generated by adding constant side wall mass injection (
V
r
w
=
-
1
) from the injection surface of the prescribed chamber. Similar experiment in [18, 27] was carried out and compared with the current mathematical models in Figure 2. The asymptotic analysis solution represents the solid line, while the dotted line represents the converged steady state numerical solution for viscous flow with wall Reynolds number R
e
w
=
5
*
10
3. The nonfilled points refer to the experimental measured data. The experimental results are recorded for flow in a channel with fully transpired wall. Three different porous surfaces using three different sizes of honeycomb cells are used, [7]. In this figure, a comparison between the analytical and the numerical solutions at the midlength of the chamber is illustrated. This comparison shows a reasonable agreement between the three approaches. The comparison between the theoretical and the computational approaches shows little bit deviation. This deviation is logically expected since the used honey combs have three different sizes. The smaller size of honey combs, the smaller deviations between the two approaches.The second set of the validations show the solution for the unsteady flow computations at nonresonance forced mode. In this set of the results, the computational results are validated by comparing with the asymptotic analysis. These results represent the asymptotic solution when the sidewall injection velocity is represented byv
r
w
(
x
,
1
,
t
)
=
-
1
+
ε
φ
r
w
(
1
-
cos
⁡
ω
t
), φ
r
w
=
cos
⁡
β
n
x. In the current work and the work by Hegab [7], the disturbance forcing frequency is taken to be ω
=
1 (nonresonance) with the amplitude of ε
=
0.4 and the eigenvalues (β
n
=
n
π
/
2, n
=
1,3
,
5,7
,
9
,
…
,
etc.). The temperature time history (T
^) at the medlength location is presented in Figure 3 for R
e
=
3
*
10
5, M
=
0.02, δ
=
20, and n
=
1.0. For the nonresonance case, the solution shows stable conditions but it is not purely harmonic. This interesting behavior shows the ability of the constructed code with the good treatment to the boundary conditions to detect the non-harmonic trend as detected by the analytical approach. A comparison of the perturbed temperature (T
^) (27) with the numerical results in Figure 3 shows good agreement between the two approaches especially for t
<
20, while the computational solution shows little less deviation than the analytical results for t
>
20 due to the nonlinearity effects:
(27)
T
^
o
x
,
t
=
∑
n
=
0
∞
2
ε
ω
γ
-
1
∫
0
1
φ
r
w
(
x
)
cos
⁡
β
n
x
ω
2
-
β
n
2
×
ω
β
n
sin
β
n
t
-
sin
ω
t
cos
⁡
β
n
x
.Generally, the coexistence of the forced and the natural modes of the system is detected in the adapted computational and analytical approaches as seen in the earlier investigation by Hegab [7]. In contrast, the computational techniques in [10] showonly pure harmonic oscillations. The source of these deficiencies may be arising from the treatment of the boundary conditions especially at the head and aft ends of the simulated chamber.


## 3. Results and Discussions

The impact of the injected wave number from the sidewall along with the frequencies of the forced mode on the generated complex acoustical fields and their interaction with the fluid are examined. The results have been categorized into the following two sets.


### 3.1. Nonresonance Solution

In the current study, a wide range of forcing frequencies away from or very close to the first fundamental modes of the prescribed chamber is considered [7], a wide range of forcing frequencies away from or very close to the first fundamental modes of the prescribed chamber is considered. In Figure 4, the irrotational and rotational fields across the whole chamber are presented when t
=
42 for ω
=
1, R
e
=
3
*
10
5, M
=
0.02, ε
=
0.4, n
=
1, 3, and 5, and λ
n
=
π
/
2. The latter one represents the Eigen values of the chamber. These interesting graphs show that the rotational fields are generated at the injection surface. As time increases the generated unsteady vorticity is transferred to fill most of the combustion chamber. Moreover, the constructed codes are able to run for enough time and count about seven vorticity wave cycles as t
=
42. It is also seen that the amplitude of the rotational field attains maximum values near the wall and found to be order of ten for the monotonic variation (n
=
1) and nonresonance frequency (
ω
=
1
). One notes that these magnitudes meet the scaling of the asymptotic solutions in [7, 9, 11, 14]. Moreover, for n
=
1, monotonic distribution of unsteady vorticity is noticed, while the variation goes in nonmonotonic axial distribution for n
>
1. For the nonresonance solution, the longitudinal variation of the vorticity at the wall ∂
Ω
w
(
x
,
t
)
/
∂
x
>
0 for λ
n
=
π
/
2 (n
=
1) is revealing that the magnitude of rotational flow fields on the injection surface increases as axial distance increases.Vorticity topography forω
=
1 (nonresonance) and ε
=
0.4, M
=
0.02, and R
e
=
3
*
10
5 at t
=
42 at several axial wave numbers (a) n
=
1, (b) n
=
3, and (c) n
=
5.
(a)
(b)
(c)The vorticity contours adjacent to the injection surface for the same parameters as in Figure4, for several eigenvalues (a) λ
n
=
π
/
2 (n
=
1), (b) 3
π
/
2 (n
=
3), and (c) 5
π
/
2 (n
=
5), are present in Figure 5. It is clearly seen that, adjacent to the combustion surface, the generated unsteady vorticity may impact the surface at downstream locations for n
=
1, at two locations for n
=
3, and many locations for n
=
5. The rotational flow fluctuates from negative to positive and vice versa at several successive times and may cause scouring on the combustion surface and, in turn, may increase the shearing stresses adjacent to the combustion surface. This important conclusion may lead to increase of the burning rates over the designed ones, especially for low wave number, since the amplitude of the unsteady vorticity decreases as the wave number increases.The vorticity contours adjacent to the injection surface for the same parameters as in Figure4, for several eigenvalues (a) λ
n
=
π
/
2 (n
=
1), (b) 3
π
/
2 (n
=
3), and (c) 5
π
/
2 (n
=
5).
(a)
(b)
(c)Fortunately, recent experimental work by Meteb [28] in 2009 was directed to study the acoustic instability in hybrid rocket system using different fuel grains and ignition systems. The experimental results verified qualitatively the current analytical and computational results for the generated complex wave patterns and the effect of the combustion surface on the acoustic fluid dynamics interaction mechanism. In this experiment the pure polyethylene (PE) fuel grain and PE with additives of aluminum (Al) powder with different percentage are used to predict the hybrid rocket motor regression rates. The flame geometry and structure for the burning of the pure PE and PE with additives of aluminum (Al) powder at 2.5%, 7.5%, and 15% of the total weight of the whole samples are illustrated in Figure 6. The effect of the Al powder on the flame length, geometry, and the rate of reaction for the burning of each sample is clearly seen. In addition, the pressure time histories at a point before the combustion chamber, in the way of oxidizer, are recorded and illustrated in Figure 7 for the pure PE and in Figures 8 and 9 for PE with additives of aluminum (Al) powder at 7.5% and 15% of the total weight, respectively. It is clearly seen that the amplitude of the pressure waves decreases as the percentage of Al fuel grain increases in the sample. This interesting behavior is similar to what we have seen when we increase the wave number of the combustion surface. As we know, most of the Al fuel grains burn away from the combustion surface consisting of many fragments. The escaping processes of Al fuel grains from the combustion surface generate many of cavities that make the combustion surface more rough than the smooth one in case of pure PE. The experimental and computational results are qualitatively in agreement. For the first time, the analytical and computational results for the current cold model and the acoustic instabilities found by Meteb [28] state an important and essential conclusion about the effect of the combustion surface morphology on the acoustic instability inside the rocket chamber and the accompanying rotational flow. Instantaneous total velocity vector map for the parameters in Figure 4 at (a) t
=
0 “steady state solution” and (b) t
=
16 is presented in Figures 10(a) and 10(b). The injection surface is located at y
=
1 on the upper part of the plot, where the flow is completely vertical, while the lower part of the plot represents the center line at y
=
0, where the flow is completely axial. The closed end is located at x
=
0 on the left hand side of the plot, while at x
=
1 it represents the exit plane on the right hand side. One may notice the smooth and undisturbed fluid flow through the entire map for the steady state solution at t
=
0. After initiating the side wall disturbances, the coexisting of acoustics and the generated rotational flow is clearly seen in the upper half part of the chamber and illustrates how the interaction between the travelling acoustic waves and the injected fluid affects the fluid flow in the chamber. In addition, at t
=
32 in Figure 11, there is an instantaneous reverse flow found at a location near the combustion surface. The complex acoustic fluid dynamics interaction structure may increase the shear stresses adjacent to the combustion surface and through the whole cavity. These results may provide the rocket designers with conceptual tools and data to be considered in case of redesigning strategies.The flame geometry for the burning of hybrid rocket systems using (a) pure PE fuel grain, (b) 97.5% PE +2.5% Al fuel grain, (c) 92.5% PE +7.5% Al fuel grain, and (d) 85% PE +15% Al fuel grain [18].
(a)
(b)
(c)
(d)Figure 7
The recorded pressure time history for a burning of hybrid motor using pure polyethylene (PE) fuel grain.Figure 8
The recorded pressure time history for a burning of hybrid motor using 7.5% aluminum (Al) additive to PE fuel grain.Figure 9
The recorded pressure time history for a burning of hybrid motor using 15% aluminum (Al) additive to PE fuel grain.Instantaneous total velocity vector map for the parameters in Figure4:  (a)  t
=
0 “steady state solution” and (b) t
=
16.
(a)
(b)Figure 11
Instantaneous total velocity vector map att
=
32 for the parameters in Figure 4.


### 3.2. Near and at Resonance Solution

A comparison of the asymptotic analysis perturbed pressure time history with the computational oneat and near resonance frequencies is presented in Figure 12. For t
<
10 the results show good agreement between the two approaches even with the asymptotic beat condition (ω
=
π
/
2
-
0.0018) which is very close to the resonance frequency. While for t
>
10 the asymptotic solution shows linear growth with time based on the linear theory for the derived wave equation, where ω
=
β
n
*
=
π
/
2 as derived by Hegab [7]. In contrast, the computational solution shows beats at the resonance condition due to the nonlinearity and the asymptotic beat solutions have smaller amplitude than the resonance one. The beats appear due to the interaction between the deriving frequency mode ω
=
π
/
2
-
0.0018 and the first eigenfunction mode (
π
/
2
). It is clearly seen that the amplitude of the vorticity at the resonance frequency is found to be ten times than the amplitudes at the nonresonance one in Figure 4. Moreover, it is found that the fluctuation of the transverse temperature gradient at the injection surface is found to be about ten times than the fluctuation of the temperature disturbances. This concludes that since the heat transfer along the sidewall is directly proportional to the transverse temperature gradient there is an unexpectedly large heat flux which exists both at the wall and within the entire chamber given the small size of the temperature distribution. The asymptotic analysis by Staab et al. [11] is used to prove this important conclusion.Figure 12
Comparison between the computational and the asymptotic solutions for the perturbed pressure time history withω
=
π
/
2 (resonance), ε
=
0.4, M
=
0.02, and R
e
=
3
*
10
5 at the half-length of the chamber.Finally, the asymptotic and computational methods are employed to study the dynamic process in an acoustically simulated solid propellant rocket chamber with transient sidewall mass addition at different wave numbers.


## 3.1. Nonresonance Solution

In the current study, a wide range of forcing frequencies away from or very close to the first fundamental modes of the prescribed chamber is considered [7], a wide range of forcing frequencies away from or very close to the first fundamental modes of the prescribed chamber is considered. In Figure 4, the irrotational and rotational fields across the whole chamber are presented when t
=
42 for ω
=
1, R
e
=
3
*
10
5, M
=
0.02, ε
=
0.4, n
=
1, 3, and 5, and λ
n
=
π
/
2. The latter one represents the Eigen values of the chamber. These interesting graphs show that the rotational fields are generated at the injection surface. As time increases the generated unsteady vorticity is transferred to fill most of the combustion chamber. Moreover, the constructed codes are able to run for enough time and count about seven vorticity wave cycles as t
=
42. It is also seen that the amplitude of the rotational field attains maximum values near the wall and found to be order of ten for the monotonic variation (n
=
1) and nonresonance frequency (
ω
=
1
). One notes that these magnitudes meet the scaling of the asymptotic solutions in [7, 9, 11, 14]. Moreover, for n
=
1, monotonic distribution of unsteady vorticity is noticed, while the variation goes in nonmonotonic axial distribution for n
>
1. For the nonresonance solution, the longitudinal variation of the vorticity at the wall ∂
Ω
w
(
x
,
t
)
/
∂
x
>
0 for λ
n
=
π
/
2 (n
=
1) is revealing that the magnitude of rotational flow fields on the injection surface increases as axial distance increases.Vorticity topography forω
=
1 (nonresonance) and ε
=
0.4, M
=
0.02, and R
e
=
3
*
10
5 at t
=
42 at several axial wave numbers (a) n
=
1, (b) n
=
3, and (c) n
=
5.
(a)
(b)
(c)The vorticity contours adjacent to the injection surface for the same parameters as in Figure4, for several eigenvalues (a) λ
n
=
π
/
2 (n
=
1), (b) 3
π
/
2 (n
=
3), and (c) 5
π
/
2 (n
=
5), are present in Figure 5. It is clearly seen that, adjacent to the combustion surface, the generated unsteady vorticity may impact the surface at downstream locations for n
=
1, at two locations for n
=
3, and many locations for n
=
5. The rotational flow fluctuates from negative to positive and vice versa at several successive times and may cause scouring on the combustion surface and, in turn, may increase the shearing stresses adjacent to the combustion surface. This important conclusion may lead to increase of the burning rates over the designed ones, especially for low wave number, since the amplitude of the unsteady vorticity decreases as the wave number increases.The vorticity contours adjacent to the injection surface for the same parameters as in Figure4, for several eigenvalues (a) λ
n
=
π
/
2 (n
=
1), (b) 3
π
/
2 (n
=
3), and (c) 5
π
/
2 (n
=
5).
(a)
(b)
(c)Fortunately, recent experimental work by Meteb [28] in 2009 was directed to study the acoustic instability in hybrid rocket system using different fuel grains and ignition systems. The experimental results verified qualitatively the current analytical and computational results for the generated complex wave patterns and the effect of the combustion surface on the acoustic fluid dynamics interaction mechanism. In this experiment the pure polyethylene (PE) fuel grain and PE with additives of aluminum (Al) powder with different percentage are used to predict the hybrid rocket motor regression rates. The flame geometry and structure for the burning of the pure PE and PE with additives of aluminum (Al) powder at 2.5%, 7.5%, and 15% of the total weight of the whole samples are illustrated in Figure 6. The effect of the Al powder on the flame length, geometry, and the rate of reaction for the burning of each sample is clearly seen. In addition, the pressure time histories at a point before the combustion chamber, in the way of oxidizer, are recorded and illustrated in Figure 7 for the pure PE and in Figures 8 and 9 for PE with additives of aluminum (Al) powder at 7.5% and 15% of the total weight, respectively. It is clearly seen that the amplitude of the pressure waves decreases as the percentage of Al fuel grain increases in the sample. This interesting behavior is similar to what we have seen when we increase the wave number of the combustion surface. As we know, most of the Al fuel grains burn away from the combustion surface consisting of many fragments. The escaping processes of Al fuel grains from the combustion surface generate many of cavities that make the combustion surface more rough than the smooth one in case of pure PE. The experimental and computational results are qualitatively in agreement. For the first time, the analytical and computational results for the current cold model and the acoustic instabilities found by Meteb [28] state an important and essential conclusion about the effect of the combustion surface morphology on the acoustic instability inside the rocket chamber and the accompanying rotational flow. Instantaneous total velocity vector map for the parameters in Figure 4 at (a) t
=
0 “steady state solution” and (b) t
=
16 is presented in Figures 10(a) and 10(b). The injection surface is located at y
=
1 on the upper part of the plot, where the flow is completely vertical, while the lower part of the plot represents the center line at y
=
0, where the flow is completely axial. The closed end is located at x
=
0 on the left hand side of the plot, while at x
=
1 it represents the exit plane on the right hand side. One may notice the smooth and undisturbed fluid flow through the entire map for the steady state solution at t
=
0. After initiating the side wall disturbances, the coexisting of acoustics and the generated rotational flow is clearly seen in the upper half part of the chamber and illustrates how the interaction between the travelling acoustic waves and the injected fluid affects the fluid flow in the chamber. In addition, at t
=
32 in Figure 11, there is an instantaneous reverse flow found at a location near the combustion surface. The complex acoustic fluid dynamics interaction structure may increase the shear stresses adjacent to the combustion surface and through the whole cavity. These results may provide the rocket designers with conceptual tools and data to be considered in case of redesigning strategies.The flame geometry for the burning of hybrid rocket systems using (a) pure PE fuel grain, (b) 97.5% PE +2.5% Al fuel grain, (c) 92.5% PE +7.5% Al fuel grain, and (d) 85% PE +15% Al fuel grain [18].
(a)
(b)
(c)
(d)Figure 7
The recorded pressure time history for a burning of hybrid motor using pure polyethylene (PE) fuel grain.Figure 8
The recorded pressure time history for a burning of hybrid motor using 7.5% aluminum (Al) additive to PE fuel grain.Figure 9
The recorded pressure time history for a burning of hybrid motor using 15% aluminum (Al) additive to PE fuel grain.Instantaneous total velocity vector map for the parameters in Figure4:  (a)  t
=
0 “steady state solution” and (b) t
=
16.
(a)
(b)Figure 11
Instantaneous total velocity vector map att
=
32 for the parameters in Figure 4.


## 3.2. Near and at Resonance Solution

A comparison of the asymptotic analysis perturbed pressure time history with the computational oneat and near resonance frequencies is presented in Figure 12. For t
<
10 the results show good agreement between the two approaches even with the asymptotic beat condition (ω
=
π
/
2
-
0.0018) which is very close to the resonance frequency. While for t
>
10 the asymptotic solution shows linear growth with time based on the linear theory for the derived wave equation, where ω
=
β
n
*
=
π
/
2 as derived by Hegab [7]. In contrast, the computational solution shows beats at the resonance condition due to the nonlinearity and the asymptotic beat solutions have smaller amplitude than the resonance one. The beats appear due to the interaction between the deriving frequency mode ω
=
π
/
2
-
0.0018 and the first eigenfunction mode (
π
/
2
). It is clearly seen that the amplitude of the vorticity at the resonance frequency is found to be ten times than the amplitudes at the nonresonance one in Figure 4. Moreover, it is found that the fluctuation of the transverse temperature gradient at the injection surface is found to be about ten times than the fluctuation of the temperature disturbances. This concludes that since the heat transfer along the sidewall is directly proportional to the transverse temperature gradient there is an unexpectedly large heat flux which exists both at the wall and within the entire chamber given the small size of the temperature distribution. The asymptotic analysis by Staab et al. [11] is used to prove this important conclusion.Figure 12
Comparison between the computational and the asymptotic solutions for the perturbed pressure time history withω
=
π
/
2 (resonance), ε
=
0.4, M
=
0.02, and R
e
=
3
*
10
5 at the half-length of the chamber.Finally, the asymptotic and computational methods are employed to study the dynamic process in an acoustically simulated solid propellant rocket chamber with transient sidewall mass addition at different wave numbers.


## 4. Conclusion Remarks

The current study emphasizes the effect of the combustion surface morphology of the solid rocket system on the acoustic instability and the accompanying vorticity and temperature dynamics generation. The rotational and irrotational fields generated in the prescribed solid rocket motor chamber with Mach number are chosen to beO(10−1) while the characteristic Reynolds number is O(105-106) at nonresonance and near and at resonance frequencies have been described. The NSCBC techniques that applied in [7, 29] are used to treat the numerical boundary conditions.The numerical solution shows that the convergence to steady state solution is achieved first and then additional transient side wall mass addition is initiated. Moreover, the result shows that the forced mode of the nonresonance and resonance frequencies generates longitudinal acoustics waves which interact again with its sources at the injection surface to generate the rotational flow fields. The generated unsteady vorticity is transferred from the injection surface toward the centerline to fill the entire chamber as time increases after more than seven wave cycles. The solution shows nonpure harmonic oscillations as suggested from the theory in [11]. Vorticity 3D graphs and contours illustrate that spatially the maximum amplitudes are found to be near the injection surface. Moreover, the largest shear stresses are found to be at the injection surface near the exit end, while the shear stresses impact the injection surface at many locations as the wave number increases. The longitudinal variation of the rotational flow fields is found to be nonmonotonic for n
>
1.Here, for the first time, the analytical and computational results for the current cold model and the recent experimental work by Meteb [28] revealed that the combustion surface morphology has a great effect on the acoustic instability inside the rocket chamber and the accompanying rotational flow. The computational approach reveals that there are unexpectedly high heat transfer and generated rotational flow fields found to be adjacent to the injection surface in both nonresonance and resonance cases. The amplitude of the latter one is found to be approximately ten times than that with the nonresonance solution. One may anticipate that these large oscillatory gradients will persist to the surface of a burning solid propellant in a real rocket chamber.More traditional computational studies of turbulence in channel flows with wall injection are reviewed by Chaouat [30]. Mean axial velocity profiles are found to be laminar until about halfway down the channel and then show signs of transition. Chaouat’s results do not show the acoustic-related effects discussed in the present work. As a result, it is difficult to know from the turbulent theory results if the source of the downstream vorticity is classical transition to turbulence, or partly the result of the kind of mechanism discussed in this paper. In our recent studies in 2011 and 2012 [31, 32] the assessment of turbulence modeling for gas flow in the chamber and in the convergent-divergent rocket nozzle for the steady state flow is studied. The detailed information and realities in the current research about the transient coexisting of the irrotational (acoustics), rotational (generated vorticity from the interaction between the injected flow and the traveling acoustic waves), and our turbulent techniques that we have adapted in our studies in 2011 and 2012 for the steady flow field in SRM chamber/nozzle geometry will be used to account all the complex flow fields in our future study.Finally, results of the asymptotic and computational approaches may provide the rocket designer with codes and information making them able to interpret some reasons behind the combustion instability inside the combustion chamber of real SRM based on the composition of the solid propellants and in turn the morphology of the combustion surface.



---

*Source: 485302-2015-01-12.xml*