WADC TECHNICAL REPORT 54-254 April, 1954 CONVECTIVE HEAT TRANSFER FROM A POROUS SURFACE MYRON TRIBUS JOHN KLEIN WRIGHT AIR DEVELOPMENT CENTER, U. S. AIR FORCE CONTRACT NO. AF 18(600)-51, E. O. NO. 462 Br-1

NOMENCLATURE (Dimensions are given in M-mass, L-length, t-time, T-temperature) C = Concentration of water vapor in air (M/L3). Cp = Heat capacity per unit mass of fluid (L2/t2T) D = Diffusion coefficient (L2/t). Ew = Rate of evaporation of water per unit area (M/L2t). G = Mass flow per unit area from surface (M/L2t). k = Thermal conductivity of fluid (ML/t3T) L = Latent heat of vaporization of water (L2/t2) m = Euler modulus (see Equation 13 and Fig. 1). p = Pressure of fluid (M/Lt2). po = Stagnation pressure of fluid stream (M/Lt2). T = Temperature of fluid (T). Th = Temperature of the supply fluid that passes through wall (T). u = Component of velocity in boundary layer parallel to surface (L/t). l = Component of velocity outside boundary layer parallel to surface (L/t). v = Component of velocity in boundary layer normal to surface (L/t). o = Velocity normal to surface at surface (L/t). x = Coordinate parallel to surface (L). y = Coordinate normal to surface (L). p = Density of fluid (M/L3). = Viscosity of fluid (M/Lt). C = Length along wedge (L). WADC TECHNICAL REPORT 54-254 iii

a = Angle in wedge (see Fig. 1). = Nondimensionalized distance normal to surface (see Equation 14). v = Kinematic viscosity of fluid (L2/t). = Stream function for flow. = Momentum boundary thickness (L) (see Equation 25). A = Operator indicating difference. WADC TECHNICAL REPORT 54-254 iv

CONVECTIVE HEAT TRANSFER FROM A POROUS SURFACE INTRODUCTION The use of porous walls in a wing leading edge to permit the outflow of hot air and thereby prevent ice formation raises some interesting technical problems. This paper is intended as a first contribution to the analytical prediction of the heat and mass transfer processes which occur in such a system. The complete problem of predicting the performance of a porous wall anti-icing system requires for its solution, that techniques be at hand for the determination of the following: Given: (a) The airfoil shape (b) The porosity of the leading edge, i.e., the air-flow pressure drop relation (c) The supply air temperature and pressure (d) The speed, altitude and angle of attack (e) The meteorological conditions, i.e., temperature, liquid water content of the cloud, drop sizes and snow content, if any Find: (a) (b) (c) (d) (e) The air outflow for each position over the airfoil The rate of impingement of moisture at each airfoil position The local surface temperatures on the airfoil The rates of evaporation of moisture The location of the transition point from laminar to turbulent boundary layer flow. For the case considered here the following assumptions are made: (1) The airflow is incompressible, laminar, and two dimensional (2) The flow of air through the porous leading edge is always directed outward (3) The fluid properties are independent of temperature (4) The effects of aerodynamic heating are neglected (5) The evaporation of the impinging water droplets is zero WADC TECHNICAL REPORT 54-254 1

