Back to Inverse Heat Conduction: Ill-Posed Problems Home
Table of Contents
SEC. 3 .1 C HAPTER 3 A PPROXIMATE M ETHODS FOR DIRECT H EAT C ONDUCTION P ROBLEMS 3.1 I NTRODUCTION In this chapter some approximate methods are given for solving direct transient heat conduction problems. These direct techniques are the first stage of solution procedures for the inverse heat conduction problems. The second stage, which involves specific algorithms for the I HCP, is developed in Chapter 4. The approximate heat conduction solution procedures o f this chapter are combined in Chapters 5 and 6 with the I HCP algorithms of Chapter 4. 3.1.1 78 heat conduction problems are those using Duhamel's theorem and Green's functions both of which yield convolution-type integrals. Since these methods are very similar, only Duhamel's theorem is given. A numerical approximation of the integral is also given. The heat-conducting body can be of arbitrary shape and of nonhomogeneous material. Although the heat transfer may be two- or three-dimensional, in this chapter there is a single forcing input which is only time dependent. A new technique which utilizes Duhamel's theorem for solving connected basic geometries is called the unsteady surface element method.! - ~ I t also involves solving an inverse problem but there are no measurement errors and the matching conditions are a t surfaces rather than at interior points o f the bodies. The solution of a nonlinear transient heat conduction problem requires an approach that discretizes the partial differential equation. Two methods for solving the nonlinear transient heat conduction equation are the finite difference and finite element (FE) methods. An alternative method is used in this chapter in which conservation o f energy is applied directly to finite control volumes; it is called the finite control volume procedure. A nonlinear problem is one for which the heat conduction equation o r a boundary condition is a nonlinear function o f temperature. O ne example of a nonlinear partial differential equation is: _~ ax (k aT)=pc aatT ax (3.1.1 ) where k a nd/or pc are functions of temperature T. I f k and pc are functions of position, x, only, then Eq. (3 .1.1) is linear. An example o f a nonlinear boundary condition is the radiation condition, a - k - TI =u£[T!,-T4(O, t)] a x .. = 0 V arious N umerical A pproaches Partial differential equations including the transient heat conduction equation can be solved using a variety of methods including exact and numerical procedures. The exact methods include the classical methods of separation of variables and Laplace transforms. Two types o f numerical procedures are discussed in this chapter. O ne is based on an integral formulation of the mathematical model and the other on a differential form o f the model. The transient heat conduction equation can be either linear o r nonlinear. F or the linear case the partial differential equation formulation can be equivalently presented by an integral equation. The advantage of the integral form is that the space dependence can be represented exactly leaving only the time dependence to be approximated. Two different approaches that employ integral equations for solving transient 79 I NTRODUCTION 3.1.2 (3.1.2) S cope o f C hapter This chapter gives a derivation of Duhamel's integral (Section 3.2), and a numerical approximation thereof. Section 3.3 discusses finite control volume solutions of the transient heat conduction equation. Both linear and nonlinear cases are included. Section 3.3.5 demonstrates that the linear finite control volume approach can yield expressions that are analogous to those obtained from the Duhamel's theorem method. This is important because the same linear I HCP algorithms can be employed using many different methods of approximating the heat conduction equation and boundary conditions; such methods include a numerical form o f Duhamel's theorem, finite differences, and finite elements. CHAP.3 80 3 .2 3.2.1 M ETHOOS FOR DIRECT HEAT C ONDUCTION PROBLEMS D UHAMEL'S T HEOREM SEC. 3 .2 D UHAMEL'S T HEOREM 81 are used to represent the heat fluxes between times of, respectively, OtOAt, AI t OAz,···, AM-I to 'M D erivation o f D uhamel's T heorem The corresponding heat fluxes are denoted q l' qz, . .. , qM o r in general Duhamel's theorem can be considered to be a result of the principle of superposition and thus it is only valid for linear cases. There are several ways of deriving it, one of which uses Laplace transforms. See Ozisik (Reference 6, p. 197) and Luikov (Reference 7, p. 344). Another derivation uses the concept o f superposition; that is the method used in Arpaci (Reference 8, p. 307) and Myers (Reference 9, p. 155). T he following derivation employs the principle of superposition and it also produces a convenient numerical approximation. Duhamel's theorem employs a "building block" solution which is used with the suPerposition principle to obtain the temperature a t any point and time. One such solution is </J(r, t) which is for the temperature rise at a point r in a heat-conducting body due to the heat flux, O, t<O q(t)= { 1, t >O (3.2.1) This is sometimes called a unit step heat flux o r even a constant unit heat flux. The thermal properties of the body are independent of temperature but can be functions of position, and the temperature distribution need not be onedimensional. F or the derivation herein the only requirement is t hat the prescribed surface heat flux is given by the product, q(s, t) = f(t)g(s) where s is a coordinate measured along the surface. For convenience, the heat flux is written as a function o f time only. The surface heat flux is approximated in the manner shown in Figure 3.1. The heat fluxes a t times (3.2.2) The initial temperature distribution in the body is the constant value of To. Then using the principle o f superposition, the temperature a t location r, at time t M , is composed o f contributions due to the heat flux components ql through qM inclusive: T(r, t M)= TO +ql[cp(r, tM -Ao)-</J(r, tM - AI)] +qz[</J(r, t M-Ad-cp(r, t M-AZ)] (3.2.3) +qM[</J(r, tM - AM- d-</J(r, tM - AM)] where </J(r, tM - AM) = </J(r, 0 )=0. Using equal time steps and times o f Aj=i.:1A, ; =0,1, . .. , M (3.2.4) A j- Ai = j.:1A- i.:1A=(j - ;).:1A=Aj-i (3 .2.5) yields: Then Eq. (3.2.3) can be written as (3.2.6) For the limiting case o f .:1A-+O a nd for a fixed time tM = t, Eq. (3.2.6) becomes the integral (3.2.7) From the relation q (t) a</J(r, t - A) a</J(r, 1 - A) aA at (3.2.8) Duhamel's theorem (sometimes called Duhamel's integral) becomes T(r, 1)=To+ f~ q(A) a</J(r~:-A) dA ~~-...4-:'-----'----:-'-_~£-.-1 1M o Al A2 A M_1 FIGURE 3.1 Approximate representation o f q (t) by M steps in heat Hux . (3.2.9) Equation (3.2.9) is a heat flux form o f Duhamel's theorem; it is a convolution because there is a product o f two functions, one of A and the other o f t -A; that is, there is folding o r convoluting o f one function with respect to the other. For a temperature-based form of Duhamel's theorem, see Eq. (2.4.1) a nd References 6 - 10. C HAP.3 82 3 .2.2 M ETHODS FOR DIRECT H EAT C ONDUCTION PROBLEMS SEC. 3.2 D UHAMEL'S THEOREM Then the T+ values at , ; = 0.1, 0.2, a nd 0.3 a re used to obtain the ~tPl values, N umerical A pproximation o f D uhamel's T heorem x 0.0\ + ~l/Io=l/I, = k 1/11 = 40 (0.003943)=9.8575E-7 C-m1/W A numerical approximation for Eq. (3.2.9) is o btained from Eq. (3.2.6). Note that <p(r, tM - A.._ d - <p(r, tM - A..) = <p(r, t M-.+ d - <p(r, t M-.)=d<p(r, t M-.) 83 X (3.2.10) ~I/I, =1/12-1/1, = k (tP; -1/Ii>=6.69725E-6 C-m 2/W a nd thus Eq. (3.2.6) c an be written as ~1/I2 =tPl-1/I2 = i (tPj - tP;) = 1.029025E-5C-m z/W M L T (r, t M )= T o+ q.d<p{r, t M -.) (3.2.11) n= 1 T he q values are evaluated at 1= 0.5, 1.5, a nd 2.5 s so that T he notation can be simplified by omitting the r dependence and writing q l=0.5E6, qz=1.5E6, q l=2 .5E6W/m z and thus the temperatures found from Eqs. (a), (b), a nd (c) a re M T M=To+ L q .d<PM-., d<Pi=<Pi+l-<P, (3.2.12) TI = 30 + 0.5E6(9 .8575E-7) = 3O.493°C n =1 where q. is best evaluated a t time ( n-t)dt as indicated by Eq. (3.2.2). I f t he actual heat flux is constant over each time step, Eq. (3 .2.12) gives an exact expression for the temperature, TM ; o n the other hand, Eq. (3.2.12) yields approximate results i fthe t rue heat flux varies over the time steps. This equation is quite important in investigation o f the I HCP because it gives a convenient expression for the temperature in terms o f t he heat flux components. The observation that the sum o f t he subscripts on the right side o f Eq. (3.2.12) is M can aid in memorizing this equation. Tz = 30 +0.5E6(6.69725E-6) + 1.5E6(9.8575E-7) = 34.827°C Tl = To + 0.5E6( 1.029025E-5) + 1.5E6(6.69725E-6) + 2.5E6(9.8575E-7) = 47.655°C T he exact temperatures are obtained from Carslaw and Jaeger (Reference 10, p. 77), T = To+(:~)(4a1)3/2 i l erfc[(41;)-I/ Z ] (3.2.13) where EXAMPLE 3.1. Calculate the temperature at 1 cm inside a thick steel plate ( k= 40 W/m-C, a = 1 0- 5 m1/s) that is initially a t 30°C with an applied heat flux o f (3.2.14) The temperatures T" T2 , a nd Tl a re obtained from Eq. (3.2.13) at times 1 = 1 ,1=2, a nd 1= 3 s o r equivalently for I ; = 0.1, 0.2, a nd 0.3. T he resulting exact values are q =4 ot, 4 o=106 W/m 2-s, t ins starting a t time zero. Find values of the temperature a t times 1,2, and 3 s. Use Eq. (3.2.12) to approximate each temperature and compare it with the exact value. Solution. T he temperatures arc denoted TJ> T2 , and Tl a nd correspond t o 1 = 1, 2, a nd 3 s; from Eq. (3 .2.12) they are given by TI = To+ql~</>o (a) Tl = To+ql~</>I +q2~</>O (b) Tl=To+q'~</>2+q2~</>' +ql~tPO (c) Hence the errors in the approximate temperatures are + 0.305,0.753, a nd 0.930, respectively. These errors are increasing, but relative to the exact temperature rises they are 0 decreasing; namely, 163%, 18.5%, and 5.6%, respectively. T he tPi values are for a semi-infinite body with a unit step increase in surface heat flux; tPi values can be found using entries in Table 1.2 for T+ defined by Eq. (1.6.18c) and letting q c= 1, T(x, t)=</>, and To=O; hence t P=(:) T+ =e~l) T+ Also the time 1 is related to , ; by Eq. (1.6.18d) which gives 10- 5 1 = -=--=0.11 2 •x 10- 4 + al 3.2.3 M atrix F orm o f D uhamel's T heorem I t is frequently advantageous to perform algebraic manipulations utilizing a matrix form o f t he model. This subsection displays a matrix form for the numerical statement o f D uhamel's theorem. An expansion o f t he constant heat flux numerical approximation o f D uhamel's integral, Eq. (3.2.12), with M replaced by 1 to M + r-1, is 1'. = To +q1d<po T l = TO+q1d<Pl +q1d<po 84 C HAP.3 M ETHODS FOR DIRECT HEAT C ONDUCTION PROBLEMS T",+ 1 = To+qIAtP",+q2AtP"'-1 + . .. +q",AtPl +q",+ ,AtPo T2 TJ AtPo AtPl AtP2 AtPo AtP, AtP",_, AtP", AtP"'-2 AtP",_, AtP"'-J AtP"'-l T"'+'_I AtP"'+,-2 AtP"'+,-J AtP",+,-4 q", (3.2.16) q"'+1 where 1 is a vector o f ones. I n a more compact form Eq. (3.2.16) is (3.2.17) T =Xq+Tol where T and q are the appropriate vectors displayed in Eq. (3.2.16) a nd X is the lower triangular matrix. X= AtPo AtP"'-1 AtP", AtP"'-2 AtP",_, AtP",- 3 AtP"'-2 (3.2.19) 1 X=[~:: AtP,-1 [q",] q = ~"'+I AtPo AtP,-2 .0. _ 1 1 ;Oq . .J AtPo AtP, AtPo AtP'-1 AtPr-2" . AtPo (3.2.18a) The X matrix is called the pulse sensitivity coefficient matrix for q. Its structure is lower triangular with AtPo's along the main diagonal. AtP,'s along the diagonal t.",+d''';''''I;o [ (3 .2.20a, b) qU+,-1 t"'I''';O AtP",+,-2 AtP"'+, - J AtP",+,-4 (3.2.18b) This is called the s tandard;{orm for the linear I HCP . This equation is n ot restricted to representing a numerical form o f Duhamel's Theorem. I t can be derived using Taylor series for the linear inverse heat conduction problem. The various components o f Eq. (3.2.19) are T'" T = : "'+1 • [ 1''''+'-1 q"'+,-I AtPo AtP, i <j T =Xq+1'lq;o + Tol AtPo AtP, AtPl O. oqi f e yv1 iJ <2 rO+GI .,.-e. ql ql q3 X . I} which shows that X i) depends only o n the difference between i a ndj for linear problems. In the sequential I HCP methods investigated in C hapter 4. the heat flux components. q l. q2 • ... , q", - I. are considered previously estimated and are denoted iiI>' ..• ii", - " Estimates of q"" q"'+I • ...• q"'+, _1 are needed. Moreover. there may be a known heat flux at x = L , a known convection condition a t x =L. o r known internal volumetric heating. F or these cases it is convenient to use the modification of Eq. (3.2.17). AtPo T", T",+, x .. =01i={AtPi-i. i~j (3.2.15) which can be written in expanded matrix form as 'Ii 85 D UHAMEL'S THEOREM just below the main one. and so on. The components of the X matrix are called sensitivity coefficients because they are equal to the first derivatives of the temperature with respect t o h eat flux components. From Eq. (3.2.15). one can verify that the ijth component o f X is T",=To+qIAtP"'-1 +q2 AtP"' - 2+ '" +q"' - IAtPl +q",AtPo T",+,_, = To+q IAtP",+,-2 + ... +Q",+,-2AtPl +q",+,_,AtPo S EC.3.2 (3.2.21) ] (3.2.22) f",+r - d,.,;,.,.I ; '" ; ,., •• -1 = 0 Equations (3.2.19) - (3.2.22) apply for the numerical convolution as well as for the finite difference methods. F or the case of known heating o r energy sources other than where q is to be estimated. X and Tlq ; o must be properly interpreted. The quantity tPi needed in X is the temperature rise at the sensor location for a unit step change in the surface heat flux a t t = 0 for the same partial differential equation and boundary conditions as the original problem except the differential equation and boundary conditions (other than at x=O) are homo- 86 C HAP.3 M ETHODS FOR DIRECT HEAT C ONDUCTION PROBLEMS geneous. An example is the problem for a hollow cylinder, ~~ r ar (r O T)=PC a T _ g(r, t ), or at k aTI ar -k - =q(t), a <r<b unknown heat flux (3 .2.23) (3.2.24a) r=G aTI -kar =h[T(b,t)-T,., (t)] (3.2.24b) r =b (3 .2.25) T (r, 0) = F(r) S EC.3.3 DIFFERENCE M ETHODS 87 I t is a ppropriate to compare Eqs. (3 .2. 17) a nd (3.2 . 19). Although the symbolism used in these two equations is very similar, the two equations are really quite different. Equation (3 .2.17) is displayed more clearly by Eq. (3 . 1.16) for which temperatures are given starting at time t I a nd continuing up to time tM + .-1' Moreover, there is only one heat input t o the body; namely, the unknown surface heat flux, and also the initial temperature is t he uniform value of To . Equation (3 .2.19) is much more general and can include the conditions for Eq . (3.2.16) o r Eq. (3.2.17 ), ~15' seKiilg 111 :;±wla:1-fl:,_o- +o . E quation (3.2.19) can also cover the case of an arbitrary temperature distribution at time t M - 1 a nd a k nown heat flux at x =L [ or o ther nonhomogeneous conditions as shown in Eqs. (3.2.23) a nd (3 .2.24)]. T he eI>(r, t) problem is the solution of the problem of aT/aqc where q (t)=qc' See Section 1.6. T he eI>(r, t) problem is ~~ r ar (r ael» = pc ael> or at _ k ael> or -k I =1 (3 .2.26) (3.2.27a) r =G ~~Lb =hel>(b, t) cf>(r, 0) = 0 (3.2.27b) (3 .2.28) Note that the only nonhomogeneous term is a t the heated surface, r =a . The el>i value is the solution of Eqs. (3.2.26)- (3.2.28) for the sensor location and time t i ' T he symbol t MI, .. =o represents the calculated temperature for the model at time tM for the estimated heat flux c omponents q l, q2, ·· ·' q M-I ' b ut qM is set equal to zero. F or example, for the model given by Eqs. (3.2.23)-(3.2.25), known values o f g(r, t) a re used in Eq. (3.2.23), known values of T,.,(t) are used in Eq. (3.2.24b), and the known initial temperature distribution, F(r), is used; the temperature at the sensor location is calculated a t time tM with these conditions and with the heat flux, q(t), equal to q I for 0 < t < I lt, q2 for I lt t o 21lt, . . . , qM - I for t M- 2 < t < t M- I , a nd q = 0 for t M- 1< t < tM' T he temperature t M+ d, .. : , . .. , =0 is similarly found a t the sensor location a nd a t t M + I with the known qi values for i =1 , 2, . .. , M-l a nd also q M=qM+I=O. F or the case g(r, t) = 0 in Eq. (3 .2.23), T,.,(t) = To in Eq. (3.2.24b) and T (r, 0) = To in Eq. (3.2.25), the t lq=o vector is simply given by M -I L qillel>M - i+ To 3.3 DIFFERENCE M ETHODS Duhamel's theorem is a powerful technique for solving a wide variety of linear heat conduction problems. Unfortunately, many situations exist in which either (I) this technique becomes cumbersome o r (2) it is n ot applicable. The most severe limitation is for nonlinear problems such as for bodies with temperaturedependent thermal properties. In many heat transfer problems, the temperature change of a body exposed to a heat flux is sufficiently large that the changes in thermal properties are appreciable. Fortunately, numerical methods can be used t o convert the nonlinear partial differential equation of heat conduction into a system of linear algebraic equations involving the temperature a t discrete locations. 3.3.1 F inite C ontrol V olume P rocedure f or C onstant P roperty P lanar G eometries Many methods exist for discretizing partial differential equations. T wo p opular ones are (1) finite differences a nd (2) finite elements. An alternative approach in which conservation o f energy is applied directly to control volumes is a dopted here. The finite control volume (FCV) procedure uses a control volume o f arbitrary and finite size as in Figure 3.2 a nd the integral form of the energy equation is applied to this control volume. Conservation of energy for a stationary solid requires that the sum o fthe rate of heat flow across the control volume boundary and the time rate of change o f energy within the control volume is equal to the rate a t which energy is produced within the control volume. In equation form, the conservation o f energy can be written as i =1 ff q.dA+ fr fff ped-Y fff e'"d-Y M -I Tlq =o= L qillel>M-i+1 + To i =1 M -I L i= 1 q illel>M-i+r-I+To = (3.2.29) A -r (3.3.1) -r where A a nd - Yare the closed surface area and volume respectively o f the control volume. q is the heat flux vector leaving the control volume surface, e is t he 88 C HAP.3 M ETHODS F OR DIRECT HEAT C ONDUCTION PROBLEMS Region of interest FIGURE 3 .2 Typical finite control volume. specific energy (energy per unit mass), and em is the energy production rate per unit volume. In the control volume shown in Figure 3.2, a p art o f the control volume boundary coincides with the boundary o f t he body, a nd t he remaining control volume boundary is completely within the body. Consequently, the surface integral o fthe h eat flux can involve both heat flux from external sources (boundary conditions) and heat conducted t o t he given control volume from all neighboring control volumes. I n its present form, Eq. (3.3.1) is too general to demonstrate the details o f the FCV procedure. Let us restrict o ur a ttention to one-dimensional p lanar p roblems and divide the body o f interest into increments o f e qual size (called elements) as shown in Figure 3.3. Each element has its own thermal properties a nd a n element boundary can also correspond t o a material interface. The control volume boundaries a re chosen t o lie a t t he midplane o f each element; this selection o f c ontrol volume boundary is somewhat arbitrary. T he choice for control volume boundaries permits two different materials t o be present within a single control volume. Each element contains two nodes, o ne a t each boundary o f the element which are identified by a (. ) i n Figure 3.3. First, the general conservation o f energy Eq. (3.3.1) is applied t o t he control volume surrounding node 1; Figure 3.40 provides an expanded view o f this control volume. T he result is x , p osition_ ~~jr r-- - I I q (t) - --• Element 2 : I ,=~~':''::'''';-=-::J'- Control , volume I boundary N ode: Node I I Element boundary 2 3 I ~~-~-~~-~-~-~-~-~~~~~ 3 Lx FIGURE 3 .3 boundaries. q~t) '=-=-==-=-=~-=--:::-~ N -l N 2 One-dimensional model for planar geometry showing element and control volume r - --, 1 q(t)+ I +--k :!' I 1 I 1 "" liT " , _ __ J 1 ( a) _ k a TI - q(t)+ !!.. ax Axl2 dt ( Ax/2 Jo p cTdx=O (3.3.2) where the energy content is given by e =cT a nd e"'=O. F or t he arbitrary interior node j a nd the backface node N (see Figure 3.4b a nd c, respectively), application o f the conservation o f energy equation yields r -----' • -+-. -+--. -k:;1 L _____ l-k :!'I . . "" I X j-T j -l a TI aTI d -k+k+ax XJ+ A x/2 ax x r A x/2 dt f ,XJ+Ax/2 pcTdx=O (3.3.3) 1 I .J j % J+T j +l (b) x r A x/2 a TI d qN(t)+k - a +-d- f ,XN pcTdx=O x X N-Ax/2 ' xN-Ax/2 (3.3 .4) T his control volume approach insures that the heat conducted o ut o f o ne control volume is c onducted into the adjacent control volume. In order t o evaluate the integrals and derivatives in Eqs. (3.3.2) - (3.3.4) a nd convert them (e) FIGURE 3 .4 Expanded view o f various con· trol volumes. (a), Control volume surrounding surface node I ; (b). control volume surrounding arbitrary interior node j ; (c). control volume surrounding backface surface node N . 89 90 C HAP.3 M ETHODS FOR DIRECT H EAT C ONDUCTION P ROBLEMS into useful computational algorithms, it is necessary to make some assumptions about the temperature profile within each element. Many assumptions can be made concerning the element temperature profile; however, the more physics exercised, the better. As an example, for a steady-state, planar geometry with specified temperatures as boundary conditions and constant k, the temperature profile is linear with position. Thus for the transient case, a linear temperature profile is assumed for the element. F or an element bounded by X j - l ~X~Xj, the assumed temperature profile is T ( x )=T. l r ( x-Xj) T. ( x-xj_d + J ( Xj- Xj_ d ' Xj) ( Xj-l - X j-l ~X~Xj (3.3.5) where 1 j= T (xj) is a node point (or nodal) temperature. Equation (3.3.5) is written in the form of a Lagrangian interpolation polynomial. The local heat flux is determined by differentiating Eq. (3.3.5) and using Fourier's law k aT ( 1j-1j- d q (x)=- -;-=-k , X j_l<X<Xj (IX X j-Xj-l (3.3.6) The linear element temperature profile produces a constant heat flux within elements; however, the heat flux changes from element to element. Assuming that the volumetric heat capacity (pc) is a constant for each element, a numerical expression for the energy storage term requires the evaluation of integrals of the form T(x)dx. I t is convenient to perform these integrals over portions of elements and then add the appropriate contribution to the control volume energy balance. The two integrals that occur are evaluated using Eq. (3.3.5) t o obtain J (3.3.7) (3.3.8) Note that Eqs. (3.3.7) and (3.3.8) are proportional to the average element temperature over h alf of the element. From the foregoing linear element temperature profiles, the finite control volume energy balance for nodes 1, j and N, are k ( T1 -T2 ) Ax d Ax - q(t)+pc 2 "dt ( iTl + iT2 )=0 k (1j-1j+d _ k(1j-l-1j)+ Ax !!-.(1T, Ax Ax pc 2 dt 4 r + pc Ax d 2dt q N(t)-k (i1j+i1j+d=0, ( TN-1-TN) Ax d I Ax +PC (4 TN-l+i TN)=0 2dt These three equations represent a system of N first-order ordinary differential equations in time. By assuming an element temperature profile and integrating out the spatial dependence, the F ey procedure yields a system of ordinary differential equations (as opposed to a single partial differential equation). Note that each of the energy storage terms in Eqs. (3.3.9)- (3.3.11) contains a contribution from a node outside the control volume; this occurs because the linear element temperature profile depends on two nodal temperatures, one of which lies outside the control volume of interest. The foregoing results are sometimes referred to as distributed capacitance effects because the thermal capacitance of a control volume may be distributed over as many as three nodes. The finite difference method is another. technique for arriving at a system of ordinary differential equations representing heat conduction within a solid. This procedure is also called lumped capacitance because all of the control volume thermal capacitance is lumped at the center node point for all interior control volumes. This is accomplished by replacing (i 1j _ I + i 1j) by 1j and ( i1j+i1j+d with 1j in Eq. (3.3.10); for surface control volumes, ( iT1 + !T2 ) becomes TI and ( iTN- 1 + iTN) becomes TN' Additional details on difference methods are available in Myers,9 Smith, I S Patankar,21 Anderson et al.,22 and S hihP Another method for the numerical solution of heat conduction problems is the finite element (FE) method. I1 - 14 Many methods can be used to develop finite element type systems of differential equations. The Galerkin-type method of weighted residuals is one of the simpler approaches and requires that the integral weighted average of the partial differential equation be zero over some domain. F or example, for a slab i L (aT aT) 2 w(x) - -ex dx=O o at O X2 (3.3.12) where L is the thickness of the slab and w(x) is a weighting function to be specified. There is an infinite number of possible weighting functions; the trick is to choose a weighting function that gives good results. A popular weighting function is the so-called tent function which is shown in Figure 3.5 : X j+l-X X j+l J j =2,3, . .. , N -1 91 DIFFERENCE M ETHODS (3.3.9) "JT,) 1 +4 SEC. 3.3 = 0, (3.3.10) (3 .3.11) - x/ all other x (3.3.13) Note the similarity between the weighting function Wj(x) and the Lagrangian interpolation polynomial in Eq. (3 .3.5), and that the weighting function is different for each node j . Since the weighting function Wj(x) is zero over a large 92 C HAP.3 M ETHOOS FOR DIRECT H EAT C ONDUCnON PROBLEMS w (x) SEC. 3.3 DIFFERENCE M ETHODS 93 All o fthe F CV, LC, a nd F E results can be p ut in a similar form : d (X q(t) d (fJT, + 1'Tz)= - A2 (T, - T2 )+ ---:.t uX p cux d - (1'TN-' dt (X + P TN)=Al ( TN-, uX qN(t) TN)- - - p cAx (3.3.19) (3 .3.21) where fJ a nd l' h ave t he following values : FIGURE 3.5 Weighting runction ror Galerkin-typc method o r weighted residuals. P LC F CV FE p ortion o f t he d omain O~x~L, Eq. (3.3.12) c an be written as "} i +1 " }-I (aT aT) 2 Wj(x) - at -(X -2 ax d x=O (3.3.14) Utilizing t he l inear element temperature profile assumption a nd E q. (3.3.13) for Wj(x), t he F E e quation for a n a rbitrary interior node j c an be written as k ( 1)-1)+1) - k ( 1)-, - 1) A x!!.-, z Ax Ax + pc 2 dt h 1)-, +"3"1) Ax d + pc T di ( i1)+l1)+d=O, j =2, 3, . .. , N -l T he o~ly difference between the F CV a nd F E results is the different weighting coefficients o n t he t hermal capacitance terms; the F CV gives a greater weighting t o t he center temperature 1). T he F E e quations for t he surface nodes are l' P+1' 1/8 1/6 1/2 1/2 1/2 o (3.3.22) E quations (3.3.19) - (3.3.21) a re restricted t o t emperature-independent thermal properties b ut t he c oncepts c an readily be extended t o T-variable cases. 3.3_2 (3.3.15) 1/2 3/8 2/6 O ther B oundary C onditions a nd M aterial I nterfaces In Section 3.3.1, a p lanar b ody with c onstant p roperties subjected t o a specified heat flux b oundary c ondition was considered. A heat flux given by a convective heat transfer coefficient a nd a t emperature difference driving potential is also a common b oundary c ondition. F or t his' l atter case t he energy balance o n t he control volume s urrounding n ode I c an b e obtained directly from Eq. (3.3.19). Ax Ax d -q(t)+pc Tdi<1 T, +1 Tl)=O (3.3.16) d (X h - (PT, + 1'Tl )= - - 2 ( T,-T2 )+--(T",,-T,) d t Ax p cAx (3.3.17) k ( T1 -T2 ) Note t hat convective b oundary c onditions involve the surface temperature. T he F CV a pproach c an be readily extended to elements o f n onuniform size and differing material properties; such a control volume is shown in Figure 3.6. L emmon a nd H eaton'6 d emonstrated t hat if the weighting function w(x) is chosen t o be t he D irac delta function Wj ~ \ - () ~ x = u(x-Xj)= {I, 0, X =X· } all o ther x Material interface I' ( e+ 1) (3.3.18) t hen t he F E results for interior nodes a re identical t o t he lumped thermal capacitance (LC). This method is also called the collocation method. Additional details o f the F E m ethod c an be found in Zienkiewicz, 1 2 N orrie a nd deVries ' 3 Becker et aI., ' 4 a nd B aker. 24 ' j-l j +1 j FIGURE 3.6 C ontrol volume for material interface. (3.3.23) 94 C HAP.3 M ETHODS FOR DIRECT HEAT C ONDUCTION PROBLEMS T he two adjacent elements are. denoted (e) a nd ( e+ 1) with each having its own temperature dependent properties (pc)le) a nd kle ). Applying an energy balance kle + l ) 1 )-1)+1 _ k le ) 1)-I-1)+(pc)e(XrXj_I)~(1'1)_1 + /11)) X j+l-x j X j-Xj_1 dt + (pc)le+ I)(Xj+ 1 - Xj) ~ (/11)+1'1)+ d=O Expanding U(l + I\l) in a Taylor series about time l, gives dU(l) u (t+l\t)=U(l)+ """dt - I\l+ . .. (3.3.24) (3.3.25) where the subscript R indicates right half o f element (e). Similarly, for the left half o f element ( e+ 1), the average temperature is f r+ I ) = 11)+*1)+ 1 96 DIFFERENCE METHODS (3.3.28) If the terms 0 (l\t 2 ) are ignored, then u(t + I\t) is given by T he element thermal conductivity is evaluated a t the mid-plane of each element; for a linear element temperature profile, the appropriate temperature for element (e) is f le)=(1)_1 +1))/2. There are several choices for the appropriate temperature a t which t o evaluate the volumetric heat capacity. The simplest procedure is to use the average element temperature to evaluate (pc)le) and assign this value to both halves o f the element. A slightly more complicated approach is to compute the average temperature of each half of an element and evaluate the volumetric heat capacity a t this temperature. F or example, the right half o f element (e) contributes t o control volume j a nd the average temperature of this half of element (e) is f~) = l[fle ) + 1)] = 1[t(1)-1 + 1)) + 1)] = i1)+*1)-1 SEC. 3 .3 (3.3.26) At this point, a word of caution is a ppropriate for the evaluation of temperaturedependent thermal properties. Most thermal properties are known only within ±5% at best. Consequently, it is not particularly fruitful to develop highly sophisticated techniques to precisely handle temperature-dependent thermal properties that are known imprecisely. Many times, the simplest approach is more than adequate. u(t + I\t) z u(t)+ f [u(t), l ]l\t (3.3.29) I f u(t.) is denoted as u ., then Eq. (3.3.28) becomes U .+ 1 = u.+ f (u., t .)l\t=u.+ ~~I.l\l (FO) (3.3.30) The algorithm given by Eq. (3.3.30) is shown in Figure 3.7; it is simply the extrapolation of the slope a t p oint (un' tn) t o (un+ I , tn+ d. As an example of the Euler method, it is applied to Eq. (3.3.19) with /1=1/2, 1 '=0 (LC). Since the subscript o n T represents node number, a superscript is used t o indicate time, 2 nl\l Ti)+~ p cux T j+ 1 = T j - 2p(Tj - (3.3.31) where p =a.l\t/l\x 2 is a F ourier number based o n the time-spac~ .grid. The Euler method for the LC formulation is also known as the explIcit method since T j+ 1 in Eq. (3.3.31) can be written as an explicit function of known temperatures only at time n ; a nother name is the forward difference method because the time derivative is replaced by a forward time difference. Explicit results are also obtained from Eqs. (3.3.20) and (3.3.21) for the LC formulation. As can be seen from Figure 3.7, the Euler method underpredicts the temperature if the true solution is concave upward. In an attempt to improve stability, a backward time difference (BO) is often used; analogous to Eq. (3.3.28), the I / 3 .3.3 N umerical T echniques f or S olving S ystems o f F irst-Order O rdinary D ifferential E quations u (t) '/ / / '- I n principle, an analytical solution of the system of energy balance Eqs. (3.3.19) (3.3.21) is possible. However, numerical techniques offer solutions t o more general problems and can be readily extended to the temperature-dependent thermal properties. Numerous numerical techniques are available. O ne o f the simpler techniques is the Euler (or forward difference, FO) method which can be developed from a Taylor series expansion. Let u(t) be a function of time that satisfies the first-order ordinary differential equation du d t = f (u, t) sIOpeU~+l I (3.3.27) FIGURE 3 .7 methods. Extrapolation of slope u~ = { (un,tn) Schematic o f the Euler (forward difference) and implicit (backward difference) 96 C HAP.3 M ETHODS FOR DIRECT H EAT C ONDUCTION PROBLEMS 2 ( 3.3.32) Neglecting terms of 0(llt 2), the backward difference integrator gives (3.3.33) The backward time difference integrator is also known as the fully implicit method because the right-hand side contains the unknown Uo + I . A central difference (CD) time derivative is known to be more accurate (under certain conditions) than either the forward o r backward difference and can be derived from the following Taylor series expansions: Ilt u" ( Ilt)2 uO+I/2=uo+u~2+2~ 2 + ... (3.3.34) , (3.3.35) U +I/1=UO +I-UO +1 O 97 DIFFERENCE M ETHODS d a. 0+ 1 0 - (RT, + T )0+1 = - 0 _ (To+l_ T~+I)+O -q dt I ' 1 Y 2 Ilx1 1 pcllx equation is , .. " Ilt U o=Uo+l-uo+lut+Uo+1 2 1-··· SEC. 3.3 Ilt U:+ 1 ( Ilt)2 2 +T! 2 _ . .. O' -d (fJ TI+yT1)0=-8' ..a. 2(T~-T~)+0' q: d uX p cux t where 8' = 1 - 0 a nd the superscript n means evaluated a t time to· Adding the foregoing two equations gives d d a. 0 - (fJT1+yT1)0+ 1 + 8' -d (fJT1+ yT1)0= - 0 A 2 (T~+ 1 - T~+ I ) dt t uX - 8' - -; (T~ - Ilx T~)+~ (Oq"+ 1 + 8'q") p cux The general integrator can be applied directly to the left-hand side ofEq. (3.3.39) with the result Tft+ 1 _ TO) ( T o+1 - TO) fJ ( 1 1+ 2 2 Ilt Y Ilt (3.3.40) Equating these two expressions and solving for Uo + 1 gives the central difference integrator (3.3.36) which is also known as the Crank-Nicolson method. von Rosenberg 1 8 describes the Crank-Nicolson procedure as successive applications of the forward and backward integrators. The forward, backward, and central difference integrators can all be combined into a general integrator (GI) as follows: U +1 = uo+ [OU~+I +(1-0)u~]llt o (3.3.37) 1 0 = 1, q "+e= OqO+ 1 +O'q" (3.3 .41) Equation (3.3.40) can be simplified to yield Ilt p cux (pO+fJ)T~+I_(pO-y)T~+ 1 = - (pO' -fJ)T~ +(pO' +y)T~ + ---"A q0+e (3.3.42) Applying the G I t o Eqs. (3.3.20) and (3.3.21) gives (3.3.43) +2(fJ - pO')Tj+(p8' +y)Tj+ I , j =2, 3, . .. , N -l forward difference o r Euler 0 =2' central difference o r Crank-Nicolson where for convenience the following definition is employed -(pO-y)Tj~: +2(pO+ f J)Tr 1 _(pO-y)TH: =(p8' + y)Tj_1 The three special cases are: 0=0, (3.3.39) (3.3.38) Ilt -(pO-y)TN~11 +(pO+fJ)TN+I =(pO' + y)TN- 1 +(fJ - pO')TN- --;- qN+6 p cux backward difference o r fully implicit 3 .3.4 G eneral F orm o f D ifference E quations f or H eat C onduction i n P lanar B ody The general integrator given by Eq. (3.3.37) for first-order equations converts an ordinary differential equation into a single difference equation. The next step is to convert the system o f ordinary differential equations, Eqs. (3.3.19) (3.3 .21), into a system o f algebraic difference equations. The procedure for applying the GI, Eq. (3.3 .37), is to evaluate the differential equation at times n + 1 and n, multiply these equations by 0 a nd ( 1- 0), respectively, and add the two equations. Using Eq. (3.3.19) as an example gives (3.3.44) Equations (3.3.42) - (3.3.44) form a linear system of N algebraic equations with N unknowns. These algebraic equations have a very special structure that is known as tridiagonal and is best visualized by writing them in matrix form for constant thermal properties, pO+fJ - (pO-y) - (pO-y) 2(pO+:P) - (pO-y) TN~\ - (pO-y) pO+fJ TN 1 + C HAP.3 98 fJ - pO' y+p(:J' y+p(:J' 2(fJ - p(:J') M ETHOOS FOR DIRECT H EAT CONDUCTION P ROBLEMS' y + p(:J' y+p(:J' [ y+p(:J' fJ - p(:J' ][ T~] T~ At S EC.3.3 Comments on Accuracy. I n order t o o btain the desired accuracy it is necessary to use a sufficient number o f nodes N and t o t ake a "small" time step. In many cases a value of N equal t o 20 is sufficient. The time steps must be a t least as small as the time steps in the measured temperatures but can be considerably smaller; it is convenient t o choose the computational time step, At, so that the time step between measurements, At m• as , is equal t o [qn+8] 0 6 i'N-I + pcAx T'N -q'N+ 8 (3.3.45) T he tridiagonal coefficient matrix has nonzero terms along only three diagonals. This tridiagonal structure allows for a very efficient specialization of the Gauss elimination procedure for solving linear systems o f algebraic equations. The details of the algorithm, known as the Thomas algorithm, are presented in Section_ .3.2. 6 T he form o f Eq. (3.3.45) shows the explicit dependence o f the new temperatures on the old temperatures and the boundary conditions. From a computational point o f view, it is more convenient t o lump the right-hand side of Eq. (3.3.45) into a constant vector. The Euler integrator (forward difference, (:J = 0) for the LC (fJ = 1/2, y = 0) is worthy o f special note. F or this case, the coefficient matrix for the unknown temperatures become diagonal; consequently, each unknown temperature T7+ I can be expressed explicitly in terms o f known temperatures ( Ti-I' Ti, Ti+d. F or either the FCY (fJ=3/8, y =I/8) o r F E (fJ=2/6, y =I/6) method in conjunction with the Euler integrator «(:J = 0), all three diagonals o f the coefficient matrix are nonzero; hence, the term explicit integrator is inappropriate for these cases. EXAMPLE 3.2. Give the difference equations for the Crank-Nicolson approximation for a plate that is heated a t x =O a nd insulated a t x =L. Let Ax = L/3 a nd use a lumped thermal capacitance. The thermal properties are independent of temperature. Solution. The difference equations can be obtained using Eq. (3.3.45) for N = 4 (see Figure 3.3); there are two boundary nodes and two interior nodes. From Eq. (3.3.22) for the lumped capacitance approximation, fJ equals 1/2 and y equals O. F rom Eq. (3.3.38) for the Crank-Nicolson approximation, 8 equals 1/2. F or the insulated boundary, q~+ 1 /2 = 0. Using these values in Eq. (3.4.45) and multiplying the first and last equations by 2 gives the four difference equations in matrix form as, [ ,+, -t ~. ['-' p/2 0 0 -p p +l - p/2 0 p 0 l -p p/2 l -p p/2 p 0 or·] oJrJ [~."'J 0 - p/2 p +l -p o T j+1 o Ti p/2 l -p T: T4 A t= Atm. .. i (3.3.46) where i is a positive integer. In order t o choose the space step, Ax, a s well as the time step, At, experience is needed for each particular case because their values are affected by many factors including the time period of interest, variation of thermal properties with temperature, and the time variation o f the surface heat flux. F or more discussion, see books o n numerical methods such as References 6, 8, 9 a nd 15. 3.3.5 S tandard F orm o f T emperature E quations f or I HCP The purpose of this section is t o demonstrate that the difference equations produced by the FD, FE, and FCY methods can be used t o o btain the same standard equation, Eq. (3.2.19), that was obtained from Duhamel's theorem. The discussion in this section is confined t o the linear I HCP b ut nearly the same approach can be used for the quasilinear analysis which is discussed in Section 6.2. T he analysis can be performed using either algebraic o r matrix notation but the latter analysis is preferred because it is more compact and general. I f n + 1 is replaced by the index M, Eq. (3.3.45) can be written as (3.3.47) where A is the square matrix o n the left of Eq. (3.3.45), B is the square matrix on the right of Eq. (3.3.45) and C = A t/ pcAx. The vectors TM a nd qM are . ~~[::l ~[f"-') (n~) T i+ 1 - p/2 99 DIFFERENCE M ETHODS p +l w herelhe subscripts relate to the space nodes and the superscripts t o time. The = gM vector represents any known energy sources due to an applied heat flux T~+I 2 At + pcAx at x =L, volumetric heating, o r convective heat transfer. I f the Crank-Nicolson method of solution is selected, (:J = 1/2 and 0 0 0 q M+8-1 = qM-I/2 0 which is the surface heat flux evaluated a t time tM - 1/ 2 = (M - 1/2)At. Observe 100 C HAP.3 M ETHODS FOR DIRECT H EAT C ONDUCTION PROBLEMS t hat q M-I/2 m eans exactly the same heat ftux component as denoted q M in the Duhamel procedure [see Eq. (3.2.2)]. Equation (3.3.47) is solved for TM by mUltiplying by A - I t o o btain TM = DTM- I + EqM + A - lgM (3.3.49a) where D =A - 'B, E =CA- ' SEC. 3.3 DIFFERENCE M ETHODS E quations (3.3.53a,b,c) are vector equations for the temperatures a t each node in the body. Temperature sensors are located a t only a few o fthese nodes, in general. F or simplicity, a single sensor is considered a nd t he space dependence is indicated by a subscript K. T hen Eqs. (3.3.53a,b,c) c an be written for the K th node as (3.3.54a) (3.3.49b,c) (Though it is n ot a g ood numerical p rocedure t o find the inverse o f A when ~lc~lating T M it is convenient here t o d o so symbolically.) Notice t hat TM ,. IS a h near funchon o f qM . Replacing M in Eq. (3.3.49a) with M + 1 gives (3.3 .50) T M+ 1 ; "DTM+ EqM+l + A - lgM+l a nd t hen using Eq. (3.3.49a) for TM yields T M+ I = D 2T M- 1 + DEqM + EqM+l + DA - lgM+ A - lgM+1 (3.3.51) (3 .3.54b) M+2 T" =tM+21,.,=,.,+1 =,".>=0+ X 31q M+ X 21q M+I + X I lqM +2 " ( 33 .54c) . where (J in q M + 9- I is c hosen t o be 1 t o simplify the notation. T he s ymbols X I I' X 21, a nd X 31 d enote sensitivity coefficients a t l ocation K a nd a re n ot functions o f t he time index M a s indicated by the corresponding matrix terms of E, D E a nd D 2E. T he sensitivity components X I I' X 21, X 31 a re given by (3.3.54d) This expression is l inear in b oth qM a nd qM + I. R epeating the process for M replaced by M + 2 in Eq. (3.3.49a) gives T M+l=D3TM-l + D1EqM + DEqM+l + EqM+l + D2A - lgM + DA - lgM+I + A - lgM+l 101 (3.3.54e) (3.3.52) which is linear in qM, qM + 1 a nd qM + 2. D ue t o t he linearity o f Eqs. (3.3.49a), (3.3.51), a nd (3.3.52), they c an be written as (3.3.54f) TM = TMlq.,=o+EqM T M+ I = TM+ Ilq"=q.,+'=o+DEqM +EqM+ 1 (3.3.53a) (3.3.53b) Equations (3.3.54) can be written in the standard matrix form with the K subscript o n T a nd Tlq=o o mitted as (3.3.55a) T M+l=tM+llq"=q"" _q.,+1=0+D1EqM + DEqM+I + EqM+l (3.3.53c) where where the temperatures evaluated a t qM, . .. , equal t o 0 a re defined by TMlq., =o=DTM-I+A - lgM T=[~:+I] , T (3.3.53d) M T M+ Ilq"=q.,+1 =0=D1T M- 1 + DA - lgM + A - lgM+ I (3.3.53e) TM+llq"=q.,+I=q.,+> =0=D3TM-l + D1A - lgM + DA - lgM+1 + A - lgM+l (3.3.55b,c) +2 • (3.3.55d) (3.3.53f) ~otice t hat. the temperature components defined by Eqs. (3.3.53d,e,f) are hnear functIOns o f only the temperature distribution a t t ime t M - 1 a nd t he known heat sources a t times t M , tM + I , a nd tM + 1 • T hey are equal t o t he corresponding temperatures given by Eqs. (3.3.35a,b,c) if the future heat ftux vectors are equal t o zero. Equations (3.3.35a,b,c) display the important characteristic o f linearity in t he t emperatures defined by Eqs. (3 .3.35d,e,f) a nd t he h eat fluxes qM, qM+I, a nd qM+2 Whereas the temperature vectors in Eqs. (3.3.47)- (3.3.53) are for the nodes a t time t , t he temperature vectors in Eqs. (3.3.55) are for different times. E quation M (3 .3.55a) is called the s tandard form for the temperature because the same equation was derived using Duhamel's theorem [see Eq. (3 .2.19)]. In Eqs. (3 .3.55b,c,d) the superscript M d enotes the same time dependence as the M subscripts in Eqs. (3 .2.20a,b) a nd (3 .2.22a). T he X 1 1, X 2 " a nd X 31 values 1 02 C HAP.3 M ETHODS FOR DIRECT HEAT C ONDUCTION PROBLEMS correspond to 17. used in the Duhamel theorem method and the numerical values of X I I and X 2 1 and llcp I , and so on, will be nearly identical if the finite difference (or F E or FCV) calculations are accurate. The symbol cpj represents the temperature rise a t time t l and node X K for the numerical solution of Eq. (3.3.47) for q M+9-1 = 1 for i =M = 1 ,2,3, . .. a nd gM = 0. Since the same standard form, Eq. (3 .2.19) o r (3 .3.55a), is obtained for T, I HCP algorithms can be developed that do not depend on the type o f numerical approximation of the transient heat conduction equation. Advantage is taken of this fact in Chapter 4 to derive general algorithms. In Chapter 5 the algorithms of Chapter 4 are implemented using Duhamel's integral, and Chapter 6 uses the finite control volume method. llcpo, Keltner, N. R. a nd Beck, J. V., Unsteady Surface Element Method, J . Heat Transfer 103, 759 764 (1981). 2. Beck, J. V. a nd Keltner, N. R., Transient Thermal Contact o f T wo Semi-Infinite Bodies O ver a C ircular Area, in Spacecraft Radiative Transfer and Temperature Control, T. E. H orton, ed., Vol. 83 o f Progress ill Astronautics and Aeronautics (1982), AIAA, 61 - 82. 3. Beck, J . V., Schisler, l . P. a nd Keltner, N. R., Simplified Laplace Transform Inversion for Unsteady Surface Element Method for Transient Conduction, A I A A Journal, 22 , 1328- 1333, (1984). 4. Myers, G. E ., T he C ritica l Time Step for Finite Element Solutions t o T wo-Dimensional Heat Conduction Transients, A SME J . Heat Transfer 1 00,120- 127, (1978). 18. von Rosenberg, D. U., Methods for the Numerical Solulion o f Partial DijJerential Equation .•, Americal Elsevier, New York, 1969. 19. Richtmyer, R. D. a nd M orton, K. W ., Difference Methods for Inirial Value Problems, 2nd ed., Interscience Publishers, New York, 1967. 20 . Beck, J. V., G reen's Function Solution for Transient Heat Conduction Problems, Int . J . Heat Mass Transfer, 2 7,1235 - 1244 (1984). 21. P atankar, S. V., Numerical Heat Transfer and Flllid Flow, Hemisphere, W ashington, 1980. 22 . Anderson, D. A., T annehill, J. C., a nd Pletcher, R. H ., Computational Fluid Mechanics and Heat Transfer, Hemisphere, Washingto n, 1984. 23 . Shih, T. M., Numerical Heat Transfer, Hemisphere, Washington, 1984. 24. Baker, A. J., Finite Element Compulational Fluid Mechanics, Hemisphere, 1983. PROBLEMS 3.1. O btain mathematical expressions for the surface temperature o f a semi-infinite body that is initially at the uniform temperature of To and then is exposed to the heat flux history given by, REFERENCES 1. 1 03 PROBLEMS Keltner, N. R. a nd Beck, J . V., Surface Temperature Measurement Errors, J . Heat Transfer, lOS, 312 - 318 (1983). qo, O <t<t l q(t)= - qo, tl < t<2t l 0, t <O a nd l t >2tl A different expression is needed for each time interval. Plot the group o f T - To ( k PC)1/2 qo tl versus t it 1 for t it I from 0 to 4. Why d o the temperatures change after t = 2t 1 even though the net energy added is zero? S. Litkouhi, B. a nd Beck, J. V., Intrinsic Thermocouple Analysis Using M ultinode Unsteady Surface Element Method, AIAA Paper N o . 83 - 1937 . (To be published in AIAA Journal.) 6. Ozisik, M. N. Heat Conduction, Wiley, New York, 1980. 7. Luikov, A. V., Analyt ical Hear Diffusion Theory, Academic Press, New York, 1968. 8. Arpaci, V. S., Conduction Hear Transfer, Addison-Wesley, Reading, MA, 1966. 9. Myers, G . E., Analytical Methods in Conduction Heat Transfer, M cGraw-HiII,"New York, 1971. 10. 11. ~ \ .) C arslaw, H. S. a nd J aeger, J. C , Conduclion o f Heat in Solids, 2nd ed ., O xford, L ondon, 19~ . Bathe, Klaus-Jurgen, Finite Element Procedures in Engineering Analysis, Prentice-Hall, Englewood Cliffs, NJ, 1982. 12. Zienkiewicz, O. C , T he Finite Element M ethod, 3rd ed., McGraw-Hili, New York, 1977 . 13 . Norrie, D. H. a nd deVries, G., An Introduclion 10 Finite Element Analysis, Academic P re ss, New York, 1978. 14. Becker, E. B., C arey, G . F., a nd O den, J . T., Finile Elements an Introduction, Vol. I , PrenticeHall, Englewood Cliffs, NJ, 1981. IS. Smith, G . D., Numerical SolUlion o f Parlial Differential Equations, O xford University Press, New York, 1965. 16. Lemmon, E. C. a nd H eaton, H. S ., Accuracy, Stability, and Oscillation C haracteristics o f F inite Element Method for Solving Heat Conduction Equation, ASME P aper N o. 69-WA/HT35 . 3 .2. T he heat flux history a t the surface of a body is 10 W/m2, q(t) = { 0 2 s<t<4s otherwise The temperature a t X l is 400 K until t =2s, 401 K a t t =4s, and 406 K a t t = 6s. The heat conduction problem is linear. Find the temperature rise a t X 1 a t times 2 and 4s for a unit step rise a t t = 0 in the surface heat flux . 3 .3. Find exact expressions for the surface temperature of a semi-infinite body that is initially a t the uniform temperature of To a nd that is subjected to a heat flux that varies in a triangular fashion with time as given by 40 t , l q(t)= ~0(2tl-t), where 40 is a constant. O <t<t l tl < t<t 2 =2t l otherwise 1 04 C HAP.3 M ETHODS FOR DIRECT HEAT C ONDUCTION PROBLEMS 105 PROBLEMS 20 w here G(x, l )=(kla)ocp(x, t)/Ol. (G(x, t) is a Green's function .) a. Is it possible t o use a numerical approximation o f Eq. (a) when G(O, t) is needed a t t = O? G ive a reason for your answer. b. D erive the numerical approximation o f Eq. (a), G ive expressions for the three time regions o f (a) 0 < 1< 1I , (b) 1I < 1< 211 a nd (c) 2 tl < t. P lot t he surface temperature in a dimensionless form versus t it I for t he r ange o f 0 t o 3. M 3.4. A h eat flux is given by T M=To+ L qft H M-ft+IL11 (b) n =1 j qG' q (l) = ,,-qG qG + ql> _ tG ( 1-1G ), t a H i == G (x, l i-1/2) Ie F or a semi-infinite body initially a t To c alculate the surface temperature a t time tl> using Eq. (3.2.12). . Use L11=0.5, IG= 1, qG= 1, q ,,=2, a nd k pc= 1. C ompare t he answer with the exact value. c. U se Eq. (b) for t he l inear h eat flux case o f E xample 3.1 a nd c ompare t he values. 3.9. S how t hat if CPi is given by 3 .5. Derive t he a pproximation o f t he Duhamel theorem given by M T M=To+ L qftaM - ft+1 .. .. 1 where q i=q a t I i a nd ; t hat t he s um ((I) 4>; , := 0 M -",-l for c"'S 0 j=(cpl')-2cpl~1 +CPj~2)~ ; =2, 3, .~ .. ,.... a nd ' f'j IS e t emperature n se a t time t i for the heat flux q = I . T he h eat flux is composed o f c onnected linear segments starting with q o=O . ( The superscript 1 denotes a linear h eat flux " building block.") T M= 1= cp\l)/M; . .1.(1)' T M + I = e - I TM + B(e-< - l)qM-",e-("+ 1)< 3 .10. 3 .6. W rite a c omputer o r p rogrammable calculator p rogram for the numerical convolution as given by Eq. (3.2.12). Allow for a t least 10 qft a nd 10 L1CPi c omponents. T he i nput is t o b e To a nd t he qft'S a nd L1CP/s. T he o utput is Assume t hat t he temperature rise for a n u nknown surfa~e h eat flux is L1 T (r, t) which is a k nown function. By t aking t he L aplace transform o f Eq. (3.2.9), d erive q (t)==~-I ~[L1T(r, I )] s 9'[q,(r, In w here s is t he L aplace transform parameter. F or TI , T2 ,···, 7 ;0' 1 Verify y our p rogram by checking Example 3.1. c p=2 ( _ k- n pc 3 .7. A semi-infinite body initially a t t he uniform temperature o f z ero is exposed a t t = 0 t o t he heat flux. )1/2 ,L1T 440(al)3/2 3a ( )1/2 k n find q(I). q = (nt)-1/2 3 .11. Let a = I , k = I, a nd L1t = 1. C alculate the temperature a t x =O using Eq. (3.2.12). Use a t least 10 time steps. T he e xact temperature is 1. C omment o n t he accuracy o f t he procedure. P artial Answers: 0.900, 0.784, 0.841, 0.865, a nd 0.880 q (l)G(x, 1- l)dl U sing the Laplace transform show t hat t he t emperature rise in a body exposed t o t he linearly varying heat flux, q (l) == 401 is e qual t o L1T(r, 1)==40 3 .8. An alternate way to write the convolution integral is i f~ q~CPM - ft c an b e written in t he s equential form th T (x, I) = To + L ft=1 (a) f~ cp(r,I')dl' where q,(r, I) is for a unit step heat flux starting a t 1 =0. 106 C HAP.3 METHODS FOR DIRECT HEAT C ONDUCTION PROBLEMS 3.12. Using the divergence theorem of Gauss, demonstrate that Eq. (3.3.1) becomes iJe iJt 3 .13. F or one-dimensional (radial) cylindrical geometries, use the steadystate temperature profile as the element temperature profile and develop the equations analogous t o Eqs. (3.3.9) - (3.3.11). Assume the control volume boundaries are a t the midpoint of the elements. 3 .14. Repeat Problem 3.13 for one-dimensional (radial) spherical geometries. 3 .15. Verify the finite element equation, Eq. (3.3.15). A x=L/4 3 .17. a. F or two nodes in a plate, one at x = 0 a nd t he other at x = L , give the two Crank-Nicolson, lumped capacitance equations. The surface at x = 0 is exposed to a heat flux a nd the surface at L is insulated. b. F or the temperature a t n ode 2, derive a n expression analogous to Eq. (3.3.54a). Give T 1'l f "=0 a s a function of T~ - 1, T 1' - 1, p, a nd A t/p cAx . c. Give X 1 1 as a function o f p a nd A t/p cAx. 3 .18. T he one-dimensional Green's function equation analogous to Duhamel's theorem is given in Problem 3.8. Write the equation in the form tM)=To+~ E" ex M q(A.)G(x, tM-A.)dA. J.'; =To+kj~1 "_1 q(A.)G(x, tu-A.)dA. Approximate q(A.) for t j _ 1 < A. < t l by t j-A. A .-t j _ 1 q(A.) = qj-I M + qj - -;;t" where qj = q (tJ, t j = i At . Derive the approximate relation, M T (x, t M)=To + I i)==~ 12 (x, M, i ) == ex k L-, J.". ' [ iql-1It(X, M, i ) i =l (a) _I G(X,IM-A.)dA. A A.t G (x, tM - A.)d..1. Show t hat E q. (a) c an be written as M T (x , t M)=To+ I j =O W hat are the a ju(x) expressions? which is exposed t o a time-variable heat flux at x =O a nd a convective boundary condition at x = L . T he properties are temperature independent. Use t he backward difference time approximation a nd lumped capacitance. T (x, where I I(x, M , V ·q+p - =e'" 3 .16. Give the difference equations in matrix form for a plate with 1 07 PROBLEMS a jM(x)qM-j