_1-12_ Karima E. Amori 1 Al-khwarizmi Engineering Journal Al-Khwarizmi Engineering Journal, Vol.3, No3,pp 1-12, (2007) Parametric study of thermal behavior of thrust chamber cooling channels Dr. Karima E. Amori Mech. Eng. Dep., University of Baghdad, Iraq (Received 6 March 2007 ; accepted 10 June 2007) Abstract A numerical investigation is adopted for two dimensional thermal analysis of rocket thrust chamber wall (RL10), employing finite difference model with iterative scheme (implemented under relaxation factor of 0.9 for convergence) to compute temperature distribution within thrust chamber wall (which is composed of Nickel and Copper layers). The analysis is conducted for different boundary conditions: only convection boundary conditions then combined radiation, convection boundary conditions also for different aspect ratio (AR) of cooling channel. The results show that Utilizing cooling channels of high aspect ratio leads to decrease in temperature variation across thrust chamber wall, while no effects on heat transferred to the coolant is indicated. The radiation has a considerable effect on the computed wall temperature values. ٌِKeywords: Thrust chamber, multi layer wall, cooling duct, CFD. Introduction The evaluating of convective heat transfer rates from high temperature combustion gases to the converging diverging nozzle of liquid fueled rocket engine is crucial for optimizing the design of thermal control systems and evaluation of thermomechanical behavior and fatigue of chamber liner to insure adequate long life operation. The liner of most current thrust chamber is fabricated from a high conductivity material, such as Copper or Copper based alloy, and closed out with Nickel, incorporating cooling channels with liquid Hydrogen as the coolant fluid as shown in fig. (1). Armstrong and Brogren conducted an experimental investigation of fatigue life of thrust chamber under cyclic thermal loading , Quentmyer identified a correlation between number of cycles to failure and temperature predicted through the application of SINDA ( System Improved Numerical Differencing Analyzer) software. Saha extended the computer program (BLIMP) for boundary layer analysis to predict a thrust chamber wall temperature by coupling this analysis with regenerative process along axial nozzle cross section passing through cooling channel . The objective of the present work is to compute the temperature distribution within material of temperature dependent conductivity by utilizing finite difference method employing different space divisions in r and θ direction with special application to thrust chamber wall of rocket thrust chamber RL10 composed of multi layers and exposed to different outside and inside boundary conditions. Dr. Kareema E. Amori /Al-khwarizmi Engineering Journal ,Vol.3, No.3, PP 1-12(2007) 2 c Copper Nickel Soot T T oT o g o 1o Cell of thrust chamber wall cooling channel Hot gas (a) (b) (c) Fig.(1). a).Rocket Thrust Chamber Nozzle. b )Cross Section Through AA. c). Layers of the Wall Theory and computational scheme The temperature distribution within thrust chamber wall can be obtained by satisfying the first law of thermodynamic combined with Fourier's law of heat conduction with the assumptions: 1-Steady state two dimensional heat transfer 2-Temperature dependent conductivity 3- No internal heat generation yields the following governing differential equation: 0)()( = ∂ ∂ ∂ ∂ + ∂ ∂ ∂ ∂ θθ T k r T r r …(1) where T : temperature (K) K : thermal conductivity (W/m.K) The interested region is divided into unequal space divisions (Δr and Δθ) in r and θ directions respectively as shown in fig.(2). A A toward i njector hot gas in coolant in hot gas out coolant out Dr. Kareema E. Amori /Al-khwarizmi Engineering Journal ,Vol.3, No.3, PP 1-12(2007) 3 R i R ci R co R o i , j r r 2 1r 2 1 R4 R4 R2 R2 R1R1 R3 R3 i , j-1 i , j+1 i -1, j i+1 , j i , j Fig.(2): Discritization of the Physical Domain in r-θ Plane Utilizing the finite difference method the temperature at each node can be computed in terms of temperature of the surrounding nodes. Hence equation (1) is approximated by this mean for a middle node (i,j) shown in fig.(3) such that [Ozisik1980]:- Fig.(3). Thermal Resistances Around Node i,j 01,,11,,1 =+++ −−++ jijijiji qqqq ….(2) where        − = − = − = − = − − − − + + + + 2 ,1, 1, ,,1 ,1 2 ,1, 1, 1 ,,1 ,1 ; 3 ; R TT q R TT q R TT q R TT q jiji ji jiji ji jiji ji jiji ji .(3) where (i,j) : indices in r and θ direction q : heat transferred in W/m R : thermal resistance (for conduction heat transfer )(*tan depthunitperareakcedisR = ) or ))(*(1 depthunitperareahR = (h is convection heat transfer coefficient hcon, or for radiation heat transfer )( 3, 22 , 3 , TTTTTTh jijijirad +++= εσ ) ε : emissivity factor = 0.9 σ :Stefan-Boltzman constant =5.669*10-8 W/m2 K4 Hence equation (2) can be rewritten such that: Dr. Kareema E. Amori /Al-khwarizmi Engineering Journal ,Vol.3, No.3, PP 1-12(2007) 4 4321 4 1, 3 ,1 2 1, 1 ,1 , 1111 RRRR R T R T R T R T T jijijiji ji +++ ++ = −−++ .……(4) where ji ji krr r krr r R RR ,21 1 ,121 1// 1 / 11 *)22( 2 *)22( 2 ∆+∆ ∆ + ∆+∆ ∆ =+= + θ θ ) 11 ( )( ,,121 1 1 jiji kkrr r R + ∆+∆ ∆ = + θ ………(5a) )2)(22( 2 )2)(22( 2 221, 2 2211, 2 // 2 / 22 rrk r rrk r R ji ji RR ∆+∆+∆ ∆ + ∆+∆+∆ ∆ = += + θθ θθ ) 11 ( ))(2/( ,,1212 2 2 jiji kkrr r R + ∆+∆∆+ ∆ = +θθ ..(5b) )22( 2 )22( 2 21, 2 21,1 2// 3 / 33 rrk r rrk r R ji ji RR ∆+∆ ∆ + ∆+∆ ∆ =+= − θ θ ) 11 ( )( ,,121 2 3 jiji kkrr r R + ∆+∆ ∆ = − θ ……….(5c) )2)(22( 2 )2)(22( 2 121, 1 1211, 1 // 4 / 44 rrk r rrk r R ji ji RR ∆−∆+∆ ∆ + ∆−∆+∆ ∆ = += − θθ θθ ) 11 ( ))(2/( ,1,211 1 4 jiji kkrr r R + ∆+∆∆− ∆ = −θθ …….(5d) where ( r ) is radius of node (i,j), and since the conductivity depends on temperature then: ( 1,,11,,1, −−++ ≠≠≠≠ jijijijiji kkkkk ). The expressions of thermal resistances for other special kind nodes are presented in appendix A: Geometry and boundary conditions Detailed dimensions of thrust chamber cross section employed in the present work. Is given in table (1) based on [Naraghi, et. al. 1987]. The left and right sides of thrust cell are assumed to be adiabatic walls, since symmetric temperature fields cross these sides (lines of symmetry). In the present study two sets of boundary conditions are analyzed: 1- A convective boundary condition applied at inner radius of thrust chamber (which exposed to hot gas), and cooling channel walls ( where heat is removed from chamber wall to liquid hydrogen). Free convection boundary condition at chamber outer radius which is exposed to outside environment where heat transfer coefficient at this surface is given by [Holman 1981]: ( ) 25.03 2 2 1, Pr)( 53.0 θ µ ρβ θ ∆ − ∆ = ∞ o ji o air conv R TTg R k h …(6) where g: gravitational acceleration m/s2 β :thermal expansion=1/T∞ (K-1) ρ :air density= 1.2 kg/m3 Pr: Prandtle No. for air=0.7 2-Radiation boundary conditions at these surfaces are considered beside the convective boundary conditions mentioned above, heat transmitted in this case is calculated as:    += −= ∞ radconv jiji hhh TThAq )( ,, …...(7) Dr. Kareema E. Amori /Al-khwarizmi Engineering Journal ,Vol.3, No.3, PP 1-12(2007) 5 The specifications of fluids surrounding the thrust cell with convective heat transfer coefficients employed are presented in table (1). Computational algorithm 1-Assume initial nodal’s temperature such that: 3/)( 1, ∞∞ ++= TTTT cgji …….(8) 2- Compute thermal conductivity of each node based on its temperature by employing linear interpolation to the data given in fig.(4). 3-Compute nodal’s temperature value which is given by       = ∑∑ == 4 1 4 1 , /1/)( k kk k kji n RRTT ….(9) for nodes located at the insulated sides (k=1→ 3). Iterative successive method with under relaxation factor (ω=0.9) has been used 4-The iterations are terminated according to the following convergence criterion. 7 ,,, 10/)( −≤− jinjiji n TTT …...(10) where n refers to new temperature value. 4-Calculate heat transferred from hot gas, to coolant, and to outside of thrust wall cell applying: )( 1 kg k kig TTAhq −= ∞ = ∑ )( 1 ck k kcc TTAhq −= ∑ = …(11) )( 1 1 ∞ = −= ∑ TTAhq k k koo Results and discussion Figure (5) presents radial temperature distribution within thrust wall cell (for cooling channel aspect ratio (AR=8.25) (which is the ratio of cooling channel height to width) when only convection boundary condition is considered ) for different selected lines (from A-A to H-H ). Highest temperature values indicated at lines located far from cooling channel, while the coldest lines are those pass through it due to cooling process. Also it is clear that temperature values are decreasing toward outside wall surface lines (AA to DD), this is agreed well with the results of (Naraghi 2001). Parametric study 1-Effect of cooling channel aspect ratio Figures (6) and (7) show temperature profile at lines AA and HH respectively for different aspect ratios. It is clearly indicated that the nodes at lines AA and HH (belong to region below cooling channel) have the highest temperature values for (AR=2.5) compared with other (AR) values while the lowest temperature is indicated for (AR=8.25) for the same nodes. This is due to increas in the cooling channel width for the same height for (AR=2.5) which leads to increase in the surface area of heat transfer to coolant. This results is agreed with [Naraghi 2001] for utilizing high aspect ratio cooling channel design. Figure (8) shows clearly that temperature profile at cooling channel wall for (AR=8.25) is preferred since no large variations exist where the subscript k refers to element indices that belong to either high temperature wall or cooling channel wall or thrust cell outer wall. 5- Verify that qg+ qc+ qo=0.0 2-Effect of radiation boundary condition In order to show the effects of radiation heat transfer on the wall temperature of cooled rocket thrust chamber, Conjugated Boundary Conditions (radiative and convective) Dr. Kareema E. Amori /Al-khwarizmi Engineering Journal ,Vol.3, No.3, PP 1-12(2007) 6 have been applied at inner surface, cooling channel wall and outside surface. Figure (9) shows the Temperature profiles for selected lines AA and EE, which indicate higher temperature values when radiation is considered. Figure (10) shows lower temperature values at coolant channel wall when combined radiation convection heat transfer is applied. It is deduced form fig.(11) that larger heat is transferred from thrust chamber wall to coolant when radiation heat transfer is considered while no effect of aspect ratio is indicated . Conclusions The effects of cooling channel aspect ratio and radiation boundary condition on thermal characteristics of thrust chamber wall are studied. It was found that Utilizing cooling channels of high aspect ratio leads to decrease in the temperature variation across thrust chamber wall, while no effects on heat transferred to the coolant. The radiation has a considerable effect on the wall temperature value. REFERENCES: Armstrong, W.H. and Brognen, E.W., (1975), " Thrust Chamber Life Prediction Volume II Plug Nozzle Centerbody and Cylinder Life Analysis", NASA CR 143822 Arya, V.K. and Arnold, S. M. , (1991), "Viscoplastic Analysis of an Experimental Cylindrical Thrust Chamber Liner", NASA TM 103287 Butler D. T. and Pindera, M-J. (2004) “Analysis of Factors Affecting the Performance of RLV Thrust Cell Liners", NASA/CR-2004-213141 Holman, J.P. ;(1981) ;”Heat Transfer”, Fifth Edition, McGraw-Hill Book Company Kacynski, K. J. (1992),"Thermal Stratification Potential in Rocket Engine Coolant Channels", NASA TM 4378 Naraghi, M. H.(1987),”A Two- Dimensional Thermal Analysis of Rocket Thrust Chambers” NASA-TM- 100191 Naraghi, M. H.;Quentmyer, R. J.;Mohr, D. H.(2001),”Effect of a Blocked Channel on the Wall Temperature of a Regeneratively Cooled Rocket Thrust Chamber” AIAA 2001-3406 Ozisik, H. H. (1980), " Heat Conduction" John Wiley & Sons,1980 Saha, H. (1976), " Rocket Engine Thrust Chamber Wall Temperature Distribution Calculation and Analysis", NASA NSG-8022 Dr. Kareema E. Amori /Al-khwarizmi Engineering Journal ,Vol.3, No.3, PP 1-12(2007) 7 A B C D E F G H HF E DCBA G Table (1):Data employed in computer program Run [Naraghi 1987] 200 400 600 800 1000 1200 1400 1600 1800 Temperature (K) 0 10 20 30 40 50 60 70 80 90 100 110 120 N ic ke l K (W /m .K ) 0 50 100 150 200 250 300 350 400 450 C up pe r k (W /m .K ) Nickel Cupper Fig. (4): Thermal Conductivity of Cupper and Nickel [Kacynski] 0.88 0.90 0.92 0.94 0.96 0.98 1.00 r/Ro 50 52 54 56 58 60 62 64 66 68 70 72 74 76 Te m pe ra tu re (K ) Line A-A Line A-A Line B-B Line C-C Line D-D Line E-E Line G-G Line h-h (a) (b) Fig.(5):a). Selected Lines Through Thrust Cell. b). Radial Temperature Distribution of Thrust Chamber Cell For (AR=8.25 ) Hot gas temperature (2655 K ) Hot gas heat transfer coefficient (W/m2K ) (290) Coolant temperature (50 K) Coolant heat transfer coefficient (W/m2K ) (15290.65) Outside temperature (294 K) No. of cooling channels (100) thrust chamber inner radius Ri (44.44 mm ) thrust chamber outer radius Ro (49.675 mm) cooling channel inner radius Rci (44.95 mm) cooling channel height H (4.2144 mm cooling channel width based on Rci (0.511 mm) Dr. Kareema E. Amori /Al-khwarizmi Engineering Journal ,Vol.3, No.3, PP 1-12(2007) 8 Fig.(6): Aspect Ratio Effect on Temperature Distribution Along Line A-A( Only Convection Boundary Condition) 0.988 0.990 0.992 0.994 0.996 0.998 1.000 r/Ro 50 52 54 56 58 60 62 64 66 Te m pe ra tu re (K ) AR=8.25 AR=6 AR=5 AR=2.5 0.894 0.895 0.896 0.897 0.898 0.899 0.900 0.901 0.902 0.903 r/Ro 72.5 73.0 73.5 74.0 74.5 75.0 Te m pe ra tu re (K ) AR=8.25 AR=6.5 AR=5 AR=2.5 ( a) (b) Fig.(7): Aspect Ratio Effect on Temperature Distribution Along Line H-H . a) Wall Region Over the Cooling Channel b) Wall Region Below the Cooling Channel 0.88 0.90 0.92 0.94 0.96 0.98 1.00 r/Ro 56 58 60 62 64 66 68 70 72 74 76 Te m pe ra tu re (K ) AR=8 AR=6.5 AR=5 AR=2.5 Dr. Kareema E. Amori /Al-khwarizmi Engineering Journal ,Vol.3, No.3, PP 1-12(2007) 9 cooling channel 39 8 7 5 4 123 6 10 12 0 1 2 3 4 5 6 7 8 9 10 11 12 Node No. 55 65 75 85 95 105 Te m pe ra tu re (K ) AR=8.25 AR=6.5 AR=5 AR=2 (a) (b) Fig.(8):a). Nodes Numbering of Cooling Channel Wall b). Effect of Aspect Ratio on Cooling Channel Wall Temperature Profile (only Convection Boundary Condition) 0.88 0.90 0.92 0.94 0.96 0.98 1.00 r/Ro 65 70 75 80 85 90 95 100 105 Te m pe ra tu re (K ) A-A Conv E-E Conv. A-A Rad. E-E Rad. Fig.(9): Temperature Distribution Along Line A-A and E-E for Radiation and convection Boundary Conditions Dr. Kareema E. Amori /Al-khwarizmi Engineering Journal ,Vol.3, No.3, PP 1-12(2007) 10 cooling channel 39 8 7 5 4 123 6 10 12 1 2 3 4 5 6 7 8 9 10 11 12 Node No. 75 80 85 90 95 100 Te m pe ra tu re (K ) Conv. and Rad. Conv. Fig.(10): a). Nodes Numbering of Cooling Channel Wall. b) Temperature Distribution Along Cooling Channel wall 0 1 2 3 4 5 6 7 8 9 10 0 20 40 60 80 100 120 140 160 180 200 220 he at (W /m ) Conv. Rad. and Conv. AR Fig.(11 ): Aspect Ratio Effect on Heat Transferred From Hot Gas Wall for Convection and Combined Convection Radiation Boundary Condition Dr. Kareema E. Amori /Al-khwarizmi Engineering Journal ,Vol.3, No.3, PP 1-12(2007) 11 Appendix Thermal resistance of nodes belong to the interfacial line between two different materials A and B as shown in fig.(a1) can be written as: ) )( 1 )( 1 (* )( )( ,1,21 1 1 AjiAji A kkrr r R + + ∆+∆ ∆ = θ .. (a1) ) )( 1 )( 1 ( )( )( ,,121 1 1 BjiBji B kkrr r R + ∆+∆ ∆ = + θ ..(a2) BA BA RR RR R )()( )()( 11 11 1 + = ………….(a3) substituting equations (a1) and (a2) in (a3) yields: ) 11 ( ))(2/( ,,1212 1 1 jiji kkrr r R + ∆+∆∆+ ∆ = +θθ θ ……..(a4) )22( 2 )22( 2 21, 2 21,1 2 // 3 / 33 rrk r rrk r R jiji RR ∆+∆ ∆ + ∆+∆ ∆ = += − θθ ) 11 ( )( ,,121 2 3 jiji kkrr r R + ∆+∆ ∆ = − θ ….(a5) )2)(22( 2 )2)(22( 2 121, 1 1211, 1 // 4 / 44 rrk r rrk r R ji ji RR ∆−∆+∆ ∆ + ∆−∆+∆ ∆ = += − θθ θθ ) 11 ( ))(2/( ,1,211 1 4 jiji kkrR r R + ∆+∆∆− ∆ = −θθ …….(a6) The thermal resistance of nodes belong to surface exposed to fluid of T∞1 as shown in fig.(a2) can be written as: )()( 111 RRR ′′+′= ……(a7) 2/)(** 1 21 2 θθ ∆+∆ = rh R …….(a8) )()( 333 RRR ′′+′= ……….(a9) )()( 444 RRR ′′+′= ……..(a10) Fig.(a1 ): Thermal Resistances Around Node (i,j) Located at interfacial line Fig.(a2 ): Thermal Resistances Around Node (i,j) Located at surface exposed to environment B A R3 R3 1 R1R1 i+1 , j i -1, j i , j+1 i , j-1 R3 R3 R1 R1R 2 R2 R4 R4 2 r 1 2r r i , j 1ooT i , j i+1 , j i -1, j i , j-1 R3 R3 R1 R1 R2 R4 R4 1 2 2r r i , j Dr. Kareema E. Amori /Al-khwarizmi Engineering Journal ,Vol.3, No.3, PP 1-12(2007) 12 راري لجدران غرفة الدفع بوجود مجاري التبریددراسة متغیرات السلوك الح كریمة اسماعیل عموري . د قسم الھندسة المیكانیكیة/ كلیة الھندسة/ جامعة بغداد الخالصة وأستخدم معام ل (بأستخدام طریقة الفروقات المحددة التكراریة ) RL10(تم أجراء تحلیل حراري ثنائي البعد لجدار غرفة الدفع لصاروخ والذي یتكون من طبقة النیكل وأخ رى (ألستحصال توزیع درجات الحرارة خالل الجدار ) لغرض تسریع الوصول الى الحل ٠.٩تخمید ال ولق یم مختلف ة , ثم نقل حرارة باألشعاع والحمل معا, نقل حرارة بالحمل فقط , أجریت الدراسة بتطبیق ظروف حدیة متعددة ). من النحاس بینت النتائج أن أستخدام مجرى تبرید ذو نس بة باعی ة عالی ة یول د توزی ع ح راري ذو تغی رات أق ل . لمجرى التبرید) AR( من النسبة الباعیة بین ت النت ائج أن ھ ال یمك ن اھم ال ت أثیر . في حین ال تتأثر كمیة الحرارة المنتقل ة ال ى م ادة التبری د بتغیی ر النس بة الباعی ة , على مدى الجدار . األشعاع الحراري