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