Under these assumptions the method provides a means for finding the air outflow for each position over the airfoil and the local surface temperature. The calculations are started at the stagnation point and can be carried along the airfoil until the point of transition from laminar to turbulent flow is reached. Beyond this point the theory is not valid. METHOD OF SOLUTION As is usual in boundary-layer theory a coordinate system is chosen so that x represents distance along the airfoil surface from the stagnation point and y represents distance perpendicular to the airfoil surface. This coordinate system is then treated as a Cartesian orthogonal system with the justification that the thickness of the boundary layer is very much less than the radius of curvature of the airfoil surface, so that the effects of this curvature on quantities inside the boundary layer are negligible. The problem may then be stated mathematically as finding a solution to the following set of partial differential equations and boundary conditions. the momentum equation, A L4a~ a _ - -Kx t La. the continuity equation, -t -O and the energy equation, e Cp X7 f -P^ T Cb (1) (2) =-T an2 au^z (3) For equations (1) and (2) the following boundary conditions obtain ( (X, c) -- i (74) OL44 -+* (4) (5) (6) rs()-> LA, (r) + constaht as and (,o) - - = 3'to (() (7) WADC TECHNICAL REPORT 54-254 2

when ul(x) is the inviscid potential flow at the edge of the boundary layer and G(x) is the mass outflow along the airfoil surface as determined from the porosity and pressure distribution. The boundary conditions for Equation (3) are as follows ('Y:)' -C-TXO <s L>_ (8) (-A( - G() CP (TC -T) (9) where TH is air supply temperature. Equation (9) expresses the condition that the heat conduction into the airstream at the wall is equal to the enthalpy loss of the air flowing through the porous wall. The solution of the system, Equations (1) through (9), is divided into two parts. First the momentum and continuity equations are solved by a method introduced by Eckert1 and extended by Eckert and Livingood.2 Then the energy equation is solved by an approximate method on a digital computer using as input information the velocity data obtained by solving the momentum and continuity equations. The results of this last computation give the temperature at a finite set of points inside the boundary layer and along the airfoil surface. SOLUTION OF MOMENTUM AND CONTINUITY EQUATIONS Although the method used here is that of Eckert and Livingood it differs in details and is simpler because of the assumption of constant fluid properties. For these reasons and for the sake of completeness the theory is developed from the basic equations. Begin by considering the case of "wedge flows", which are those flows where the velocity of the fluid just outside the boundary layer follows the law u,3) - (10) where m and c are constants and: represents distance along the wedge from the apex. (Note that 5 is not to be' confused with x which is reserved for the arbitrary airfoil shape.) If the wedge half-angle is a, the parameter m is related to this angle and the local flow conditions by dC, De = 3 t5 (11)' 4' ", ^~~~~~~~~~(~ WADC TECHNICAL REPORT 54-254 5

and is called the Euler modulus. The following figure illustrates these quantities. Fig. 1. In wedge flows the velocity profiles in the boundary layer do not depend on the distance along the wedge ~ and normal to the surface y separately, but depend only on the local velocity Ul( ) and a parameter \ =F - Ho^ JL YV (12) where v is kinematic viscosity. This change of variable for the special boundary conditions imposed by wedge flows permits the partial differential equations of momentum and continuity to be converted into total differential equations and effect a solution by numerical means. If the wedge is porous and there is an outflow at the wedge surface, this outflow must vary in a prescribed manner along the wedge in order to make the boundary conditions dependent only n and not ~ or y. From the Bernoulli equation where Po is the stagnation pressure I, P (13) z e a differentiating with respect to x obtains _x (14) Using Equation (14), Equations (1) and (2) are rewritten as follows Uf^ au 3rrc c -^eg,1% (15) WADC TECHNICAL REPORT 54-254 4

+ - = Now define a new variable ~ = 4(~,y) such that r = T - I In order to have velocity profile a function of il alone * must have the i r - (I-am)'" (vu) j)G. f(E) Substitution of t and its derivatives into the momentum equation yields where subscripts indicate partial differentiation. Note that the contin equation is automatically satisfied by this choice of r. From the relationship between r and f('r) it is a straightforwa: cedure to compute the following "a,' - -. (~) - cA) -t ) ~,,'His a * s e12 -')( <) +Xt, ~'v L, = ^/ ^ t*( — y' f-i S^izz o54+^l1-/ ~. The boundary conditions are converted to the new variable as follows U so (0) 0 U w ^.o -= v-to (S);l~) = Jo (S) ( n ) L (Lot5) (16) (16.1) (16.2) f(orm (163.5) (17) uity rd pro(18) (19) (20) (21) (22) (23) WADC TECHNICAL REPORT 54-254 5

It is necessary, therefore, for the outflow v(S) to have the value Vr0( ) = - U. ( )(Iw^y)( U. S, 4( (24) ( fo) (24) where f(O) is constant. Hence, different outflow velocities correspond to different values of f(O), but the outflow distribution must satisfy Equation (24). The remaining boundary conditions are as 3 0 - 90 D'. e I,'o' ~l) /'YL() -o Thus, once the functions f(rl) are known for various values of m, the boundarylayer flow distribution is easily calculated from Equations (18) and (21). The momentum boundary thickness is defined by 00 > -= (, w aU (25) which in terms of TR can be written as (i+ ( )^ /= i - 1F (0 A r) (26) Once the differential Equation (23) has been solved numerically, the function F(m,f(O)) may be computed. A table of values of [2/(1 + m)]l/2 F(m,f(O)) appears on page 66 of Eckert and Livingood.2 The fundamental postulate of Eckert which permits a "patching" together of wedge solutions to determine the flow over an arbitrary body is now introduced. This postulate may be stated as follows: If at a particular point x on a laminar boundary layer the velocity just outside the boundary layer is ul(x), the velocity gradient is dul/dx, the outflow is v0(x), and the momentum thickness is,;(x), and if a wedge is found which at some point f (not necessarily equal to x), the following relations are valid., (SU) LL,(Y ) d% aT d and>(:) =/,(x) then the rate of growth of the momentum thickness along the body and along the wedge are equal, i.e., WADC TECHNICAL REPORT 54-254 6

The law for the growth of the boundary layer on a wedge is found by differentiating Equation (26). This procedure leads to as _ /v, J Thus, if at a given point x on an arbitrary airfoil v(x),,A(x), ul(x),n and dul/dx are known, the problem is to find the proper value of m and of 5 to put into this equation to find d4/d = d /dx. The equations obtained for computation are as follows ( )L (- w U) ( (27) T = _ a M' (28) (28) (, ( ) (A, ) (29),,.,, -t du- (30) H, dS Eliminate ( by substituting Equation (30) into Equations (27) through (29) to obtain __ -) -, (31) _ 1 - 8 d 1 A slight regrouping of terms yields xi F(,( b f ru, +to)) (32) <<x - s ~a, __ (55) (^/^. —' Y(yd4 o (34) To start the calculation, begin at the stagnation point. Near this point the flow corresponds to "stagnation flow" or m = 1. Hence, in this region a graph of Ul(x) will be a straight line given by a formula Ul(x) = cx and from a plot of Ul(x) from which c may then be determined. Substituting m = 1 and U1 = cx into Equations (32) through (34) results in WADC TECHNICAL REPORT 54-254 7

,3 =^ F(l,~) )( ~~C) (35) ado) s h- V(36) These represent the starting values of the calculation and Equations (355) and (36) are used until dul/dx begins to change from its initial value. Then from Equation (54) a new value of f(O) is computed. This value when substituted into Equation (32) allows a new value of m to be determined. Equation (33) is rewritten in the form ar' i9d4, aX (37) to determine the increment in 9 once m has been found. It is clear that the accuracy depends on the size of Ax and a convenient criterion for judging the proper value of Ax is a comparison of the results obtained for one value of Ax with those obtained with one-half that value. If the difference in the results are small this value of Ax is probably satisfactory, otherwise it is necessary to continue halving the value of Ax until good agreement is obtained between two successive computations. Then the values of J and m at the point x + Ax can be computed and starting with these values the procedure is repeated and continued in this manner until) and m are known as functions of x. Finally, once the variation in m is known Equations (18) and (21) (with t computed from Equation (30)) are used to determine u and v as functions of x and y. This completes the solution of the momentum and continuity equations. SOLUTION OF ENERGY EQUATION An approximate solution to the energy partial differential equation is obtained by replacing it with an equivalent difference equation, which is solved on an automatic digital computer. This procedure requires that a grid of points be chosen at which a solution is desired which raises problems as to the proper choice of a mesh size and a proper ratio of Ax to Ay in the grid. A compromise must be effected between the demands of accuracy which require a small mesh size and the desire to keep the computations at a minimum, which requires that the solutions be obtained at as few points as possible. A choice of a proper ratio of Ax to Ay is necessary to insure the convergence of the solution of the difference equation, i.e., small variations in the results should tend to decrease as the solution proceeds along the grid instead of increasing as may happen with the choice of an improper ratio. Since the theory on WADC TECHNICAL REPORT 54-254 8

choosing an optimum grid size is meager, it has been found necessary to use a trial-and-error procedure on a problem whose solution is known, so that a comparison may be made. A useful case for comparison is that of a porous flat plate (m = O) whose face is parallel to the direction of the free stream. A table of values for this case is given on page 40 of Brown and Donoughe.3 The energy equation may be rewritten as follows ar A_ I CT _ S^r 7 ase ecp A~ll w) aa' ads Ar) aat (38) with boundary condition at y = 0 -T(hL, = A k/ a aI I + r'' e cp -cr(Y., T ) a 0 H fOr (39) A grid of points is now introduced, illustrated in Fig. 2, A X X x x x (i-, j)-,-~x ( L) t i X x AX. Fig. 2. Then rewrite Equation (38) in difference form using forward differences ^- Ti, e ii _ I AX, ecp Ai.,j Ti05t (40) - L; - i.-..,,j Ao Tij is computed from the formula'TT - - -71 -.I'4.. 4. a I. (41) WADC TECHNICAL REPORT 54-254 9

Substituting Equation (40) into Equation (41) obtains Tj =C r as -U T.-^. T 7! r cC- -r ~ j -l- + (42) - CLUE ( 7T -. 4 -' where C1 and C2 are constants given by C, -a- (43) and CL - (44) Note that Equation (42) shows that the temperature at any point (i,j) in the ith column depends only on three values in the (i-l)st column, namely at (i-1, j-l), (i-1, j), and (i-l, j+l). Hence, this equation may be used to compute all points, except those on the airfoil surface (j = 0), provided an initial set of values is given. For the boundary points Equation (39) can be rewritten in difference form -, - It, (45) This equation may be solved for TiO,0 -7 TC3-; 1 "0 TH (46) C3 + Lo where c3 is a constant given by cS= ek -s (47) kC1 Equation (46) shows that the temperature at a boundary point depends on the temperature only at the point immediately above that boundary point in the same column. Hence, the procedure is to compute the temperature at the interior points of a given column using Equation (42) and then Equation (46) to compute the temperature at the boundary before moving on to the next column. The choice of the initial values to start the computation is somewhat arbitrary. One restriction is placed on these values by the condition that at the edge of the boundary layer the temperature is equal to To and has this value for all points above the boundary layer. Another restriction is given by the physical requirement that the temperature in the boundary layer lie WADC TECHNICAL REPORT 54-254

between TH and Too, since aerodynamic heating is neglected. A possible choice of starting values would be to give the point at the surface the value TH and change the values in a linear fashion until TO is reached at the edge of the boundary layer. If the mesh has been properly chosen the inaccuracies introduced in the initial values will disappear as the computation is carried along. Note that in any column it is necessary only to carry the computations out to a value of j which is equal to the value of j in the preceding column at which the temperature first equalled T,. All points above this automatically have the value To as can be seen by examination of Equation (42). Hence, the number of points at which a computation is required increases by one each time the value of i increases by one. In programming this problem for a digital computer it is necessary to devise a method of computing uil, and vii,j. One possibility is to store the values of f(Tr) and fl(i) for a set of values of tl corresponding to points inside the boundary layer, compute Tr for the point in question, and then use an interpolation procedure to find these functions at this point. Equations (18) and (21) (with t replaced by x) then give the values of u and v. The temperature may be computed directly from the appropriate expressions to complete the solution of the problem. Extension of Problem In order to extend the problem to cover the effects of evaporation of the impinging water droplets, it is necessary to add another partial differential equation to the basic set, Equations (1) through (3). This is the diffusion equation, which for the assumptions made here may be written in the form U^a- + v = (48) where C is concentration of water vapor, and with the boundary conditions Cma,) - c g -s5 (49) -D (a EW (50) C(,9) = a function of T(x,O), i.e., the vapor pressure of water at the surface is taken from tables of the saturation pressure, (51) where E is the rate of evaporation of water per unit area. In addition the boundary condition in Equation (9) must be modified to include the effects of evaporation as follows WADC TECHNICAL REPORT 54-254 11

n(,),o 3 = G(c) (H -T) + L ( a ) (52) where the last term represents the energy absorbed by the vaporizing water. If it is assumed that the effect of the evaporation on the boundarylayer velocity distribution is negligible, then the momentum and continuity equations may again be solved separately from the energy and diffusion equations by the method described above. It would then be necessary to devise an approximate procedure to solve the energy and diffusion equations together, since they are connected through Equations (51) and (52). Through the use of a difference equation technique it should be possible to develop such a procedure for use on a digital computer, although this has not yet been attempted. CONCLUDING REMARKS The methods described in this report have been used to set up the problem of finding the temperature distribution above a porous flat plate for solution by the MIDAC, an automatic digital computer located at the Willow Run Research Laboratory of the University of Michigan. Funds were not available to make the computation, but preliminary results indicate that it will be possible to make a comparison with the known solution for this case in order to check the feasibility of the procedure. If the techniques are successful the extension to an arbitrary airfoil is a straightforward process as the only change in the computer code will be in the methods used for computing u and v. Thus, the basic logic and most of the detailed steps in the machine program can be used directly in the more complicated case. REFERENCES 1. Eckert, E.R.G., "Calculation of the Heat Transfer in the Boundary Layer of Bodies in Flow", VDI Forschung 416, Translation, University of California, Berkeley, April 16, 1948. 2. Eckert, and Livingood, "Method for Calculation of Heat Transfer in Laminar Region of Air Flow around Cylinders of Arbitrary Cross Section (Including Large Temperature Differences and Transpiration Cooling)", NACA TN2733, June, 1952. WADC TECHNICAL REPORT 54-254 12

3. Brown, W. B., and Donoughe, P. L., "Tables of Exact Laminar-Boundary-Layer Solutions when the Wall is Porous and Fluid Properties are Variable", NACA TN2479, September, 1951. WADC TECHNICAL REPORT 54-254 13