Lecture 7: Convective Heat Transfer: Reynolds Analogy

16.512, Rocket Propulsion Prof. Manuel Martinez-Sanchez Lecture 7: Convective Heat Transfer: Reynolds Analogy Heat Transfer in Rocket Nozzles General Heat transfer to walls can affect a rocket in at least two ways: (a) Reducing the performance. This tends to be a 1-3% effect on Isp only, and is therefore secondary. (b) Creating great difficulties in the design of hot-side structures that have to survive heat fluxes in the 107 − 108 w / m2 range. The principal modes of heat transfer to nozzle and combustor walls are convection and radiation. Of these, convection dominates, and radiation tends to be important only for particle-laden flows from solid propellant rockets. Convective Heat Transfer We will review here the compressible 2D boundary layer equations in order to extract information on wall heat transfer. The governing equations are (in the B.L. approximation) Continuity ∂ ( ρu ) ∂x + ∂ ( ρv) ∂y =0 ∂u ∂u ∂p ∂τxy ∂ ⎛ ∂u ⎞ + ρv + = = ⎜µ ⎟ ∂x ∂y ∂x ∂y ∂y ⎝ ∂y ⎠ X-Momentum ρu Y-Momentum ∂P =0 ∂y (1) (2) (3) ρu Total enthalpy ∂ht ∂h ∂ ∂ ⎛ ∂T ⎞ uτxy + + ρv t = ⎜k ⎟ ∂x ∂y ∂y ∂y ⎝ ∂y ⎠ ( ) (4) u2 is the specific total enthalpy, and µ is the viscosity. For a 2 laminar flow, µ = µ ( T ) is a fluid property. Rocket boundary layers are almost always where ht = h + turbulent, and µ is then the “turbulent viscosity”, where momentum transport is effected by the random motion of turbulent “eddies”. If these eddies have a velocity scale u' and a length scale l ' , we have, in order-of-magnitude. µ turb. ∼ ρ u'l' (5) where u' is some fraction of the local u, and l ' tends to be of the order of the wall distance y. The important points about (5) are (a) µ turb. µ , mostly because l ' mean free path and (b) µ turb. is proportional to density (whereas µ is not, because the m.f.p. is inversely proportional to ρ ). Similarly, the last term on the right in the energy balance, representing the convergence of heat flux, contains the “turbulent thermal conductivity” K ∼ ρ cpu'l' . Once again, we notice that K is here proportional to density. We also note that the “turbulent Prandtl number” Pr = µ t cp kt ∼ 1 (from the orders of magnitude) It is of some interest to note the origin and composition of the viscous term in equation (4). If we collect the dot products u i τ t around a fluid element as shown (in B.L. approximation), we obtain the term ∂ uτxy ∂y ( ) as written in (4). This can be expanded as 2 ∂τxy ⎛ ∂u ⎞ ∂ ∂u ∂ ⎛ ∂u ⎞ + τxy =u uτxy = u ⎜µ ⎟ + µ⎜ ⎟ ∂y ∂y ∂y ∂y ⎝ ∂y ⎠ ⎝ ∂y ⎠ ( ) (6) The 1st term in (6) is just the velocity times the viscous net force per unit volume, so it is the part of the total viscous work that goes to accelerate the local flow. The second term in (6) is positive definite, and it is the rate of dissipation of energy into heat due to viscous effects. We will return later to this heating effect. Approximate Analysis Let us manipulate the right hand side of equation (4): ∂ ⎛ ∂u ⎞ ∂ ⎛ ∂T ⎞ ∂ ⎡ ⎛ ∂u K ∂T ⎞ ⎤ + ⎢µ ⎜ u ⎜ uµ ⎟+ ⎜K ⎟= ⎟⎥ ∂y ⎝ ∂y ⎠ ∂y ⎝ ∂y ⎠ ∂y ⎣ ⎝ ∂y µ ∂y ⎠ ⎦ and, since ∂h ∂T = cp , this yields ∂y ∂y ⎡ ⎛ u2 ⎞⎤ ∂ ⎢ ⎜∂ 2 1 ∂h ⎟ ⎥ µ + ∂y ⎢ ⎜⎜ ∂y Pr ∂y ⎟⎟ ⎥ ⎢⎣ ⎝ ⎠ ⎥⎦ µcp ⎞ ⎛ ⎜⎜ Pr ≡ ⎟ k ⎟⎠ ⎝ We note here that, both for laminar and turbulent flows, Pr is a constant, independent of P and T to a good approximation. In fact, as we noted before, it is also of order unity ( ∼ 0.9 for turbulent flows). So, the RHS of the energy equation becomes ∂ ∂y ⎡ ∂ ⎛ h u2 ⎞ ⎤ ⎢µ ⎜ + ⎟⎥ 2 ⎠⎟ ⎥⎦ ⎢⎣ ∂y ⎜⎝ Pr (7) If we made the approximation Pr = 1 , then this would reduce ∂ ⎛ ∂ht ⎞ u2 . If in addition, we made ⎜µ ⎟ with ht = h + ∂y ⎝ ∂y ⎠ 2 ∂P the “flat plate” approximation 0 , then the pair of equations ∂x (1), (4) would become further to ∂u ∂u ∂ ⎛ ∂u ⎞ + ρv = ⎜µ ⎟ ∂x ∂y ∂y ⎝ ∂y ⎠ ∂h ∂h ∂ ⎛ ∂ht ρu t + ρ v t = ⎜µ ∂x ∂y ∂y ⎝ ∂y ρu ⎫ ⎪ ⎪ ⎬ ⎞⎪ ⎟⎪ ⎠⎭ approximations (8) These are identical equations for u and ht . The same equation would also govern the linearly transformed variables u= u ; ue where the ( )e h= ht − htw (9) hte − htw subscript denotes the value of a variable in the local “external” flow (just outside the boundary layer). Both u and h satisfy identical boundary conditions: uw = hw = 0; ue = he = 1 (10) and, as noted, identical governing equations. We conclude that, under the assumption ∂P ⎛ ⎞ ⎜ Pr = 1, ∂x = 0 ⎟ , ⎝ ⎠ ht − hw u = hte − hw ue (11) u2w = hw . This similarity relation between velocity 2 and total enthalpy profiles is known as Crocco’s analogy. where we also noticed htw = hw + Approximate heat flux at the wall We are interested in the magnitude of the wall heat flux ⎛ ∂T ⎞ qw = ⎜ K ⎟ ⎝ ∂y ⎠w ⎛ K ∂h ⎞ ⎛ K ∂h qw = ⎜ µ t ⎟ =⎜ ⎜ cp ∂y ⎟ ⎜ ⎝ ⎠w ⎝ µcp ∂y (12) ⎞ ⎟ ⎟ ⎠w 0, since uw = 0 ⎛ ∂h ⎞ ⎛ ∂h ⎞ ⎛ ∂h ⎞ u ⎞ ∂ ⎛ where we used ⎜ t ⎟ = ⎜⎜ h + ⎟⎟ = ⎜ ⎟ + ⎜ u ⎟ 2 ⎠ ⎝ ∂y ⎠w ∂y ⎝ ⎝ ∂y ⎠w ⎝ ∂y ⎠w w 2 The group K 1 should be set equal to unity, for consistency with the stated = µcp Pr approximations. Thus ⎛ ∂h ⎞ qw = ⎜ µ t ⎟ ⎝ ∂y ⎠w Use now equation (11): ht − hw ⎛ ∂u ⎞ ⎛ ∂ ⎡ u ⎤⎞ qw = ⎜ µ hw + hte − hw ⎟ = e ⎢ ⎥ ⎜µ ⎟ ⎜ ∂y ue ⎦ ⎟⎠ ue ⎝ ∂y ⎠w ⎣ ⎝ w ( ) ⎛ ∂u ⎞ and notice that ⎜ µ ⎟ is the wall shear stress, τw . So ⎝ ∂y ⎠w qw = hte − hw ue (13) τw which is also called Reynolds analogy. A more compact form of this can be written in terms of the Friction Coefficient cf ≡ τw (14) 1 ρ u2 2 e e and the Stanton number St = ( qw ρeue hte − hw ) (15) with the result (from (13)) St = cf 2 (16) One important point can be made about the result (13): The heat flux to the wall is driven by the enthalpy (or temperature) difference between Total external and Wall values, not between static values. This can be nonintuitive. Consider the situation near the exit of a highly expanded space nozzle, where the bulk temperature Te may have dropped to, say, 300K due to the strong expansion from a chamber temperature of, say, 3000K. The wall could be made of Tungsten so as to be able to sustain relatively high temperature and cool itself by radiation to space, so Tw could be, say, 1500K. Is the nozzle wall being heated or cooled by the 300K gas? The answer is that it is being heated, because Tte = Tc = 3000K, while Tw = 1500 < Tt . Are we violating the 2nd Principle of e Thermodynamics. Read on. Simplified Profiles, Across the Boundary Layer To better understand this situation, let us return to Crocco’s analogy (equation 11) 2 and write ht = h + u , and solve for h: 2 ( h = hte − hw ) uu − e u2 2 (17) This is a quadratic relationship between h and u. For low subsonic flows h ht , so the last term is not strong, and the relationship becomes linear in the limit. The relationship between slopes at the wall flows from (17): 0 ⎛ dh ⎞ ⎛ ∂u ⎞ 1 ⎛ du ⎞ ⎜ ⎟ = hte − hw ⎜ ⎟ − ⎜u ⎟ dy u dy ⎝ ⎠w ⎠w ⎝ ∂y ⎠w e ⎝ ( or ) ⎛ hte − hw ⎛ dh ⎞ ⎜ du ⎟ = ⎜⎜ u ⎝ ⎠w ⎝ e ⎞ ⎟ ⎟ ⎠ (18) We can use (17) and (18) to sketch h vs. u across the boundary layer. For a case with he > hw , this looks like Note that whenever ue > hte − hw there develops an intermediate temperature maximum. But in any case, the wall slope is as if the line were coming from h t e , not from he . The case when he < hw is more revealing even: Now the wall slope is seen to be positive (heat into the wall), despite he < hw (as long as ht > hw ) e So, the quadratic portion of the Crocco relationship is responsible for the extra wall heat; this can in turn be traced to viscous dissipation, which accumulates in the boundary layer and elevates its temperature, so that the wall is heated even when the outside temperature is low (as long as the flow has high speed). Modification for Pr ≠ 1 We leave for now the issue of the non zero pressure gradient, except to note that it introduces small modifications down to the throat. The deviations of Pr from unity are small, and, for gases Pr < 1 (~ 0.9 for turbulent flow). This breaks the perfect balance between dissipation and conduction responsible for Crocco’s analogy, in the sense of favoring conduction of the dissipated heat. As a second consequence, the temperature overshoot is reduced, and so is the wall slope of T and the heat flux to the wall. The direct effect of higher conduction ( Pr < 1 ) is accounted for approximately by modifying Reynolds analogy to St = cf 2Pr 0.6 (19) The secondary effect (reduced overshoot) is accounted for by replacing the driving enthalpy difference ht − hw by haw − hw , where haw is the “Adiabatic-wall enthalpy”, e defined as haw = he + r u2e 2 ;r (20) 0.9 (turbulent) and r is the “Recovery factor”. With these changes, the heat flux is now qw = ρeue (haw − hw ) cf 2Pr 0.6 (21) The Bartz heat flux formula A very crude, but surprisingly effective representation for the friction factor cf is that supplied by the well-studied case of fully developed turbulent flow in a pipe. cf = 0.046 R 0.2 e ; Re = ρeueD (22) µe where R e is the Reynolds number based on diameter D, and µe is the laminar viscosity. Putting also h = cp T + constant, equation (21) now gives qw = ρeuecp ( Taw 0.2 0.023 ⎛ µe ⎞ − Tw ) ⎜ ⎟ Pr 0.6 ⎜⎝ ρeueD ⎟⎠ = 0.026 0.2 D ( ρeue )0.8 µ0.2 e cp ( Taw − Tw ) (23) 0.026 It is common practice to define a heat transfer “gas-side film coefficient”, hg (not an enthalpy!) by hg ≡ qw Taw − Tw (24) And, so far, we have hg = 0.026 D0.2 ( ρeue )0.8 µ0.2 e cp (25) At this point we note that the formulation so far has ignored the strong variations of ρ and µ across the boundary layer since these quantities depend on temperature as ρ ∼ 1 (at P=constant) ; µ ∼ T w ( w T 0.6 ) (26) A commonly used approach to including these variations is to replace ρe and µe in equation (25) by their values at some intermediate temperature : ⎛< T >⎞ µe → µe ⎜⎜ ⎟⎟ ⎝ Te ⎠ Te ; ρe → ρe w (27) and can be evaluated by several empirical rules. For Mach numbers not much higher than 1, we can simply use Te + Tw 2 (28) Making the replacements of (27) in equation (25), we obtain hg = 0.026 D0.2 ( ρeue ) 0.8 0.8 − 0.2w ⎛ Te ⎞ ⎜ ⎟ ⎝< T >⎠ µ0.2 e cp (29) which is one form of Bartz’ formula. A more useful form follows from the continuity equation: i P At m ρeue = = c , with A c* A R g Tc Γ (γ) , and where A is the local cross-section, and A t the throat cross-section. Substituting 2 A ⎛D ⎞ in (29), and using t = ⎜ t ⎟ , the final form is A ⎝D⎠ 0.8 hg = 0.026 ⎛ Pc ⎞ ⎜ ⎟ D0.2 ⎝ c *⎠ t 1.8 ⎛ Dt ⎞ ⎜ ⎟ ⎝D⎠ ⎛ Te ⎞ cp µ0.2 ⎟ e ⎜ ⎝< T >⎠ 0.8 − 0.2w (30) Several important trends and observations can be made now: ⎛ 1 (a) Smaller throat diameter leads to larger heat flux ⎜ ∼ 0.2 ⎜ D t ⎝ straight from the Reynolds no. dependence of cf . ( ⎞ ⎟ . This comes ⎟ ⎠ ) (b) Heat flux is almost linear in chamber pressure ∼ Pc0.8 . This limits the feasibility of high chamber pressures, which are otherwise very desirable. ⎛ ⎛ D ⎞1.8 ⎞ (c) Maximum heat flux occurs at the throat ⎜ ∼ ⎜ t ⎟ ⎟ . One critical design ⎜ ⎝D⎠ ⎟ ⎝ ⎠ consideration is therefore the thermal integrity of the throat structure. (d) Lighter gasses lead to higher heat fluxes, through the combined effects of cp 1 ⎞ ⎛ and c * ⎜ hg ∼ 0.6 ⎟ M ⎠ ⎝ 0.8 − 0.2w 0.68 ⎛ T ⎞ ⎛ Te ⎞ (e) The factor ⎜ e ⎟ is greater than unity. This ⎜ ⎟ ⎝< T >⎠ ⎝< T >⎠ enhancement of heat flux follows mainly from the fact that the gas in the boundary layer is mostly cooler than in the core, hence denser. We showed before that the turbulent heat conductivity is proportional to density. Example Consider the Space Shuttle Main Engine (SSME), which is a Hydrogen-Oxygen rocket with (roughly) these characteristics: 2.2 × 107 Pa Pc = 220 atm Tc = 3600 K M = 15g / mol r 1.25 c* = R g Tc Γ (γ) 2600 m / s cp = γ R γ −1 M µe 3 × 10−5 Kg / m / s Tthroat = 2800 J / Kg / K 2 Tc γ +1 3200 K = Te Tw = 1000 K We calculate then Te + Tw 3200 + 1000 = = 2100 K 2 2 < T >= 0.8 − 0.2w ⎛ Te ⎞ ⎜ ⎟ ⎝< T >⎠ 0.61 ⎛ 3000 ⎞ ⎜ 2100 ⎟ ⎝ ⎠ 1.3 (at the throat) and so, using equation (30), 160, 000 w / m2 / K hg and qw ( 160, 000 Tawt − 1000 ) 1 ( Taw )t qw γ −1 2⎞ ⎛ = Tt ⎜1 + r µ ⎟ = 1.057 × 3200 2 ⎝ ⎠ 0.9 3400 K (slightly less than Tc ) 160, 000 × 2400 = 3.8 × 108 W / m2 This is a very high level of heat flux. To visualize the implications, suppose this qw had to be transmitted through a thin metal plate (thickness δ , thermal conductivity k). ∆T where ∆T is the temperature drop through the metal . δ As an initial guess, suppose the metal were stainless steel (K 20 W / m / K ) , and One would have qw = K δ =1mm. Then ∆T = qw δ 3.8 × 108 × 10−3 = = 19, 000 K !! K 20 Obviously, this is unacceptable. Try using Copper instead, with K 400 W / m / K (twenty times better). This gives ∆T =950K, still not acceptable (copper would be very soft then). The plate would have to be thinner and made of copper. Not an easy problem. More rationally The St or hg should depend on x, distance from start of nozzle, since the B.L. is still developing (not fully developed). In addition, there should be some accounting for • • • acceleration property variation through B.L. cylindrical geometry The article by Rubsin and Inonye (ch. 8 in Rosenhow and Hartnett’s Handbook of Heat Transfer, McGraw-Hill, 1973) gives a general formula for turbulent B.L. In an cylinder, with acceleration: St ( x ) = s= A (and hg = ρeuecpSt ) n ⎛ρu x ⎞ s ⎜ e e eff ⎟ Fc1−n FRnθ ⎝ µe ⎠ cf 2 ch 1 found walls. A⎫ ⎬ = constants, depending on Reynolds no. based on mom. th. n⎭ R eθ > 4000 , A = 0.0131, n = 1 7 R eθ < 4000 , A = 0.0293 , n = 1 5 Fc ⎫⎪ ⎬ = Factors for property variability. Can take several nearly equivalent forms. A FR θ ⎪ ⎭ simple one from Eckert, is Fc = FRθ = ρe ρ (< Τ >) µe µ (< Τ >) Taw = Te + r = <Τ > Te ⎛ Te ⎞ ⎜ ⎟ ⎝< Τ >⎠ w (w ⎫ ⎪ ⎪ T T ⎪< Τ > = 0.28 + 0.50 w + 0.22 aw ⎬ Te Te ⎪ Te ⎪ 0.6 ) ⎪ ⎭ u2e γ −1 2⎞ ⎛ = Te ⎜ 1 + r Me ⎟ 2 2 ⎝ ⎠ r 0.9 (recovery factor) The “effective distance” xeff is related to the actual distance x through an integral (accounting for memory of past acceleration) xeff ( x ) = ∫ where f = z= x 0 f ( x ') f (x) dx ' ( ρeue zRµne haw − hw hte − hw ) 1 1−n n Fc FR1 −n θ ⎛ u2 ⎜ haw = he + r e ⎜ 2 ⎝ ⎞ ⎟ ⎟ ⎠ R=R(x)= body radius at x. For a quick estimate of R eθ , we can simplify further to the flat-plate case, in which dθ c f , = dx 2 ⎧0.0128 R1 4 eθ cf ⎪ with =⎨ 2 ⎪0.0065 R1e 6 θ ⎩ dR eθ dR ex or (R (R eθ eθ ) , and with > 4000 ) < 4000 ⎧ 0.0128 ⎧4 5 4 ⎪ 14 ⎪⎪ 5 R eθ = 0.0128 R ex ⎪ R eθ =⎨ →⎨ → R eθ ⎪ 0.0065 ⎪ 6 R 7 6 = 0.0065R ex ⎪ R1 6 ⎪⎩ 7 e eθ ⎩ R eθ R ex ⎧ 0.0366 ⎪ 0.2 θ ⎪ R ex = =⎨ x ⎪ 0.0152 ⎪ R1 7 ex ⎩ dR eθ dθ = dx dR ex ⎧0.0366 R 4 5 ex ⎪ ⎪ ⎨ ⎪ ⎪0.0152 R 6e 7 x ⎩ R eθ < 4000 (R ex < 1.99 × 101 R eθ > 4000 ) For large rockets, R ex throat tends to be ∼ 107 − 108 , so the high R e formulas should be better, despite the common use of Bartz’s formulae, which are based on the low R e formulation. Fortunately, differences tend to be small, and are marked often by other uncertainties (surface films, fluid properties). Example and Comparisons: Consider nozzle R = x tan x + x= R ± R 2 − R 2t 2 tan α R c − R 2c − R 2t with origin at x = xc = and R 2t 1 4 tan α x 2 tan α Rc = 1.5 , α = 15o Rt and going through throat at x = xt = Rt 2 tan α Using γ=1.25 and the R eθ > 4000 option, we find ⎛ xeff ⎜⎜ ⎝ Rt ⎞ = ⎟⎟ ⎠throat where xt Rt ∫ xc Rt 5 M 12 ⎛ 1.125 ⎜⎜ 2 ⎝ 1 + 0.125M R x = tan15o + Rt Rt ⎛ x ⎜ ⎝ Rt 1 ⎛ 1 + 0.125M2 and 1 ⎜ ⎜ 1.125 M 2⎝ 2.25 ⎞ ⎟⎟ ⎠ = 1.875 ⎞ ⎟⎟ ⎠ ⎛ 0.6979 ⎜⎜ 2 ⎝ 0.6515 + 0.0464 M 1 R 4 ⇒ (x) R ⎞ t o ⎟ tan15 ⎠ R ⇒ M (x) Rt ⎛x ⎞ The integration gives ⎜⎜ eff ⎟⎟ = 1.0892 ⎝ R t ⎠throat Compared to x t − xc = 1.153 Rt ⎛ ⎞ xc = 0.713 ⎟⎟ ⎜⎜ and Rt ⎝ ⎠ 0.9 ⎞ ⎟⎟ ⎠ ⎛ x d ⎜⎜ ⎝ Rt ⎞ ⎟⎟ ⎠ Since xeff appears to the 1 power, the memory/acceleration effect (up to the 7 throat) is insignificant. The throat St is then ( Taw (0.9) = 3263 ) (using r=1, Tw = 1000 K , Tc = 3300 K, (St )throat Te 2 = 3300 = 2933K ) t 2.25 0.0131 = 1 ⎛ ρ u ( x t − xc ) ⎞ 7 ⎜⎜ ⎟⎟ µ ⎝ ⎠throat 0.771 ⎛< T >⎞ ⎜ ⎟ ⎝ Te ⎠throat (3263) (0.6952) 1000 3300 = 0.28 + 0.50 + 0.22 = 0.6979 Te 2933 2933 Take Pc = 2 × 107 N / m2 , M= 25 g/mol * c = R g Tc = Γ ⇒ ( ρu ) t = 8.314 3300 0.025 2.25 = 1592 m / s ⎛ 2 ⎞ 0.5 1.25 ⎜ ⎟ ⎝ 2.25 ⎠ Pc * c = 12560 Kg / s / m2 xt 1 = = 1.866 Rt 2 tan15o xc 1.5 − 1.52 − 1 = = 0.7128 Rt 2 tan15o and, with R t = 0.3 m , 0.6 µ ⎛ T ⎞ 6.8 × 10−5 ⎜ ⎟ ⎝ 3000 ⎠ ⇒ µthroat = 6.70 × 10−5 Kg m sec One gets, (St )throat (0.00124 using xt instead of xt − xc ) = 0.00133 Using the R eθ < 4000 option (small rockets) (St )throat = 0.0293 0.68 0.2 ⎛< T >⎞ ⎛ ρ uxeff ⎞ ⎜ ⎟ ⎜ ⎟ ⎝ µ ⎠throat ⎝ Te ⎠throat and using again xeff = xt − xc , etc, we get ( St )throat = 0.00102 (0.000933 using xt ) For comparison, the “fully developed pipe flows” formulation would give St = 0.2 Rg ρ ue cp = 0.026 ⎛ c* ⎞ ⎜⎜ ⎟⎟ D0.2 ⎝ Pc ⎠ t 0.8 − 0.2w ⎛ Te ⎞ µ0.2 ⎟ e ⎜ ⎝< T >⎠ ⎛ At ⎞ ⎜ ⎟ ⎝ A ⎠ 0.9 1 at throat (St )throat = 0.000958 This is close to the R θ < 4000 results above (and, indeed, the coefficients are for R θ < 4000 ). But this appears coincidental, based on the fact that for most nozzles, ∆x ∼ R t .