cen58933_ch05 - cen58933_ch05.qxd 9/4/2002 11:41 AM Page...

Info iconThis preview shows page 1. Sign up to view the full content.

View Full Document Right Arrow Icon
This is the end of the preview. Sign up to access the rest of the document.

Unformatted text preview: cen58933_ch05.qxd 9/4/2002 11:41 AM Page 265 CHAPTER NUMERICAL METHODS I N H E AT C O N D U C T I O N o far we have mostly considered relatively simple heat conduction problems involving simple geometries with simple boundary conditions because only such simple problems can be solved analytically. But many problems encountered in practice involve complicated geometries with complex boundary conditions or variable properties and cannot be solved analytically. In such cases, sufficiently accurate approximate solutions can be obtained by computers using a numerical method. Analytical solution methods such as those presented in Chapter 2 are based on solving the governing differential equation together with the boundary conditions. They result in solution functions for the temperature at every point in the medium. Numerical methods, on the other hand, are based on replacing the differential equation by a set of n algebraic equations for the unknown temperatures at n selected points in the medium, and the simultaneous solution of these equations results in the temperature values at those discrete points. There are several ways of obtaining the numerical formulation of a heat conduction problem, such as the finite difference method, the finite element method, the boundary element method, and the energy balance (or control volume) method. Each method has its own advantages and disadvantages, and each is used in practice. In this chapter we will use primarily the energy balance approach since it is based on the familiar energy balances on control volumes instead of heavy mathematical formulations, and thus it gives a better physical feel for the problem. Besides, it results in the same set of algebraic equations as the finite difference method. In this chapter, the numerical formulation and solution of heat conduction problems are demonstrated for both steady and transient cases in various geometries. S 5 CONTENTS 5–1 Why Numerical Methods 266 5–2 Finite Difference Formulation of Differential Equations 269 5–3 One-Dimensional Steady Heat Conduction 272 5–4 Two-Dimensional Steady Heat Conduction 282 5–5 Transient Heat Conduction 291 Topic of Special Interest: Controlling the Numerical Error 309 265 cen58933_ch05.qxd 9/4/2002 11:41 AM Page 266 266 HEAT TRANSFER 5–1 · g0 1d — — r 2 dT + — = 0 — r 2 dr dr k () r0 Sphere 0 dT(0) ——– = 0 dr I WHY NUMERICAL METHODS? The ready availability of high-speed computers and easy-to-use powerful software packages has had a major impact on engineering education and practice in recent years. Engineers in the past had to rely on analytical skills to solve significant engineering problems, and thus they had to undergo a rigorous training in mathematics. Today’s engineers, on the other hand, have access to a tremendous amount of computation power under their fingertips, and they mostly need to understand the physical nature of the problem and interpret the results. But they also need to understand how calculations are performed by the computers to develop an awareness of the processes involved and the limitations, while avoiding any possible pitfalls. In Chapter 2 we solved various heat conduction problems in various geometries in a systematic but highly mathematical manner by (1) deriving the governing differential equation by performing an energy balance on a differential volume element, (2) expressing the boundary conditions in the proper mathematical form, and (3) solving the differential equation and applying the boundary conditions to determine the integration constants. This resulted in a solution function for the temperature distribution in the medium, and the solution obtained in this manner is called the analytical solution of the problem. For example, the mathematical formulation of one-dimensional steady heat conduction in a sphere of radius r0 whose outer surface is maintained at a uni· form temperature of T1 with uniform heat generation at a rate of g0 was expressed as (Fig. 5–1) T(r0) = T1 1 d 2 dT r dr r 2 dr dT(0) 0 dr Solution: · g0 T(r) = T1 + — (r02 – r 2) 6k · 4 π g0r 3 · dT ——— Q(r) = –kA — = dr 3 FIGURE 5–1 The analytical solution of a problem requires solving the governing differential equation and applying the boundary conditions. · g0 k 0 and T(r0) T1 (5-1) whose (analytical) solution is T(r) T1 · g0 2 (r 6k 0 r 2) (5-2) This is certainly a very desirable form of solution since the temperature at any point within the sphere can be determined simply by substituting the r-coordinate of the point into the analytical solution function above. The analytical solution of a problem is also referred to as the exact solution since it satisfies the differential equation and the boundary conditions. This can be verified by substituting the solution function into the differential equation and the boundary conditions. Further, the rate of heat flow at any location within the sphere or its surface can be determined by taking the derivative of the solution function T(r) and substituting it into Fourier’s law as · Q (r) kA dT dr k(4 r 2) · g0r 3k · 4 g0r 3 3 (5-3) The analysis above did not require any mathematical sophistication beyond the level of simple integration, and you are probably wondering why anyone would ask for something else. After all, the solutions obtained are exact and cen58933_ch05.qxd 9/4/2002 11:41 AM Page 267 267 CHAPTER 5 easy to use. Besides, they are instructive since they show clearly the functional dependence of temperature and heat transfer on the independent variable r. Well, there are several reasons for searching for alternative solution methods. T ,h k = constant No radiation T ,h h, T h = constant T = constant FIGURE 5–2 Analytical solution methods are limited to simplified problems in simple geometries. An oval-shaped body 2 Better Modeling We mentioned earlier that analytical solutions are exact solutions since they do not involve any approximations. But this statement needs some clarification. Distinction should be made between an actual real-world problem and the mathematical model that is an idealized representation of it. The solutions we get are the solutions of mathematical models, and the degree of applicability of these solutions to the actual physical problems depends on the accuracy of the model. An “approximate” solution of a realistic model of a physical problem is usually more accurate than the “exact” solution of a crude mathematical model (Fig. 5–3). When attempting to get an analytical solution to a physical problem, there is always the tendency to oversimplify the problem to make the mathematical model sufficiently simple to warrant an analytical solution. Therefore, it is common practice to ignore any effects that cause mathematical complications such as nonlinearities in the differential equation or the boundary conditions. So it comes as no surprise that nonlinearities such as temperature dependence of thermal conductivity and the radiation boundary conditions are seldom considered in analytical solutions. A mathematical model intended for a numerical solution is likely to represent the actual problem better. Therefore, the numerical solution of engineering problems has now become the norm rather than the exception even when analytical solutions are available. No radiation Long cylinder 1 Limitations Analytical solution methods are limited to highly simplified problems in simple geometries (Fig. 5–2). The geometry must be such that its entire surface can be described mathematically in a coordinate system by setting the variables equal to constants. That is, it must fit into a coordinate system perfectly with nothing sticking out or in. In the case of one-dimensional heat conduction in a solid sphere of radius r0, for example, the entire outer surface can be described by r r0. Likewise, the surfaces of a finite solid cylinder of radius r0 and height H can be described by r r0 for the side surface and z 0 and z H for the bottom and top surfaces, respectively. Even minor complications in geometry can make an analytical solution impossible. For example, a spherical object with an extrusion like a handle at some location is impossible to handle analytically since the boundary conditions in this case cannot be expressed in any familiar coordinate system. Even in simple geometries, heat transfer problems cannot be solved analytically if the thermal conditions are not sufficiently simple. For example, the consideration of the variation of thermal conductivity with temperature, the variation of the heat transfer coefficient over the surface, or the radiation heat transfer on the surfaces can make it impossible to obtain an analytical solution. Therefore, analytical solutions are limited to problems that are simple or can be simplified with reasonable approximations. h, T Simplified model Realistic model A sphere Exact (analytical) solution of model, but crude solution of actual problem Approximate (numerical) solution of model, but accurate solution of actual problem FIGURE 5–3 The approximate numerical solution of a real-world problem may be more accurate than the exact (analytical) solution of an oversimplified model of that problem. cen58933_ch05.qxd 9/4/2002 11:41 AM Page 268 268 HEAT TRANSFER z 3 Flexibility Engineering problems often require extensive parametric studies to understand the influence of some variables on the solution in order to choose the right set of variables and to answer some “what-if ” questions. This is an iterative process that is extremely tedious and time-consuming if done by hand. Computers and numerical methods are ideally suited for such calculations, and a wide range of related problems can be solved by minor modifications in the code or input variables. Today it is almost unthinkable to perform any significant optimization studies in engineering without the power and flexibility of computers and numerical methods. L T T(r, z) 0 r0 T0 4 Complications r Analytical solution: T(r, z) – T J0(λnr) sinh λn(L – z) ————– = ∑ ———— –————— λnJ1( λnr0) sinh (λnL) T0 – T n=1 where λn’s are roots of J0(λnr0) = 0 FIGURE 5–4 Some analytical solutions are very complex and difficult to use. Some problems can be solved analytically, but the solution procedure is so complex and the resulting solution expressions so complicated that it is not worth all that effort. With the exception of steady one-dimensional or transient lumped system problems, all heat conduction problems result in partial differential equations. Solving such equations usually requires mathematical sophistication beyond that acquired at the undergraduate level, such as orthogonality, eigenvalues, Fourier and Laplace transforms, Bessel and Legendre functions, and infinite series. In such cases, the evaluation of the solution, which often involves double or triple summations of infinite series at a specified point, is a challenge in itself (Fig. 5–4). Therefore, even when the solutions are available in some handbooks, they are intimidating enough to scare prospective users away. 5 Human Nature FIGURE 5–5 The ready availability of highpowered computers with sophisticated software packages has made numerical solution the norm rather than the exception. As human beings, we like to sit back and make wishes, and we like our wishes to come true without much effort. The invention of TV remote controls made us feel like kings in our homes since the commands we give in our comfortable chairs by pressing buttons are immediately carried out by the obedient TV sets. After all, what good is cable TV without a remote control. We certainly would love to continue being the king in our little cubicle in the engineering office by solving problems at the press of a button on a computer (until they invent a remote control for the computers, of course). Well, this might have been a fantasy yesterday, but it is a reality today. Practically all engineering offices today are equipped with high-powered computers with sophisticated software packages, with impressive presentation-style colorful output in graphical and tabular form (Fig. 5–5). Besides, the results are as accurate as the analytical results for all practical purposes. The computers have certainly changed the way engineering is practiced. The discussions above should not lead you to believe that analytical solutions are unnecessary and that they should be discarded from the engineering curriculum. On the contrary, insight to the physical phenomena and engineering wisdom is gained primarily through analysis. The “feel” that engineers develop during the analysis of simple but fundamental problems serves as an invaluable tool when interpreting a huge pile of results obtained from a computer when solving a complex problem. A simple analysis by hand for a limiting case can be used to check if the results are in the proper range. Also, cen58933_ch05.qxd 9/4/2002 11:41 AM Page 269 269 CHAPTER 5 nothing can take the place of getting “ball park” results on a piece of paper during preliminary discussions. The calculators made the basic arithmetic operations by hand a thing of the past, but they did not eliminate the need for instructing grade school children how to add or multiply. In this chapter, you will learn how to formulate and solve heat transfer problems numerically using one or more approaches. In your professional life, you will probably solve the heat transfer problems you come across using a professional software package, and you are highly unlikely to write your own programs to solve such problems. (Besides, people will be highly skeptical of the results obtained using your own program instead of using a wellestablished commercial software package that has stood the test of time.) The insight you will gain in this chapter by formulating and solving some heat transfer problems will help you better understand the available software packages and be an informed and responsible user. 5–2 I FINITE DIFFERENCE FORMULATION OF DIFFERENTIAL EQUATIONS The numerical methods for solving differential equations are based on replacing the differential equations by algebraic equations. In the case of the popular finite difference method, this is done by replacing the derivatives by differences. Below we will demonstrate this with both first- and second-order derivatives. But first we give a motivational example. Consider a man who deposits his money in the amount of A0 $100 in a savings account at an annual interest rate of 18 percent, and let us try to determine the amount of money he will have after one year if interest is compounded continuously (or instantaneously). In the case of simple interest, the money will earn $18 interest, and the man will have 100 100 0.18 $118.00 in his account after one year. But in the case of compounding, the interest earned during a compounding period will also earn interest for the remaining part of the year, and the year-end balance will be greater than $118. For example, if the money is compounded twice a year, the balance will be 100 100 (0.18/2) $109 after six months, and 109 109 (0.18/2) $118.81 at the end of the year. We could also determine the balance A directly from A A0(1 i)n ($100)(1 0.09)2 $118.81 (5-4) where i is the interest rate for the compounding period and n is the number of periods. Using the same formula, the year-end balance is determined for monthly, daily, hourly, minutely, and even secondly compounding, and the results are given in Table 5–1. Note that in the case of daily compounding, the year-end balance will be $119.72, which is $1.72 more than the simple interest case. (So it is no wonder that the credit card companies usually charge interest compounded daily when determining the balance.) Also note that compounding at smaller time intervals, even at the end of each second, does not change the result, and we suspect that instantaneous compounding using “differential” time intervals dt will give the same result. This suspicion is confirmed by obtaining the differential TABLE 5–1 Year-end balance of a $100 account earning interest at an annual rate of 18 percent for various compounding periods Compounding Period Number of Year-End Periods, n Balance 1 year 1 $118.00 6 months 2 118.81 1 month 12 119.56 1 week 52 119.68 1 day 365 119.72 1 hour 8760 119.72 1 minute 525,600 119.72 1 second 31,536,000 119.72 Instantaneous 119.72 cen58933_ch05.qxd 9/4/2002 11:41 AM Page 270 270 HEAT TRANSFER equation dA/dt stitution yields iA for the balance A, whose solution is A A ($100)exp(0.18 1) A0 exp(it). Sub- $119.72 which is identical to the result for daily compounding. Therefore, replacing a differential time interval dt by a finite time interval of t 1 day gave the same result, which leads us into believing that reasonably accurate results can be obtained by replacing differential quantities by sufficiently small differences. Next, we develop the finite difference formulation of heat conduction problems by replacing the derivatives in the differential equations by differences. In the following section we will do it using the energy balance method, which does not require any knowledge of differential equations. Derivatives are the building blocks of differential equations, and thus we first give a brief review of derivatives. Consider a function f that depends on x, as shown in Figure 5–6. The first derivative of f(x) at a point is equivalent to the slope of a line tangent to the curve at that point and is defined as f (x) df(x) dx f (x + ∆ x) lim x→0 f x lim f(x x→0 x) x f(x) (5-5) ∆f f (x) which is the ratio of the increment f of the function to the increment x of the independent variable as x → 0. If we don’t take the indicated limit, we will have the following approximate relation for the derivative: ∆x Tangent line x x + ∆x FIGURE 5–6 The derivative of a function at a point represents the slope of the function at that point. Plane wall Tm + 1 Tm Tm – 1 01 2 L m m +1… M M –1 m– 1 m+ 1 – – 2 2 … m –1 x FIGURE 5–7 Schematic of the nodes and the nodal temperatures used in the development of the finite difference formulation of heat transfer in a plane wall. f(x x) x f(x) (5-6) This approximate expression of the derivative in terms of differences is the finite difference form of the first derivative. The equation above can also be obtained by writing the Taylor series expansion of the function f about the point x, f(x T (x) 0 df(x) dx x x) f(x) x df(x) dx 1 2 d 2f(x) x 2 dx2 ··· (5-7) and neglecting all the terms in the expansion except the first two. The first term neglected is proportional to x2, and thus the error involved in each step of this approximation is also proportional to x2. However, the commutative error involved after M steps in the direction of length L is proportional to x since M x2 (L/ x) x2 L x. Therefore, the smaller the x, the smaller the error, and thus the more accurate the approximation. Now consider steady one-dimensional heat transfer in a plane wall of thickness L with heat generation. The wall is subdivided into M sections of equal thickness x L/M in the x-direction, separated by planes passing through M 1 points 0, 1, 2, . . . , m 1, m, m 1, . . . , M called nodes or nodal points, as shown in Figure 5–7. The x-coordinate of any point m is simply xm m x, and the temperature at that point is simply T(xm) Tm. The heat conduction equation involves the second derivatives of temperature with respect to the space variables, such as d 2T/dx2, and the finite difference formulation is based on replacing the second derivatives by appropriate cen58933_ch05.qxd 9/4/2002 11:41 AM Page 271 271 CHAPTER 5 differences. But we need to start the process with first derivatives. Using 1 Eq. 5–6, the first derivative of temperature dT/dx at the midpoints m 2 and 1 m 2 of the sections surrounding the node m can be expressed as dT dx Tm Tm x 1 2 m 1 dT dx and Tm m 1 2 Tm 1 (5-8) x Noting that the second derivative is simply the derivative of the first derivative, the second derivative of temperature at node m can be expressed as dT dx 2 dT dx2 m dT dx 1 2 Tm 1 2 m x Tm 1 Tm x 1 x x m Tm Tm 1 2Tm x2 Tm 1 (5-9) which is the finite difference representation of the second derivative at a general internal node m. Note that the second derivative of temperature at a node m is expressed in terms of the temperatures at node m and its two neighboring nodes. Then the differential equation · g k d 2T dx 2 0 (5-10) which is the governing equation for steady one-dimensional heat transfer in a plane wall with heat generation and constant thermal conductivity, can be expressed in the finite difference form as (Fig. 5–8) Tm 1 2T m x2 Tm 1 · gm k 0, m 1, 2, 3, . . . , M 1 (5-11) · where gm is the rate of heat generation per unit volume at node m. If the surface temperatures T0 and TM are specified, the application of this equation to each of the M 1 interior nodes results in M 1 equations for the determination of M 1 unknown temperatures at the interior nodes. Solving these equations simultaneously gives the temperature values at the nodes. If the temperatures at the outer surfaces are not known, then we need to obtain two more equations in a similar manner using the specified boundary conditions. Then the unknown temperatures at M 1 nodes are determined by solving the resulting system of M 1 equations in M 1 unknowns simultaneously. Note that the boundary conditions have no effect on the finite difference formulation of interior nodes of the medium. This is not surprising since the control volume used in the development of the formulation does not involve any part of the boundary. You may recall that the boundary conditions had no effect on the differential equation of heat conduction in the medium either. The finite difference formulation above can easily be extended to two- or three-dimensional heat transfer problems by replacing each second derivative by a difference equation in that direction. For example, the finite difference formulation for steady two-dimensional heat conduction in a region with Plane wall Differential equation: · g d 2T — —– + =0 2 k dx Valid at every point Finite difference equation: · Tm – 1 – 2Tm + Tm + 1 gm —–—————— + — = 0 — k ∆x2 Valid at discrete points ∆x FIGURE 5–8 The differential equation is valid at every point of a medium, whereas the finite difference equation is valid at discrete points (the nodes) only. cen58933_ch05.qxd 9/4/2002 11:41 AM Page 272 272 HEAT TRANSFER heat generation and constant thermal conductivity can be expressed in rectangular coordinates as (Fig. 5–9) Tm m, n + 1 n+1 ∆y m – 1, n ∆y n n–1 m, n m–1 m m+1 5–3 · Qcond, left Volume element of node m · gm · Qcond, right A general interior node 0 01 2 1, n Tm, n 1 2Tm, n y2 Tm, n 1 · gm, n k 0 (5-12) for m 1, 2, 3, . . . , M 1 and n 1, 2, 3, . . . , N 1 at any interior node (m, n). Note that a rectangular region that is divided into M equal subregions in the x-direction and N equal subregions in the y-direction has a total of (M 1)(N 1) nodes, and Eq. 5–12 can be used to obtain the finite difference equations at (M 1)(N 1) of these nodes (i.e., all nodes except those at the boundaries). The finite difference formulation is given above to demonstrate how difference equations are obtained from differential equations. However, we will use the energy balance approach in the following sections to obtain the numerical formulation because it is more intuitive and can handle boundary conditions more easily. Besides, the energy balance approach does not require having the differential equation before the analysis. FIGURE 5–9 Finite difference mesh for twodimensional conduction in rectangular coordinates. Plane wall Tm m + 1, n ∆x ∆x x 2Tm, n x2 m, n – 1 y 1, n m–1 m m+1 ∆x ∆x ∆x FIGURE 5–10 The nodal points and volume elements for the finite difference formulation of one-dimensional conduction in a plane wall. L M x I ONE-DIMENSIONAL STEADY HEAT CONDUCTION In this section we will develop the finite difference formulation of heat conduction in a plane wall using the energy balance approach and discuss how to solve the resulting equations. The energy balance method is based on subdividing the medium into a sufficient number of volume elements and then applying an energy balance on each element. This is done by first selecting the nodal points (or nodes) at which the temperatures are to be determined and then forming elements (or control volumes) over the nodes by drawing lines through the midpoints between the nodes. This way, the interior nodes remain at the middle of the elements, and the properties at the node such as the temperature and the rate of heat generation represent the average properties of the element. Sometimes it is convenient to think of temperature as varying linearly between the nodes, especially when expressing heat conduction between the elements using Fourier ’s law. To demonstrate the approach, again consider steady one-dimensional heat · transfer in a plane wall of thickness L with heat generation g(x) and constant conductivity k. The wall is now subdivided into M equal regions of thickness x L/M in the x-direction, and the divisions between the regions are selected as the nodes. Therefore, we have M 1 nodes labeled 0, 1, 2, . . . , m 1, m, m 1, . . . , M, as shown in Figure 5–10. The x-coordinate of any node m is simply xm m x, and the temperature at that point is T(xm) Tm. Elements are formed by drawing vertical lines through the midpoints between the nodes. Note that all interior elements represented by interior nodes are full-size elements (they have a thickness of x), whereas the two elements at the boundaries are half-sized. To obtain a general difference equation for the interior nodes, consider the element represented by node m and the two neighboring nodes m 1 and m 1. Assuming the heat conduction to be into the element on all surfaces, an energy balance on the element can be expressed as cen58933_ch05.qxd 9/4/2002 11:41 AM Page 273 273 CHAPTER 5 Rate of heat conduction at the left surface Rate of heat conduction at the right surface Rate of heat generation inside the element Rate of change of the energy content of the element or · Q cond, left · Q cond, right Eelement t · Gelement 0 (5-13) since the energy content of a medium (or any part of it) does not change under steady conditions and thus Eelement 0. The rate of heat generation within the element can be expressed as · Gelement · gmVelement · gm A x (5-14) · where gm is the rate of heat generation per unit volume in W/m3 evaluated at node m and treated as a constant for the entire element, and A is heat transfer area, which is simply the inner (or outer) surface area of the wall. Recall that when temperature varies linearly, the steady rate of heat conduction across a plane wall of thickness L can be expressed as · Q cond kA T L (5-15) where T is the temperature change across the wall and the direction of heat transfer is from the high temperature side to the low temperature. In the case of a plane wall with heat generation, the variation of temperature is not linear and thus the relation above is not applicable. However, the variation of temperature between the nodes can be approximated as being linear in the determination of heat conduction across a thin layer of thickness x between two nodes (Fig. 5–11). Obviously the smaller the distance x between two nodes, the more accurate is this approximation. (In fact, such approximations are the reason for classifying the numerical methods as approximate solution methods. In the limiting case of x approaching zero, the formulation becomes exact and we obtain a differential equation.) Noting that the direction of heat transfer on both surfaces of the element is assumed to be toward the node m, the rate of heat conduction at the left and right surfaces can be expressed as · Q cond, left kA Tm Tm 1 · Q cond, right and x kA Tm Tm 1 x (5-16) Tm Tm 1 x kA Tm Tm 1 x · gm A x 0 (5-17) which simplifies to Tm 1 2Tm x2 Tm 1 · gm k 0, m 1, 2, 3, . . . , M 1 (5-18) Volume element k Tm + 1 Tm Linear Linear ∆x m–1 Tm – 1 – Tm k A ————– ∆x A Substituting Eqs. 5–14 and 5–16 into Eq. 5–13 gives kA Tm – 1 ∆x m m+1 Tm + 1 – Tm k A ————– ∆x A FIGURE 5–11 In finite difference formulation, the temperature is assumed to vary linearly between the nodes. cen58933_ch05.qxd 9/4/2002 11:41 AM Page 274 274 HEAT TRANSFER · g2 A∆ x T1 – T2 k A ——— ∆x T2 – T3 k A ——— ∆x 1 2 3 Volume element of node 2 T2 – T3 · T1 – T2 k A ——— – k A ——— + g2 A∆ x = 0 ∆x ∆x or · T1 – 2T2 + T3 + g2 A∆ x 2 / k = 0 (a) Assuming heat transfer to be out of the volume element at the right surface. · g2 A∆ x T1 – T2 k A ——— ∆x T3 – T2 k A ——— ∆x 1 2 which is identical to the difference equation (Eq. 5–11) obtained earlier. Again, this equation is applicable to each of the M 1 interior nodes, and its application gives M 1 equations for the determination of temperatures at M 1 nodes. The two additional equations needed to solve for the M 1 unknown nodal temperatures are obtained by applying the energy balance on the two elements at the boundaries (unless, of course, the boundary temperatures are specified). You are probably thinking that if heat is conducted into the element from both sides, as assumed in the formulation, the temperature of the medium will have to rise and thus heat conduction cannot be steady. Perhaps a more realistic approach would be to assume the heat conduction to be into the element on the left side and out of the element on the right side. If you repeat the formulation using this assumption, you will again obtain the same result since the heat conduction term on the right side in this case will involve Tm Tm 1 instead of Tm 1 Tm, which is subtracted instead of being added. Therefore, the assumed direction of heat conduction at the surfaces of the volume elements has no effect on the formulation, as shown in Figure 5–12. (Besides, the actual direction of heat transfer is usually not known.) However, it is convenient to assume heat conduction to be into the element at all surfaces and not worry about the sign of the conduction terms. Then all temperature differences in conduction relations are expressed as the temperature of the neighboring node minus the temperature of the node under consideration, and all conduction terms are added. 3 Volume element of node 2 T3 – T2 · T1 – T2 k A ——— + k A ——— + g2 A∆ x = 0 ∆x ∆x or · T1 – 2T2 + T3 + g2 A∆ x 2 / k = 0 (b) Assuming heat transfer to be into the volume element at all surfaces. FIGURE 5–12 The assumed direction of heat transfer at surfaces of a volume element has no effect on the finite difference formulation. Boundary Conditions Above we have developed a general relation for obtaining the finite difference equation for each interior node of a plane wall. This relation is not applicable to the nodes on the boundaries, however, since it requires the presence of nodes on both sides of the node under consideration, and a boundary node does not have a neighboring node on at least one side. Therefore, we need to obtain the finite difference equations of boundary nodes separately. This is best done by applying an energy balance on the volume elements of boundary nodes. Boundary conditions most commonly encountered in practice are the specified temperature, specified heat flux, convection, and radiation boundary conditions, and here we develop the finite difference formulations for them for the case of steady one-dimensional heat conduction in a plane wall of thickness L as an example. The node number at the left surface at x 0 is 0, and at the right surface at x L it is M. Note that the width of the volume element for either boundary node is x/2. The specified temperature boundary condition is the simplest boundary condition to deal with. For one-dimensional heat transfer through a plane wall of thickness L, the specified temperature boundary conditions on both the left and right surfaces can be expressed as (Fig. 5–13) T(0) T(L) T0 TM Specified value Specified value (5-19) where T0 and Tm are the specified temperatures at surfaces at x 0 and x L, respectively. Therefore, the specified temperature boundary conditions are cen58933_ch05.qxd 9/4/2002 11:41 AM Page 275 275 CHAPTER 5 incorporated by simply assigning the given surface temperatures to the boundary nodes. We do not need to write an energy balance in this case unless we decide to determine the rate of heat transfer into or out of the medium after the temperatures at the interior nodes are determined. When other boundary conditions such as the specified heat flux, convection, radiation, or combined convection and radiation conditions are specified at a boundary, the finite difference equation for the node at that boundary is obtained by writing an energy balance on the volume element at that boundary. The energy balance is again expressed as · Gelement · Q 0 Plane wall 35°C 0 82°C 0 1 for heat transfer under steady conditions. Again we assume all heat transfer to be into the volume element from all surfaces for convenience in formulation, except for specified heat flux since its direction is already specified. Specified heat flux is taken to be a positive quantity if into the medium and a negative quantity if out of the medium. Then the finite difference formulation at the node m 0 (at the left boundary where x 0) of a plane wall of thickness L during steady one-dimensional heat conduction can be expressed as (Fig. 5–14) kA T0 T1 · g0(A x/2) x 0 T1 T0 · g0(A x/2) x · Special case: Insulated Boundary (q0 kA T1 T0 0 (5-22) 0) · g0(A x/2) x 0 (5-23) 2. Convection Boundary Condition hA(T T0) kA T1 T0 x · g0(A x/2) FIGURE 5–13 Finite difference formulation of specified temperature boundary conditions on both surfaces of a plane wall. ∆x — – 2 · g0 1. Specified Heat Flux Boundary Condition kA M T0 = 35°C TM = 82°C (5-21) where A x/2 is the volume of the volume element (note that the boundary ele· ment has half thickness), g0 is the rate of heat generation per unit volume (in 3 W/m ) at x 0, and A is the heat transfer area, which is constant for a plane wall. Note that we have x in the denominator of the second term instead of x/2. This is because the ratio in that term involves the temperature difference between nodes 0 and 1, and thus we must use the distance between those two nodes, which is x. The finite difference form of various boundary conditions can be obtained · from Eq. 5–21 by replacing Q left surface by a suitable expression. Next this is done for various boundary conditions at the left boundary. · q 0A L … (5-20) all sides · Q left surface 2 0 (5-24) Volume element of node 0 T1 – T0 k A ——— ∆x · Qleft surface 0 0 1 ∆x 2 … L x ∆x T1 – T0 · ∆ x · – Qleft surface + k A ——— + g0 A — = 0 ∆x 2 FIGURE 5–14 Schematic for the finite difference formulation of the left boundary node of a plane wall. cen58933_ch05.qxd 9/4/2002 11:41 AM Page 276 276 HEAT TRANSFER 3. Radiation Boundary Condition Tsurr 4 A(Tsurr ∆x — – 2 ε T1 – T0 k A ——— ∆x hA(T – T0) 0 0 1 2 ∆x A … hA(T L x T0 · g0(A x/2) x 0 (5-25) · q 0A Medium B kB kAA Tm + 1 – Tm kB A ————– ∆x m–1 m ∆x x m+1 ∆x ∆x ∆x —– —– 22 kA T1 T0 x · g0(A x/2) 0 (5-26) T0) kA T1 T0 x · g0(A x/2) 0 (5-27) hA(T T0) 4 A(Tsurr T04) kA T1 T0 x · g0(A x/2) 0 (5-28) 6. Interface Boundary Condition Two different solid media A and B are assumed to be in perfect contact, and thus at the same temperature at the interface at node m (Fig. 5–16). Subscripts A and B indicate properties of media A and B, respectively. Interface Tm – 1 – Tm kA A ————– ∆x T04) 5. Combined Convection, Radiation, and Heat Flux Boundary Condition FIGURE 5–15 Schematic for the finite difference formulation of combined convection and radiation on the left boundary of a plane wall. · · gA,m gB,m 4 A(Tsurr T0) hcombined A(T T1 – T0 · ∆ x + kA ——— + g0 A —– = 0 ∆x 2 A T1 or ∆x hA(T – T0) + εσ A(T 4 – T 4 ) surr 0 Medium A kA kA 4. Combined Convection and Radiation Boundary Condition (Fig. 5–15) · g0 εσ A(T 4 – T 4 ) surr 0 T04) A Tm + 1 – Tm Tm – 1 – Tm kA A ————– + kB A ————– ∆x ∆x ∆x · ∆x · A —– + g A —– = 0 + gA, m B, m 2 2 FIGURE 5–16 Schematic for the finite difference formulation of the interface boundary condition for two mediums A and B that are in perfect thermal contact. Tm Tm 1 x kBA Tm Tm 1 x · gA, m(A x/2) · gB, m(A x/2) 0 (5-29) · In these relations, q0 is the specified heat flux in W/m2, h is the convection coefficient, hcombined is the combined convection and radiation coefficient, T is the temperature of the surrounding medium, Tsurr is the temperature of the surrounding surfaces, is the emissivity of the surface, and is the Stefan– Boltzman constant. The relations above can also be used for node M on the right boundary by replacing the subscript “0” by “M” and the subscript “1” by “M 1”. Note that absolute temperatures must be used in radiation heat transfer calculations, and all temperatures should be expressed in K or R when a boundary condition involves radiation to avoid mistakes. We usually try to avoid the radiation boundary condition even in numerical solutions since it causes the finite difference equations to be nonlinear, which are more difficult to solve. Treating Insulated Boundary Nodes as Interior Nodes: The Mirror Image Concept One way of obtaining the finite difference formulation of a node on an insulated boundary is to treat insulation as “zero” heat flux and to write an energy balance, as done in Eq. 5–23. Another and more practical way is to treat the node on an insulated boundary as an interior node. Conceptually this is done cen58933_ch05.qxd 9/4/2002 11:41 AM Page 277 277 CHAPTER 5 by replacing the insulation on the boundary by a mirror and considering the reflection of the medium as its extension (Fig. 5–17). This way the node next to the boundary node appears on both sides of the boundary node because of symmetry, converting it into an interior node. Then using the general formula (Eq. 5–18) for an interior node, which involves the sum of the temperatures of the adjoining nodes minus twice the node temperature, the finite difference formulation of a node m 0 on an insulated boundary of a plane wall can be expressed as Tm 1 2Tm x2 Tm 1 · gm k 0 → T1 2T0 x2 T1 · g0 k 0 EXAMPLE 5–1 01 x 2 Mirror Equivalent interior node Mirror image (5-30) which is equivalent to Eq. 5–23 obtained by the energy balance approach. The mirror image approach can also be used for problems that possess thermal symmetry by replacing the plane of symmetry by a mirror. Alternately, we can replace the plane of symmetry by insulation and consider only half of the medium in the solution. The solution in the other half of the medium is simply the mirror image of the solution obtained. Insulated boundary node Insulation x 2 1 01 x 2 FIGURE 5–17 A node on an insulated boundary can be treated as an interior node by replacing the insulation by a mirror. Steady Heat Conduction in a Large Uranium Plate Consider a large uranium plate of thickness L 4 cm and thermal conductivity k 28 W/m · °C in which heat is generated uniformly at a constant rate of · g 5 106 W/m3. One side of the plate is maintained at 0°C by iced water while the other side is subjected to convection to an environment at T 30°C with a heat transfer coefficient of h 45 W/m2 · °C, as shown in Figure 5–18. Considering a total of three equally spaced nodes in the medium, two at the boundaries and one at the middle, estimate the exposed surface temperature of the plate under steady conditions using the finite difference approach. SOLUTION A uranium plate is subjected to specified temperature on one side and convection on the other. The unknown surface temperature of the plate is to be determined numerically using three equally spaced nodes. Assumptions 1 Heat transfer through the wall is steady since there is no indication of any change with time. 2 Heat transfer is one-dimensional since the plate is large relative to its thickness. 3 Thermal conductivity is constant. 4 Radiation heat transfer is negligible. Properties The thermal conductivity is given to be k 28 W/m · °C. Analysis The number of nodes is specified to be M 3, and they are chosen to be at the two surfaces of the plate and the midpoint, as shown in the figure. Then the nodal spacing x becomes x L M 1 0.04 m 31 0.02 m We number the nodes 0, 1, and 2. The temperature at node 0 is given to be T0 0°C, and the temperatures at nodes 1 and 2 are to be determined. This problem involves only two unknown nodal temperatures, and thus we need to have only two equations to determine them uniquely. These equations are obtained by applying the finite difference method to nodes 1 and 2. Uranium plate 0°C 0 h T k = 28 W/m·°C · g = 5 × 106 W/m3 L 0 1 2 x FIGURE 5–18 Schematic for Example 5–1. cen58933_ch05.qxd 9/4/2002 11:41 AM Page 278 278 HEAT TRANSFER Node 1 is an interior node, and the finite difference formulation at that node is obtained directly from Eq. 5–18 by setting m 1: T0 2T1 x2 T2 · g1 k 0 → 0 2T1 x2 T2 · g1 k → 2T1 0 T2 · g1 x2 k (1) Node 2 is a boundary node subjected to convection, and the finite difference formulation at that node is obtained by writing an energy balance on the volume element of thickness x/2 at that boundary by assuming heat transfer to be into the medium at all sides: hA(T T2) kA T1 T2 x · g2(A x/2) 0 Canceling the heat transfer area A and rearranging give 1 T1 hx T2 k · g2 x2 2k hx T k (2) Equations (1) and (2) form a system of two equations in two unknowns T1 and T2. Substituting the given quantities and simplifying gives T1 2T1 T2 1.032T2 71.43 36.68 (in °C) (in °C) This is a system of two algebraic equations in two unknowns and can be solved easily by the elimination method. Solving the first equation for T1 and substituting into the second equation result in an equation in T2 whose solution is T2 136.1°C This is the temperature of the surface exposed to convection, which is the desired result. Substitution of this result into the first equation gives T1 103.8°C, which is the temperature at the middle of the plate. h T Plate 0 1 2 cm 2 x 2 cm Finite difference solution: T2 = 136.1°C Exact solution: T2 = 136.0°C FIGURE 5–19 Despite being approximate in nature, highly accurate results can be obtained by numerical methods. Discussion The purpose of this example is to demonstrate the use of the finite difference method with minimal calculations, and the accuracy of the result was not a major concern. But you might still be wondering how accurate the result obtained above is. After all, we used a mesh of only three nodes for the entire plate, which seems to be rather crude. This problem can be solved analytically as described in Chapter 2, and the analytical (exact) solution can be shown to be T(x) · 0.5ghL2/k hL · gL k Th x · gx2 2k Substituting the given quantities, the temperature of the exposed surface of the plate at x L 0.04 m is determined to be 136.0°C, which is almost identical to the result obtained here with the approximate finite difference method (Fig. 5–19). Therefore, highly accurate results can be obtained with numerical methods by using a limited number of nodes. cen58933_ch05.qxd 9/4/2002 11:41 AM Page 279 279 CHAPTER 5 EXAMPLE 5–2 Heat Transfer from Triangular Fins Consider an aluminum alloy fin (k 180 W/m · °C) of triangular cross section with length L 5 cm, base thickness b 1 cm, and very large width w in the direction normal to the plane of paper, as shown in Figure 5–20. The base of the fin is maintained at a temperature of T0 200°C. The fin is losing heat to the surrounding medium at T 25°C with a heat transfer coefficient of h 15 W/m2 · °C. Using the finite difference method with six equally spaced nodes along the fin in the x-direction, determine (a) the temperatures at the nodes, (b) the rate of heat transfer from the fin for w 1 m, and (c) the fin efficiency. h, T T0 b 0 SOLUTION A long triangular fin attached to a surface is considered. The nodal temperatures, the rate of heat transfer, and the fin efficiency are to be determined numerically using six equally spaced nodes. Triangular fin w 1 2 0.05 m 61 L M 1 0 → kAleft Tm Tm 1 kAright x all sides Tm 0.01 m Tm 1 x hAconv(T Tm) 0 Note that heat transfer areas are different for each node in this case, and using geometrical relations, they can be expressed as Aleft Aright Aconv (Height Width)@m (Height Width)@m 2 Length Width 1 2 1 2 2w[L (m 2w[L (m 2w( x/cos ) 1/2) x]tan 1/2) x]tan Substituting, 2kw[L 2kw[L 1) 2 (m (m x]tan 1) 2 Tm x]tan Tm 1 x Tm Tm 1 x h 5 2w x (T cos Tm) x L ∆x ——– cos θ [L – (m + 1 )∆ x]tan θ – 2 θ The temperature at node 0 is given to be T0 200°C, and the temperatures at the remaining five nodes are to be determined. Therefore, we need to have five equations to determine them uniquely. Nodes 1, 2, 3, and 4 are interior nodes, and the finite difference formulation for a general interior node m is obtained by applying an energy balance on the volume element of this node. Noting that heat transfer is steady and there is no heat generation in the fin and assuming heat transfer to be into the medium at all sides, the energy balance can be expressed as · Q 4 180 W/m · °C. Analysis (a) The number of nodes in the fin is specified to be M 6, and their location is as shown in the figure. Then the nodal spacing x becomes x 3 ∆x Assumptions 1 Heat transfer is steady since there is no indication of any change with time. 2 The temperature along the fin varies in the x direction only. 3 Thermal conductivity is constant. 4 Radiation heat transfer is negligible. Properties The thermal conductivity is given to be k b/ 2 tan θ = —– L θ 0 m–1 m m+1 ∆x L – (m – 1 )∆ x – 2 1 )∆ x]tan θ [L – (m – – 2 FIGURE 5–20 Schematic for Example 5–2 and the volume element of a general interior node of the fin. cen58933_ch05.qxd 9/4/2002 11:41 AM Page 280 280 HEAT TRANSFER Dividing each term by 2kwL tan / x gives 1 x (Tm L 1 ) 2 (m 1 Tm) 1 (m 1 ) 2 x (Tm L 1 h( x)2 (T kL sin Tm) Tm) 0 Note that b/2 L tan Also, sin 5.71° (5.5 0.5 cm 5 cm 0.1 → tan 10.1 5.71° 0.0995. Then the substitution of known quantities gives m)Tm 1 (10.00838 2m)Tm (4.5 m)Tm 1 0.209 Now substituting 1, 2, 3, and 4 for m results in these finite difference equations for the interior nodes: m m m m kAleft T4 5 ∆x —– 2 FIGURE 5–21 Schematic of the volume element of node 5 at the tip of a triangular fin. T5 x θ 4 ∆x —– 2 8.00838T1 3.5T2 900.209 3.5T1 6.00838T2 2.5T3 0.209 2.5T2 4.00838T3 1.5T4 0.209 1.5T3 2.00838T4 0.5T5 0.209 (1) (2) (3) (4) The finite difference equation for the boundary node 5 is obtained by writing an energy balance on the volume element of length x/2 at that boundary, again by assuming heat transfer to be into the medium at all sides (Fig. 5–21): ∆x/2 —–— cos θ ∆x —– tan θ 2 1: 2: 3: 4: hAconv (T T5) 0 where Aleft 2w x tan 2 and Aconv 2w x/2 cos Canceling w in all terms and substituting the known quantities gives T4 1.00838T5 0.209 (5) Equations (1) through (5) form a linear system of five algebraic equations in five unknowns. Solving them simultaneously using an equation solver gives T1 T4 198.6°C, 194.3°C, T2 T5 197.1°C, 192.9°C T3 195.7°C, which is the desired solution for the nodal temperatures. (b) The total rate of heat transfer from the fin is simply the sum of the heat transfer from each volume element to the ambient, and for w 1 m it is determined from cen58933_ch05.qxd 9/4/2002 11:41 AM Page 281 281 CHAPTER 5 5 · Q fin 5 · Q element, m m0 hAconv, m(Tm T) m0 Noting that the heat transfer surface area is w x/cos for the boundary nodes 0 and 5, and twice as large for the interior nodes 1, 2, 3, and 4, we have · Q fin h h wx [(T0 T ) 2(T1 T ) cos 2(T4 T ) (T5 T )] wx [T0 cos 2(T1 (15 W/m2 · °C) T2 T3 2(T2 T4) (1 m)(0.01 m) [200 cos 5.71° T5 2 T) 2(T3 T) 10T ] 785.7 192.9 10 25] 258.4 W (c) If the entire fin were at the base temperature of T0 of heat transfer from the fin for w 1 m would be · Q max 200°C, the total rate hAfin, total (T0 T ) h(2wL/cos )(T0 T ) (15 W/m2 · °C)[2(1 m)(0.05 m)/cos5.71°](200 263.8 W 25)°C Then the fin efficiency is determined from fin · Q fin · Q max 258.4 W 263.8 W 0.98 which is less than 1, as expected. We could also determine the fin efficiency in this case from the proper fin efficiency curve in Chapter 3, which is based on the analytical solution. We would read 0.98 for the fin efficiency, which is identical to the value determined above numerically. The finite difference formulation of steady heat conduction problems usually results in a system of N algebraic equations in N unknown nodal temperatures that need to be solved simultaneously. When N is small (such as 2 or 3), we can use the elementary elimination method to eliminate all unknowns except one and then solve for that unknown (see Example 5–1). The other unknowns are then determined by back substitution. When N is large, which is usually the case, the elimination method is not practical and we need to use a more systematic approach that can be adapted to computers. There are numerous systematic approaches available in the literature, and they are broadly classified as direct and iterative methods. The direct methods are based on a fixed number of well-defined steps that result in the solution in a systematic manner. The iterative methods, on the other hand, are based on an initial guess for the solution that is refined by iteration until a specified convergence criterion is satisfied (Fig. 5–22). The direct methods usually require a large amount of computer memory and computation time, Direct methods: Solve in a systematic manner following a series of well-defined steps. Iterative methods: Start with an initial guess for the solution, and iterate until solution converges. FIGURE 5–22 Two general categories of solution methods for solving systems of algebraic equations. cen58933_ch05.qxd 9/4/2002 11:41 AM Page 282 282 HEAT TRANSFER and they are more suitable for systems with a relatively small number of equations. The computer memory requirements for iterative methods are minimal, and thus they are usually preferred for large systems. The convergence of iterative methods to the desired solution, however, may pose a problem. 5–4 y … N ∆y ∆y Node (m, n) … n+1 n n–1 2 1 0 ∆x ∆x … 012 … m m–1 m+1 M x FIGURE 5–23 The nodal network for the finite difference formulation of twodimensional conduction in rectangular coordinates. m, n + 1 n+1 ∆y Volume element m – 1, n n · gm, n m, n m + 1, n Rate of heat conduction at the left, top, right, and bottom surfaces m, n – 1 ∆x y m–1 ∆x m TWO-DIMENSIONAL STEADY HEAT CONDUCTION In Section 5–3 we considered one-dimensional heat conduction and assumed heat conduction in other directions to be negligible. Many heat transfer problems encountered in practice can be approximated as being one-dimensional, but this is not always the case. Sometimes we need to consider heat transfer in other directions as well when the variation of temperature in other directions is significant. In this section we will consider the numerical formulation and solution of two-dimensional steady heat conduction in rectangular coordinates using the finite difference method. The approach presented below can be extended to three-dimensional cases. Consider a rectangular region in which heat conduction is significant in the x- and y-directions. Now divide the x-y plane of the region into a rectangular mesh of nodal points spaced x and y apart in the x- and y-directions, respectively, as shown in Figure 5–23, and consider a unit depth of z 1 in the z-direction. Our goal is to determine the temperatures at the nodes, and it is convenient to number the nodes and describe their position by the numbers instead of actual coordinates. A logical numbering scheme for two-dimensional problems is the double subscript notation (m, n) where m 0, 1, 2, . . . , M is the node count in the x-direction and n 0, 1, 2, . . . , N is the node count in the y-direction. The coordinates of the node (m, n) are simply x m x and y n y, and the temperature at the node (m, n) is denoted by Tm, n. Now consider a volume element of size x y 1 centered about a gen· eral interior node (m, n) in a region in which heat is generated at a rate of g and the thermal conductivity k is constant, as shown in Figure 5–24. Again assuming the direction of heat conduction to be toward the node under consideration at all surfaces, the energy balance on the volume element can be expressed as ∆y n–1 I Rate of heat generation inside the element Rate of change of the energy content of the element or m+1 x FIGURE 5–24 The volume element of a general interior node (m, n) for twodimensional conduction in rectangular coordinates. · Q cond, left · Q cond, top · Q cond, right · Q cond, bottom · Gelement Eelement t 0 (5-31) for the steady case. Again assuming the temperatures between the adjacent nodes to vary linearly and noting that the heat transfer area is y1 y in the x-direction and Ay x1 x in the y-direction, Ax the energy balance relation above becomes cen58933_ch05.qxd 9/4/2002 11:41 AM Page 283 283 CHAPTER 5 ky Tm Tm, n 1, n kx x Tm, n kx Dividing each term by x Tm 1, n 2Tm, n Tm Tm, n 1 y Tm, n ky Tm, n 1 y Tm Tm, n 1, n x · gm, n xy 0 (5-32) y and simplifying gives Tm, n 1, n 2Tm, n 1 2 Tm, n · gm, n k 1 2 x y 0 (5-33) for m 1, 2, 3, . . . , M 1 and n 1, 2, 3, . . . , N 1. This equation is identical to Eq. 5–12 obtained earlier by replacing the derivatives in the differential equation by differences for an interior node (m, n). Again a rectangular region M equally spaced nodes in the x-direction and N equally spaced nodes in the y-direction has a total of (M 1)(N 1) nodes, and Eq. 5–33 can be used to obtain the finite difference equations at all interior nodes. In finite difference analysis, usually a square mesh is used for simplicity (except when the magnitudes of temperature gradients in the x- and y-directions are very different), and thus x and y are taken to be the same. Then x y l, and the relation above simplifies to Tm 1, n Tm 1, n Tm, n 1 Tm, n 1 4Tm, n · gm, nl 2 k 0 (5-34) That is, the finite difference formulation of an interior node is obtained by adding the temperatures of the four nearest neighbors of the node, subtracting four times the temperature of the node itself, and adding the heat generation term. It can also be expressed in this form, which is easy to remember: Tleft Ttop Tright Tbottom 4Tnode · g nodel 2 k 0 (5-35) When there is no heat generation in the medium, the finite difference equation for an interior node further simplifies to Tnode (Tleft Ttop Tright Tbottom)/4, which has the interesting interpretation that the temperature of each interior node is the arithmetic average of the temperatures of the four neighboring nodes. This statement is also true for the three-dimensional problems except that the interior nodes in that case will have six neighboring nodes instead of four. Boundary Nodes The development of finite difference formulation of boundary nodes in two(or three-) dimensional problems is similar to the development in the onedimensional case discussed earlier. Again, the region is partitioned between the nodes by forming volume elements around the nodes, and an energy balance is written for each boundary node. Various boundary conditions can be handled as discussed for a plane wall, except that the volume elements in the two-dimensional case involve heat transfer in the y-direction as well as the x-direction. Insulated surfaces can still be viewed as “mirrors, ” and the cen58933_ch05.qxd 9/4/2002 11:41 AM Page 284 284 HEAT TRANSFER Volume element of node 2 h, T 1 · Qtop 2 3 · Qleft ∆y mirror image concept can be used to treat nodes on insulated boundaries as interior nodes. For heat transfer under steady conditions, the basic equation to keep in mind when writing an energy balance on a volume element is (Fig. 5–25) Boundary subjected to convection · Q · Qright · Qbottom · g2V2 · · · · Qleft + Qtop + Qright + Qbottom + —— = 0 k FIGURE 5–25 The finite difference formulation of a boundary node is obtained by writing an energy balance on its volume element. EXAMPLE 5–3 Convection h, T y 1 2 3 4 5 6 7 8 10 11 12 13 14 15 ∆ x = ∆y = l ∆y · 9 qR ∆y ∆x ∆x 90°C ∆x x ∆x h, T h, T 2 1 2 3 5 4 (a) Node 1 (5-36) Steady Two-Dimensional Heat Conduction in L-Bars Consider steady heat transfer in an L-shaped solid body whose cross section is given in Figure 5–26. Heat transfer in the direction normal to the plane of the paper is negligible, and thus heat transfer in the body is two-dimensional. The thermal conductivity of the body is k 15 W/m · °C, and heat is generated in · the body at a rate of g 2 106 W/m3. The left surface of the body is insulated, and the bottom surface is maintained at a uniform temperature of 90°C. The entire top surface is subjected to convection to ambient air at T 25°C with a convection coefficient of h 80 W/m2 · °C, and the right surface is sub· jected to heat flux at a uniform rate of q R 5000 W/m2. The nodal network of the problem consists of 15 equally spaced nodes with x y 1.2 cm, as shown in the figure. Five of the nodes are at the bottom surface, and thus their temperatures are known. Obtain the finite difference equations at the remaining nine nodes and determine the nodal temperatures by solving them. ∆x FIGURE 5–26 Schematic for Example 5–3 and the nodal network (the boundaries of volume elements of the nodes are indicated by dashed lines). 1 0 whether the problem is one-, two-, or three-dimensional. Again we assume, for convenience in formulation, all heat transfer to be into the volume element from all surfaces except for specified heat flux, whose direction is already specified. This is demonstrated in Example 5–3 for various boundary conditions. 4 ∆x · gVelement all sides (b) Node 2 FIGURE 5–27 Schematics for energy balances on the volume elements of nodes 1 and 2. SOLUTION Heat transfer in a long L-shaped solid bar with specified boundary conditions is considered. The nine unknown nodal temperatures are to be determined with the finite difference method. Assumptions 1 Heat transfer is steady and two-dimensional, as stated. 2 Thermal conductivity is constant. 3 Heat generation is uniform. 4 Radiation heat transfer is negligible. Properties The thermal conductivity is given to be k 15 W/m · °C. Analysis We observe that all nodes are boundary nodes except node 5, which is an interior node. Therefore, we will have to rely on energy balances to obtain the finite difference equations. But first we form the volume elements by partitioning the region among the nodes equitably by drawing dashed lines between the nodes. If we consider the volume element represented by an interior node to be full size (i.e., x y 1), then the element represented by a regular boundary node such as node 2 becomes half size (i.e., x y/2 1), and a corner node such as node 1 is quarter size (i.e., x/2 y/2 1). Keeping Eq. 5–36 in mind for the energy balance, the finite difference equations for each of the nine nodes are obtained as follows: (a) Node 1. The volume element of this corner node is insulated on the left and subjected to convection at the top and to conduction at the right and bottom surfaces. An energy balance on this element gives [Fig. 5–27a] cen58933_ch05.qxd 9/4/2002 11:41 AM Page 285 285 CHAPTER 5 0 h Taking x x (T 2 y T1) y T2 T1 2 x k · xy g1 22 x T4 T1 2 y k 0 l, it simplifies to –2 hl T1 k T2 hl T k T4 · g1l 2 2k (b) Node 2. The volume element of this boundary node is subjected to convection at the top and to conduction at the right, bottom, and left surfaces. An energy balance on this element gives [Fig. 5–27b] h x(T T2) Taking x k y y T3 T2 2 x kx T5 T2 y k y T1 T2 2 x y · g2 x 2 0 l, it simplifies to T1 4 2hl T2 k T3 · g2l 2 k 2hl T k 2T5 (c) Node 3. The volume element of this corner node is subjected to convection at the top and right surfaces and to conduction at the bottom and left surfaces. An energy balance on this element gives [Fig. 5–28a] h x 2 Taking x y (T 2 y T3) k x T6 T3 2 y k y T2 T3 2 x · xy g3 22 0 3 (5) 4 5 h, T l, it simplifies to T2 2hl T3 k 2 2hl T k T6 · g3l 2 2k (d ) Node 4. This node is on the insulated boundary and can be treated as an interior node by replacing the insulation by a mirror. This puts a reflected image of node 5 to the left of node 4. Noting that x y l, the general interior node relation for the steady two-dimensional case (Eq. 5–35) gives [Fig. 5–28b] T5 or, noting that T10 2 1 Mirror h, T T1 T5 T10 4T4 · g4l 2 k (a) Node 3 0 T2 3 2 4T4 2T5 90 · g4l 2 k 4 T6 T11 4T5 · g5l 2 k 5 6 h, T 7 6 5 (e) Node 5. This is an interior node, and noting that x y l, the finite difference formulation of this node is obtained directly from Eq. 5–35 to be [Fig. 5–29a] T4 (b) Node 4 FIGURE 5–28 Schematics for energy balances on the volume elements of nodes 3 and 4. 90° C, T1 10 6 0 12 11 (a) Node 5 (b) Node 6 FIGURE 5–29 Schematics for energy balances on the volume elements of nodes 5 and 6. cen58933_ch05.qxd 9/4/2002 11:41 AM Page 286 286 HEAT TRANSFER or, noting that T11 90°C, T2 T4 4T5 T6 90 · g5l 2 k (f ) Node 6. The volume element of this inner corner node is subjected to convection at the L-shaped exposed surface and to conduction at other surfaces. An energy balance on this element gives [Fig. 5–29b] y (T 2 x 2 h ky Taking x y T3 6 7 8 T5 T6 x l and noting that T12 2hl T6 k 6 2T5 9 h x(T · qR T7) k 13 k 90°C, it simplifies to T7 180 · 3g6l 2 2k 2hl T k (g) Node 7. The volume element of this boundary node is subjected to convection at the top and to conduction at the right, bottom, and left surfaces. An energy balance on this element gives [Fig. 5–30a] h, T h, T T12 T6 y T7 T6 kx 2 x y x T3 T6 ·3xy 0 k g6 2 4 y T6) 15 FIGURE 5–30 Schematics for energy balances on the volume elements of nodes 7 and 9. Taking x y k y T6 2 T13 T7 y T8 T7 kx 2 x y T7 y · g7 x 0 2 x l and noting that T13 T6 2hl T7 k 4 90°C, it simplifies to T8 180 · g7l 2 k 2hl T k (h) Node 8. This node is identical to Node 7, and the finite difference formulation of this node can be obtained from that of Node 7 by shifting the node numbers by 1 (i.e., replacing subscript m by m 1). It gives T7 2hl T8 k 4 T9 180 · g8l 2 k 2hl T k (i ) Node 9. The volume element of this corner node is subjected to convection at the top surface, to heat flux at the right surface, and to conduction at the bottom and left surfaces. An energy balance on this element gives [Fig. 5–30b] h x (T 2 Taking x T9) y ·y qR 2 k x T15 T9 2 y l and noting that T15 T8 2 hl T9 k 90 k y T8 T9 2 x · xy g9 22 90°C, it simplifies to · qRl k hl T k · g9l 2 2k 0 cen58933_ch05.qxd 9/4/2002 11:41 AM Page 287 287 CHAPTER 5 This completes the development of finite difference formulation for this problem. Substituting the given quantities, the system of nine equations for the determination of nine unknown nodal temperatures becomes T1 T3 –2.064T1 T2 T4 4.128T2 T3 2T5 T2 2.128T3 T6 T1 4T4 2T5 T2 T4 4T5 T6 2T5 6.128T6 T7 T6 4.128T7 T8 T7 4.128T8 T9 T8 2.064T9 11.2 22.4 12.8 109.2 109.2 212.0 202.4 202.4 105.2 which is a system of nine algebraic equations with nine unknowns. Using an equation solver, its solution is determined to be T1 T4 T7 112.1°C 109.4°C 97.3°C T2 T5 T8 110.8°C 108.1°C 96.3°C T3 T6 T9 106.6°C 103.2°C 97.6°C Note that the temperature is the highest at node 1 and the lowest at node 8. This is consistent with our expectations since node 1 is the farthest away from the bottom surface, which is maintained at 90°C and has one side insulated, and node 8 has the largest exposed area relative to its volume while being close to the surface at 90°C. Irregular Boundaries In problems with simple geometries, we can fill the entire region using simple volume elements such as strips for a plane wall and rectangular elements for two-dimensional conduction in a rectangular region. We can also use cylindrical or spherical shell elements to cover the cylindrical and spherical bodies entirely. However, many geometries encountered in practice such as turbine blades or engine blocks do not have simple shapes, and it is difficult to fill such geometries having irregular boundaries with simple volume elements. A practical way of dealing with such geometries is to replace the irregular geometry by a series of simple volume elements, as shown in Figure 5–31. This simple approach is often satisfactory for practical purposes, especially when the nodes are closely spaced near the boundary. More sophisticated approaches are available for handling irregular boundaries, and they are commonly incorporated into the commercial software packages. EXAMPLE 5–4 Actual boundary Approximation Heat Loss through Chimneys Hot combustion gases of a furnace are flowing through a square chimney made of concrete (k 1.4 W/m · °C). The flow section of the chimney is 20 cm 20 cm, and the thickness of the wall is 20 cm. The average temperature of the FIGURE 5–31 Approximating an irregular boundary with a rectangular mesh. cen58933_ch05.qxd 9/4/2002 11:41 AM Page 288 288 HEAT TRANSFER hot gases in the chimney is Ti 300°C, and the average convection heat transfer coefficient inside the chimney is hi 70 W/m2 · °C. The chimney is losing heat from its outer surface to the ambient air at To 20°C by convection with a heat transfer coefficient of ho 21 W/m2 · °C and to the sky by radiation. The emissivity of the outer surface of the wall is 0.9, and the effective sky temperature is estimated to be 260 K. Using the finite difference method with x y 10 cm and taking full advantage of symmetry, determine the temperatures at the nodal points of a cross section and the rate of heat loss for a 1-m-long section of the chimney. SOLUTION Heat transfer through a square chimney is considered. The nodal temperatures and the rate of heat loss per unit length are to be determined with the finite difference method. Assumptions 1 Heat transfer is steady since there is no indication of change with time. 2 Heat transfer through the chimney is two-dimensional since the height of the chimney is large relative to its cross section, and thus heat conduction through the chimney in the axial direction is negligible. It is tempting to simplify the problem further by considering heat transfer in each wall to be one-dimensional, which would be the case if the walls were thin and thus the corner effects were negligible. This assumption cannot be justified in this case since the walls are very thick and the corner sections constitute a considerable portion of the chimney structure. 3 Thermal conductivity is constant. Symmetry lines (Equivalent to insulation) Properties 0.9. h1 T1 1 2 3 4 5 6 7 8 h0, T0 Tsky 9 Representative section of chimney FIGURE 5–32 Schematic of the chimney discussed in Example 5–4 and the nodal network for a representative section. h, T 1 h, T 2 1 2 (a) Node 1 (a) Node 1. On the inner boundary, subjected to convection, Figure 5–33a hi x (T 2i T1) k y T2 T1 2 x k x T3 T1 2 y 4 (b) Node 2 FIGURE 5–33 Schematics for energy balances on the volume elements of nodes 1 and 2. 1.4 W/m · °C and Analysis The cross section of the chimney is given in Figure 5–32. The most striking aspect of this problem is the apparent symmetry about the horizontal and vertical lines passing through the midpoint of the chimney as well as the diagonal axes, as indicated on the figure. Therefore, we need to consider only one-eighth of the geometry in the solution whose nodal network consists of nine equally spaced nodes. No heat can cross a symmetry line, and thus symmetry lines can be treated as insulated surfaces and thus “mirrors” in the finite difference formulation. Then the nodes in the middle of the symmetry lines can be treated as interior nodes by using mirror images. Six of the nodes are boundary nodes, so we will have to write energy balances to obtain their finite difference formulations. First we partition the region among the nodes equitably by drawing dashed lines between the nodes through the middle. Then the region around a node surrounded by the boundary or the dashed lines represents the volume element of the node. Considering a unit depth and using the energy balance approach for the boundary nodes (again assuming all heat transfer into the volume element for convenience) and the formula for the interior nodes, the finite difference equations for the nine nodes are determined as follows: 0 3 The properties of chimney are given to be k Taking x y l, it simplifies to –2 hi l T1 k T2 T3 hi l T ki 0 0 cen58933_ch05.qxd 9/4/2002 11:42 AM Page 289 289 CHAPTER 5 (b) Node 2. On the inner boundary, subjected to convection, Figure 5–33b k Taking x y T1 T2 2 x y x (T 2i hi T2) 0 kx T2 (4) hi l T ki 2T4 3 4 5 (8) 7 8 9 Mirror T1 T2 T4 T4 T5 T8 T6 T7 T8 4T3 4T4 4T5 Mirror FIGURE 5–34 Converting the boundary nodes 3 and 5 on symmetry lines to interior nodes by using mirror images. (c) Nodes 3, 4, and 5. (Interior nodes, Fig. 5–34) Node 3: T4 Node 4: T3 Node 5: T4 Mirror images Mirror image hi l T2 k 3 (4) 0 y l, it simplifies to T1 2 6 T4 1 0 0 0 (d ) Node 6. (On the outer boundary, subjected to convection and radiation) 0 ho Taking x y T2 x T3 T 6 2 y k x (T 2o k y T7 T6 2 x x4 (T 2 sky T6) T64) 0 l, it simplifies to 2 T3 ho l T6 k ho l T ko l k 4 (Tsky T64) (e) Node 7. (On the outer boundary, subjected to convection and radiation, Fig. 5–35) k Taking x 2T4 y T6 y T6 T7 T4 T 7 kx 2 x y ho x(To T7) y T8 T7 2 x 4 x(Tsky T74) 4 k 0 6 l, it simplifies to 4 2ho l T7 k 2ho l T ko T8 2 l k 4 (Tsky T74) (f ) Node 8. Same as Node 7, except shift the node numbers up by 1 (replace 4 by 5, 6 by 7, 7 by 8, and 8 by 9 in the last relation) 2T5 T7 4 2ho l T8 k 2ho l T ko T9 2 l k 4 (Tsky T84) (g) Node 9. (On the outer boundary, subjected to convection and radiation, Fig. 5–35) k y T8 T9 2 x 0 ho Insulation x (T 2o T9) x4 (T 2 sky T94) 0 7 8 9 h, T Tsky FIGURE 5–35 Schematics for energy balances on the volume elements of nodes 7 and 9. cen58933_ch05.qxd 9/4/2002 11:42 AM Page 290 290 HEAT TRANSFER Taking x y l, it simplifies to T8 1 ho l T9 k ho l T ko l k 4 (Tsky T94) This problem involves radiation, which requires the use of absolute temperature, and thus all temperatures should be expressed in Kelvin. Alternately, we could use °C for all temperatures provided that the four temperatures in the radiation terms are expressed in the form (T 273)4. Substituting the given quantities, the system of nine equations for the determination of nine unknown nodal temperatures in a form suitable for use with the Gauss-Seidel iteration method becomes T1 T2 T3 T4 T5 T6 T7 T8 T9 (T2 (T1 (T1 (T2 (2T4 (T2 (2T4 (2T5 (T8 T3 2865)/7 2T4 2865)/8 2T4 T6)/4 T3 T5 T7)/4 2T8)/4 T3 456.2 0.3645 10 9 T64)/3.5 T6 T8 912.4 0.729 10 9 T74)/7 T7 T9 912.4 0.729 10 9 T84)/7 456.2 0.3645 10 9 T94)/2.5 which is a system of nonlinear equations. Using an equation solver, its solution is determined to be 23 40 55 60 55 40 T1 T4 T7 40 Temperature, °C 55 60 55 40 89 152 89 138 152 256 273 273 256 138 273 152 138 256 273 256 138 89 138 152 138 89 545.7 K 411.2 K 328.1 K 272.6°C 138.0°C 54.9°C T2 T5 T8 529.2 K 362.1 K 313.1 K 40 55 60 55 FIGURE 5–36 The variation of temperature in the chimney. 40 T3 T6 T9 425.2 K 332.9 K 296.5 K 152.1°C 59.7°C 23.4°C 23 40 55 The variation of temperature in the chimney is shown in Figure 5–36. Note that the temperatures are highest at the inner wall (but less than 300°C) and lowest at the outer wall (but more that 260 K), as expected. The average temperature at the outer surface of the chimney weighed by the surface area is 60 Twall, out 55 40 0.5 332.9 (0.5T6 T7 (0.5 1 328.1 313.1 3 23 256.1°C 89.0°C 39.9°C T8 0.5T9) 1 0.5) 0.5 296.5 318.6 K 23 Then the rate of heat loss through the 1-m-long section of the chimney can be determined approximately from cen58933_ch05.qxd 9/4/2002 11:42 AM Page 291 291 CHAPTER 5 · Q chimney 4 4 ho Ao (Twall, out To) Ao (Twall, out Tsky) 2 (21 W/m · K)[4 (0.6 m)(1 m)](318.6 293)K 0.9(5.67 10 8 W/m2 · K4) [4 (0.6 m)(1 m)](318.6 K)4 (260 K)4] 1291 702 1993 W We could also determine the heat transfer by finding the average temperature of the inner wall, which is (272.6 256.1)/2 264.4°C, and applying Newton’s law of cooling at that surface: · Q chimney hi Ai (Ti Twall, in) (70 W/m2 · K)[4 (0.2 m)(1 m)](300 264.4)°C 1994 W The difference between the two results is due to the approximate nature of the numerical analysis. Discussion We used a relatively crude numerical model to solve this problem to keep the complexities at a manageable level. The accuracy of the solution obtained can be improved by using a finer mesh and thus a greater number of nodes. Also, when radiation is involved, it is more accurate (but more laborious) to determine the heat losses for each node and add them up instead of using the average temperature. 5–5 I TRANSIENT HEAT CONDUCTION So far in this chapter we have applied the finite difference method to steady heat transfer problems. In this section we extend the method to solve transient problems. We applied the finite difference method to steady problems by discretizing the problem in the space variables and solving for temperatures at discrete points called the nodes. The solution obtained is valid for any time since under steady conditions the temperatures do not change with time. In transient problems, however, the temperatures change with time as well as position, and thus the finite difference solution of transient problems requires discretization in time in addition to discretization in space, as shown in Figure 5–37. This is done by selecting a suitable time step t and solving for the unknown nodal temperatures repeatedly for each t until the solution at the desired time is obtained. For example, consider a hot metal object that is taken out of the oven at an initial temperature of Ti at time t 0 and is allowed to cool in ambient air. If a time step of t 5 min is chosen, the determination of the temperature distribution in the metal piece after 3 h requires the determination of the temperatures 3 60/5 36 times, or in 36 time steps. Therefore, the computation time of this problem will be 36 times that of a steady problem. Choosing a smaller t will increase the accuracy of the solution, but it will also increase the computation time. In transient problems, the superscript i is used as the index or counter of time steps, with i 0 corresponding to the specified initial condition. In the case of the hot metal piece discussed above, i 1 corresponds to t1 t 5 min, i 2 corresponds to t 2 t 10 min, and a general t i i i Tm +1 Tm +1 Tm +1 –1 +1 i+1 ∆t i i Tm – 1 Tm i 1 0 ∆t ∆x 0 1 ∆x m–1 i Tm +1 ∆x m m+1 x FIGURE 5–37 Finite difference formulation of timedependent problems involves discrete points in time as well as space. cen58933_ch05.qxd 9/4/2002 11:42 AM Page 292 292 HEAT TRANSFER i time step i corresponds to ti i t. The notation Tm is used to represent the temperature at the node m at time step i. The formulation of transient heat conduction problems differs from that of steady ones in that the transient problems involve an additional term representing the change in the energy content of the medium with time. This additional term appears as a first derivative of temperature with respect to time in the differential equation, and as a change in the internal energy content during t in the energy balance formulation. The nodes and the volume elements in transient problems are selected as they are in the steady case, and, again assuming all heat transfer is into the element for convenience, the energy balance on a volume element during a time interval t can be expressed as Heat transferred into the volume element from all of its surfaces during t Heat generated within the volume element during t The change in the energy content of the volume element during t or · Q t t · G element Eelement (5-37) All sides · where the rate of heat transfer Q normally consists of conduction terms for interior nodes, but may involve convection, heat flux, and radiation for boundary nodes. Velement C T, where is density and C is Noting that Eelement mC T the specific heat of the element, dividing the earlier relation by t gives · Q · G element All sides Eelement t Velement C T t (5-38) or, for any node m in the medium and its volume element, · Q All sides Volume element (can be any shape) Node m ρ = density V = volume ρV = mass C = specific heat ∆ T = temperature change i i ∆U = ρVC∆T = ρVC(Tm+ 1 – Tm ) FIGURE 5–38 The change in the energy content of the volume element of a node during a time interval t. · G element Velement C i Tm 1 i Tm t (5-39) i i where Tm and Tm 1 are the temperatures of node m at times ti i t and ti 1 i i Tm represents the temperature change (i 1) t, respectively, and Tm 1 of the node during the time interval t between the time steps i and i 1 (Fig. 5–38). i i Note that the ratio (Tm 1 Tm)/ t is simply the finite difference approximation of the partial derivative T/ t that appears in the differential equations of transient problems. Therefore, we would obtain the same result for the finite difference formulation if we followed a strict mathematical approach instead of the energy balance approach used above. Also note that the finite difference formulations of steady and transient problems differ by the single term on the right side of the equal sign, and the format of that term remains the same in all coordinate systems regardless of whether heat transfer is one-, i i two-, or three-dimensional. For the special case of Tm 1 Tm (i.e., no change in temperature with time), the formulation reduces to that of steady case, as expected. cen58933_ch05.qxd 9/4/2002 11:42 AM Page 293 293 CHAPTER 5 The nodal temperatures in transient problems normally change during each time step, and you may be wondering whether to use temperatures at the previous time step i or the new time step i 1 for the terms on the left side of Eq. 5–39. Well, both are reasonable approaches and both are used in practice. The finite difference approach is called the explicit method in the first case and the implicit method in the second case, and they are expressed in the general form as (Fig. 5–39) · Qi Explicit method: · G ielement Velement C i Tm 1 i Tm t All sides · Qi Implicit method: 1 ·i 1 G element Velement C i Tm 1 i Tm t All sides If expressed at i + 1: Implicit method · · ∑ Q + Gelement All sides i i Tm+ 1 – Tm = ρVelementC ————– ∆t If expressed at i: Explicit method FIGURE 5–39 The formulation of explicit and implicit methods differs at the time step (previous or new) at which the heat transfer and heat generation terms are expressed. (5-40) (5-41) It appears that the time derivative is expressed in forward difference form in the explicit case and backward difference form in the implicit case. Of course, it is also possible to mix the two fundamental formulations of Eqs. 5–40 and 5–41 and come up with more elaborate formulations, but such formulations offer little insight and are beyond the scope of this text. Note that both formulations are simply expressions between the nodal temperatures before and i after a time interval and are based on determining the new temperatures Tm 1 i using the previous temperatures Tm. The explicit and implicit formulations given here are quite general and can be used in any coordinate system regardless of the dimension of heat transfer. The volume elements in multidimensional cases simply have more surfaces and thus involve more terms in the summation. The explicit and implicit methods have their advantages and disadvantages, and one method is not necessarily better than the other one. Next you will see that the explicit method is easy to implement but imposes a limit on the allowable time step to avoid instabilities in the solution, and the implicit method requires the nodal temperatures to be solved simultaneously for each time step but imposes no limit on the magnitude of the time step. We will limit the discussion to one- and two-dimensional cases to keep the complexities at a manageable level, but the analysis can readily be extended to three-dimensional cases and other coordinate systems. A Plane wall · gm Volume element of node m i Tm+ 1 i Tm kA Transient Heat Conduction in a Plane Wall Consider transient one-dimensional heat conduction in a plane wall of thick· ness L with heat generation g(x, t) that may vary with time and position and constant conductivity k with a mesh size of x L/M and nodes 0, 1, 2, . . . , M in the x-direction, as shown in Figure 5–40. Noting that the volume element of a general interior node m involves heat conduction from two sides and the volume of the element is Velement A x, the transient finite difference formulation for an interior node can be expressed on the basis of Eq. 5–39 as kA Tm Tm 1 x kA Tm Tm 1 x · gm A x A xC i Tm 1 i Tm t (5-42) i i Tm – 1 – Tm ————– i i Tm + 1 – Tm kA ————– ∆x ∆x ∆x 01 2 ∆x m –1 m m +1 ∆x M –1 Mx FIGURE 5–40 The nodal points and volume elements for the transient finite difference formulation of one-dimensional conduction in a plane wall. cen58933_ch05.qxd 9/4/2002 11:42 AM Page 294 294 HEAT TRANSFER Canceling the surface area A and multiplying by x/k, it simplifies to Tm 2Tm 1 Tm · gm x2 k 1 x2 i (T tm 1 i Tm) (5-43) where k/ C is the thermal diffusivity of the wall material. We now define a dimensionless mesh Fourier number as αt x2 (5-44) Then Eq. 5–43 reduces to Tm 1 2Tm Tm · gm x2 k 1 i Tm 1 i Tm (5-45) Note that the left side of this equation is simply the finite difference formulation of the problem for the steady case. This is not surprising since the formui i Tm. Also, we are still not lation must reduce to the steady case for Tm 1 committed to explicit or implicit formulation since we did not indicate the time step on the left side of the equation. We now obtain the explicit finite difference formulation by expressing the left side at time step i as i Tm i 2Tm 1 i Tm ·i gm x2 k 1 i Tm 1 i Tm (explicit) i This equation can be solved explicitly for the new temperature Tm the name explicit method) to give i Tm 1 i (Tm i Tm 1) 1 1 ·i gm x2 k i 2 ) Tm (1 (5-46) (and thus (5-47) for all interior nodes m 1, 2, 3, . . . , M 1 in a plane wall. Expressing the left side of Eq. 5–45 at time step i 1 instead of i would give the implicit finite difference formulation as ∆– x — 2 A i hA(T – T0 ) ∆– x ρA — C 2 i Tm T0i + 1 – T0i ————– i 2Tm 1 i Tm ·i gm 1 1 1 x2 i Tm 1 i Tm k (implicit) (5-48) ∆t which can be rearranged as · g0 i T1 – T0i kA ——— ∆x i Tm ∆x ∆x 0 1 1 1 2 … FIGURE 5–41 Schematic for the explicit finite difference formulation of the convection condition at the left boundary of a plane wall. L x 1 1 (1 i 2 ) Tm 1 i Tm 1 1 ·i gm 1 k x2 i Tm 0 (5-49) The application of either the explicit or the implicit formulation to each of the M 1 interior nodes gives M 1 equations. The remaining two equations are obtained by applying the same method to the two boundary nodes unless, of course, the boundary temperatures are specified as constants (invariant with time). For example, the formulation of the convection boundary condition at the left boundary (node 0) for the explicit case can be expressed as (Fig. 5–41) hA(T T0i) kA T1i T0i x x · gi0 A 2 i1 A x T0 C 2 T0i t (5-50) cen58933_ch05.qxd 9/4/2002 11:42 AM Page 295 295 CHAPTER 5 which simplifies to T0i 1 1 2 2 hx i T0 k 2 T1i 2 hx T k · g i0 x 2 k (5-51) Note that in the case of no heat generation and 0.5, the explicit i finite difference formulation for a general interior node reduces to Tm 1 i i (Tm 1 Tm 1)/2, which has the interesting interpretation that the temperature of an interior node at the new time step is simply the average of the temperatures of its neighboring nodes at the previous time step. Once the formulation (explicit or implicit) is complete and the initial condition is specified, the solution of a transient problem is obtained by marching in time using a step size of t as follows: select a suitable time step t and determine the nodal temperatures from the initial condition. Taking the initial i i temperatures as the previous solution Tm at t 0, obtain the new solution Tm 1 at all nodes at time t t using the transient finite difference relations. Now i using the solution just obtained at t t as the previous solution Tm, obtain i1 the new solution Tm at t 2 t using the same relations. Repeat the process until the solution at the desired time is obtained. Stability Criterion for Explicit Method: Limitation on t The explicit method is easy to use, but it suffers from an undesirable feature that severely restricts its utility: the explicit method is not unconditionally stable, and the largest permissible value of the time step t is limited by the stability criterion. If the time step t is not sufficiently small, the solutions obtained by the explicit method may oscillate wildly and diverge from the actual solution. To avoid such divergent oscillations in nodal temperatures, the value of t must be maintained below a certain upper limit established by the stability criterion. It can be shown mathematically or by a physical argument based on the second law of thermodynamics that the stability criterion is sati i isfied if the coefficients of all Tm in the Tm 1 expressions (called the primary coefficients) are greater than or equal to zero for all nodes m (Fig. 5–42). Of i course, all the terms involving Tm for a particular node must be grouped together before this criterion is applied. Different equations for different nodes may result in different restrictions on the size of the time step t, and the criterion that is most restrictive should be used in the solution of the problem. A practical approach is to identify the equation with the smallest primary coefficient since it is the most restrictive and to determine the allowable values of t by applying the stability criterion to that equation only. A t value obtained this way will also satisfy the stability criterion for all other equations in the system. For example, in the case of transient one-dimensional heat conduction in a plane wall with specified surface temperatures, the explicit finite difference equations for all the nodes (which are interior nodes) are obtained from i i Eq. 5–47. The coefficient of Tm in the Tm 1 expression is 1 2 , which is independent of the node number m, and thus the stability criterion for all nodes in this case is 1 2 0 or t x2 1 2 interior nodes, one-dimensional heat transfer in rectangular coordinates (5-52) Explicit formulation: T0i 1 a0T0i ··· T1i 1 a1T1i ··· M i Tm 1 i amTm ··· M i TM 1 i aMTM ··· Stability criterion: am 0, m 0, 1, 2, . . . m, . . . M FIGURE 5–42 The stability criterion of the explicit method requires all primary coefficients to be positive or zero. cen58933_ch05.qxd 9/4/2002 11:42 AM Page 296 296 HEAT TRANSFER When the material of the medium and thus its thermal diffusivity is known and the value of the mesh size x is specified, the largest allowable value of the time step t can be determined from this relation. For example, in the case of a brick wall ( 0.45 10 6 m2/s) with a mesh size of x 0.01 m, the upper limit of the time step is t (0.01 m)2 2(0.45 10 6 m2/s) 1 x2 2 111s 1.85 min The boundary nodes involving convection and/or radiation are more restrictive than the interior nodes and thus require smaller time steps. Therefore, the most restrictive boundary node should be used in the determination of the maximum allowable time step t when a transient problem is solved with the explicit method. To gain a better understanding of the stability criterion, consider the explicit finite difference formulation for an interior node of a plane wall (Eq. 5–47) for the case of no heat generation, i Tm 80°C 50°C 50°C 20°C m–1 m m–1 m+1 m m+1 Time step: i + 1 Time step: i FIGURE 5–43 The violation of the stability criterion in the explicit method may result in the violation of the second law of thermodynamics and thus divergence of solution. 0°C k = 28 W/m·°C · g = 5 × 106 W/m3 h T α = 12.5 × 10 – 6 m2/s ∆x ∆x 0 0 1 L 2 Tinitial = 200°C FIGURE 5–44 Schematic for Example 5–5. x i (Tm 1 i Tm 1) (1 i 2 )Tm i i Assume that at some time step i the temperatures Tm 1 and Tm 1 are equal but i i i i less than Tm (say, Tm 1 Tm 1 50°C and Tm 80°C). At the next time step, we expect the temperature of node m to be between the two values (say, 70°C). However, if the value of exceeds 0.5 (say, 1), the temperature of node m at the next time step will be less than the temperature of the neighboring nodes (it will be 20°C), which is physically impossible and violates the second law of thermodynamics (Fig. 5–43). Requiring the new temperature of node m to remain above the temperature of the neighboring nodes is equivalent to requiring the value of to remain below 0.5. The implicit method is unconditionally stable, and thus we can use any time step we please with that method (of course, the smaller the time step, the better the accuracy of the solution). The disadvantage of the implicit method is that it results in a set of equations that must be solved simultaneously for each time step. Both methods are used in practice. EXAMPLE 5–5 Uranium plate 1 Transient Heat Conduction in a Large Uranium Plate Consider a large uranium plate of thickness L 4 cm, thermal conductivity k 28 W/m · °C, and thermal diffusivity 12.5 10 6 m2/s that is initially at a uniform temperature of 200°C. Heat is generated uniformly in the plate at a · constant rate of g 5 106 W/m3. At time t 0, one side of the plate is brought into contact with iced water and is maintained at 0°C at all times, while the other side is subjected to convection to an environment at T 30°C with a heat transfer coefficient of h 45 W/m2 · °C, as shown in Figure 5–44. Considering a total of three equally spaced nodes in the medium, two at the boundaries and one at the middle, estimate the exposed surface temperature of the plate 2.5 min after the start of cooling using (a) the explicit method and (b) the implicit method. cen58933_ch05.qxd 9/4/2002 11:42 AM Page 297 297 CHAPTER 5 SOLUTION We have solved this problem in Example 5–1 for the steady case, and here we repeat it for the transient case to demonstrate the application of the transient finite difference methods. Again we assume one-dimensional heat transfer in rectangular coordinates and constant thermal conductivity. The number of nodes is specified to be M 3, and they are chosen to be at the two surfaces of the plate and at the middle, as shown in the figure. Then the nodal spacing x becomes 0.04 m 31 L x M 1 0.02 m We number the nodes as 0, 1, and 2. The temperature at node 0 is given to be T0 0°C at all times, and the temperatures at nodes 1 and 2 are to be determined. This problem involves only two unknown nodal temperatures, and thus we need to have only two equations to determine them uniquely. These equations are obtained by applying the finite difference method to nodes 1 and 2. (a) Node 1 is an interior node, and the explicit finite difference formulation at that node is obtained directly from Eq. 5–47 by setting m 1: T1i 1 T2i) (T0 · g1 x2 k 2 ) T1i (1 (1) Node 2 is a boundary node subjected to convection, and the finite difference formulation at that node is obtained by writing an energy balance on the volume element of thickness x/2 at that boundary by assuming heat transfer to be into the medium at all sides (Fig. 5–45): T2i) hA(T kA T1i T2i x i1 x · g2 A 2 A x T2 C 2 A T2i i i T1 – T2 kA ——— ∆x x Dividing by kA/2 x and using the definitions of thermal diffusivity the dimensionless mesh Fourier number t/( x)2 gives k/ C and T2i) i which can be solved for T2 T2i 1 1 2 1 2 2(T1i T2i) · g2 x2 k T2i 1 T2i to give hx i T2 k 2T1i 2 · g2 x2 k hx T k (2) Note that we did not use the superscript i for quantities that do not change with time. Next we need to determine the upper limit of the time step t from the i stability criterion, which requires the coefficient of T1 in Equation 1 and the coi efficient of T2 in the second equation to be greater than or equal to zero. The i coefficient of T2 is smaller in this case, and thus the stability criterion for this problem can be expressed as 1 2 2 hx k 0 → 2(1 1 → h x/k) t 2 (1 x2 h x/k) i hA(T – T2 ) i T2 + 1 i T2 0 2h x (T k · g2 Volume element of node 2 0 1 ∆– x — 2 2 x FIGURE 5–45 Schematic for the explicit finite difference formulation of the convection condition at the right boundary of a plane wall. cen58933_ch05.qxd 9/4/2002 11:42 AM Page 298 298 HEAT TRANSFER since t/( x)2. Substituting the given quantities, the maximum allowable value of the time step is determined to be t 2(12.5 10 6 m2/s)[1 (0.02 m)2 (45 W/m2 · °C)(0.02 m)/28 W/m · °C] 15.5 s Therefore, any time step less than 15.5 s can be used to solve this problem. For convenience, let us choose the time step to be t 15 s. Then the mesh Fourier number becomes αt ( x)2 (12.5 10 6 m2/s)(15 s) (0.02 m)2 0.46875 (for t 15 s) Substituting this value of and other given quantities, the explicit finite difference equations (1) and (2) developed here reduce to T1i T2i 1 1 0.0625T1i 0.9375T1i 0.46875T2i 33.482 0.032366T2i 34.386 The initial temperature of the medium at t 0 and i 0 is given to be 200°C 0 0 1 throughout, and thus T1 200°C. Then the nodal temperatures at T1 T2 1 and T2 at t t 15 s are determined from these equations to be T11 T21 2 2 Similarly, the nodal temperatures T1 and T2 at t determined to be TABLE 5–2 The variation of the nodal temperatures in Example 5–5 with time obtained by the explicit method Time Step, i 0 1 2 3 4 5 6 7 8 9 10 20 30 40 Time, s 0 15 30 45 60 75 90 105 120 135 150 300 450 600 0.0625T10 0.46875T20 33.482 0.0625 200 0.46875 200 33.482 139.7°C 0.9375T10 0.032366T20 34.386 0.9375 200 0.032366 200 34.386 228.4°C Node Temperature, °C i T1 200.0 139.7 149.3 123.8 125.6 114.6 114.3 109.5 108.9 106.7 106.3 103.8 103.7 103.7 i T2 200.0 228.4 172.8 179.9 156.3 157.1 146.9 146.3 141.8 141.1 139.0 136.1 136.0 136.0 T12 T22 2t 2 15 30 s are 0.0625T11 0.46875T21 33.482 0.0625 139.7 0.46875 228.4 33.482 149.3°C 0.9375T11 0.032366T21 34.386 0.9375 139.7 0.032366 228.4 34.386 172.8°C Continuing in the same manner, the temperatures at nodes 1 and 2 are determined for i 1, 2, 3, 4, 5, . . . , 50 and are given in Table 5–2. Therefore, the temperature at the exposed boundary surface 2.5 min after the start of cooling is 2.5 TL min T210 139.0°C (b) Node 1 is an interior node, and the implicit finite difference formulation at that node is obtained directly from Eq. 5–49 by setting m 1: T0 (1 2 ) T1i 1 T2i 1 · g0 x2 k T1i 0 (3) Node 2 is a boundary node subjected to convection, and the implicit finite difference formulation at that node can be obtained from this formulation by expressing the left side of the equation at time step i 1 instead of i as cen58933_ch05.qxd 9/4/2002 11:42 AM Page 299 299 CHAPTER 5 2h x (T k T2i 1) 2(T1i 1 T2i 1) hx i T2 k 1 2 · g2 x2 k T2i 1 T2i which can be rearranged as 2 T1i 1 1 2 2 · g2 x2 k hx T k T2i 0 (4) Again we did not use the superscript i or i 1 for quantities that do not change with time. The implicit method imposes no limit on the time step, and thus we can choose any value we want. However, we will again choose t 15 s, and thus 0.46875, to make a comparison with part (a) possible. Substituting this value of and other given quantities, the two implicit finite difference equations developed here reduce to –1.9375T1i 1 0.46875T2i 0.9375T1i 1 1.9676T2i 0 Again T1 and for i 1 1 T1i T2i 33.482 34.386 0 0 0 T2 200°C at t 0 and i 0 because of the initial condition, 0, these two equations reduce to 1.9375T11 0.46875T21 0.9375T11 1.9676T21 200 200 33.482 34.386 1 1 The unknown nodal temperatures T1 and T2 at t solving these two equations simultaneously to be T11 Similarly, for i 168.8°C T21 and t 0 0 15 s are determined by 199.6°C TABLE 5–3 1, these equations reduce to 1.9375T12 0.46875T22 0.9375T12 1.9676T22 168.8 199.6 33.482 34.386 0 0 2 2 The unknown nodal temperatures T1 and T2 at t t 2 15 determined by solving these two equations simultaneously to be T12 150.5°C and T22 30 s are 190.6°C Continuing in this manner, the temperatures at nodes 1 and 2 are determined for i 2, 3, 4, 5, . . . , 40 and are listed in Table 5–3, and the temperature at the exposed boundary surface (node 2) 2.5 min after the start of cooling is obtained to be 2.5 TL min T210 143.9°C which is close to the result obtained by the explicit method. Note that either method could be used to obtain satisfactory results to transient problems, except, perhaps, for the first few time steps. The implicit method is preferred when it is desirable to use large time steps, and the explicit method is preferred when one wishes to avoid the simultaneous solution of a system of algebraic equations. The variation of the nodal temperatures in Example 5–5 with time obtained by the implicit method Time Step, i Time, s 0 1 2 3 4 5 6 7 8 9 10 20 30 40 0 15 30 45 60 75 90 105 120 135 150 300 450 600 Node Temperature, °C i T1 i T2 200.0 168.8 150.5 138.6 130.3 124.1 119.5 115.9 113.2 111.0 109.4 104.2 103.8 103.8 200.0 199.6 190.6 180.4 171.2 163.6 157.6 152.8 149.0 146.1 143.9 136.7 136.1 136.1 cen58933_ch05.qxd 9/4/2002 11:42 AM Page 300 300 HEAT TRANSFER EXAMPLE 5–6 South Warm air Sun’s rays Trombe wall Heat loss Heat gain Vent Cool air Glazing FIGURE 5–46 Schematic of a Trombe wall (Example 5–6). TABLE 5–4 The hourly variation of monthly average ambient temperature and solar heat flux incident on a vertical surface for January in Reno, Nevada Time of Day 7 10 1 4 7 10 1 4 Ambient Solar Temperature, Radiation, °F Btu/h · ft2 AM–10 AM AM–1 PM PM–4 PM PM–7 PM PM–10 PM PM–1 AM AM–4 AM AM–7 AM 33 43 45 37 32 27 26 25 114 242 178 0 0 0 0 0 Solar Energy Storage in Trombe Walls Dark painted thick masonry walls called Trombe walls are commonly used on south sides of passive solar homes to absorb solar energy, store it during the day, and release it to the house during the night (Fig. 5–46). The idea was proposed by E. L. Morse of Massachusetts in 1881 and is named after Professor Felix Trombe of France, who used it extensively in his designs in the 1970s. Usually a single or double layer of glazing is placed outside the wall and transmits most of the solar energy while blocking heat losses from the exposed surface of the wall to the outside. Also, air vents are commonly installed at the bottom and top of the Trombe walls so that the house air enters the parallel flow channel between the Trombe wall and the glazing, rises as it is heated, and enters the room through the top vent. Consider a house in Reno, Nevada, whose south wall consists of a 1-ft-thick Trombe wall whose thermal conductivity is k 0.40 Btu/h · ft · °F and whose thermal diffusivity is 4.78 10 6 ft2/s. The variation of the ambient tem· perature Tout and the solar heat flux q solar incident on a south-facing vertical surface throughout the day for a typical day in January is given in Table 5–4 in 3-h intervals. The Trombe wall has single glazing with an absorptivity-transmissivity product of 0.77 (that is, 77 percent of the solar energy incident is absorbed by the exposed surface of the Trombe wall), and the average combined heat transfer coefficient for heat loss from the Trombe wall to the ambient is determined to be hout 0.7 Btu/h · ft2 · °F. The interior of the house is maintained at Tin 70°F at all times, and the heat transfer coefficient at the interior surface of the Trombe wall is hin 1.8 Btu/h · ft2 · °F. Also, the vents on the Trombe wall are kept closed, and thus the only heat transfer between the air in the house and the Trombe wall is through the interior surface of the wall. Assuming the temperature of the Trombe wall to vary linearly between 70°F at the interior surface and 30°F at the exterior surface at 7 AM and using the explicit finite difference method with a uniform nodal spacing of x 0.2 ft, determine the temperature distribution along the thickness of the Trombe wall after 12, 24, 36, and 48 h. Also, determine the net amount of heat transferred to the house from the Trombe wall during the first day and the second day. Assume the wall is 10 ft high and 25 ft long. SOLUTION The passive solar heating of a house through a Trombe wall is considered. The temperature distribution in the wall in 12-h intervals and the amount of heat transfer during the first and second days are to be determined. Assumptions 1 Heat transfer is one-dimensional since the exposed surface of the wall is large relative to its thickness. 2 Thermal conductivity is constant. 3 The heat transfer coefficients are constant. Properties The wall properties are given to be k 0.40 Btu/h · ft · °F, 4.78 10 6 ft2/s, and 0.77. Analysis The nodal spacing is given to be x 0.2 ft, and thus the total number of nodes along the Trombe wall is M L x 1 1 ft 0.2 ft 1 6 We number the nodes as 0, 1, 2, 3, 4, and 5, with node 0 on the interior surface of the Trombe wall and node 5 on the exterior surface, as shown in Figure 5–47. Nodes 1 through 4 are interior nodes, and the explicit finite difference formulations of these nodes are obtained directly from Eq. 5–47 to be cen58933_ch05.qxd 9/4/2002 11:42 AM Page 301 301 CHAPTER 5 Node 1 (m Node 2 (m Node 3 (m Node 4 (m T1i T2i T3i T4i 1): 2): 3): 4): 1 (T0i (T1i (T2i (T3i 1 1 1 T2i) T3i) T4i) T5i) (1 (1 (1 (1 )T1i )T2i )T3i )T4i 2 2 2 2 Trombe wall (2) k = 0.40 Btu/h·ft·°F α = 4.78 × 10 – 6 ft2/s (3) (4) The interior surface is subjected to convection, and thus the explicit formulation of node 0 can be obtained directly from Eq. 5–51 to be T0i 1 1 hin x i T0 2 k 2 Substituting the quantities hin, into this equation gives T0i 1 · Aq isolar T5i) 0 (2T1i 126.0) (5) kA T4i T5i i1 A x x T5 C 2 T5i t (5-53) which simplifies to T5i 1 1 2 hout x i T5 k 2 2 T4i 2 hout x i Tout k 2 ·i qsolar x (5-54) k where t/ x2 is the dimensionless mesh Fourier number. Note that we kept the superscript i for quantities that vary with time. Substituting the quantities hout, x, k, and , which do not change with time, into this equation gives T5i 1 2.70 ) T5i (1 (2T4i i 0.70Tout 0.770q·isolar) (6) · where the unit of q isolar is Btu/h · ft2. Next we need to determine the upper limit of the time step t from the stability criterion since we are using the explicit method. This requires the identification of the smallest primary coefficient in the system. We know that the boundary nodes are more restrictive than the interior nodes, and thus we examine the formulations of the boundary nodes 0 and 5 only. The smallest and thus i the most restrictive primary coefficient in this case is the coefficient of T0 in the formulation of node 0 since 1 3.8 1 2.7 , and thus the stability criterion for this problem can be expressed as 1 3.80 0 → αx x2 1 3.80 Substituting the given quantities, the maximum allowable value of the time step is determined to be t x2 3.80α 3.80 30°F ∆ x = 0.2 ft The exterior surface of the Trombe wall is subjected to convection as well as to heat flux. The explicit finite difference formulation at that boundary is obtained by writing an energy balance on the volume element represented by node 5, i hout A(Tout hout, Tout 0 1 2 3 4 5L x x, k, and Tin, which do not change with time, 3.80 ) T0i (1 2 Initial temperature distribution at 7 AM (t = 0) 70°F hin, Tin hin x 2 Tin k T1i · qsolar (1) (0.2 ft)2 (4.78 10 6 ft2/s) 2202 s FIGURE 5–47 The nodal network for the Trombe wall discussed in Example 5–6. cen58933_ch05.qxd 9/4/2002 11:42 AM Page 302 302 HEAT TRANSFER Therefore, any time step less than 2202 s can be used to solve this problem. For convenience, let us choose the time step to be t 900 s 15 min. Then the mesh Fourier number becomes αt ( x)2 (4.78 10 6 ft2/s)(900 s) (0.2 ft)2 0.10755 (for t 15 min) Initially (at 7 AM or t 0), the temperature of the wall is said to vary linearly between 70°F at node 0 and 30°F at node 5. Noting that there are five nodal spacings of equal length, the temperature change between two neighboring nodes is (70 30)°F/5 8°F. Therefore, the initial nodal temperatures are T00 T30 70°F, 46°F, Then the nodal temperatures at t from these equations to be T01 T11 T21 T31 Temperature T41 °F 170 T51 150 1st day 2nd day 130 7 PM 110 1 AM 90 1 PM 70 50 Initial temperature 30 0 7 AM 0.2 0.4 0.6 0.8 Distance along the Trombe wall FIGURE 5–48 The variation of temperatures in the Trombe wall discussed in Example 5–6. 1 ft T10 T40 62°F, 38°F, t T20 T50 54°F, 30°F 15 min (at 7:15 AM) are determined (1 3.80 ) T00 (2T10 126.0) (1 3.80 0.10755) 70 0.10755(2 62 126.0) 68.3° F (T00 T20) (1 2 ) T10 0.10755(70 54) (1 2 0.10755)62 62°F (T10 T30) (1 2 ) T20 0.10755(62 46) (1 2 0.10755)54 54°F (T20 T40) (1 2 ) T30 0.10755(54 38) (1 2 0.10755)46 46°F (T30 T50) (1 2 ) T40 0.10755(46 30) (1 2 0.10755)38 38°F 0 · (1 2.70 ) T50 (2T40 0.70Tout 0.770q 0 ) solar (1 2.70 0.10755)30 0.10755(2 38 0.70 33 0.770 41.4°F 114) Note that the inner surface temperature of the Trombe wall dropped by 1.7°F and the outer surface temperature rose by 11.4°F during the first time step while the temperatures at the interior nodes remained the same. This is typical of transient problems in mediums that involve no heat generation. The nodal temperatures at the following time steps are determined similarly with the help of a computer. Note that the data for ambient temperature and the incident solar radiation change every 3 hours, which corresponds to 12 time steps, and this must be reflected in the computer program. For example, the value of · · · q isolar must be taken to be q isolar 75 for i 1–12, q isolar 242 for i 13–24, ·i ·i q solar 178 for i 25–36, and q solar 0 for i 37–96. The results after 6, 12, 18, 24, 30, 36, 42, and 48 h are given in Table 5–5 and are plotted in Figure 5–48 for the first day. Note that the interior temperature of the Trombe wall drops in early morning hours, but then rises as the solar energy absorbed by the exterior surface diffuses through the wall. The exterior surface temperature of the Trombe wall rises from 30 to 142°F in just 6 h because of the solar energy absorbed, but then drops to 53°F by next morning as a result of heat loss at night. Therefore, it may be worthwhile to cover the outer surface at night to minimize the heat losses. cen58933_ch05.qxd 9/4/2002 11:42 AM Page 303 303 CHAPTER 5 TABLE 5–5 The temperatures at the nodes of a Trombe wall at various times Time 0 6 12 18 24 30 36 42 48 h h h h h h h h h (7 (1 (7 (1 (7 (1 (7 (1 (7 Nodal Temperatures, °F Time Step, i AM) PM) PM) AM) AM) PM) PM) AM) AM) T0 T1 T2 T3 T4 T5 0 24 48 72 96 120 144 168 192 70.0 65.3 71.6 73.3 71.2 70.3 75.4 75.8 73.0 62.0 61.7 74.2 75.9 71.9 71.1 81.1 80.7 75.1 54.0 61.5 80.4 77.4 70.9 74.3 89.4 83.5 72.2 46.0 69.7 88.4 76.3 67.7 84.2 98.2 83.0 66.0 38.0 94.1 91.7 71.2 61.7 108.3 101.0 77.4 66.0 30.0 142.0 82.4 61.2 53.0 153.2 89.7 66.2 56.3 The rate of heat transfer from the Trombe wall to the interior of the house during each time step is determined from Newton’s law using the average temperature at the inner surface of the wall (node 0) as QiTrombe wall · Q iTrombe wall t hin A(T0i Tin) t hin A[(T0i T0i 1)/2 Therefore, the amount of heat transfer during the first time step (i during the first 15-min period is Q1 Trombe wall hin A[(T01 T00)/2 Tin] t (1.8 Btu/h · ft2 · °F)(10 25 ft2)[(68.3 95.6 Btu 70)/2 Tin] t 1) or 70°F](0.25 h) The negative sign indicates that heat is transferred to the Trombe wall from the air in the house, which represents a heat loss. Then the total heat transfer during a specified time period is determined by adding the heat transfer amounts for each time step as I QTrombe wall i1 I · Q iTrombe wall hin A[(T0i T0i 1)/2 Tin] t (5-55) i1 where I is the total number of time intervals in the specified time period. In this case I 48 for 12 h, 96 for 24 h, and so on. Following the approach described here using a computer, the amount of heat transfer between the Trombe wall and the interior of the house is determined to be QTrombe wall QTrombe wall QTrombe wall QTrombe wall 17, 048 Btu after 12 h 2483 Btu after 24 h 5610 Btu after 36 h 34, 400 Btu after 48 h ( 17, 078 Btu during the first 12 h) (14, 565 Btu during the second 12 h) (8093 Btu during the third 12 h) (28, 790 Btu during the fourth 12 h) Therefore, the house loses 2483 Btu through the Trombe wall the first day as a result of the low start-up temperature but delivers a total of 36,883 Btu of heat to the house the second day. It can be shown that the Trombe wall will deliver even more heat to the house during the third day since it will start the day at a higher average temperature. cen58933_ch05.qxd 9/4/2002 11:42 AM Page 304 304 HEAT TRANSFER Two-Dimensional Transient Heat Conduction m, n + 1 n+1 ∆y Volume element m – 1, n n · gm, n m, n m + 1, n ∆y m, n – 1 n–1 ∆x y m–1 ∆x m m+1 x FIGURE 5–49 The volume element of a general interior node (m, n) for twodimensional transient conduction in rectangular coordinates. Consider a rectangular region in which heat conduction is significant in the x- and y-directions, and consider a unit depth of z 1 in the z-direction. · Heat may be generated in the medium at a rate of g(x, y, t), which may vary with time and position, with the thermal conductivity k of the medium assumed to be constant. Now divide the x-y-plane of the region into a rectangular mesh of nodal points spaced x and y apart in the x- and y-directions, respectively, and consider a general interior node (m, n) whose coordinates are x m x and y n y, as shown in Figure 5–49. Noting that the volume element centered about the general interior node (m, n) involves heat conduction from four sides (right, left, top, and bottom) and the volume of the element is x y1 x y, the transient finite difference formulation for Velement a general interior node can be expressed on the basis of Eq. 5–39 as ky Tm kx Tm, n 1, n kx x Tm, n Tm, n 1 Tm 1, n Tm 1, n Tm, n y 1 Tm, n 1 ky y · gm, n x y y Taking a square mesh ( x simplifying, Tm, n x yC Tm Tm, n 1, n x i Tm 1 i Tm (5-56) t l ) and dividing each term by k gives after Tm, n 1 4Tm, n · gm, nl 2 k i Tm 1 i Tm (5-57) where again k/ C is the thermal diffusivity of the material and t/l 2 is the dimensionless mesh Fourier number. It can also be expressed in terms of the temperatures at the neighboring nodes in the following easy-to-remember form: Tleft Ttop Tright Tbottom 4Tnode · gnodel 2 k i1 Tnode i Tnode (5-58) Again the left side of this equation is simply the finite difference formulation of the problem for the steady case, as expected. Also, we are still not committed to explicit or implicit formulation since we did not indicate the time step on the left side of the equation. We now obtain the explicit finite difference formulation by expressing the left side at time step i as i Tleft i Ttop i Tright i Tbottom i 4Tnode ·i gnodel 2 k i1 Tnode i Tnode (5-59) Expressing the left side at time step i 1 instead of i would give the implicit formulation. This equation can be solved explicitly for the new temperature i1 Tnode to give i1 Tnode i (Tleft i Ttop i Tright i Tbottom) (1 i 4 ) Tnode ·i gnodel 2 k for all interior nodes (m, n) where m 1, 2, 3, . . . , M 1 and n 3, . . . , N 1 in the medium. In the case of no heat generation and (5-60) 1, 2, the 1 , 4 cen58933_ch05.qxd 9/4/2002 11:42 AM Page 305 305 CHAPTER 5 explicit finite difference formulation for a general interior node reduces to i1 i i i i Tnode (Tleft Ttop Tright Tbottom)/4, which has the interpretation that the temperature of an interior node at the new time step is simply the average of the temperatures of its neighboring nodes at the previous time step (Fig. 5–50). i i The stability criterion that requires the coefficient of Tm in the Tm 1 expression to be greater than or equal to zero for all nodes is equally valid for twoor three-dimensional cases and severely limits the size of the time step t that can be used with the explicit method. In the case of transient two-dimensional i i heat transfer in rectangular coordinates, the coefficient of Tm in the Tm 1 expression is 1 4 , and thus the stability criterion for all interior nodes in this case is 1 4 0, or t l2 1 4 (interior nodes, two-dimensional heat transfer in rectangular coordinates) Time step i: 30°C i Tm 20°C 40°C Node m 10°C Time step i + 1: (5-61) i Tm+ 1 25°C Node m where x y l. When the material of the medium and thus its thermal diffusivity are known and the value of the mesh size l is specified, the largest allowable value of the time step t can be determined from the relation above. Again the boundary nodes involving convection and/or radiation are more restrictive than the interior nodes and thus require smaller time steps. Therefore, the most restrictive boundary node should be used in the determination of the maximum allowable time step t when a transient problem is solved with the explicit method. The application of Eq. 5–60 to each of the (M 1) (N 1) interior nodes gives (M 1) (N 1) equations. The remaining equations are obtained by applying the method to the boundary nodes unless, of course, the boundary temperatures are specified as being constant. The development of the transient finite difference formulation of boundary nodes in two- (or three-) dimensional problems is similar to the development in the one-dimensional case discussed earlier. Again the region is partitioned between the nodes by forming volume elements around the nodes, and an energy balance is written for each boundary node on the basis of Eq. 5–39. This is illustrated in Example 5–7. EXAMPLE 5–7 FIGURE 5–50 In the case of no heat generation 1 and , the temperature of an 4 interior node at the new time step is the average of the temperatures of its neighboring nodes at the previous time step. Transient Two-Dimensional Heat Conduction in L-Bars Consider two-dimensional transient heat transfer in an L-shaped solid body that is initially at a uniform temperature of 90°C and whose cross section is given in Figure 5–51. The thermal conductivity and diffusivity of the body are k 15 W/m · °C and 3.2 10 6 m2/s, respectively, and heat is generated in · the body at a rate of g 2 106 W/m3. The left surface of the body is insulated, and the bottom surface is maintained at a uniform temperature of 90°C at all times. At time t 0, the entire top surface is subjected to convection to ambient air at T 25°C with a convection coefficient of h 80 W/m2 · °C, · and the right surface is subjected to heat flux at a uniform rate of q R 5000 2 W/m . The nodal network of the problem consists of 15 equally spaced nodes with x y 1.2 cm, as shown in the figure. Five of the nodes are at the bottom surface, and thus their temperatures are known. Using the explicit method, determine the temperature at the top corner (node 3) of the body after 1, 3, 5, 10, and 60 min. Convection h, T y 1 2 3 4 5 6 7 8 10 11 12 13 14 15 ∆ x = ∆y = l ∆y · 9 qR ∆y ∆x ∆x 90°C ∆x x ∆x ∆x FIGURE 5–51 Schematic and nodal network for Example 5–7. cen58933_ch05.qxd 9/4/2002 11:42 AM Page 306 306 HEAT TRANSFER h, T 1 SOLUTION This is a transient two-dimensional heat transfer problem in rectangular coordinates, and it was solved in Example 5–3 for the steady case. Therefore, the solution of this transient problem should approach the solution for the steady case when the time is sufficiently large. The thermal conductivity and heat generation rate are given to be constants. We observe that all nodes are boundary nodes except node 5, which is an interior node. Therefore, we will have to rely on energy balances to obtain the finite difference equations. The region is partitioned among the nodes equitably as shown in the figure, and the explicit finite difference equations are determined on the basis of the energy balance for the transient case expressed as h, T 2 1 2 3 · Qi · G ielement Velement C i Tm 1 t All sides 5 4 (a) Node 1 (b) Node 2 FIGURE 5–52 Schematics for energy balances on the volume elements of nodes 1 and 2. i Tm · · The quantities h, T , g , and q R do not change with time, and thus we do not need to use the superscript i for them. Also, the energy balance expressions are simplified using the definitions of thermal diffusivity k/ C and the dimensionless mesh Fourier number t/l2, where x y l. (a) Node 1. (Boundary node subjected to convection and insulation, Fig. 5–52a) h x (T 2 T1i) y T2i T1i 2 x k i k i x T4 T1 2 y · xy g1 22 i x y T1 C 22 1 T1i t Dividing by k/4 and simplifying, 2hl (T k T1i) 2(T2i i which can be solved for T1 T1i 1 1 4 1 T1i) 2(T4i · g1l 2 k T1i) T1i 1 T1i to give 2 hl T1i k 2 T2i hl T k T4i · g1l 2 2k (b) Node 2. (Boundary node subjected to convection, Fig. 5–52b) T2i) h x(T y T3i T2i 2 x y T1i T2i k 2 x k kx · g2 T5i y y x 2 i Dividing by k/2, simplifying, and solving for T2 T2i 1 1 4 2 hl T2i k T1i T2i 1 T3i x y T2i C 2 1 T2i t gives 2T5i 2hl T k · g2l 2 k cen58933_ch05.qxd 9/4/2002 11:42 AM Page 307 307 CHAPTER 5 (c) Node 3. (Boundary node subjected to convection on two sides, Fig. 5–53a) i i k x T6 T3 2 y y T2i T3i 2 x · xy g3 22 y (T 2 x 2 h T3i) k i Dividing by k/4, simplifying, and solving for T3 T3i 1 1 4 hl T3i k 4 2 3 (5) 4 5 h, T i x y T3 22 1 1 T3i t 10 6 gives T4i 2 1 Mirror h, T (a) Node 3 T6i 2 FIGURE 5–53 Schematics for energy balances on the volume elements of nodes 3 and 4. · g3l 2 2k hl T k (b) Node 4 (d ) Node 4. (On the insulated boundary, and can be treated as an interior node, Fig. 5–53b). Noting that T10 90°C, Eq. 5–60 gives T4i 1 4 ) T4i (1 2T5i T1i · g4l 2 k 90 (e) Node 5. (Interior node, Fig. 5–54a). Noting that T11 gives T5i 1 4 ) T5i (1 T2i T4i T6i 90°C, Eq. 5–60 · g5l 2 k 90 (f ) Node 6. (Boundary node subjected to convection on two sides, Fig. 5–54b) h y (T 2 x 2 T6i) i k i x T3 T6 2 y y 2 T7i T6i kx x ·3xy g6 4 i T12 T6i y 3 x y T6i C 4 ky 1 T5i T6i x 3 2 4 5 6 T6i 1 1 4 3 4 2T3i 1 T6i t gives 12 2T7i 4 90 4 FIGURE 5–54 Schematics for energy balances on the volume elements of nodes 5 and 6. hl T k 3 · g6l 2 k (g) Node 7. (Boundary node subjected to convection, Fig. 5–55) T7i) h x(T k i T13 T7i y T8i T7i kx 2 x y i i y y T6 T7 · k g7 x 2 2 x i Dividing by k/2, simplifying, and solving for T7 T7i 1 1 4 2 hl T7i k T6i T8i 1 h, T h, T 6 7 8 9 · qR x y T7i C 2 1 T7i t gives 2 (b) Node 6 (a) Node 5 hl Ti 3k 3 4T5i 7 6 5 11 i Dividing by 3k/4, simplifying, and solving for T6 h, T 90 2hl T k · g7l 2 k 13 15 FIGURE 5–55 Schematics for energy balances on the volume elements of nodes 7 and 9. cen58933_ch05.qxd 9/4/2002 11:42 AM Page 308 308 HEAT TRANSFER (h) Node 8. This node is identical to node 7, and the finite difference formulation of this node can be obtained from that of node 7 by shifting the node numbers by 1 (i.e., replacing subscript m by subscript m 1). It gives T8i 1 1 4 hl T8i k 2 T7i T9i 2 · g8l 2 k 2hl T k 90 (i ) Node 9. (Boundary node subjected to convection on two sides, Fig. 5–55) h i i y x T15 T9 q·R k 2 2 y k y T8i T9i · xy g9 2 22 x x (T 2 T9i) i Dividing by k/4, simplifying, and solving for T9 T9i 1 1 4 2 hl T9i k 2 T8i i x y T9 C 22 1 1 T9i t gives · qR l k 90 hl T k · g9l 2 2k This completes the finite difference formulation of the problem. Next we need to determine the upper limit of the time step t from the stability criterion, i i which requires the coefficient of Tm in the Tm 1 expression (the primary coefficient) to be greater than or equal to zero for all nodes. The smallest primary coi efficient in the nine equations here is the coefficient of T3 in the expression, and thus the stability criterion for this problem can be expressed as 1 4 hl k 4 → 0 1 4(1 hl/k) → t l2 4 (1 hl/k) since t /l 2. Substituting the given quantities, the maximum allowable value of the time step is determined to be t 4(3.2 10 6 2 m /s)[1 (0.012 m)2 (80 W/m2 · °C)(0.012 m)/(15 W/m · °C)] 10.6 s Therefore, any time step less than 10.6 s can be used to solve this problem. For convenience, let us choose the time step to be t 10 s. Then the mesh Fourier number becomes t l (3.2 2 10 6 m2/s)(10 s) (0.012 m)2 0.222 (for t 10 s) Substituting this value of and other given quantities, the developed transient finite difference equations simplify to T1i T2i T3i T4i T5i 1 1 1 1 1 0.0836T1i 0.444(T2i T4i 11.2) 0.0836T2i 0.222(T1i T3i 2T5i 22.4) 0.0552T3i 0.444(T2i T6i 12.8) 0.112T4i 0.222(T1i 2T5i 109.2) 0.112T1i 0.222(T2i T4i T6i 109.2) cen58933_ch05.qxd 9/4/2002 11:42 AM Page 309 309 CHAPTER 5 T6i T7i T8i T9i 1 1 1 1 0.0931T6i 0.0836T7i 0.0836T8i 0.0836T9i 0.074(2T3i 4T5i 2T7i 0.222(T6i T8i 202.4) 0.222(T7i T9i 202.4) 0.444(T8i 105.2) 424) Using the specified initial condition as the solution at time t 0 (for i 0), sweeping through these nine equations will give the solution at intervals of 10 s. The solution at the upper corner node (node 3) is determined to be 100.2, 105.9, 106.5, 106.6, and 106.6°C at 1, 3, 5, 10, and 60 min, respectively. Note that the last three solutions are practically identical to the solution for the steady case obtained in Example 5–3. This indicates that steady conditions are reached in the medium after about 5 min. TOPIC OF SPECIAL INTEREST Controlling the Numerical Error A comparison of the numerical results with the exact results for temperature distribution in a cylinder would show that the results obtained by a numerical method are approximate, and they may or may not be sufficiently close to the exact (true) solution values. The difference between a numerical solution and the exact solution is the error involved in the numerical solution, and it is primarily due to two sources: • The discretization error (also called the truncation or formulation error), which is caused by the approximations used in the formulation of the numerical method. • The round-off error, which is caused by the computer ’s use of a limited number of significant digits and continuously rounding (or chopping) off the digits it cannot retain. Below we discuss both types of errors. T(xm, t) Discretization Error The discretization error involved in numerical methods is due to replacing the derivatives by differences in each step, or the actual temperature distribution between two adjacent nodes by a straight line segment. Consider the variation of the solution of a transient heat transfer problem with time at a specified nodal point. Both the numerical and actual (exact) solutions coincide at the beginning of the first time step, as expected, but the numerical solution deviates from the exact solution as the time t increases. The difference between the two solutions at t t is due to the approximation at the first time step only and is called the local discretization error. One would expect the situation to get worse with each step since the second step uses the erroneous result of the first step as its starting point and adds a second local discretization error on top of it, as shown in Figure 5–56. The accumulation of the local discretization errors continues with the increasing number of time steps, and the total discretization error at any Local error Actual solution T(x0, t) Global error T3 T2 T0 T1 Step 1 t0 Numerical solution Step 2 t1 Step 3 t2 t3 Time FIGURE 5–56 The local and global discretization errors of the finite difference method at the third time step at a specified nodal point. cen58933_ch05.qxd 9/4/2002 11:42 AM Page 310 310 HEAT TRANSFER step is called the global or accumulated discretization error. Note that the local and global discretization errors are identical for the first time step. The global discretization error usually increases with the increasing number of steps, but the opposite may occur when the solution function changes direction frequently, giving rise to local discretization errors of opposite signs, which tend to cancel each other. To have an idea about the magnitude of the local discretization error, consider the Taylor series expansion of the temperature at a specified nodal point m about time ti, T(xm, ti t) T(xm, ti) t T(xm, ti) t 12 t 2 2 T(xm, ti) t2 ··· (5-62) The finite difference formulation of the time derivative at the same nodal point is expressed as T(xm, ti) t T(xm, ti t) t T(xm, ti) i Tm 1 i Tm t (5-63) or T(xm, ti t) T(xm, ti) t T(xm, ti) t (5-64) which resembles the Taylor series expansion terminated after the first two terms. Therefore, the third and later terms in the Taylor series expansion represent the error involved in the finite difference approximation. For a sufficiently small time step, these terms decay rapidly as the order of derivative increases, and their contributions become smaller and smaller. The first term neglected in the Taylor series expansion is proportional to t 2, and thus the local discretization error of this approximation, which is the error involved in each step, is also proportional to t 2. The local discretization error is the formulation error associated with a single step and gives an idea about the accuracy of the method used. However, the solution results obtained at every step except the first one involve the accumulated error up to that point, and the local error alone does not have much significance. What we really need to know is the global discretization error. At the worst case, the accumulated discretization error (t0/ t)( t)2 t0 t, after I time steps during a time period t0 is i( t)2 which is proportional to t. Thus, we conclude that the local discretization error is proportional to the square of the step size t 2 while the global discretization error is proportional to the step size t itself. Therefore, the smaller the mesh size (or the size of the time step in transient problems), the smaller the error, and thus the more accurate is the approximation. For example, halving the step size will reduce the global discretization error by half. It should be clear from the discussions above that the discretization error can be minimized by decreasing the step size in space or time as much as possible. The discretization error approaches zero as the difference quantities such as x and t approach the differential quantities such as dx and dt. cen58933_ch05.qxd 9/4/2002 11:42 AM Page 311 311 CHAPTER 5 Round-off Error If we had a computer that could retain an infinite number of digits for all numbers, the difference between the exact solution and the approximate (numerical) solution at any point would entirely be due to discretization error. But we know that every computer (or calculator) represents numbers using a finite number of significant digits. The default value of the number of significant digits for many computers is 7, which is referred to as single precision. But the user may perform the calculations using 15 significant digits for the numbers, if he or she wishes, which is referred to as double precision. Of course, performing calculations in double precision will require more computer memory and a longer execution time. In single precision mode with seven significant digits, a computer will register the number 44444.666666 as 44444.67 or 44444.66, depending on the method of rounding the computer uses. In the first case, the excess digits are said to be rounded to the closest integer, whereas in the second case they are said to be chopped off. Therefore, the numbers a 44444.12345 and b 44444.12032 are equivalent for a computer that performs calculations using seven significant digits. Such a computer would give a b 0 instead of the true value 0.00313. The error due to retaining a limited number of digits during calculations is called the round-off error. This error is random in nature and there is no easy and systematic way of predicting it. It depends on the number of calculations, the method of rounding off, the type of computer, and even the sequence of calculations. In algebra you learned that a b c a c b, which seems quite reasonable. But this is not necessarily true for calculations performed with a computer, as demonstrated in Figure 5–57. Note that changing the sequence of calculations results in an error of 30.8 percent in just two operations. Considering that any significant problem involves thousands or even millions of such operations performed in sequence, we realize that the accumulated round-off error has the potential to cause serious error without giving any warning signs. Experienced programmers are very much aware of this danger, and they structure their programs to prevent any buildup of the round-off error. For example, it is much safer to multiply a number by 10 than to add it 10 times. Also, it is much safer to start any addition process with the smallest numbers and continue with larger numbers. This rule is particularly important when evaluating series with a large number of terms with alternating signs. The round-off error is proportional to the number of computations performed during the solution. In the finite difference method, the number of calculations increases as the mesh size or the time step size decreases. Halving the mesh or time step size, for example, will double the number of calculations and thus the accumulated round-off error. Controlling the Error in Numerical Methods The total error in any result obtained by a numerical method is the sum of the discretization error, which decreases with decreasing step size, and the round-off error, which increases with decreasing step size, as shown in Figure 5–58. Therefore, decreasing the step size too much in order to get more Given: a b c D E Find: 7777777 7777776 0.4444432 a a b c c b Solution: D 7777777 7777776 0.4444432 1 0.4444432 1.444443 (Correct result) E 7777777 0.4444432 7777776 7777777 7777776 1.000000 (In error by 30.8%) FIGURE 5–57 A simple arithmetic operation performed with a computer in single precision using seven significant digits, which results in 30.8 percent error when the order of operation is reversed. Error Total error Discretization error Round-off error Optimum step size Step size FIGURE 5–58 As the mesh or time step size decreases, the discretization error decreases but the round-off error increases. cen58933_ch05.qxd 9/4/2002 11:42 AM Page 312 312 HEAT TRANSFER accurate results may actually backfire and give less accurate results because of a faster increase in the round-off error. We should be careful not to let round-off error get out of control by avoiding a large number of computations with very small numbers. In practice, we will not know the exact solution of the problem, and thus we will not be able to determine the magnitude of the error involved in the numerical method. Knowing that the global discretization error is proportional to the step size is not much help either since there is no easy way of determining the value of the proportionality constant. Besides, the global discretization error alone is meaningless without a true estimate of the round-off error. Therefore, we recommend the following practical procedures to assess the accuracy of the results obtained by a numerical method. • Start the calculations with a reasonable mesh size x (and time step size t for transient problems) based on experience. Then repeat the calculations using a mesh size of x/2. If the results obtained by halving the mesh size do not differ significantly from the results obtained with the full mesh size, we conclude that the discretization error is at an acceptable level. But if the difference is larger than we can accept, then we have to repeat the calculations using a mesh size x/4 or even a smaller one at regions of high temperature gradients. We continue in this manner until halving the mesh size does not cause any significant change in the results, which indicates that the discretization error is reduced to an acceptable level. • Repeat the calculations using double precision holding the mesh size (and the size of the time step in transient problems) constant. If the changes are not significant, we conclude that the round-off error is not a problem. But if the changes are too large to accept, then we may try reducing the total number of calculations by increasing the mesh size or changing the order of computations. But if the increased mesh size gives unacceptable discretization errors, then we may have to find a reasonable compromise. It should always be kept in mind that the results obtained by any numerical method may not reflect any trouble spots in certain problems that require special consideration such as hot spots or areas of high temperature gradients. The results that seem quite reasonable overall may be in considerable error at certain locations. This is another reason for always repeating the calculations at least twice with different mesh sizes before accepting them as the solution of the problem. Most commercial software packages have built-in routines that vary the mesh size as necessary to obtain highly accurate solutions. But it is a good engineering practice to be aware of any potential pitfalls of numerical methods and to examine the results obtained with a critical eye. SUMMARY Analytical solution methods are limited to highly simplified problems in simple geometries, and it is often necessary to use a numerical method to solve real world problems with complicated geometries or nonuniform thermal conditions. The cen58933_ch05.qxd 9/4/2002 11:42 AM Page 313 313 CHAPTER 5 numerical finite difference method is based on replacing derivatives by differences, and the finite difference formulation of a heat transfer problem is obtained by selecting a sufficient number of points in the region, called the nodal points or nodes, and writing energy balances on the volume elements centered about the nodes. For steady heat transfer, the energy balance on a volume element can be expressed in general as · Q · gVelement 0 also be used to solve a system of equations simultaneously at the press of a button. The finite difference formulation of transient heat conduction problems is based on an energy balance that also accounts for the variation of the energy content of the volume element during a time interval t. The heat transfer and heat generation terms are expressed at the previous time step i in the explicit method, and at the new time step i 1 in the implicit method. For a general node m, the finite difference formulations are expressed as All sides whether the problem is one-, two-, or three-dimensional. For convenience in formulation, we always assume all heat transfer to be into the volume element from all surfaces toward the node under consideration, except for specified heat flux whose direction is already specified. The finite difference formulations for a general interior node under steady conditions are expressed for some geometries as follows: One-dimensional steady conduction in a plane wall: Twodimensional steady Ttop T conduction left in rectangular coordinates: Tm 2Tm ( x)2 1 Tright Tbottom Tm 1 4Tnode · gm k 0 · gnodel 2 k 0 where x is the nodal spacing for the plane wall and x y l is the nodal spacing for the two-dimensional case. Insulated boundaries can be viewed as mirrors in formulation, and thus the nodes on insulated boundaries can be treated as interior nodes by using mirror images. The finite difference formulation at node 0 at the left boundary of a plane wall for steady one-dimensional heat conduction can be expressed as · Q left surface kA T1 T0 x · g0(A x/2) 0 · where A x/2 is the volume of the volume, g0 is the rate of heat generation per unit volume at x 0, and A is the heat transfer area. The form of the first term depends on the boundary condition at x 0 (convection, radiation, specified heat flux, etc.). The finite difference formulation of heat conduction problems usually results in a system of N algebraic equations in N unknown nodal temperatures that need to be solved simultaneously. There are numerous systematic approaches available in the literature. Several widely available equation solvers can Explicit method: Implicit method: · G ielement · Qi Velement C i Tm 1 t All sides · Qi ·i 1 G element 1 Velement C i Tm i Tm 1 All sides i Tm t i i where Tm and Tm 1 are the temperatures of node m at times i i ti i t and ti 1 (i 1) t, respectively, and Tm 1 Tm represents the temperature change of the node during the time interval t between the time steps i and i 1. The explicit and implicit formulations given here are quite general and can be used in any coordinate system regardless of heat transfer being one-, two-, or three-dimensional. The explicit formulation of a general interior node for oneand two-dimensional heat transfer in rectangular coordinates can be expressed as Onei Tm 1 dimensional case: Twoi1 Tnode dimensional case: i (Tm i (Tleft (1 1 i Tm 1) (1 i Tright i Ttop i 4 ) Tnode i 2 ) Tm ·i gm x2 k i Tbottom) ·i gnodel 2 k where t x2 is the dimensionless mesh Fourier number and k/ C is the thermal diffusivity of the medium. The implicit method is inherently stable, and any value of t can be used with that method as the time step. The largest value of the time step t in the explicit method is limited by the stai bility criterion, expressed as: the coefficients of all Tm in the i1 Tm expressions (called the primary coefficients) must be greater than or equal to zero for all nodes m. The maximum value of t is determined by applying the stability criterion to the equation with the smallest primary coefficient since it is the cen58933_ch05.qxd 9/4/2002 11:42 AM Page 314 314 HEAT TRANSFER most restrictive. For problems with specified temperatures or heat fluxes at all the boundaries, the stability criterion can be 1 expressed as for one-dimensional problems and 2 the two-dimensional problems in rectangular coordinates. 1 4 for REFERENCES AND SUGGESTED READING 1. D. A. Anderson, J. C. Tannehill, and R. H. Pletcher. Computational Fluid Mechanics and Heat Transfer. New York: Hemisphere, 1984. 7. W. J. Minkowycz, E. M. Sparrow, G. E. Schneider, and R. H. Pletcher. Handbook of Numerical Heat Transfer. New York: John Wiley & Sons, 1988. 2. C. A. Brebbia. The Boundary Element Method for Engineers. New York: Halsted Press, 1978. 8. G. E. Myers. Analytical Methods in Conduction Heat Transfer. New York: McGraw-Hill, 1971. 3. G. E. Forsythe and W. R. Wasow. Finite Difference Methods for Partial Differential Equations. New York: John Wiley & Sons, 1960. 9. D. H. Norrie and G. DeVries. An Introduction to Finite Element Analysis. New York: Academic Press, 1978. 4. B. Gebhart. Heat Conduction and Mass Diffusion. New York: McGraw-Hill, 1993. 5. K. H. Huebner and E. A. Thornton. The Finite Element Method for Engineers. 2nd ed. New York: John Wiley & Sons, 1982. 10. M. N. Özisik. Finite Difference Methods in Heat Transfer. ¸ Boca Raton, FL: CRC Press, 1994. 11. S. V. Patankhar. Numerical Heat Transfer and Fluid Flow. New York: Hemisphere, 1980. 12. T. M. Shih. Numerical Heat Transfer. New York: Hemisphere, 1984. 6. Y. Jaluria and K. E. Torrance. Computational Heat Transfer. New York: Hemisphere, 1986. PROBLEMS* Why Numerical Methods? 5–1C What are the limitations of the analytical solution methods? 5–2C How do numerical solution methods differ from analytical ones? What are the advantages and disadvantages of numerical and analytical methods? 5–3C What is the basis of the energy balance method? How does it differ from the formal finite difference method? For a specified nodal network, will these two methods result in the same or a different set of equations? 5–4C Consider a heat conduction problem that can be solved both analytically, by solving the governing differential equation and applying the boundary conditions, and numerically, by a software package available on your computer. Which *Problems designated by a “C” are concept questions, and students are encouraged to answer them all. Problems designated by an “E” are in English units, and the SI users can ignore them. Problems with an EES-CD icon are solved using EES, and complete solutions together with parametric studies are included on the enclosed CD. Problems with a computer-EES icon are comprehensive in nature, and are intended to be solved with a computer, preferably using the EES software that accompanies this text. approach would you use to solve this problem? Explain your reasoning. 5–5C Two engineers are to solve an actual heat transfer problem in a manufacturing facility. Engineer A makes the necessary simplifying assumptions and solves the problem analytically, while engineer B solves it numerically using a powerful software package. Engineer A claims he solved the problem exactly and thus his results are better, while engineer B claims that he used a more realistic model and thus his results are better. To resolve the dispute, you are asked to solve the problem experimentally in a lab. Which engineer do you think the experiments will prove right? Explain. Finite Difference Formulation of Differential Equations 5–6C Define these terms used in the finite difference formulation: node, nodal network, volume element, nodal spacing, and difference equation. 5–7 Consider three consecutive nodes n 1, n, and n 1 in a plane wall. Using the finite difference form of the first derivative at the midpoints, show that the finite difference form of the second derivative can be expressed as Tn 1 2Tn x2 Tn 1 0 cen58933_ch05.qxd 9/4/2002 11:42 AM Page 315 315 CHAPTER 5 T(x) nodes 0, 1, 2, 3, 4, and 5 with a uniform nodal spacing of x. Using the finite difference form of the first derivative (not the energy balance approach), obtain the finite difference formulation of the boundary nodes for the case of insulation at the left boundary (node 0) and radiation at the right boundary (node 5) with an emissivity of and surrounding temperature of Tsurr. Tn + 1 Tn One-Dimensional Steady Heat Conduction Tn – 1 5–11C Explain how the finite difference form of a heat conduction problem is obtained by the energy balance method. 5–12C In the energy balance formulation of the finite difference method, it is recommended that all heat transfer at the boundaries of the volume element be assumed to be into the volume element even for steady heat conduction. Is this a valid recommendation even though it seems to violate the conservation of energy principle? ∆x ∆x n–1 n n+1 x FIGURE P5–7 5–8 The finite difference formulation of steady twodimensional heat conduction in a medium with heat generation and constant thermal conductivity is given by Tm 1, n 2Tm, n Tm 1, n Tm, n 1 2Tm, n 2 Tm, n 1 2 x y · gm, n k 0 in rectangular coordinates. Modify this relation for the threedimensional case. 5–9 Consider steady one-dimensional heat conduction in a plane wall with variable heat generation and constant thermal conductivity. The nodal network of the medium consists of nodes 0, 1, 2, 3, and 4 with a uniform nodal spacing of x. Using the finite difference form of the first derivative (not the energy balance approach), obtain the finite difference formulation of the boundary nodes for the case of uniform heat flux q·0 at the left boundary (node 0) and convection at the right boundary (node 4) with a convection coefficient of h and an ambient temperature of T . 5–10 Consider steady one-dimensional heat conduction in a plane wall with variable heat generation and constant thermal conductivity. The nodal network of the medium consists of 5–13C How is an insulated boundary handled in the finite difference formulation of a problem? How does a symmetry line differ from an insulated boundary in the finite difference formulation? 5–14C How can a node on an insulated boundary be treated as an interior node in the finite difference formulation of a plane wall? Explain. 5–15C Consider a medium in which the finite difference formulation of a general interior node is given in its simplest form as Tm (a) (b) (c) (d) (e) 2Tm x2 1 Tm 1 · gm k 0 Is heat transfer in this medium steady or transient? Is heat transfer one-, two-, or three-dimensional? Is there heat generation in the medium? Is the nodal spacing constant or variable? Is the thermal conductivity of the medium constant or variable? 5–16 Consider steady heat conduction in a plane wall whose left surface (node 0) is maintained at 30°C while the right surface (node 8) is subjected to a heat flux of 800 W/m2. Express the finite difference formulation of the boundary nodes 0 and 8 30°C Insulation Tsurr Radiation 1 2 FIGURE P5–10 3 4 ∆x ε ∆x 0 W – 800 —2 m No heat generation · g(x) 01 2 3 5 FIGURE P5–16 4 5 6 78 cen58933_ch05.qxd 9/4/2002 11:42 AM Page 316 316 HEAT TRANSFER for the case of no heat generation. Also obtain the finite difference formulation for the rate of heat transfer at the left boundary. 5–17 Consider steady heat conduction in a plane wall with variable heat generation and constant thermal conductivity. The nodal network of the medium consists of nodes 0, 1, 2, 3, and 4 with a uniform nodal spacing of x. Using the energy balance approach, obtain the finite difference formulation of the boundary nodes for the case of uniform heat flux q·0 at the left boundary (node 0) and convection at the right boundary (node 4) with a convection coefficient of h and an ambient temperature of T . 5–18 Consider steady one-dimensional heat conduction in a plane wall with variable heat generation and constant thermal conductivity. The nodal network of the medium consists of nodes 0, 1, 2, 3, 4, and 5 with a uniform nodal spacing of x. Using the energy balance approach, obtain the finite difference formulation of the boundary nodes for the case of insulation at the left boundary (node 0) and radiation at the right boundary (node 5) with an emissivity of and surrounding temperature of Tsurr. 5–19 Consider steady one-dimensional heat conduction in a plane wall with variable heat generation and constant thermal conductivity. The nodal network of the medium consists of nodes 0, 1, 2, 3, 4, and 5 with a uniform nodal spacing of x. The temperature at the right boundary (node 5) is specified. Using the energy balance approach, obtain the finite difference formulation of the boundary node 0 on the left boundary for the case of combined convection, radiation, and heat flux at the left boundary with an emissivity of , convection coefficient of h, ambient temperature of T , surrounding temperature of Tsurr, and uniform heat flux of q·0. Also, obtain the finite difference formulation for the rate of heat transfer at the right boundary. · q0 Tsurr ε Convection T5 · g(x) Radiation ∆x 0 1 2 3 45 boundary (node 0) and radiation at the right boundary (node 2) with an emissivity of and surrounding temperature of Tsurr. 5–21 Consider steady one-dimensional heat conduction in a plane wall with variable heat generation and variable thermal conductivity. The nodal network of the medium consists of nodes 0, 1, and 2 with a uniform nodal spacing of x. Using the energy balance approach, obtain the finite difference formulation of this problem for the case of specified heat flux q·0 to the wall and convection at the left boundary (node 0) with a convection coefficient of h and ambient temperature of T , and radiation at the right boundary (node 2) with an emissivity of and surrounding surface temperature of Tsurr. Tsurr · q0 · g(x) k(T) Convection h T Radiation ε ∆x 0 1 2 FIGURE P5–21 5–22 Consider steady one-dimensional heat conduction in a pin fin of constant diameter D with constant thermal conductivity. The fin is losing heat by convection to the ambient air at T with a heat transfer coefficient of h. The nodal network of the fin consists of nodes 0 (at the base), 1 (in the middle), and 2 (at the fin tip) with a uniform nodal spacing of x. Using the energy balance approach, obtain the finite difference formulation of this problem to determine T1 and T2 for the case of specified temperature at the fin base and negligible heat transfer at the fin tip. All temperatures are in °C. 5–23 Consider steady one-dimensional heat conduction in a pin fin of constant diameter D with constant thermal conductivity. The fin is losing heat by convection to the ambient air at T with a convection coefficient of h, and by radiation to the surrounding surfaces at an average temperature of Tsurr. h, T Tsurr FIGURE P5–19 5–20 Consider steady one-dimensional heat conduction in a composite plane wall consisting of two layers A and B in perfect contact at the interface. The wall involves no heat generation. The nodal network of the medium consists of nodes 0, 1 (at the interface), and 2 with a uniform nodal spacing of x. Using the energy balance approach, obtain the finite difference formulation of this problem for the case of insulation at the left T0 Radiation ε 0 ∆x 1 h, T Convection FIGURE P5–23 D 2 cen58933_ch05.qxd 9/4/2002 11:42 AM Page 317 317 CHAPTER 5 The nodal network of the fin consists of nodes 0 (at the base), 1 (in the middle), and 2 (at the fin tip) with a uniform nodal spacing of x. Using the energy balance approach, obtain the finite difference formulation of this problem to determine T1 and T2 for the case of specified temperature at the fin base and negligible heat transfer at the fin tip. All temperatures are in °C. 5–24 Consider a large uranium plate of thickness 5 cm and thermal conductivity k 28 W/m · °C in which heat is gener· ated uniformly at a constant rate of g 6 105 W/m3. One side of the plate is insulated while the other side is subjected to convection to an environment at 30°C with a heat transfer coefficient of h 60 W/m2 · °C. Considering six equally spaced nodes with a nodal spacing of 1 cm, (a) obtain the finite difference formulation of this problem and (b) determine the nodal temperatures under steady conditions by solving those equations. 5–25 Consider an aluminum alloy fin (k 180 W/m · °C) of triangular cross section whose length is L 5 cm, base thickness is b 1 cm, and width w in the direction normal to the plane of paper is very large. The base of the fin is maintained at a temperature of T0 180°C. The fin is losing heat by con25°C with a heat transfer vection to the ambient air at T coefficient of h 25 W/m2 · °C and by radiation to the surrounding surfaces at an average temperature of Tsurr 290 K. Using the finite difference method with six equally spaced nodes along the fin in the x-direction, determine (a) the temperatures at the nodes and (b) the rate of heat transfer from the fin for w 1 m. Take the emissivity of the fin surface to be 0.9 and assume steady one-dimensional heat transfer in the fin. from 100°C to 200°C. Plot the fin tip temperature and the rate of heat transfer as a function of the fin base temperature, and discuss the results. 5–27 Consider a large plane wall of thickness L 0.4 m, thermal conductivity k 2.3 W/m · °C, and surface area A 20 m2. The left side of the wall is maintained at a constant temperature of 80°C, while the right side loses heat by con15°C with a heat transvection to the surrounding air at T fer coefficient of h 24 W/m2 · °C. Assuming steady onedimensional heat transfer and taking the nodal spacing to be 10 cm, (a) obtain the finite difference formulation for all nodes, (b) determine the nodal temperatures by solving those equations, and (c) evaluate the rate of heat transfer through the wall. 5–28 Consider the base plate of a 800-W household iron having a thickness of L 0.6 cm, base area of A 160 cm2, and thermal conductivity of k 20 W/m · °C. The inner surface of the base plate is subjected to uniform heat flux generated by the resistance heaters inside. When steady operating conditions are reached, the outer surface temperature of the plate is measured to be 85°C. Disregarding any heat loss through the upper part of the iron and taking the nodal spacing to be 0.2 cm, (a) obtain the finite difference formulation for the nodes and (b) determine the inner surface temperature of the plate by Answer: (b) 100°C solving those equations. Insulation Resistance heater, 800 W ∆ x = 0.2 cm Triangular fin w 0 h, T 1 2 3 x 160 cm2 T0 FIGURE P5–28 b 0 b/2 tan θ = —– L θ 1 2 3 ∆x 4 5 x L FIGURE P5–25 5–26 85°C Base plate Reconsider Problem 5–25. Using EES (or other) software, investigate the effect of the fin base temperature on the fin tip temperature and the rate of heat transfer from the fin. Let the temperature at the fin base vary 5–29 Consider a large plane wall of thickness L 0.3 m, thermal conductivity k 2.5 W/m · °C, and surface area A 12 m2. The left side of the wall is subjected to a heat flux · of q 0 700 W/m2 while the temperature at that surface is measured to be T0 60°C. Assuming steady one-dimensional heat transfer and taking the nodal spacing to be 6 cm, (a) obtain the finite difference formulation for the six nodes and (b) determine the temperature of the other surface of the wall by solving those equations. 5–30E A large steel plate having a thickness of L 5 in., thermal conductivity of k 7.2 Btu/h · ft · °F, and an emissivity of 0.6 is lying on the ground. The exposed surface of cen58933_ch05.qxd 9/4/2002 11:42 AM Page 318 318 HEAT TRANSFER Tsurr = 295 K Tsky Radiation ε 6 T = 25°C h = 13 W/m2·°C Convection h, T 0 1 2 Plate 3 4 5 1 in. 5 4 18 cm 3 2 1 3 cm 0 0.6 ft 6 Soil 95° 7 8 FIGURE P5–32 9 10 5–34 x FIGURE P5–30E the plate exchanges heat by convection with the ambient air 80°F with an average heat transfer coefficient of at T h 3.5 Btu/h · ft2 · °F as well as by radiation with the open sky at an equivalent sky temperature of Tsky 510 R. The ground temperature below a certain depth (say, 3 ft) is not affected by the weather conditions outside and remains fairly constant at 50°F at that location. The thermal conductivity of the soil can be taken to be ksoil 0.49 Btu/h · ft · °F, and the steel plate can be assumed to be in perfect contact with the ground. Assuming steady one-dimensional heat transfer and taking the nodal spacings to be 1 in. in the plate and 0.6 ft in the ground, (a) obtain the finite difference formulation for all 11 nodes shown in Figure P5–30E and (b) determine the top and bottom surface temperatures of the plate by solving those equations. 5–31E Repeat Problem 5–30E by disregarding radiation heat transfer from the upper surface. Answers: (b) 78.7°F, 78.4°F 5–32 Consider a stainless steel spoon (k 15.1 W/m · C, 0.6) that is partially immersed in boiling water at 95°C in a kitchen at 25°C. The handle of the spoon has a cross section of about 0.2 cm 1 cm and extends 18 cm in the air from the free surface of the water. The spoon loses heat by convection to the ambient air with an average heat transfer coefficient of h 13 W/m2 · °C as well as by radiation to the surrounding surfaces at an average temperature of Tsurr 295 K. Assuming steady one-dimensional heat transfer along the spoon and taking the nodal spacing to be 3 cm, (a) obtain the finite difference formulation for all nodes, (b) determine the temperature of the tip of the spoon by solving those equations, and (c) determine the rate of heat transfer from the exposed surfaces of the spoon. 5–33 Repeat Problem 5–32 using a nodal spacing of 1.5 cm. Reconsider Problem 5–33. Using EES (or other) software, investigate the effects of the thermal conductivity and the emissivity of the spoon material on the temperature at the spoon tip and the rate of heat transfer from the exposed surfaces of the spoon. Let the thermal conductivity vary from 10 W/m · °C to 400 W/m · °C, and the emissivity from 0.1 to 1.0. Plot the spoon tip temperature and the heat transfer rate as functions of thermal conductivity and emissivity, and discuss the results. 5–35 One side of a 2-m-high and 3-m-wide vertical plate at 130°C is to be cooled by attaching aluminum fins (k 237 W/m · °C) of rectangular profile in an environment at 35°C. The fins are 2 cm long, 0.3 cm thick, and 0.4 cm apart. The heat transfer coefficient between the fins and the surrounding air for combined convection and radiation is estimated to be 30 W/m2 · °C. Assuming steady one-dimensional heat transfer along the fin and taking the nodal spacing to be 0.5 cm, determine (a) the finite difference formulation of this problem, (b) the nodal temperatures along the fin by solving these equations, (c) the rate of heat transfer from a single fin, 130°C T = 35°C 0.4 cm 0.5 cm 0.3 cm 0 1 2 3 x 4 3m 2 cm FIGURE P5–35 cen58933_ch05.qxd 9/4/2002 11:42 AM Page 319 319 CHAPTER 5 and (d) the rate of heat transfer from the entire finned surface of the plate. 10 cm 9.2 cm 5–36 A hot surface at 100°C is to be cooled by attaching 3-cm-long, 0.25-cm-diameter aluminum pin fins (k 237 W/m · °C) with a center-to-center distance of 0.6 cm. The temperature of the surrounding medium is 30°C, and the combined heat transfer coefficient on the surfaces is 35 W/m2 · °C. Assuming steady one-dimensional heat transfer along the fin and taking the nodal spacing to be 0.5 cm, determine (a) the finite difference formulation of this problem, (b) the nodal temperatures along the fin by solving these equations, (c) the rate of heat transfer from a single fin, and (d) the rate of heat transfer from a 1-m 1-m section of the plate. Tair = 8°C 1 cm 1 cm 20 cm Steam 200°C 3 cm FIGURE P5–38 0.6 cm 5–39 Reconsider Problem 5–38. Using EES (or other) software, investigate the effects of the steam temperature and the outer heat transfer coefficient on the flange tip temperature and the rate of heat transfer from the exposed surfaces of the flange. Let the steam temperature vary from 150°C to 300°C and the heat transfer coefficient from 15 W/m2 · °C to 60 W/m2 · °C. Plot the flange tip temperature and the heat transfer rate as functions of steam temperature and heat transfer coefficient, and discuss the results. 0.25 cm 100°C 5–40 0.5 cm 0 1 2 3 4 5 6 x Using EES (or other) software, solve these systems of algebraic equations. (a) 3x1 x2 3x3 x1 2x2 x3 2x1 x2 x3 FIGURE P5–36 5–37 Repeat Problem 5–36 using copper fins (k W/m · °C) instead of aluminum ones. 386 (b) 4x1 Answers: (b) 98.6°C, 97.5°C, 96.7°C, 96.0°C, 95.7°C, 95.5°C 5–38 Two 3-m-long and 0.4-cm-thick cast iron (k 52 W/m · °C, 0.8) steam pipes of outer diameter 10 cm are connected to each other through two 1-cm-thick flanges of outer diameter 20 cm, as shown in the figure. The steam flows inside the pipe at an average temperature of 200°C with a heat transfer coefficient of 180 W/m2 · °C. The outer surface of the pipe is exposed to convection with ambient air at 8°C with a heat transfer coefficient of 25 W/m2 · °C as well as radiation with the surrounding surfaces at an average temperature of Tsurr 290 K. Assuming steady one-dimensional heat conduction along the flanges and taking the nodal spacing to be 1 cm along the flange (a) obtain the finite difference formulation for all nodes, (b) determine the temperature at the tip of the flange by solving those equations, and (c) determine the rate of heat transfer from the exposed surfaces of the flange. x3 Answers: (a) x1 1.62 5–41 (a) 2x2 0.5x3 2 x3 x2 x3 1 x1 x2 x3 2, x2 3, x3 2 11.964 3 1, (b) x1 2.33, x2 2.29, Using EES (or other) software, solve these systems of algebraic equations. 3x1 2x1 (b) 0 3 2 2x2 x3 x4 x1 2x2 x4 x2 3x3 x4 3x2 x3 4x4 6 3x1 x2 2 x2 3x2 1 2x1 x4 2 8 2x3 2x3 4x3 3 2 6 6.293 12 cen58933_ch05.qxd 9/4/2002 11:42 AM Page 320 320 HEAT TRANSFER 5–42 Using EES (or other) software, solve these systems of algebraic equations. 4x1 x2 2x3 x4 x1 3x2 x3 4x4 x1 2x2 5x4 2x2 4x3 3x4 (a) (b) x4 2 2x1 x2 1 4x2 x1 3x1 2x3 x4 2x2 2x4 3 x4 5x3 2 x2 8x4 3 Insulation 6 1 · g = 6 × 106 W/ m3 260 5 1 3 10 15 5 cm 5 cm (a) (b) (c) (d) (e) Tright Tbottom 4Tnode · gnodel 2 k 0 Is heat transfer in this medium steady or transient? Is heat transfer one-, two-, or three-dimensional? Is there heat generation in the medium? Is the nodal spacing constant or variable? Is the thermal conductivity of the medium constant or variable? 5–44C Consider a medium in which the finite difference formulation of a general interior node is given in its simplest form as Tnode (a) (b) (c) (d) (e) (Tleft Ttop Tright 2 325 Convection T = 20°C, h = 50 W/ m2 · °C 5–43C Consider a medium in which the finite difference formulation of a general interior node is given in its simplest form as Ttop 1 200°C Two-Dimensional Steady Heat Conduction Tleft 290 240 350°C 305 3 5 Tbottom)/4 Is heat transfer in this medium steady or transient? Is heat transfer one-, two-, or three-dimensional? Is there heat generation in the medium? Is the nodal spacing constant or variable? Is the thermal conductivity of the medium constant or variable? 5–45C What is an irregular boundary? What is a practical way of handling irregular boundary surfaces with the finite difference method? 5–46 Consider steady two-dimensional heat transfer in a long solid body whose cross section is given in the figure. The temperatures at the selected nodes and the thermal conditions at the boundaries are as shown. The thermal conductivity of the body is k 45 W/m · °C, and heat is generated in the body uni· formly at a rate of g 6 106 W/m3. Using the finite difference method with a mesh size of x y 5.0 cm, determine (a) the temperatures at nodes 1, 2, and 3 and (b) the rate of heat loss from the bottom surface through a 1-m-long section of the body. FIGURE P5–46 5–47 Consider steady two-dimensional heat transfer in a long solid body whose cross section is given in the figure. The measured temperatures at selected points of the outer surfaces are as shown. The thermal conductivity of the body is k 45 W/m · °C, and there is no heat generation. Using the finite difference method with a mesh size of x y 2.0 cm, determine the temperatures at the indicated points in the medium. Hint: Take advantage of symmetry. 150 180 200 150°C 1 180 2 3 180 180 4 5 6 200 200 7 8 9 180 180 2 cm 2 cm 150 180 200 180 150 FIGURE P5–47 5–48 Consider steady two-dimensional heat transfer in a long solid bar whose cross section is given in the figure. The measured temperatures at selected points of the outer surfaces are as shown. The thermal conductivity of the body is k 20 W/m · °C, and there is no heat generation. Using the finite difference method with a mesh size of x y 1.0 cm, determine the temperatures at the indicated points in the medium. Answers: T1 185°C, T2 T3 T4 190°C 5–49 Starting with an energy balance on a volume element, obtain the steady two-dimensional finite difference equation for a general interior node in rectangular coordinates for T(x, y) for the case of variable thermal conductivity and uniform heat generation. cen58933_ch05.qxd 9/4/2002 11:42 AM Page 321 321 CHAPTER 5 200°C W/m3. Plot the temperatures at nodes 1 and 3, and the rate of heat loss as functions of thermal conductivity and heat generation rate, and discuss the results. 1 2 3 4 180 5–52 Consider steady two-dimensional heat transfer in a long solid bar whose cross section is given in the figure. The measured temperatures at selected points on the outer surfaces are as shown. The thermal conductivity of the body is k 20 W/m · °C, and there is no heat generation. Using the finite difference method with a mesh size of x y 1.0 cm, determine the temperatures at the indicated points in the medium. Hint: Take advantage of symmetry. 200 Insulation (a) 120 100 120 1 100°C Answers: (b) T1 143°C, T2 T4 100°C 2 120 x 120 y 150 3 140 1 200 FIGURE P5–48 100 4 6 300 (a) Insulation 100 100 1 100 2 100 3 100°C 4 1 cm 1 cm 200 200 200 200 200°C (b) · g = 107 W/ m3 FIGURE P5–52 1 2 120 120 3 4 150 150 0.1 m 0.1 m Convection 200 200 1 Reconsider Problem 5–50. Using EES (or other) software, investigate the effects of the thermal conductivity and the heat generation rate on the temperatures at nodes 1 and 3, and the rate of heat loss from the top surface. Let the thermal conductivity vary from 10 W/m · °C to 400 W/m · °C and the heat generation rate from 105 W/m3 to 108 · qL 2 4 h, T 3 5 6 1.5 cm 1.5 cm 120°C FIGURE P5–53 7 8 Insulation 200 5–53 Consider steady two-dimensional heat transfer in an L-shaped solid body whose cross section is given in the figure. The thermal conductivity of the body is k 45 W/m · °C, and · heat is generated in the body at a rate of g 5 106 W/m3. FIGURE P5–50 5–51 250 5 Insulation 100°C 100 3 1 cm 300 5–50 Consider steady two-dimensional heat transfer in a long solid body whose cross section is given in the figure. The temperatures at the selected nodes and the thermal conditions on the boundaries are as shown. The thermal conductivity of the body is k 180 W/m · °C, and heat is generated in the body · uniformly at a rate of g 107 W/m3. Using the finite difference method with a mesh size of x y 10 cm, determine (a) the temperatures at nodes 1, 2, 3, and 4 and (b) the rate of heat loss from the top surface through a 1-m-long section of the body. 200 2 250 (b) 200 1 cm Insulation 100 150 4 140 136°C T3 Insulation 180 Insulation 150 cen58933_ch05.qxd 9/4/2002 11:43 AM Page 322 322 HEAT TRANSFER The right surface of the body is insulated, and the bottom surface is maintained at a uniform temperature of 120°C. The entire top surface is subjected to convection with ambient air 30°C with a heat transfer coefficient of h 55 at T W/m2 · °C, and the left surface is subjected to heat flux at a uniform rate of q·L 8000 W/m2. The nodal network of the problem consists of 13 equally spaced nodes with x y 1.5 cm. Five of the nodes are at the bottom surface and thus their temperatures are known. (a) Obtain the finite difference equations at the remaining eight nodes and (b) determine the nodal temperatures by solving those equations. 5–54E Consider steady two-dimensional heat transfer in a long solid bar of square cross section in which heat is gener· ated uniformly at a rate of g 0.19 105 Btu/h · ft3. The cross section of the bar is 0.4 ft 0.4 ft in size, and its thermal conductivity is k 16 Btu/h · ft · °F. All four sides of the bar are 70°F with subjected to convection with the ambient air at T a heat transfer coefficient of h 7.9 Btu/h · ft2 · °F. Using the finite difference method with a mesh size of x y 0.2 ft, determine (a) the temperatures at the nine nodes and (b) the rate of heat loss from the bar through a 1-ft-long section. Answer: (b) 3040 Btu/h h, T 1 2 3 · g 4 5 7 h, T 6 8 h, T 9 h, T FIGURE P5–54E 5–55 Hot combustion gases of a furnace are flowing through a concrete chimney (k 1.4 W/m · °C) of rectangular cross Tsky section. The flow section of the chimney is 20 cm 40 cm, and the thickness of the wall is 10 cm. The average temperature of the hot gases in the chimney is Ti 280°C, and the average convection heat transfer coefficient inside the chimney is hi 75 W/m2 · °C. The chimney is losing heat from its outer surface 15°C by convection with a heat to the ambient air at To 18 W/m2 · °C and to the sky by transfer coefficient of ho radiation. The emissivity of the outer surface of the wall is 0.9, and the effective sky temperature is estimated to be 250 K. Using the finite difference method with x y 10 cm and taking full advantage of symmetry, (a) obtain the finite difference formulation of this problem for steady twodimensional heat transfer, (b) determine the temperatures at the nodal points of a cross section, and (c) evaluate the rate of heat loss for a 1-m-long section of the chimney. 5–56 Repeat Problem 5–55 by disregarding radiation heat transfer from the outer surfaces of the chimney. 5–57 Reconsider Problem 5–55. Using EES (or other) software, investigate the effects of hot-gas temperature and the outer surface emissivity on the temperatures at the outer corner of the wall and the middle of the inner surface of the right wall, and the rate of heat loss. Let the temperature of the hot gases vary from 200°C to 400°C and the emissivity from 0.1 to 1.0. Plot the temperatures and the rate of heat loss as functions of the temperature of the hot gases and the emissivity, and discuss the results. 5–58 Consider a long concrete dam (k 0.6 W/m · °C, 0.7 m2/s) of triangular cross section whose s exposed surface is subjected to solar heat flux of q·s 800 W/m2 and to convection and radiation to the environment at 25°C with a combined heat transfer coefficient of 30 W/m2 · °C. The 2-m-high vertical section of the dam is subjected to convection by water at 15°C with a heat transfer coefficient of 150 W/m2 · °C, and heat transfer through the 2-m-long base is considered to be negligible. Using the finite difference method with a mesh size of x y 1 m and assuming steady two-dimensional heat transfer, determine the temperature of the top, middle, and bottom of the exposed surAnswers: 21.3°C, 43.2°C, 43.6°C face of the dam. ε To, ho Chimney 1 · qs 10 cm h T Hot gases 20 cm Ti, hi 1m 2 3 Water 10 cm 1m 4 10 cm FIGURE P5–55 40 cm 1m 5 1m 6 10 cm FIGURE P5–58 cen58933_ch05.qxd 9/4/2002 11:43 AM Page 323 323 CHAPTER 5 5–59E Consider steady two-dimensional heat transfer in a V-grooved solid body whose cross section is given in the figure. The top surfaces of the groove are maintained at 32°F while the bottom surface is maintained at 212°F. The side surfaces of the groove are insulated. Using the finite difference method with a mesh size of x y 1 ft and taking advantage of symmetry, determine the temperatures at the middle of the insulated surfaces. 32°F 2 4 3 5 9 7 6 10 8 Insulation Insulation 1 11 1 ft 1 ft dimensional heat transfer and (b) determine the unknown nodal temperatures by solving those equations. Answers: (b) 85.7°C, 86.4°C, 87.6°C 5–62 Consider a 5-m-long constantan block (k 23 W/m · °C) 30 cm high and 50 cm wide. The block is completely submerged in iced water at 0°C that is well stirred, and the heat transfer coefficient is so high that the temperatures on both sides of the block can be taken to be 0°C. The bottom surface of the bar is covered with a low-conductivity material so that heat transfer through the bottom surface is negligible. The top surface of the block is heated uniformly by a 6-kW resistance heater. Using the finite difference method with a mesh size of x y 10 cm and taking advantage of symmetry, (a) obtain the finite difference formulation of this problem for steady two-dimensional heat transfer, (b) determine the unknown nodal temperatures by solving those equations, and (c) determine the rate of heat transfer from the block to the iced water. 6 kW heater Insulation 212°F FIGURE P5–59E Reconsider Problem 5–59E. Using EES (or other) software, investigate the effects of the temperatures at the top and bottom surfaces on the temperature in the middle of the insulated surface. Let the temperatures at the top and bottom surfaces vary from 32°F to 212°F. Plot the temperature in the middle of the insulated surface as functions of the temperatures at the top and bottom surfaces, and discuss the results. 0°C 5–60 5–61 Consider a long solid bar whose thermal conductivity is k 12 W/m · °C and whose cross section is given in the figure. The top surface of the bar is maintained at 50°C while the bottom surface is maintained at 120°C. The left surface is insulated and the remaining three surfaces are subjected to 25°C with a heat transfer convection with ambient air at T coefficient of h 30 W/m2 · °C. Using the finite difference method with a mesh size of x y 10 cm, (a) obtain the finite difference formulation of this problem for steady two50°C Insulation Convection h, T 10 cm 10 cm 10 cm 10 cm Insulation FIGURE P5–62 Transient Heat Conduction 5–63C How does the finite difference formulation of a transient heat conduction problem differ from that of a steady heat conduction problem? What does the term i i A xC(Tm 1 Tm)/ t represent in the transient finite difference formulation? 5–64C What are the two basic methods of solution of transient problems based on finite differencing? How do heat transfer terms in the energy balance formulation differ in the two methods? 5–65C The explicit finite difference formulation of a general interior node for transient heat conduction in a plane wall is given by i Tm 120°C FIGURE P5–61 0°C 1 i 2Tm i Tm 1 · g im x2 k T im 1 i Tm Obtain the finite difference formulation for the steady case by simplifying the relation above. cen58933_ch05.qxd 9/4/2002 11:43 AM Page 324 324 HEAT TRANSFER 5–66C The explicit finite difference formulation of a general interior node for transient two-dimensional heat conduction is given by i1 Tnode i (Tleft (1 i Ttop 4 i Tright i )Tnode · g(x, t ) · q0 i Tbottom) ·i gnodel 2 k 5–67C Is there any limitation on the size of the time step t in the solution of transient heat conduction problems using (a) the explicit method and (b) the implicit method? 5–68C Express the general stability criterion for the explicit method of solution of transient heat conduction problems. 5–69C Consider transient one-dimensional heat conduction in a plane wall that is to be solved by the explicit method. If both sides of the wall are at specified temperatures, express the stability criterion for this problem in its simplest form. 5–70C Consider transient one-dimensional heat conduction in a plane wall that is to be solved by the explicit method. If both sides of the wall are subjected to specified heat flux, express the stability criterion for this problem in its simplest form. 5–71C Consider transient two-dimensional heat conduction in a rectangular region that is to be solved by the explicit method. If all boundaries of the region are either insulated or at specified temperatures, express the stability criterion for this problem in its simplest form. 5–72C The implicit method is unconditionally stable and thus any value of time step t can be used in the solution of transient heat conduction problems. To minimize the computation time, someone suggests using a very large value of t since there is no danger of instability. Do you agree with this suggestion? Explain. 5–73 Consider transient heat conduction in a plane wall whose left surface (node 0) is maintained at 50°C while the right surface (node 6) is subjected to a solar heat flux of 600 W/m2. The wall is initially at a uniform temperature of 50°C. Express the explicit finite difference formulation of the boundary nodes 0 and 6 for the case of no heat generation. Also, obtain the finite difference formulation for the total amount of heat transfer at the left boundary during the first three time steps. 5–74 Consider transient heat conduction in a plane wall with variable heat generation and constant thermal conductivity. The nodal network of the medium consists of nodes 0, 1, 2, 3, and 4 with a uniform nodal spacing of x. The wall is initially at a specified temperature. Using the energy balance approach, obtain the explicit finite difference formulation of the boundary · nodes for the case of uniform heat flux q 0 at the left boundary ∆x 0 Obtain the finite difference formulation for the steady case by simplifying the relation above. h, T 1 2 3 x 4 FIGURE P5–74 (node 0) and convection at the right boundary (node 4) with a convection coefficient of h and an ambient temperature of T . Do not simplify. 5–75 tion. Repeat Problem 5–74 for the case of implicit formula- 5–76 Consider transient heat conduction in a plane wall with variable heat generation and constant thermal conductivity. The nodal network of the medium consists of nodes 0, 1, 2, 3, 4, and 5 with a uniform nodal spacing of x. The wall is initially at a specified temperature. Using the energy balance approach, obtain the explicit finite difference formulation of the boundary nodes for the case of insulation at the left boundary (node 0) and radiation at the right boundary (node 5) with an emissivity of and surrounding temperature of Tsurr. 5–77 Consider transient heat conduction in a plane wall with variable heat generation and constant thermal conductivity. The nodal network of the medium consists of nodes 0, 1, 2, 3, and 4 with a uniform nodal spacing of x. The wall is initially at a specified temperature. The temperature at the right boundary (node 4) is specified. Using the energy balance approach, obtain the explicit finite difference formulation of the boundary node 0 for the case of combined convection, radiation, and heat flux at the left boundary with an emissivity of , convection coefficient of h, ambient temperature of T , surrounding temper· ature of Tsurr, and uniform heat flux of q 0 toward the wall. Also, obtain the finite difference formulation for the total amount of heat transfer at the right boundary for the first 20 time steps. Tsurr Radiation TL · g(x, t ) ∆x · q0 0 h, T Convection FIGURE P5–77 L 1 2 3 4 x cen58933_ch05.qxd 9/4/2002 11:43 AM Page 325 325 CHAPTER 5 5–78 Starting with an energy balance on a volume element, obtain the two-dimensional transient explicit finite difference equation for a general interior node in rectangular coordinates for T(x, y, t) for the case of constant thermal conductivity and no heat generation. 5–79 Starting with an energy balance on a volume element, obtain the two-dimensional transient implicit finite difference equation for a general interior node in rectangular coordinates for T(x, y, t) for the case of constant thermal conductivity and no heat generation. 5–80 Starting with an energy balance on a disk volume element, derive the one-dimensional transient explicit finite difference equation for a general interior node for T(z, t) in a cylinder whose side surface is insulated for the case of constant thermal conductivity with uniform heat generation. 5–81 Consider one-dimensional transient heat conduction in a composite plane wall that consists of two layers A and B with perfect contact at the interface. The wall involves no heat generation and initially is at a specified temperature. The nodal network of the medium consists of nodes 0, 1 (at the interface), and 2 with a uniform nodal spacing of x. Using the energy balance approach, obtain the explicit finite difference formulation of this problem for the case of insulation at the left boundary (node 0) and radiation at the right boundary (node 2) with an emissivity of and surrounding temperature of Tsurr. Tsurr Insulation A ε B Radiation ∆x 0 ∆x 1 2 5–84 Consider a large uranium plate of thickness L 8 cm, thermal conductivity k 28 W/m · °C, and thermal diffusivity 12.5 10 6 m2/s that is initially at a uniform temperature of 100°C. Heat is generated uniformly in the plate at a constant · rate of g 106 W/m3. At time t 0, the left side of the plate is insulated while the other side is subjected to convection with 20°C with a heat transfer coefficient of an environment at T h 35 W/m2 · °C. Using the explicit finite difference approach with a uniform nodal spacing of x 2 cm, determine (a) the temperature distribution in the plate after 5 min and (b) how long it will take for steady conditions to be reached in the plate. 5–85 Reconsider Problem 5–84. Using EES (or other) software, investigate the effect of the cooling time on the temperatures of the left and right sides of the plate. Let the time vary from 5 min to 60 min. Plot the temperatures at the left and right surfaces as a function of time, and discuss the results. 5–86 Consider a house whose south wall consists of a 30-cmthick Trombe wall whose thermal conductivity is k 0.70 W/m · °C and whose thermal diffusivity is 0.44 10 6 m2/s. The variations of the ambient temperature Tout and the · solar heat flux q solar incident on a south-facing vertical surface throughout the day for a typical day in February are given in the table in 3-h intervals. The Trombe wall has single glazing with an absorptivity-transmissivity product of 0.76 (that is, 76 percent of the solar energy incident is absorbed by the exposed surface of the Trombe wall), and the average combined heat transfer coefficient for heat loss from the Trombe wall to the ambient is determined to be hout 3.4 W/m2 · °C. 20°C at all The interior of the house is maintained at Tin times, and the heat transfer coefficient at the interior surface of the Trombe wall is hin 9.1 W/m2 · °C. Also, the vents on the Trombe wall are kept closed, and thus the only heat transfer between the air in the house and the Trombe wall is through the Interface FIGURE P5–81 5–82 Consider transient one-dimensional heat conduction in a pin fin of constant diameter D with constant thermal conductivity. The fin is losing heat by convection to the ambient air at T with a heat transfer coefficient of h and by radiation to the surrounding surfaces at an average temperature of Tsurr. The nodal network of the fin consists of nodes 0 (at the base), 1 (in the middle), and 2 (at the fin tip) with a uniform nodal spacing of x. Using the energy balance approach, obtain the explicit finite difference formulation of this problem for the case of a specified temperature at the fin base and negligible heat transfer at the fin tip. 5–83 Repeat Problem 5–82 for the case of implicit formulation. Sun’s rays Trombe wall Heat gain Heat loss Glazing FIGURE P5–86 cen58933_ch05.qxd 9/4/2002 11:43 AM Page 326 326 HEAT TRANSFER 5–88 TABLE P5–86 The hourly variations of the monthly average ambient temperature and solar heat flux incident on a vertical surface Ambient Temperature, °C Time of Day 7 10 1 4 7 10 1 4 AM–10 AM Solar Insolation, W/m2 0 4 6 1 2 3 4 4 AM–1 PM PM–4 PM PM–7 PM PM–10 PM PM–1 AM AM–4 AM AM–7 AM 375 750 580 95 0 0 0 0 Reconsider Problem 5–87. Using EES (or other) software, plot the temperature at the top corner as a function of heating time varies from 2 min to 30 min, and discuss the results. 5–89 Consider a long solid bar (k 28 W/m · °C and 12 10 6 m2/s) of square cross section that is initially at a uniform temperature of 20°C. The cross section of the bar is 20 cm 20 cm in size, and heat is generated in it uniformly at · a rate of g 8 105 W/m3. All four sides of the bar are sub30°C with jected to convection to the ambient air at T a heat transfer coefficient of h 45 W/m2 · °C. Using the explicit finite difference method with a mesh size of x y 10 cm, determine the centerline temperature of the bar (a) after 10 min and (b) after steady conditions are established. h, T 1 interior surface of the wall. Assuming the temperature of the Trombe wall to vary linearly between 20°C at the interior surface and 0°C at the exterior surface at 7 AM and using the explicit finite difference method with a uniform nodal spacing of x 5 cm, determine the temperature distribution along the thickness of the Trombe wall after 6, 12, 18, 24, 30, 36, 42, and 48 hours and plot the results. Also, determine the net amount of heat transferred to the house from the Trombe wall during the first day if the wall is 2.8 m high and 7 m long. 5–87 Consider two-dimensional transient heat transfer in an L-shaped solid bar that is initially at a uniform temperature of 140°C and whose cross section is given in the figure. The thermal conductivity and diffusivity of the body are k 15 W/m · °C and 3.2 10 6 m2/s, respectively, and heat is · generated in the body at a rate of g 2 107 W/m3. The right surface of the body is insulated, and the bottom surface is maintained at a uniform temperature of 140°C at all times. At time t 0, the entire top surface is subjected to convection 25°C with a heat transfer coefficient with ambient air at T of h 80 W/m2 · °C, and the left surface is subjected to · uniform heat flux at a rate of q L 8000 W/m2. The nodal network of the problem consists of 13 equally spaced nodes with x y 1.5 cm. Using the explicit method, determine the temperature at the top corner (node 3) of the body after 2, 5, and 30 min. Convection · qL 2 4 h, T 3 5 6 1.5 cm 1.5 cm 7 5 10 cm 10 cm 7 8 6 h, T 9 h, T FIGURE P5–89 5–90E Consider a house whose windows are made of 0.375-in.-thick glass (k 0.48 Btu/h · ft · °F and 4.2 10 6 ft2/s). Initially, the entire house, including the walls and the windows, is at the outdoor temperature of To 35°F. It is observed that the windows are fogged because the indoor temperature is below the dew-point temperature of 54°F. Now the heater is turned on and the air temperature in the house is raised to Ti 72°F at a rate of 2°F rise per minute. The heat transfer coefficients at the inner and outer surfaces of the wall can be taken to be hi 1.2 and ho 2.6 Btu/h · ft2 · °F, respectively, and the outdoor temperature can be assumed to remain constant. Using the explicit finite difference method with a mesh size of x 0.125 in., determine how long it will take Window glass ∆x 12 House Fog 140°C FIGURE P5–87 4 h, T hi 8 3 · g Ti Insulation 1 2 FIGURE P5–90E 34 0.375 in To ho Outdoors cen58933_ch05.qxd 9/4/2002 11:43 AM Page 327 327 CHAPTER 5 for the fog on the windows to clear up (i.e., for the inner surface temperature of the window glass to reach 54°F). 5–91 A common annoyance in cars in winter months is the formation of fog on the glass surfaces that blocks the view. A practical way of solving this problem is to blow hot air or to attach electric resistance heaters to the inner surfaces. Consider the rear window of a car that consists of a 0.4-cm-thick glass (k 0.84 W/m · °C and 0.39 10 6 m2/s). Strip heater wires of negligible thickness are attached to the inner surface of the glass, 4 cm apart. Each wire generates heat at a rate of 10 W/m length. Initially the entire car, including its windows, 3°C. The heat transfer is at the outdoor temperature of To coefficients at the inner and outer surfaces of the glass can be taken to be hi 6 and ho 20 W/m2 · °C, respectively. Using the explicit finite difference method with a mesh size of x 0.2 cm along the thickness and y 1 cm in the direction normal to the heater wires, determine the temperature distribution throughout the glass 15 min after the strip heaters are turned on. Also, determine the temperature distribution when steady conditions are reached. Thermal symmetry line Inner surface Outer surface Glass Heater 10 W/ m 0.2 cm 1 cm Thermal symmetry line Tsky Radiation Convection ho, To ε 15 cm ε Convection hi, Ti Radiation Ti FIGURE P5–93 are maintained at a constant temperature of 20°C during the night, and the emissivity of both surfaces of the concrete roof is 0.9. Considering both radiation and convection heat transfers and using the explicit finite difference method with a time step of t 5 min and a mesh size of x 3 cm, determine the temperatures of the inner and outer surfaces of the roof at 6 AM. Also, determine the average rate of heat transfer through the roof during that night. 5–94 Consider a refrigerator whose outer dimensions are 1.80 m 0.8 m 0.7 m. The walls of the refrigerator are constructed of 3-cm-thick urethane insulation (k 0.026 W/m · ° C and 0.36 10 6 m2/s) sandwiched between two layers of sheet metal with negligible thickness. The refrigerated space is maintained at 3°C and the average heat transfer coefficients at the inner and outer surfaces of the wall are 6 W/m2 · °C and 9 W/m2 · °C, respectively. Heat transfer through the bottom surface of the refrigerator is negligible. The kitchen temperature remains constant at about 25°C. Initially, the refrigerator contains 15 kg of food items at an average specific heat of 3.6 kJ/kg · °C. Now a malfunction occurs and the refrigerator stops running for 6 h as a result. Assuming the FIGURE P5–91 5–92 Refrigerator Air Repeat Problem 5–91 using the implicit method with a time step of 1 min. 5–93 The roof of a house consists of a 15-cm-thick concrete slab (k 1.4 W/m · °C and 0.69 10 6 m2/s) that is 20 m wide and 20 m long. One evening at 6 PM, the slab is observed to be at a uniform temperature of 18°C. The average ambient air and the night sky temperatures for the entire night are predicted to be 6°C and 260 K, respectively. The convection heat transfer coefficients at the inner and outer surfaces of the roof can be taken to be hi 5 and ho 12 W/m2 · °C, respectively. The house and the interior surfaces of the walls and the floor Concrete roof 3°C Heat 25°C ho hi Food 15 kg FIGURE P5–94 cen58933_ch05.qxd 9/4/2002 11:43 AM Page 328 328 HEAT TRANSFER temperature of the contents of the refrigerator, including the air inside, rises uniformly during this period, predict the temperature inside the refrigerator after 6 h when the repairman arrives. Use the explicit finite difference method with a time step of t 1 min and a mesh size of x 1 cm and disregard corner effects (i.e., assume one-dimensional heat transfer in the walls). 5–95 Reconsider Problem 5–94. Using EES (or other) software, plot the temperature inside the refrigerator as a function of heating time as time varies from 1 h to 10 h, and discuss the results. Special Topic: Controlling the Numerical Error coordinates for T(x, y, z, t) for the case of constant thermal conductivity and no heat generation. 5–108 Consider steady one-dimensional heat conduction in a plane wall with variable heat generation and constant thermal conductivity. The nodal network of the medium consists of nodes 0, 1, 2, and 3 with a uniform nodal spacing of x. The temperature at the left boundary (node 0) is specified. Using the energy balance approach, obtain the finite difference formulation of boundary node 3 at the right boundary for the case of combined convection and radiation with an emissivity of , convection coefficient of h, ambient temperature of T , and surrounding temperature of Tsurr. Also, obtain the finite difference formulation for the rate of heat transfer at the left boundary. 5–96C Why do the results obtained using a numerical method differ from the exact results obtained analytically? What are the causes of this difference? 5–97C What is the cause of the discretization error? How does the global discretization error differ from the local discretization error? Tsurr Radiation 5–98C Can the global (accumulated) discretization error be less than the local error during a step? Explain. 5–99C How is the finite difference formulation for the first derivative related to the Taylor series expansion of the solution function? 5–100C Explain why the local discretization error of the finite difference method is proportional to the square of the step size. Also explain why the global discretization error is proportional to the step size itself. 5–101C What causes the round-off error? What kind of calculations are most susceptible to round-off error? 5–102C What happens to the discretization and the roundoff errors as the step size is decreased? Suggest some practical ways of reducing the round-off error. 5–103C 5–104C What is a practical way of checking if the round-off error has been significant in calculations? 5–105C What is a practical way of checking if the discretization error has been significant in calculations? ε · g(x) T0 h, T Convection ∆x 0 1 2 3 FIGURE P5–108 5–109 Consider one-dimensional transient heat conduction in a plane wall with variable heat generation and variable thermal conductivity. The nodal network of the medium consists of nodes 0, 1, and 2 with a uniform nodal spacing of x. Using the energy balance approach, obtain the explicit finite difference formulation of this problem for the case of specified heat · flux q 0 and convection at the left boundary (node 0) with a convection coefficient of h and ambient temperature of T , and radiation at the right boundary (node 2) with an emissivity of and surrounding temperature of Tsurr. 5–110 Repeat Problem 5–109 for the case of implicit formulation. 5–111 Consider steady one-dimensional heat conduction in a pin fin of constant diameter D with constant thermal conductivity. The fin is losing heat by convection with the ambient air Review Problems 5–106 Starting with an energy balance on the volume element, obtain the steady three-dimensional finite difference equation for a general interior node in rectangular coordinates for T(x, y, z) for the case of constant thermal conductivity and uniform heat generation. 5–107 Starting with an energy balance on the volume element, obtain the three-dimensional transient explicit finite difference equation for a general interior node in rectangular Convection h, T ε 0 1 ∆x FIGURE P5–111 Tsurr Radiation 2 cen58933_ch05.qxd 9/4/2002 11:43 AM Page 329 329 CHAPTER 5 at T (in °C) with a convection coefficient of h, and by radiation to the surrounding surfaces at an average temperature of Tsurr (in K). The nodal network of the fin consists of nodes 0 (at the base), 1 (in the middle), and 2 (at the fin tip) with a uniform nodal spacing of x. Using the energy balance approach, obtain the finite difference formulation of this problem for the case of a specified temperature at the fin base and convection and radiation heat transfer at the fin tip. 5–112 Starting with an energy balance on the volume element, obtain the two-dimensional transient explicit finite difference equation for a general interior node in rectangular coordinates for T(x, y, t) for the case of constant thermal conductivity and uniform heat generation. 5–113 Starting with an energy balance on a disk volume element, derive the one-dimensional transient implicit finite difference equation for a general interior node for T(z, t) in a cylinder whose side surface is subjected to convection with a convection coefficient of h and an ambient temperature of T for the case of constant thermal conductivity with uniform heat generation. 5–114E The roof of a house consists of a 5-in.-thick concrete slab (k 0.81 Btu/h · ft · °F and 7.4 10 6 ft2/s) that is 45 ft wide and 55 ft long. One evening at 6 PM, the slab is observed to be at a uniform temperature of 70°F. The ambient air temperature is predicted to be at about 50°F from 6 PM to 10 PM, 42°F from 10 PM to 2 AM, and 38°F from 2 AM to 6 AM, while the night sky temperature is expected to be about 445 R for the entire night. The convection heat transfer coefficients at the inner and outer surfaces of the roof can be taken to be hi 0.9 and ho 2.1 Btu/h · ft2 · °F, respectively. The house and the interior surfaces of the walls and the floor are maintained at a constant temperature of 70°F during the night, and the emissivity of both surfaces of the concrete roof is 0.9. Considering both radiation and convection heat transfers and using the explicit finite difference method with a mesh size of x 1 in. and a time step of t 5 min, determine the temperatures of the inner and outer surfaces of the roof at 6 AM. Also, determine the average rate of heat transfer through the roof during that night. 5–115 Solar radiation incident on a large body of clean water (k 0.61 W/m · °C and 0.15 10 6 m2/s) such as a lake, a river, or a pond is mostly absorbed by water, and the amount of absorption varies with depth. For solar radiation incident at a 45° angle on a 1-m-deep large pond whose bottom surface is black (zero reflectivity), for example, 2.8 percent of the solar energy is reflected back to the atmosphere, 37.9 percent is absorbed by the bottom surface, and the remaining 59.3 percent is absorbed by the water body. If the pond is considered to be four layers of equal thickness (0.25 m in this case), it can be shown that 47.3 percent of the incident solar energy is absorbed by the top layer, 6.1 percent by the upper mid layer, 3.6 percent by the lower mid layer, and 2.4 percent by the bottom layer [for more information see Çengel and Özisik, Solar En¸ ergy, 33, no. 6 (1984), pp. 581–591]. The radiation absorbed by the water can be treated conveniently as heat generation in the heat transfer analysis of the pond. Consider a large 1-m-deep pond that is initially at a uniform temperature of 15°C throughout. Solar energy is incident on the pond surface at 45° at an average rate of 500 W/m2 for a period of 4 h. Assuming no convection currents in the water and using the explicit finite difference method with a mesh size of x 0.25 m and a time step of t 15 min, determine the temperature distribution in the pond under the most favorable conditions (i.e., no heat losses from the top or bottom surfaces of the pond). The solar energy absorbed by the bottom surface of the pond can be treated as a heat flux to the water at that surface in this case. Sun Tsky Radiation Convection ho, To ε Solar radiation Concrete roof 45° 0 1 ε Radiation Convection hi, Ti 2 3 Ti 4 Top layer Solar pond Upper mid layer Lower mid layer Bottom layer x FIGURE P5–114E · qs, W/ m2 FIGURE P5–115 Black cen58933_ch05.qxd 9/4/2002 11:43 AM Page 330 330 HEAT TRANSFER 5–116 Reconsider Problem 5–115. The absorption of solar radiation in that case can be expressed more accurately as a fourth-degree polynomial as · g(x) · q s(0.859 3.415x 6.704x2 6.339x3 2.278x4), W/m3 · where q s is the solar flux incident on the surface of the pond in 2 W/m and x is the distance from the free surface of the pond in m. Solve Problem 5–115 using this relation for the absorption of solar radiation. 5–117 A hot surface at 120°C is to be cooled by attaching 8 cm long, 0.8 cm in diameter aluminum pin fins (k 237 W/m · °C and 97.1 10 6 m2/s) to it with a center-tocenter distance of 1.6 cm. The temperature of the surrounding medium is 15°C, and the heat transfer coefficient on the surfaces is 35 W/m2 · °C. Initially, the fins are at a uniform temperature of 30°C, and at time t 0, the temperature of the hot surface is raised to 120°C. Assuming one-dimensional heat conduction along the fin and taking the nodal spacing to be x 2 cm and a time step to be t 0.5 s, determine the nodal temperatures after 5 min by using the explicit finite difference method. Also, determine how long it will take for steady conditions to be reached. T0 Convection h, T 2 cm 0 1 2 3 4 FIGURE P5–117 5–118E Consider a large plane wall of thickness L 0.3 ft and thermal conductivity k 1.2 Btu/h · ft · °F in space. The wall is covered with a material having an emissivity of 0.80 and a solar absorptivity of s 0.45. The inner surface of the wall is maintained at 520 R at all times, while the outer surface is exposed to solar radiation that is incident at a · rate of q s 300 Btu/h · ft2. The outer surface is also losing heat αs T0 by radiation to deep space at 0 R. Using a uniform nodal spacing of x 0.1 ft, (a) obtain the finite difference formulation for steady one-dimensional heat conduction and (b) determine the nodal temperatures by solving those equations. Answers: (b) 522 R, 525 R, 527 R 5–119 Frozen food items can be defrosted by simply leaving them on the counter, but it takes too long. The process can be speeded up considerably for flat items such as steaks by placing them on a large piece of highly conducting metal, called the defrosting plate, which serves as a fin. The increased surface area enhances heat transfer and thus reduces the defrosting time. Consider two 1.5-cm-thick frozen steaks at 18°C that resemble a 15-cm-diameter circular object when placed next to each other. The steaks are now placed on a 1-cm-thick blackanodized circular aluminum defrosting plate (k 237 0.90) whose outer W/m · °C, 97.1 10 6 m2/s, and diameter is 30 cm. The properties of the frozen steaks are 1.55 kJ/kg · °C, k 1.40 W/m · °C, 970 kg/m3, Cp 0.95, and the heat of fusion is 0.93 10 6 m2/s, and hif 187 kJ/kg. The steaks can be considered to be defrosted when their average temperature is 0°C and all of the ice in the steaks is melted. Initially, the defrosting plate is at the room temperature of 20°C, and the wooden countertop it is placed on can be treated as insulation. Also, the surrounding surfaces can be taken to be at the same temperature as the ambient air, and the convection heat transfer coefficient for all exposed surfaces can be taken to be 12 W/m2 · °C. Heat transfer from the lateral surfaces of the steaks and the defrosting plate can be neglected. Assuming one-dimensional heat conduction in both the steaks and the defrosting plate and using the explicit finite difference method, determine how long it will take to defrost the steaks. Use four nodes with a nodal spacing of x 0.5 cm for the steaks, and three nodes with a nodal spacing of r 3.75 cm for the exposed portion of the defrosting plate. Also, use a time step of t 5 s. Hint: First, determine the total amount of heat transfer needed to defrost the steaks, and then determine how long it will take to transfer that much heat. Symmetry line 20°C Heat transfer · qs 6 ∆x 0 1 2 3 ε 5 1 2 3 4 Tsurr Defrosting plate Radiation FIGURE P5–118E Steaks FIGURE P5–119 cen58933_ch05.qxd 9/4/2002 11:43 AM Page 331 331 CHAPTER 5 5–120 Repeat Problem 5–119 for a copper defrosting plate using a time step of t 3 s. Design and Essay Problems 5–121 Write a two-page essay on the finite element method, and explain why it is used in most commercial engineering software packages. Also explain how it compares to the finite difference method. 5–122 Numerous professional software packages are available in the market for performing heat transfer analysis, and they are widely advertised in professional magazines such as the Mechanical Engineering magazine published by the American Society of Mechanical Engineers (ASME). Your company decides to purchase such a software package and asks you to prepare a report on the available packages, their costs, capabilities, ease of use, and compatibility with the available hardware, and other software as well as the reputation of the software company, their history, financial health, customer support, training, and future prospects, among other things. After a preliminary investigation, select the top three packages and prepare a full report on them. 5–123 Design a defrosting plate to speed up defrosting of flat food items such as frozen steaks and packaged vegetables and evaluate its performance using the finite difference method (see Prob. 5–119). Compare your design to the defrosting plates currently available on the market. The plate must perform well, and it must be suitable for purchase and use as a household utensil, durable, easy to clean, easy to manufacture, and affordable. The frozen food is expected to be at an initial temperature of 18°C at the beginning of the thawing process and 0°C at the end with all the ice melted. Specify the material, shape, size, and thickness of the proposed plate. Justify your recommendations by calculations. Take the ambient and surrounding surface temperatures to be 20°C and the convection heat transfer coefficient to be 15 W/m2 · °C in your analysis. For a typical case, determine the defrosting time with and without the plate. 5–124 Design a fire-resistant safety box whose outer dimensions are 0.5 m 0.5 m 0.5 m that will protect its combustible contents from fire which may last up to 2 h. Assume the box will be exposed to an environment at an average temperature of 700°C with a combined heat transfer coefficient of 70 W/m2 · °C and the temperature inside the box must be below 150°C at the end of 2 h. The cavity of the box must be as large as possible while meeting the design constraints, and the insulation material selected must withstand the high temperatures to which it will be exposed. Cost, durability, and strength are also important considerations in the selection of insulation materials. cen58933_ch05.qxd 9/4/2002 11:43 AM Page 332 ...
View Full Document

Ask a homework question - tutors are online