Convective Heat Transfer in the Reusable Solid Rocket Motor of the Space
Transportation System
Rashid A. Ahmad
Motor Performance Department, Gas Dynamics Section, M/S 252
Science and Engineering
ATK Thiokol Propulsion, A Division of ATK Aerospace Company Inc., Brigham City, Utah
84302
This simulation involved a two-dimensional axisymmetric model of a full motor initial grain
of the Reusable Solid Rocket Motor (RSRM) of the Space Transportation System (STS). It was
conducted with CFD commercial code FLUENT.® This analysis was performed to: a) maintain
continuity with most related previous analyses, b) serve as a non-vectored baseline for any three-
dimensional vectored nozzles, c) provide a relatively simple application and checkout for various
CFD solution schemes, grid sensitivity studies, turbulence modeling and heat transfer, and d)
calculate nozzle convective heat transfer coefficients. The accuracy of the present results and the
selection of the numerical schemes and turbulence models were based on matching the rocket
ballistic predictions of mass flow rate, head end pressure, vacuum thrust and specific impulse,
and measured chamber pressure drop. Matching these ballistic predictions was found to be good.
This study was limited to convective heat transfer and the results compared favorably with
existing theory. On the other hand, qualitative comparison with backed-out data of the ratio of
the convective heat transfer coefficient to the specific heat at constant pressure was made in a
relative manner. This backed-out data was devised to match nozzle erosion that was a result of
heat transfer (convective, radiative and conductive), chemical (transpirating), and mechanical
(shear and particle impingement forces) effects combined.
Presented as Paper AIAA-200 1-3585, AIAA/ASME/SAE/ASEE 37th Joint Propulsion
Conference, Salt Lake City, UT, July 8-1 1, 2001.
*Sr. Principal Engineer, Associate Fellow AIAA.
®2002 ATK Thiokol Propulsion, a Division of ATK Aerospace Company Inc., Published by the
American Institute of Aeronautics and Astronautics
Nomenclature
A area, m" (ft )
a empirical constant (Eq. 3)
Cf friction coefficient
Cp , Cv specific heats at constant pressure and volume, respectively, J/kg-K (Btu/lbm-R)
COD coefficient of determination of a correlation
d diameter, m (in.)
2 2
g acceleration due to gravity, m/s (ft/s )
h convective heat transfer coefficient, W/m'-K (Btu/hr-ft"-R)
/ turbulence intensity, «'/«„ (%)
K an acceleration parameter
k gas thermal conductivity, W/m-K (Btu/hr-ft-R)
/ mixing length constant usually known as VonKarman length scale, m (ft)
fn mass flow rate, kg/s (Ibm/s)
MW molecular weight, kg/kgmole (Ib/lbmole)
M Mach number
n burning rate pressure exponent
Po , P total and static pressure, respectively, Pa (psia)
Pr Prandtl number
Prt turbulent Prandtl number
q wall heat flux, W/m" (Btu/hr-ft")
Re(x) local axial Reynolds number based on distance along the nozzle wall, pjx) «<» (x) Xw /
jloix)
Red(x) local axial Reynolds number based on diameter of the nozzle, pjx) m„ (x) d(x) / /ij(x)
SK recovery factor
^ propellant bum rate, m/s (in./s)
r radius, m (in.)
S source term (Eq. I)
To , T total and static temperature, respectively, K (R)
r"^ non-dimensional temperature in wall coordinates
u, V axial and radial velocity components, m/s (ft/s)
M„ velocity magnitude of u and v at the motor centerline, m/s, ft/s
u^ non-dimensional velocity in wall coordinates
u ' velocity fluctuation, m/s (ft/s)
Ui shear velocity, m/s (ft/s)
w propellant weight flow rate, N/s (Ibf/s)
X, r axial and radial axes, m (in.)
(x), (r) function of axial and radial directions, respectively
y' non-dimensional distance from the wall in wall units
a augmentation/de-augmentation factor in propellant mass flux
y ratio of specific heats, Q/Cv
K, e turbulence kinetic energy {mis (ft^/s^)) and its rate of dissipation (m7s /s (ft7s7s)),
respectively
ji, V dynamic (N-s/m^ (Ibm/ft-s)) and kinematic (m7s (ft"/s)) viscosity, respectively
p gas density, kg/m (Ibm/ft )
Xw wall shear stress, N/m" (Ibf/ff)
Subscripts
aw adiabatic wall or recovery temperature
e at nozzle exit
/ face
g gas
h hydraulic
o chamber conditions
p propellant
ref reference
s static
t turbulent
vac vacuum
w wall
<x> conditions at motor centerline
Superscripts
* throat conditions
flux
- mean and time averaged
Introduction
A computational fluid dynamics (CFD) simulation involved a two-dimensional axisymmetric
model of a full motor initial grain of the Resusable Solid Rocket Motor (RSRM) was conducted
with CFD commercial code FLUENT.'' ^ This analysis was performed to a) Maintain continuity
with the most related previous analyses ^'^ of this motor, b) Serve as a non-vectored baseline for
any three-dimensional vectored nozzles; c) Provide a relatively simple application and checkout
for various CFD solution schemes, grid sensitivity studies, turbulence modeling/heat transfer,
etc.; and d) Calculate the nozzle convective heat transfer coefficients. Nozzle heat transfer is of
interest because the convective heat transfer coefficient to the specific heat at constant pressure
ratio (h/Cp) is usually used as input in the Charring and Material Ablation (CMA) code for
nozzle erosion predictions.
The following are the four most related studies of this motor. First, Golafshani and Loh
conducted a time-dependent, axisymmetric numerical solution of the Navier-Stokes to analyze
the viscous coupled gas-particle non-reacting flow in solid rocket motors. The solution assumed
laminar internal flow. Second, Loh and Chwalowski used particles of diameters of 1 to 100 (im
in converging-diverging nozzles and a mass loading of 28.8%. Acceleration between Ig and 3g
had a minimal effect on the particles' behavior in the nozzle. Third, Whitesides et al. conducted
a two-dimensional axisymmetric two-phase flow analysis in the RSRM. The overall objective
was to determine the structure of the flow field in the recirculation region underneath the
submerged nozzle nose and to define the gas flow and particle impingement environments along
the surface of the aft case dome insulation. It was concluded that particles were impacting the
area undemeath the nozzle nose and forming a sheet of molten aluminum oxide or slag. The
sheet flows afterwards, along the undemeath nozzle nose surface as is the direction of the near
surface velocity vector during the last half of motor bum. This slag layer is then sheared from
the nozzle cowl/boot ring surface and impacts the aft dome case insulation at the location of
severe erosion. Fourth, Laubacher^ conducted two-dimensional axisymmetric analysis to
compute chamber pressure drop in the RSRM. The walls were assumed to be adiabatic and
utilized the standard K-e turbulence model and a coupled solver using in-house code, SHARP.
These studies were conducted for the RSRM and no attempt was made to calculate the
convective heat transfer coefficients. The fifth study' has the details of this study. It involved 2-
D axisynunetric and 3-D vectored nozzles. Furthermore, it involved two-phase flow where slag
concentration, accretion rates, and particle trajectories were calculated. At this time, it is not
feasible how to augment convective heat transfer by effects of particle impingement and
trajectories.
Related convective heat transfer studies are given in Refs. (9-16). When considering
convective heat transfer in solid rocket motors, surface temperatures and heat fluxes are high and
very difficult to measure. Ablative materials are used to dissipate and inhibit heat transfer by
erosion and transpiration. For the lack of reliable thermal conditions, the nozzle wall was usually
assumed to be adiabatic in CFD calculations. On the other hand, CFD calculations (velocity,
density, pressure, temperature, viscosity, etc.) and geometry enable someone to calculate heat
transfer.^ It is usually estimated^ using three well-known methods. They are the modified
Reynolds' analogy'^ for laminar flow over a flat plate, Dittus-Boelter correlation for fully-
developed turbulent pipe flow, and the Bartz^° correlation for nozzle flows. Bartz extended the
well-known Dittus-Boelter correlation''' for turbulent pipe flow to account for mass flux and
variations in velocity and temperature. Back et al."' '^"'^ conducted analytical and experimental
convective heat transfer studies in the Jet Propulsion Laboratory (JPL) nozzle. Moretti and
Kays'" conducted experimental convective heat transfer to an essentially constant property
turbulent boundary layer in a two-dimensional channel for various rates of free-stream
acceleration. Back et al. ''^'^ and Moretti and Kays'^ found that acceleration causes a depression
in heat transfer rate below what would be predicted assuming a boundary-layer structure such as
obtained for constant free-stream velocity. They attributed it to re-laminarization of the turbulent
boundary layer. They further state, it is by no means obvious that the same acceleration
parameter applicable to an axisymmetric flow. Wang focused on the capability of general-
purpose CFD codes in predicting convective heat transfer coefficients between a fluid and a solid
surface. Effects of various parameters such as grid resolution, turbulence models and near-wall
treatments, as well as numerical schemes on the accuracy of predicted convective heat transfer
were studied. Test cases included flat plate, pipe flow, JPL nozzle, and impinging jets.
The attributes of this two-dimensional (2D) axisymmetric analysis are:
• Significant effort was made to assess grid sensitivity and grid consistency with turbulence
models.
• Verifying flow/thermal solution quality represented by vs. u* vs. y'^ and i^ vs. y* against the
usual velocity and thermal laws of the wall for incompressible turbulent flows, respectively.
• Calculations of nozzle convective heat transfer, including the assessment of turbulent
boundary layer re-laminarization.
Discussion of Modeling Approach
Governing equations, geometry parameters, operating and predicted ballistic conditions, gas
thermophysical properties, grid density and turbulence modeling, boundary conditions,
computational schemes and numerical convergence (residuals) are discussed next.
Governing Equations: The numerical studies considered the solution of the Navier-Stokes
equations, energy equation, the turbulence kinetic energy with its rate of dissipation equations,
mass transfer and the necessary constitutive equations (ideal gas law, power law for gas thermal
conductivity and viscosity, etc.). The general governing equation was
v.(py^-r^v^)=5, (1)
and the mass conservation equation
V . (p v) = (2)
where ^can be velocity components (u, v, w), enthalpy (0 and turbulence quantities (k, e); Tis
an exchange coefficient for ^, 5^ is a source term for ^per unit volume.
Geometry Parameters: Table 1 gives a summary of geometry parameters. They are the
normalized chamber and exit radii along with chamber and exit area ratios.
Operating and Predicted Ballistic Conditions: Table 1 gives a summary of the ballistic
prediction parameters for the RSRM that include head end pressure and FSM-9 measured
chamber pressure drop. '^ The accuracy and the selection of the schemes, models and results are
based on matching the above ballistic prediction parameters in addition to mass flow rate, and
18
vacuum thrust and specific impulse.
Gas Thermophysical Properties: The total head-end pressure given in Table 1 along with
20
propellant formulation was used as input to the NASA-Lewis program or the ODE module of
the SPP code^' to obtain chamber gas temperature (To), dynamic viscosity (jj,,), specific heat at
constant pressure (Cp), thermal conductivity (ko), and molecular weight. They are given in Table
1. The local gas dynamic viscosity and thermal conductivity were calculated as a function of
local temperature as given at the bottom of Table 1 .
Grid Density and Turbulence Modeling: Table 2 gives the turbulence models used with the
desired values for wall y"" and the pertinent results. A coarse and fine grids with quadrilateral
cells used in this simulation and are summarized in Table 2. They were designed, solved, and
iterated on to give the desired values for y'^ so that consistency with the turbulence models was
achieved as given in Table 2. All the grids were generated by using GRIDGEN " and made
orthogonal and smoothed from one domain to another.
Boundary Conditions: Boundary conditions are applied at the propellant surface, nozzle exit,
and walls and discussed as follows:
(^^ At the prnppllant xurfac.e: Mass flux was calculated as a function of the local static pressure
as
where
(3b)
a = ^-f
[Pj
In addition, uniform chamber temperature, flow direction that was normal to the propellant
surface, an assumed turbulence intensity, (/), and hydraulic diameter {dh) were specified. The
augmentation factor, a, was used as 1 for the propellant except in the head end fin region, where
it was increased to 4.528 to account for the three-dimensional fins modeled in two-dimensional
axi-symmetric analysis.
( T\ At exit: A supersonic boundary condition was utilized where the quantities (P, T, u, v, k, e)
were calculated from cells upstream of the exit. The exit pressure, temperature, turbulence
intensity, and exit hydraulic diameter were specified to start the calculation. The exit pressure
and temperature were updated as the solution proceeded.
(^^ At wall: Three wall boundary conditions are used and discussed as follows:
(a) Velocity wall boundary condition was assumed to be no slip condition.
(b) Thermal wall boundary conditions were assumed for the submerged and the converging-
diverging part of the nozzle walls. The submerged wall was assumed to be isothermal at 2938.5
K (5289.3 R). On the other hand, the nozzle wall was assumed to be non-isothermal. The
surface temperature was taken from Ref. 23 and curve-fitted by this author using
TableCurve2D^'* as
T{x) =
a + c x+ e x^ + g x^ + i X* + k x^
\ + bx + dx^+fx^+hx* + jx^
(4)
where x was taken along the nozzle surface and where the coefficients are given as follows: a =
2789.03, b = - 4.61, c = -13307.21, J = 11.96, e = 32905.45,/= - 6.33, g = - 20712.33, h = -
0.916,7 = 1 174.06 and )fc = 815.51. The correlation coefficient was calculated to be r^ = 0.997. A
user defined function (UDF) was used to compile the specified surface temperature profile.
Assuming this surface temperature profile enables the calculation of the heat flux that in turn
enables the calculation of the convective heat transfer coefficient depending on the assumption of
a reference temperature. Calculated centerline and recovery temperatures are also shown and
will be discussed shortly.
(c) Near wall turbulence treatment: Two equation turbulence models were used. They
are the standard k-e and RNG K-e. Near wall treatment involved the standard wall functions and
two-layer zonal models as described in Table 2. The standard wall functions are given in terms
of the non-dimensional velocity in wall coordinates (m"^), non-dimensional distance from wall in
wall coordinates (y'^) and the frictional velocity («r) are defined' as
, u(r) , M, {r^. - r)
A/2
u =
u.
> y
, u^
P[r)
P.
(5a)
for velocity and
10
r =
(5b)
-lexit plane
for temperature at the exit plane.
In the two-layer zonal model, ' " the wall functions are not used. Rather, the turbulent
boundary layer is divided into two layers distinguished by a wall-distance-based turbulent
Reynolds number {Re, = p k'^ y / /i). The first layer is adjacent to the wall and is called a
viscous sublayer (viscosity (v) is much larger than eddy diffusivity for momentum (em)) where
the low standard k-£ turbulence model (Lx)w-/?e) is used. Only the k equation is solved in the
viscosity affected region while e is computed from a length scale correlation. The second layer is
fully turbulent (v « ew) and where the high standard K-e turbulence model (High-i?e) is used.
This model requires finer mesh resolution and therefore larger CPU and memory resources.
Computational Schemes: The segregated solver in the commercial code Fluent was used.
Differencing schemes utilized were 1^* and 2"'' order Upwind, Power law, and Quick schemes.
The 1^' Upwind scheme was used to start the problem and the higher order schemes to obtain the
final results. The 2"^* order Upwind and Quick schemes were found to give similar results in
terms of mass flow rate and mass imbalance, head-end pressure, chamber pressure drop, and
maximum Mach number at the nozzle exit.
Numerical Convergence (Residuals): Numerical convergence was achieved by satisfying four
requirements in the following sequence. First, the residual error diminished as the number of
iterations was increased. Second, the profiles of the variables ceased to change, at least
qualitatively. Third, a first monitor on the total pressure at the propellant surface until the
average total pressure ceased to change. Fourth, a second monitor on the mass imbalance
11
between the inlet (propellant surface) and the outlet (nozzle exit) mass flow rates until it reached
a small value (10'^ - 10'^ kg/s).
Results
Motor chamber pressure drop, turbulence results, convective heat transfer and acceleration
parameter discussed next.
Motor Chamber Pressure Drop: References 21, 26 and 27 give the chamber pressure as 6.30
MPa (913.85 psia), 6.24 MPa (905 psia), and 6.27 MPa (910 psia), respectively. The" total
pressure used in this study was 6.28 MPa (910.78 psia).
Figure 1 shows the geometry considered and the static pressure distribution in the whole
motor at 1 s post ignition. Figure 2 shows the submerged cavity and its location at prior to
ignition relative to the nozzle. Figure 3 shows the local axial static pressure along the centerline
of the RSRM chamber. The axial coordinate, x, started at the head end and ended at the nozzle
entrance and was measured along the centerline. The interest was to match the calculated motor
chamber pressure drop against measured data from static tests of QM-7 and QM-8 (Ref. 27) and
FSM-9 (Ref. 19). For the cases with the default value for C^ (Model 1 and 2a as given in Table
2), the viscosity ratio (effective (laminar plus turbulent) to laminar) was limited to 10,000 to
match the total pressure of 6.28 MPa (910.78 psia). The chamber pressure drop was calculated to
be 1.38 MPa (200 psia). In Model 2b of Table 2, the viscosity ratio was not limited, but C^ was
reduced from the default value of 0.09 to 0.055." Variable C^ was suggested by modelers and is
well substantiated by experimental evidence." For example, C^ is found to be around 0.09 in the
inertial sublayer of equilibrium boundary layers, and 0.05 in a strong homogeneous shear flow.
This was applied in Model 2b of Table 2. The motor chamber pressure drop was calculated to be
1.23 MPa (178 psia). A better agreement in chamber pressure drop (1.17 MPa (170 psia)) was
12
achieved on a finer grid using Model 2c of Table 2. On the other hand, the calculated values of
y'^ were in the range between 0.5 and 8 and with some scattering. The scattering is because of
the high values of aspect ratio of cells. In addition, it is suggested ' to locate the first grid point
adjacent to the wall at j* = 1. Therefore, the results in this study are reported as calculated on a
coarse grid with y'^ consistent with the standard wall functions. The results of all turbulence
models used will be given when compared. Otherwise, the results of Model 2b of Table 2 will be
given.
Also shown are the measurements of the local static pressure along QM-7 and QM-8 test
motors. In addition, chamber pressure drop was measured to be 1.14 MPa (165 psia) in FSM-9
static test.'' These measurements were made between the head end and in the vicinity of joint 2
in the submerged region. The above agreement is acceptable.
Turbulence Quantities: Two turbulence models were utilized in this study and given in Table
2b. The calculated values of y'^ on the fine grid^ were in the range between 0.5 and 8 and with
some scattering. The scattering is because of the high values of aspect ratio of cells and the
desired wall y'*' were not fully achieved. The size of this motor is too large and axial refinements
were not pursued because they are cost and time prohibitive. The height of the first cell adjacent
to the nozzle wall can be estimated as:
* o ^ ^ 2v . (6)
and shown in Fig. 4. To obtain y'*' < 1 would require grid refinement in the radial direction.
Additional refinement in the axial direcdon is required to keep the cell aspect ratio acceptable
numerically and to avoid cells of negative areas. Thus large number of cells (in millions) are
required for this large motor. It was found that large cell aspect rafio would generate
13
questionable scattering in wall y* and thus heat transfer coefficient. On the other hand,
calculated chamber pressure drop calculated on the finer grid compared more favorably than the
coarse grid. Since the interest here is the calculation of the convective heat transfer, the results
are mainly given on the coarse grid unless stated otherwise.
Figure 5 shows the calculated wall y* along the converging-diverging part of the nozzle. For
the nozzle wall values (j*, h, etc.) the axial coordinate, Xw, started at the beginning of the
convergent section of the nozzle and ended at the nozzle exit and was measured along the nozzle
wall. All the models used give similar wall y'^ profiles and values as shown in Fig. 5. Since the
values were between 25 and 130, the first and second turbulence models used were consistent
with grid resolution.
Calculated flow quality of the present results is shown in terms ofy* vs. u* and y'^ vs. t'^ and
in comparison with the velocity and thermal laws of the wall for incompressible turbulent flow.
Velocity Law of the Wall: Figure 6 shows the present results at the nozzle exit in comparison
with the famous incompressible turbulent velocity law of the wall for an external boundary layer.
28
i.e. the Spalding's profile:
y* =u* + exp(- K B)
I A . ^ [ku) [ku)
Q.X^[ku )-\-KU -^ -
^"^ ' 2 6
(7a)
where k and B were taken^* as 0.4 and 5.5, respectively. The main assumptions made in the
above law were incompressible, negligible stream-wise advection, no axial pressure gradient, and
no transpiration.
The first flow cell from the present study was located at y^ of 29. Very good agreement with
Spalding's profile^^ was achieved around y"" = 25. Furthermore, at y* = 100, the incompressible
velocity law of the wall yielded 16.8 for u* . The present calculations for compressible turbulent
14
flow using Model 1, Model 2a, and Model 2b (Table 2) yielded 16.1, 16.3, and 15.8,
respectively. They are depressed by -4.17%, -2.98%, and -5.95%, respectively.
28
The corresponding compressible turbulent velocity law of the wall is a lot more involved and
is unwarranted. This above comparison is partially in agreement with conclusions made by
White"^ regarding the compressible turbulent velocity law of the wall. His conclusions were:
1) The effect of increasing high Mach number (compressibility, y = 9fuT )/{2Cp Tw)), for a
given TJTaw depress u* below the incompressible turbulent velocity law of the. wall.
This is the case in Model 1 and Model 2a, but partially with Model 2b.
2) On the other hand, a cold wall (heat flux, fi=qw Vw/Ty,kw Wr) tends to raise u* above the
incompressible turbulent velocity law of the wall. This is the case partially in Model 2b.
Temperature Law of the Wall: On the other hand, neither compressible turbulent thermal law
of the wall exist nor conclusions were made by White." To our dismay, again we need to
compare with incompressible turbulent flow. Similarly, Fig. 7 shows the present results at the
nozzle exit in comparison with the temperature law of the wall. The incompressible turbulent
25
temperature law of the wall for an external boundary layer was taken as
r=Pr.y^ ; /<13.2 (7b)
t^= 13.2 Pr +— ^ In
Pr. , r v" ^
13.2
; />13.2
where Pr, and / were taken as 0.85 and 0.41, respectively. Once again, the main assumptions
made in the above law were incompressible, negligible stream-wise advection, no axial pressure
gradient, and no transpiration.
Only a qualitative agreement with the thermal law of the wall is achieved. Similarly, at y"" =
100, the incompressible thermal law of the wall yielded 10.9 for r"". The present calculations
15
using Model 1, Model 2a, and Model 2b (Table 2) yielded 8.5, 8.1, and 6, respectively. They are
depressed by -22.02%, -25.69%, and -^M.95%, respectively. The data shown in Fig. 7 from the
present study is correlated by the present author to give the following:
r = [0.1525 + 1.6181n(/)];25</<1000,C6>D = 0.80 (7c)
from Model 1 , and
f = [- 0.4793 + 1 .709 111(3;" )] ; 25 < y" < 1000, COD = 0.83 (7d)
from Model 2a, and
r = [3.965 + 0.4292 \n[y'^;29 <y< 1000, COD = 0.63 (7e)
from Model 2b of Table 2. The "COD" is the coefficient of determination of a correlation. The
higher COD, the higher the quality of the correlation. These correlations apply to compressible
turbulent flow.
Convective Heat Transfer: The ratio of the convective heat transfer coefficient to the specific
heat at constant pressure (h/Cp) was calculated by eight different methods. The eight methods
will be given, shown and finally discussed. The first two methods used finite volume CFD
(Methods 1 and 2), followed by four approximate methods (Methods 3, 4, 5 and 6), followed by
an integral method (Method 7) and then using backed-out data from measurements (Method 8).
In Method 1, it was calculated internally using the calculated heat flux based on the difference
between the local normalized specified wall temperature given by Eq. (4) and shown in Fig. 8
and the chamber temperature used as a reference temperature, i.e.,
(8a)
c^ ~cjr-r(x)]
and is shown in Fig. 9.
16
In Method 2, it was calculated using the recovery temperature and defined as
(9a)
h^jx) _ q^, {x)
and shown in Fig. 9. The recovery temperature" is given as
r^ , (X) = r (X) + 9t, (X) [r - r (x)] (9b)
and where 9i i{x) = [Pr,(x)]"^ was used for turbulent flow and was calculated locally. The
temperature at the edge of the boundary conditions has been replaced by the local axial static
temperature along the motor centerline and shown in Fig. 8. The recovery temperature (Taw. i)
was calculated at the motor centerline and shown in Fig. 8 to be less than the chamber
temperature. The chamber (7^) and recovery temperature {Taw. 2) were matched only by taking Pr
as 1 which yields 1 for 9^2 which corresponds to an ideal situation. Also shown the gas
temperature {T g, w ) in the cell adjacent and along the nozzle wall. Good agreement between
Methods 1 and 2 is shown.
The four approximate methods (Methods 3, 4, 5 and 6) are discussed next. The third method
used turbulent flow over a flat plate which is
^^^^ = n moA fd^^^mOs
c,
= 0.0296 [Re(x)]"'' Pr(x)
1/3
k(x)
(10)
V 'J
and is not shown in Fig. 9, because it overestimates the heat transfer and does not have the
correct profile. It appears approximating this large nozzle as a flat plate is not plausible. The
maximum occurred at the start of the nozzle and decreased along the wall as expected for flow
over a flat plate.
The fourth method used the modified Reynolds ' analogy for laminar flow over a flat plate
which is,
17
!}M^Lc,ix)mx)Pr(xf'
[kix)'
X
f ' 1
(11a)
where the local skin friction is calculated as
Cf(x) =
K(x)
M^Miv^i^-
(11b)
and is not shown in Fig. 9. Again, this method overestimates the heat transfer but has the correct
profile. In addition, Methods 3 and 4 are not shown to reduce clutter.
17
The fifth method used the Dittus-Boelter correlation for fully-developed turbulent pipe flow
as
h,{x)
= 0.023 [Re/x)rPr(x)''
k{x)
d{x)
'2'
v^^;
(12)
and is shown in Fig. 9. The exponent n for Pr was taken as 0.3 for cooling {Tw < TJ) and 0.4 for
heating {Tw > TJ) per Fig. 8.
10
The sixth method used the Bartz correlation
Kix)
0.026
[dT
^^0. c^ ^ ,
Pr
0,6
/)
p.sTi^i'T
V ^ J
k'^ J
/ ^. ^o^
A{x)
(7{x)
kix)
dix)
( , ^
A^";
(13a)
and shown in Fig. 9. It was based on a similar correlation of fully developed turbulent pipe flow
heat transfer, as given previously. The terms in the bracket are constant for a single nozzle. The
subscript "0", on the second term in the bracket, signifies properties are to be evaluated at
chamber conditions. Prandtl number at chamber conditions was evaluated using
Ay
Pr =
(13b)
9^-5
and was calculated to be 0.843. The ratio of specific heats is calculated from Cp and the gas
constant. Using lower values for Pr (= 0.5) as in this study, would overestimate h/Cp as
calculated by Eq. (13a). The characteristic exhaust velocity was calculated as
irRTT
(13c)
c =
f
r
Y + l
V(r-i)
The first term outside the bracket in Eq. (13a) is a function of nozzle local area (A{x) = {n/4)
d^(x)). The second term outside the bracket, aix), is a dimensionless factor accounting for
variations of density and viscosity across the boundary layer. It is shown in Bartz to be
1 (13d)
a(x) = -
~2~ '■"
2 T
5 f
J
l + ^^^.v.'W
\/5
V
where co is taken to be 0.6. This dimensionless parameter was calculated and found to vary
between 1.03 and 1.16.
29
The seventh method used the Turbulent Boundary Layer (TBL) code and is shown Fig. 9.
TBL method solves, simultaneously the integral momentum and energy equations for thin axi-
symmetric boundary layers.
The eighth method used backed-out {hJCj) data" that is used in CMA code to match the
measured nozzle erosion profile.
The following are observations and comparisons between the CFD results (Methods 1 and 2)
and the other methods (Methods 3 to 8) as:
(1) They all have similar profiles. They all increase to a maximum in the vicinity of the
throat and then decrease along the nozzle wall with the exception of the correlation for the
turbulent flow over a flat plate (Method 3). Method 3 shows the maximum to occur at the start of
the nozzle (analogous to the leading edge of a flat plate) and decreases along the nozzle wall as
expected. It was found that, along the last 40% of the nozzle length, the flat plate correlation
19
becomes reasonable. This is because curvature effects become less important along the nozzle.
(2) Methods 1 and 2 (CFD), 7 (TBL), and 8 (backed-out) show the maximum to occur
upstream of the nozzle throat. On the other hand, Methods 5 {Dittus-Boelter) and 6 (Bartz) show
the maximum to occur at the nozzle throat. The latter group of methods is heavily dependent on
the nozzle area ratio. For completeness. Method 4 (the modified Reynolds' analogy) also
predicts the maximum to occur at the throat.
(3) The profiles obtained from all Methods (with the exception of Method 3 (flat plate)i were
analogous to the y* distribution shown in Fig. 5. This confirms the dependency of heat transfer
on the turbulence models as well as the value of >'"^.
(4) At the nozzle throat: From the data used Fig. 9, the following values of {h/Cp) are taken
as 4, 4.32, 6, 12.09, 4.12, 3.89, 4.29, and 3.18 kg/m'-s for Methods 1 through 8, respectively.
The values from Methods 1 through 8 differ Method 1 by 0.0%, 8.0%, 50.0%, 202.3%, 3.1%, -
2.8%, 7.3% and -20.5%, respectively. This comparison is shown in Fig. 10. Results of Method
4 fall out of the chart. Similarly, the values from Methods 1 through 7 differ from the backed-out
data"^ (Method 8) by 25.8%, 35.8%, 88.7%, 280.2%, 29.6%, 22.3%, 34.9%, respectively. This
comparison is also shown in Fig. 10. Results of Method 3 fall out of the chart.
The comparison between Method 1 and Method 8 (backed-out) is within 26% and may seem to
be in disagreement. The present results (Methods 1 and 2) are restricted to convective heat
transfer. On the other hand. Method 8 (backed-out) was devised to match nozzle erosion that
was a result of heat transfer (convective, radiative and conductive), chemical, and mechanical
(particle impingement) effects combined. Therefore, comparing the present results and Method 8
in a relative manner would be a sound approach. In the converging section of the nozzle
(excluding the nose-tip), all methods, over-predict the ratio /z/Q. The heat transfer in the
20
convergence section may be inhibited by chemical effects. The nozzle wall is an ablating
(transpiring) surface and was not modeled in this study. The boundary layer adjacent to the
nozzle wall would be oxidizer-poor that inhibit chemical reactions and heat transfer. At the
nozzle throat, Method 1 over-predicts Method 8 by 26%. It is to be noted that 12 data points
used along the nozzle wall in Method 8. While, in the present study, there are 320 data used
along the nozzle wall. All the methods over-predict the heat transfer has great implications. One
concludes some factors that inhibit heat transfer. Some factors such as a thick boundary, layer
and a transpiring surface would inhibit heat transfer. Furthermore, the condensed phase may
flow as a film that would insulate the nozzle and inhibit heat transfer. In the diverging section of
the nozzle, less agreement is obtained with Method 8 in the last 40% of the nozzle length. This
disagreement may be attributed to eroded exit cone because of particle impingement based on
observations made in post static and flight tests of this motor." Based on this author's
experience and observations, some conclusions have to be made through heuristic arguments
alone.
(5) hi the vicinity downstream of the nozzle throat, the present results have a sudden
decrease-increase-decrease. This sudden drop was attributed to the large drop in the specified
surface temperature (Fig. 8). In the vicinity of the throat, the surface temperature dropped by 590
K (1062 R) within 0.5 m (1.64 ft). This was easily verified by specifying uniform surface
temperature at 2000 K (3600 R) and the sudden drop was not calculated. The heat transfer
decreases when turbulent flow re-lamin arize. Furthermore, this sudden drop was detected
experimentally in the Jet Propulsion Laboratory (JPL) nozzle. ' An isothermal wall
boundary condition was used and was operated at a maximum chamber pressure of 1.72 MPa
(250 psia). The JPL nozzle has a throat radius of about 2 cm (0.8 in) that amounts to about 3
21
percent of the throat size in this study. This is the reason for considering the possibility of
turbulent boundary layer re-laminarization discussed in terms of acceleration parameter in the
next section.
Acceleration Parameter: The acceleration parameter, Ki, was calculated from" based on
external boundary layer as
U fdu ^
K =
ydx J
(14a)
> 3x10"'
pSx) u^ {x)
and is shown in Fig. 1 1 . The calculated values of Ki were smaller than the transition value of
3x10'^. Therefore, re-laminarization of the turbulent boundary layer did not occur. The above
30
acceleration parameter was translated by Coon and Perkins into terms more pertment to tube
flow as
K =
4 q. (X)
fn(r^^ . (14b)
Mix)
C
> 1.5j10
-6
P. (x) M„ (x) d(x) T
The maximum calculated values of Ki and K2 were smaller than the transition values given above
as 3x10'^ and 1.5x10'^, respectively. Therefore, re-laminarization of the turbulent boundary layer
did not occur.
Summary and Conclusions
The following summary and conclusions have been reached:
• Two turbulence models and two types of wall treatment were used in this study. Based on
economy and execution time, RNG k-£ with standard wall functions gives reasonable results in
relative comparison with the backed-out data."
• Convective heat transfer coefficients have been calculated using two methods. They have
been compared with approximate methods, Bartz correlation, turbulent boundary layer (TBL)
22
theory code, and backed-out data based on post static and flight tests.
• The interdependency of the convective heat transfer and wall >>"*" has been shown.
Therefore, consistency between a grid and a turbulence model cannot be over-emphasized.
• The accuracy and the selection of the schemes and results are based on matching the RSRM
ballistic predictions of mass flow rate, maximum head end pressure, and vacuum thrust and
specific impulse and measured chamber pressure drop and was found to be good.
• Convective heat transfer coefficients were calculated and compared favorably with ejusting
theory. Qualitative comparison with backed-out data of the ratio of the convective heat transfer
coefficient to the specific heat at constant pressure was made in a relative manner. Backed-out
data was devised to match nozzle erosion that was a result of heat transfer (convective, radiative
and conductive), chemical, and mechanical (shear particle impingement forces) effects
combined.
References
^FLUENT 5 Solver Training Notes, Huent Inc., TRN-1998-006, Dec. 1998.
-FLUENT 5 User's Guide, Fluent Inc., July 1998.
^Golafshani, M. and Loh, H.T., "Computation of Two-Phase Viscous Flow in Solid Rocket
Motors Using a Flux-Split Eulerian-Lagrangian Technique," AIAA Paper 89-2785,
AIAA/ASME/SAE/ASEE 25th Joint Propulsion Conference, Monerey, CA, July 10-12, 1989.
Yoh, H.T. and Chwalowski, P., "One and Two-Phase Converging-Diverging Nozzle Flow
Study," AIAA Paper 95-0084, AIAA/SAE/ASME/ASEE 33rd Aerospace Sciences Meeting,
Reno, NV, Jan. 9-12, 1995.
^Whitesides, R.H., Dill, R.A. and Purinton, D.C., "Application of Two-Phase CFD Analysis to
the Evaluation of Asbestos-Free Insulation in the RSRM," AIAA-97-2861,
23
AIAA/ASME/SAE/ASEE 33rd Joint Propulsion Conference, Seattle, WA, July 6-9, 1997.
^Laubacher, B.A., "Internal Flow Analysis of Large L/D Solid Rocket Motors," AIAA Paper
2000-3803, AIAA/ASME/SAE/ASEE 36th Joint Propulsion Conference, Huntsville, AL, July
16-19, 2000.
^Ahmad, R.A., "Internal Flow Simulation of Enhanced Performance Solid Rocket Booster for
the Space Transportation System," AIAA-2001-3585, AIAA/ASME/SAE/ASEE 37th Joint
Propulsion Conference, Salt Lake City, UT, July 8-11, 2001.
^McBride, B.J., Reno, M.A., and Gordon, S., "CMA- Charring Material Thermal Response and
Ablation Program, SEEOl, Version 1.6," Computer Software Management and Information
Center (COSMIC), The University of Georgia Office of Computing and Information Services,
Athens, Georgia, May 1994.
^Ahmad, R.A., "Discharge Coefficients and Heat Transfer for Axisymmetric Supersonic
Nozzles," Heat Transfer Engineering, Quarterly International, Vol. 22, No. 6, Dec. 2001.
'°Bartz, D.R., "A Simple Equation for Rapid Estimation of Rocket Nozzle Convective Heat
Transfer Coefficients," Jet Propulsion, Vol. 27, No. 1, pp. 49-51, Dec. 1957.
'^Back, L.H., Massier, P.F., and Gier, H.L., "Convective Heat Transfer in a Convergent-
Divergent Nozzle," Intern. Journal. Heat Mass Transfer, Vol. 7, pp. 549-568, May 1964.
^^Moretti, P.M. and Kays, W.M., "Heat Transfer to a Turbulent Boundary Layer With Varying
Free-Stream Velocity and Varying Surface Temperature - An Experimental Study," Intern.
Journal. Heat Mass Transfer, Vol. 8, pp. 1187-1202, 1965.
'^Back, L.H., Massier, P.F., and Cuffel, R.F., "Some Observations on Reduction of Turbulent
Boundary-Layer Heat Transfer in Nozzles," A/AA J., Vol. 4, No. 12, Dec. 1966 pp. 2226-2229.
''^ack, L.H., Massier. P.F., and Cuffel, R.F., "Flow Phenomena and Convective Heat Transfer
24
in a Conical Supersonic Nozzle," Journal of Spacecraft and Rockets, Vol. 4, No. 8, Aug. 1967,
pp. 1040-1047.
'^Back, L., H., Cuffel, R.F., and Massier, P.F., "Laminarization of a Turbulent Boundary Layer
in Nozzle Flow - Boundary Layer and Heat Transfer Measurements with Wall Cooling," Journal
of Heat Transfer, Vol. 92, Series C, No. 3, Aug. 1970, pp. 333-344.
^^ang, Q., "On the Prediction of Convective Heat Transfer Coefficients Using General-
Purpose CFD Codes," AL\A Paper 2001-0361, AL\A 39* Aerospace Sciences Meeting, Reno,
NV, Jan. 8-11,2001.
'^Incropera, F.P. and Dewitt, D.P., Introduction to Heat Transfer, Wiley, New York, 1985.
'^Wiesenberg, B., "Design Data Book for Space Shuttle Reusable Solid Rocket Motor,"
Thiokol Corp., TWR-16881, Rev. A., July 1997.
'^Manuel, L.J., "Final Report for CTP-0179 - Static Test Plan for FSM-9," Thiokol Propulsion,
TWR-76760, Sept. 2001.
^°McBride, B.J., Reno, M.A., and Gordon, S., "CET93 and CETPC: An Interim Updated
Version of the NASA Lewis Computer Program for Calculating Complex Chemical Equilibria
With Applications," NASA Technical Memorandum 4557, 1994.
^'Nickerson, G.R., Coats, D.E., Dang, A.L., Dunn, S.S., Berker, D.R., Hermsen, R.L., and
L^amberty, J.T., "The Solid Propellant Rocket Motor Performance Prediction Computer Program
(SPP)," Version 6.0, Vol. 1-4, AFAL-TR-87-078, Dec. 1987.
"GRIDGEN User Manual Version 13" by Pointwise, Inc. 1998.
-■^Maw, J.F., "AeroAThermal Analysis of the RSRM Nozzle," TWR-17219, Rev. D, Thiokol
Corp., Feb. 20, 1997, p. 22.
''^TableCurve® 2D V 5.0/3D V 3.0 automated surface fitting software, Jandel Scientific, San
25
Rafael, CA, 1993.
"^Kays, W.M. and Crawford, M.E., Convective Heat and Mass Transfer, 3' ed., McGraw-Hill,
New York, Chaps. 11 and 13, 1993.
^^A user's Guide for Computer program No. SCB02 - An Internal Ballistics Program for
Segmented Solid Propellant Rocket Motors," Thiokol Propulsion, 1982.
"Gruet, L., "QM-7 and QM-8 Transient Axial Pressure Calculated from Strain Gauge Data,"
TWR-63695, Thiokol Corporation, 30 March 1992.
^Vhite, P.M., Viscous Fluid Flow, 2"'^ ed., McGraw-Hill, New York, Chaps. 6 and 7, 1991.
^^Elliot, D.G., Bartz, D.R., and Silver, S., "Calculation of Turbulent Boundary Layer Growth
and Heat Transfer in Axi-symmetric Nozzles," National Aeronautics and Space Administration
(NASA), Jet Propulsion Laboratory (JPL) Technical Report No. 32-387, Feb., 1963.
^°Coon, C.W. and Perkins, H.C., "Transition From the Turbulent to the Laminar Regime for
Internal Convective Flow with Large Property Variation," Journal of Heat Transfer, Vol. 92,
Series C, No. 3, Aug. 1970, pp. 506-512.
26
Table 1: Summary of geometry,
operating ballistic conditions, ballistic
predictions, and gas thermophysical
properties.
Geometry
rn/r*
1.50
rjrl
2.63
AdAL
2.065
A./4'
7.72
Ballistic Operating Conditions (TWR-
16881. Rev. A fRef 18))
P n ("M Pa. psia)
6.28 (910.78)
Chamber Pressure Drop (rSM-9 Static
Test (Ref. 19))
AP, (MPa. psia)
1.138.165
Gas Thermophysical Properties
r»(K.R)
3419.2.6164.56
(lo (kg/m-s, Ibm/ft-s)
9.25x10"-
6.22x10-
Q g/kg-K. Btu/lbm-R)
1966.54.0.47
kn rW/m-K. Btu/hr-ft-R)
0.397. 0.229
MWCkg/kgmole)
28.46.
// = /^ (T/To/ , k = ko (T/Tof where ^ =
0.66687 as calculated from NASA-Lewis
Program (Ref. 20).
27
Table 2 CFD calculated pertinent results for the RSRM at 1 sec burn time
Parameters
Coarse Grid (105,850 cells)
Fine Grid
(376,100
cells)
Model 1
Model 2a
Model 2b
Model 2c
Viscous Model
Standard k-p.
RNG' K-F
RNG'k-f
RNG K-F
Near Wall Treatment
Standard
Wall
Functions
Standard
Wall
Functions
Standard
Wall
Functions
Two-
Layer
Zonal
Model
Desired wall y'^
20 < j^ <
100
20 < y^ <
100
20<y^<
100
/<1
C,,
0.09
0.09
0.055
0.055-0.06
111/ lii
10^
10^
10^
10%10'
P s. max. CFD (M Pa, psia)
6.29,911.90
6.30,913.62
6.34,
919.14
6.28,
911.16
P 0. max. CFD (M Pa, psia)
6.29,911.99
6.30,913.71
6.45,
934.92
6.28,
911.26
% difference in Pn (Using Table 1)
0.13%
0.32%
2.65%
0.05%
M tt. max
3.55
3.59
3.58
3.61
% difference in mass flow rate
-1.01%
-0.92%
-0.71%
-0.11%
AP fM Pa. psia)
1.38.200
1.38,200
1.23. 178
1.17. 170
% difference in AP (Using Table 1)
21.21%
21.21%
7.88%
3.03%
Iterations used
2000
2000
4000
5000
^•iWG (renormalization group theory), % difference = (calculated / predicted - 1) x 100
28
2.0
1.5
1.0
0.5
Propellant
0.0
Nozzle
67.0 67.5 68.0 68.5 69.0 69.5 70.0
xir'
Fig. 2 RSRM cavity at prior to ignition time.
30
I
6.34«-f06
S.Q&»*Q6
3.S2«-»C6
2.56fl+06
i/aoa-toe
3.638+04
Max P
Min P
Enlarged
Contours of Static Prsssurs (pascal]
Mar OS, 2002
FLUENT 6.0 (axi, segregatfld, mgkfl)
Fig. 1 Geometry and pressure distribution in the RSRM at 1 s post ignition time, 0, < Pj
(M Pa, psia) < 6.34, 919.54; number of levels = 10; 0, < ^. (M Pa, psia) < 0.634, 91.95.
29
- - - -Model 2b (Table 2) — - -Model 2c (Table 2)
■ QM-7 Test Motor [27] O QM-8 Test Motor [27]
▲ FSM-9[19]
1 nc
■ ■
■
1 -.
■ .
A-
0.95
?; 0.9-
/sn\ ■
o:
o ^ ■
0.85 -
o . V\ ■ ■
3k ^
0.8
^\ — \
\o o;
^**^»
I
15 30 45 60 75
xir'
Fig. 3 Local axial static pressure along the
centerline of the RSRM at 1 s post ignition.
31
Fig. 4 Estimated height of cell adjacent to the
wall
32
•r/r*. Nozzle radius —
Model 2a (Table 2) —
— Model 1 (Table 2)
- -Model 2b (Table 2)
Fig. 5 Local wall y'*' along the converging-
diverging section of the RSRM nozzle at 1 s
post ignition.
33
25
20
15
10
■ Model 1 (Table 2)
Model 2a (TaWe2)
- - - 'Model 2b (Table 2)
— - -VelocityLawofthe Wall [28]
*
^:^^ —
^^'b^
.y
I
;— - •- ^ y* = 26
10
100
1000
Fig. 6 Comparison of the velocity law of the
wall with the RSRM nozzle calculations at
exit.
34
15
10
■Model 1 (Table 2)
- Model 2a (Table2)
- -Model 2b (Table 2)
- -Thermal Law of the Wall [25]
10 100
1000
Fig, 7 Comparison of the temperature law of
the wall with the RSRM nozzle calculations
at exit.
35
2£
15
0£
Taw,2 fT„
"— ' ' — Taw, I /To
r/r
-T-ir/'T-
w ' ■■ ow
* e,w ' *
Tw/To
Tcl/To
-J — I — 1 — I — ^-
1.00
0.91
0.82
0.73
0.64
0.55
-^ 1 1 : ^ : 1 1 : ^ f : : : 1 , : ^ ^ : : , 1 L_
67 68 69 70 71 72 73 74 75 76 77
x/r*
0.46
Fig. 8 Local specified surface temperature and calculated temperatures in the RSRM at 1 s
post ignition.
36
^— — r / r*, Nozzle Radius
- - -CFD, Method 2
— Bartz [10] - Method 6
O Backed-out [23] - Method 8
CFD, Method 1
— - -Dittus-Boelter[17]- Methods
■ TBL Code [29] - Method 7
M Throat Location
Fig. 9 Local convective heat transfer on the nozzle wall of the RSRM at 1 s post ignition
(Methods 3 and 4 are not shown).
37
O Difference from Method 1 D Difference from Method 8
♦ Method 1 ■ Method 8
50
40
I 20^
10
-^
TT
Q 0^
^' -10^
-20
-30
^
-s-
D
n-
-<x
XT
xr
^-
— — ! 1 1 ; : ! '•
123456789
Method
Fig. 10 Comparison of heat transfer at the
throat among all methods
38
•r/r* , Nozzle Radius
Based on External Boundary Layer
>■ - - - Based on Internal Boundary Layer
3.6E-08
2.4E-08 .S ^
1.2E-08 1 I
< £
l.OE-13
Fig. 11 An acceleration Parameter (K) based
on external and internal boundary layers
along the nozzle centerline of the RSRM at 1
s post ignition.
39