VIEWS: 0 PAGES: 17 CATEGORY: Technology POSTED ON: 1/13/2010 Public Domain
Session 8 Heat Conduction by Laplace Transform, Duhamel, Green Function and Goodman Approximate Integral Methods 1 Introduction There is a variety of heat conduction problems which although linear, are not readily solved using the separation of variables method. Fortunately, a number of methods not based on separation of variables are available for the solution of these heat conduction problems. Four such methods are brieﬂy described below, namely: the Laplace transform method, the Duhamel method, the Green’s function method and Goodman’s approximate integral method. In each case, examples are used to illustrate the application of the method. 2 Laplace Transform Method The Laplace transform method makes use of the Laplace transform to convert the original heat equation into an ordinary diﬀerential equation which is more readily solved. The so- lution of the transformed problem must then be inverted in order to obtain the solution to the original problem. 2.1 Deﬁnitions Consider a real function F (t). The Laplace transform of F (t), L[F (t)] = F (s) is deﬁned as ¯ ∞ L[F (t)] = F (s) = ¯ e−st F (t )dt t =0 The inverse transform (to recover F (t) from L[F (t)]) is 1 γ+i∞ F (t) = est F (s)ds ¯ 2πi s=γ−i∞ 1 ¯ where γ is large enough that all the singularities of F (s) lie to the left of the imaginary axis. For the transform and inverse to exist, F (t) must be at least piecewise continuous, of exponential order and tn F (t) must be bounded as t → 0+ . Laplace transforms have the following properties: Linearity. ¯ ¯ L[c1 F (t) + c2 G(t)] = c1 F (s) + c2 G(s) Transforms of Derivatives. L[F (t)] = sF (s) − F (0) ¯ L[F (t)] = s2 F (s) − sF (0) − F (0) ¯ L[F (t)] = s3 F (s) − s2 F (0) − sF (0) − F (0) ¯ L[F (n) (t)] = sn F (s) − sn−1 F (0) − sn−2 F (1) (0) − sn−3 F (2) (0) − ... − F ¯ (n−1) (0) Transforms of Integrals. Let g(t) = 0t F (τ )dτ (i.e. g (t) = F (t)). t 1¯ L[ F (τ )] = F (s) 0 s t τ2 1 ¯ L[ F (τ1 )dτ1 dτ2 ] = F (s) 0 0 s2 t τn 1 ¯ L[ ... F (τ1 )dτ1 ...dτn ] = F (s) 0 0 sn Scale change. Let a > 0 be a real number. 1¯ s L[F (at)] = F ( ) a a 1 L[F ( t)] = aF (as) ¯ a Shift. 2 L[e±at F (t)] = F (s ¯ a) Transform of translated functions. The translated function is F (t − a) t > a U (t − a)F (t − a) = 0 t<a where 1 t>a U (t − a) = 0 t<a The Laplace transform of the translated function is L[U (t − a)F (t − a)] = e−as F (s) ¯ Also 1 L[U (t − a)] = e−as s Transform of the δ function. L[δ(t)] = 1 Transform of convolution. The convolution of two functions f (t) and g(t) is t f ∗g = f (t − τ )g(τ )dτ = g ∗ f 0 The Laplace transform is ¯ g L[f ∗ g] = f (s)¯(s) Derivatives. ¯ If F (s) is the Laplace transform of F (t) then ¯ dF = F (s) = L[(−t)F (t)] ¯ ds dn F ¯ = F (s) = L[(−t)n F (t)] ¯ dsn 3 Integrals. ∞ F (t) ¯ F (s )ds = L[ ] s t Functions of More than One Variable. Let f (x, t) be a function of the two independent variables x and t. The Laplace transform ¯ of f (x, t), L[f (x, t)] = f (x, s) is deﬁned as ∞ ¯ f (x, s) = e−st f (x, t)dt 0 With this, the following useful formulae are obtained ∂f (x, t) ¯ df (x, s) L[ ]= ∂x dx ∂ 2 f (x, t) d2 f (x, s) ¯ L[ ]= ∂x2 dx2 ∂f (x, t) ¯ L[ ] = sf (x, s) − f (x, 0) ∂t ∂ 2 f (x, t) ∂f (x, 0) L[ 2 ] = s2 f (x, s) − sf (x, 0) − ¯ ∂t ∂t 2.2 Inversion of Laplace Transforms The inverse Laplace transform of L[F (t)] is 1 γ+i∞ F (t) = est F (s)ds ¯ 2πi s=γ−i∞ The integration must be carried out in the complex domain and requires a basic understand- ing of analytic functions of a complex variable and residue calculus. Fortunately, inversion formulae for many transforms have been computed and are readily available in tabular form. A short table of useful transforms is included below. 4 F (t) ¯ F (s) 1 1 s 1 t √ 2 s√ t √π 2 s3 1 √ π t s n n! t sn+1 tn−1 1 (n−1)! sn n at n! t e (s−a)n+1 k sin kt s2 +k2 s cos kt s2 +k2 k sinh kt s2 −k2 s cosh kt s2 −k2 −as ¯ F (t − a); t ≥ a 0; otherwise e F (s) 2 √ k −k √ 2 πt3 e 4t e−k s √ k 1 −k s erf c( 2√t ) s e 2 √ k k √ ae−k s −eak ea t erf c(a t + √ )+ s t erf c( 2√t ) √ s(a+ s) 2 2 √ ∞ e−λn t/a J0 (λn x/a) J0 (ix √ s) 1−2 n=1 λn J1 (λn ) ; λn roots of J0 (λn ) = 0 sJ0 (ia s) 2.3 Applications Consider the problem of ﬁnding T (x, t) in 0 ≤ x < ∞ satisfying ∂2T 1 ∂T 2 = ∂x α ∂t and subject to T (0, t) = f (t) T (x → ∞, t) = 0 and T (x, 0) = 0 Taking Laplace transforms one gets d2 T (x, s) ¯ s¯ 2 − T (x, s) = 0 dx α 5 ¯ ¯ T (0, s) = f (s) ¯ T (x → ∞, s) = 0 The solution of this problem is √ T (x, s) = f (s)e−x s/α = f (s)¯(x, s) = L[f (t) ∗ g(x, t)] ¯ ¯ ¯ g Inversion then produces t T (x, t) = f (t) ∗ g(x, t) = f (τ )g(x, t − τ )dτ 0 ¯ Inversion of g (x, s) to get g(x, t) ﬁnally gives x t f (τ ) x2 T (x, t) = √ exp[− ]dτ 4πα τ =0 (t − τ )3/2 4α(t − τ ) If f (t) = T0 = constant, the solution is x T (x, t) = T0 erf c( √ ) 2 αt Consider now the transient problem of quenching a solid sphere (radius b) initially at T0 by maintaining its surface at zero. The Laplace transform method can be used to ﬁnd and expression for T (r, t). The transformed problem is 1 d2 (rT ) ¯ s¯ T0 − T =− r dr2 α α the solution of which is √ T0 T0 b sinh(r sα) ¯ T (r, s) = − √ s rs sinh(b sα) To invert this the trigonometric functions are expressed as asymptotic series and then in- verted term by term. The ﬁnal result is ∞ bT0 b(2n + 1) − r b(2n + 1) + r T (r, t) = T0 − {erf c( √ ) − erf c( √ )} r n=0 2 αt 2 αt 3 Duhamel’s Method Duhamel’s method is based on on Duhamel’s theorem which allows expressing the solution of a problem subjected to time dependent boundary conditions in terms of an integral of the solution corresponding to time independent boundary conditions. 6 3.1 Deﬁnitions Many practical transient heat conduction problems involve time dependent boundary condi- tions. Duhamel’s theorem allows the analysis of transient heat conduction problems involving time varying thermal conditions at boundaries. Consider the following problem 2 1 1 ∂T (r, t) T (r, t) + g(r, t) = k α ∂t subject to ∂T ki + hi T = fi (r, t) ∂ni on boundaries Si for i = 1, 2, ...N and for t > 0, as well as T (r, 0) = F (r) Solution by separation of variables is not possible because of the time dependent nonhomo- geneous terms g(r, t) and fi (r, t). Consider instead an auxiliary problem obtained by introducing a new parameter τ , un- related to t. Find Φ(r, t, τ ) such that 2 1 1 ∂Φ(r, t, τ ) Φ(r, t, τ ) + g(r, τ ) = k α ∂t subject to ∂Φ ki + hi Φ = fi (r, τ ) ∂ni on boundaries Si for i = 1, 2, ...N and for t > 0, as well as Φ(r, 0, τ ) = F (r) Since here g(r, τ ) and fi (r, τ ) do not depend on time, the problem can be solved by separation of variables. Duhamel’s theorem states that once Φ(r, t, τ ) has been determined, T (r, t) can be found by t ∂ T (r, t) = F (r) + Φ(r, t − τ, τ )dτ τ =0 ∂t If the nonhomogeneity occurs only in the boundary condition at one boundary and the initial temperature is F (r) = 0, the appropriate form of Duhamel’s theorem is t ∂ T (r, t) = f (τ ) Φ(r, t − τ )dτ τ =0 ∂t 7 3.2 Discontinuities Sometimes the nonhomogeneous term at a boundary f (t) changes discontinuously with time. The integral in Duhamel’s theorem must then be done separatedly integrating by parts. For a total of N discontinuities at times over the time range 0 < t < τN , the temperature T (r, t) for τN −1 < t < τN is given by t N −1 df (τ ) T (r, t) = Φ(r, t − τ ) dτ + Φ(r, t − τj )∆fj τ =0 dτ j=0 where ∆fj = f + (τj ) − f − (τj ) is the step change in the boundary condition at t = τj . If the values of f are constant between the steps changes, Duhamel’s theorem becomes N −1 T (r, t) = Φ(r, t − j∆t)∆fj j=0 for the time interval (N − 1)∆t < t < N ∆t. 3.3 Applications Consider the following 1D transient problem in a slab of thickness L. Find T (x, t) satisfying ∂ 2T 1 ∂T 2 = ∂x α ∂t subject to T (0, t) = 0 bt; 0 < t < τ1 T (L, t) = f (t) = 0; t > τ1 and T (x, 0) = 0 The appropriate auxiliary problem here is, ﬁnd Φ(x, t) ∂ 2Φ 1 ∂Φ 2 = ∂x α ∂t subject to Φ(0, t) = 0 8 Φ(L, t) = 1 and Φ(x, 0) = 0 The desired function Φ(x, t − τ ) is obtained from the solution to the above problem by replacing t by t − τ , i.e. x 2 ∞ −αβm t 1 2 Φ(x, t) = + e sin βm x L L m=1 βm where βm = mπ/L. So, for t < τ1 t df (t) T (x, t) = Φ(x, t − τ ) dτ = τ =0 dτ x 2 ∞ (−1)m 2 =b t+b 3 (1 − e−αβm t ) sin βm x L L m=1 αβm And for t > τ1 τ1 df t df T (x, t) = Φ(x, t − τ ) dτ + Φ(x, t − τ ) dτ + Φ(x, t − τ1 )∆f1 τ =0 dτ τ1 dτ since df /dτ = b when t < τ1 , df /dτ = 0 for t > τ1 and ∆f1 = −bτ1 , so τ1 x 2 ∞ −αβm (t−τ ) (−1)m 2 T (x, t) = {+ e sin βm x}bdτ τ =0 L L m=1 βm x 2 ∞ −αβm (t−τ1 ) (−1)m 2 −bτ1 [ + e sin βm x] L L m=1 βm Consider now the transient 1D problem in a solid cylinder (radius b) with internal heat generation g(t). Initially and at the boundary the temperature is zero. The auxiliary problem is ∂ 2 Φ 1 ∂Φ 1 1 ∂Φ + + = ∂r2 r ∂r k α ∂t subject to Φ(b, t) = 0 Φ(r, 0) = 0 9 The solution to this problem is b2 − r 2 2 ∞ −αβm t J0 (βm r) 2 Φ(r, t) = − e 3 4k bk m=1 βm J1 (βm b) where βm are the roots of J0 (βm b) = 0. From Duhamel’s theorem, the desired solution is t ∂Φ T (r, t) = g(τ ) dτ = τ =0 ∂t 2α ∞ −αβm t J0 (βm r) 2 t 2 = e g(τ )eαβm τ dτ bk m=1 βm J1 (βm b) τ =0 4 Green’s Function Method The method of Green’s functions is based on the notion of fundamental solutions of the heat equation described before. Speciﬁcally, Green’s functions are solutions of heat conduction problems involving instantaneous point, line or planar sources of heat. Solutions to other heat conduction problems can then be expressed as special integrals of the Green’s functions. 4.1 Deﬁnitions Consider the problem of a 3D body (volume V ) initially at zero temperature. Suddenly, a thermal explosion at a point r instantaneously releases a unit of energy at time τ while the entire outer surface S of the body undergoes convective heat exchange with the surrounding environment at zero temperature. The thermal explosion is an instantaneous point source of thermal energy which then diﬀuses throughout the body. The mathematical formulation of this problem consists of ﬁnding the function G(r, t|r , τ ) satisfying the nonhomogenoeus transient heat equation in V , 2 1 1 ∂G G + δ(r − r )δ(t − τ ) = k α ∂t subject to ∂G k + hG = 0 ∂n on the boundary S. Here, δ(x) is Dirac’s delta function (which is zero everywhere except at the origin, where it is inﬁnite). The function G(r, t|r , τ ) represents the temperature at location r and time t resulting from an instantaneous point source of heat releasing a unit of thermal energy at location r and time τ and it is called a Green’s function. Consider now the related problem of a 3D body of volume V , initially at T (r, 0) = F (r) inside which thermal energy is generated at a volumetric rate g(r, t) while its outer surface(s) 10 Si exchange heat by convection with a medium at T∞ . The solution to this problem can be expressed in terms of the Green’s function above as T (r, t) = G(r, t|r , τ )|τ =0 F (r )dV + V α t dτ G(r, t|r , τ )g(r , τ )dV + k τ =0 V t N 1 α dτ G(r, t|r , τ )|r =ri hi T∞ dSi τ =0 i=1 Si ki The terms on the RHS above are: i) the contribution due to the initial condition, (ii) the contribution due to the internal energy generation, and (iii) the contribution due to the nonhomogeneity in the boundary conditions. Line sources and plane sources can be similarly deﬁned respectively for two- and one- dimensional systems and nonhomogeneous transient problems in 2D and 1D are the solvable in terms of appropriate Green’s functions. For instance, for a 1D system of extent L, T (x, t) = G(x, t|x , τ )|τ =0 F (x )dx + L α t dτ G(x, t|x , τ )g(x , τ )dx + k τ =0 L t 2 1 α dτ G(x, t|x , τ )|x =xi hi T∞ τ =0 i=1 ki Here, the Green’s function G(x, t|x , τ ) represents the temperature at x, t resulting from the instantaneous release of a unit of energy at time τ , by a planar source located at x . 4.2 Point, Line and Surface Sources Point, line and surface sources can be instantaneous or continuous. Subscript and superscript i signs are used to make the distinction explicit. For instance gp denotes an instantaneous point i source of heat. In Cartesian coordinates, the relationship between gp and g(x, y, z, t) is i gp δ(x − x )δ(y − y )δ(z − z )δ(t − τ ) = g(x, y, z, t) 4.3 Determination of Green’s Functions Once the Green function associated with a particular problem is known, the actual temper- ature can be calculated simply by integration. So, one needs techniques for Green’s function determination. Consider the following homogenoeus problem in a region R. Find T (r, t) satisfying 2 1 ∂T T (r, t) = α ∂t 11 subject to ∂T + Hi T = 0 ∂ni on all surfaces Si and T (r, t) = F (r) Many special cases of this problem have already been solved by separation of variables. In all instances one can symbolically write the solution as T (r, t) = K(r, r , t)F (r )dv R where K(r, r , t) is called the kernel of the integration. Now, applying the Green’s function approach to this problem yields K(r, r , t) = G(r, t|r , 0) so the kernel is the Green’s function in this case. The associated nonhomogeneous prob- lem can also be solved once G(r, t|r , 0) is known by simply replacing t by t − τ to obtain G(r, t|r , τ ). Consider as an example the 1D transient heat conduction in a cylinder. 1 ∂ ∂T 1 1 ∂T (r ) + g(r, t) = r ∂r ∂r k α ∂t subject to T (b, t) = f (t) and T (r, 0) = F (r) The associated homogeneous problem with homogeneous boundary condition is readily solved by separation of variables and the solution is b ∞ 2 2 1 ψ(r, t) = [ e−αβm t 2 2 2 r J0 (βm r)J0 (βm r )]F (r )dr = r =0 b2 m=1 J1 (βm b) b r G(r, t|r , 0)F (r )dr r =0 where βm are the roots of J0 (βm b) = 0. Therefore, the Green’s function with τ = 0 is ∞ 2 2 1 G(r, t|r , τ )|τ =0 = 2 e−αβm t 2 2 2 r J0 (βm r)J0 (βm r ) b m=1 J1 (βm b) 12 and replacing t by t − τ , the desired Green’s function is ∞ 2 2 1 G(r, t|r , τ ) = 2 e−αβm (t−τ ) 2 2 2 r J0 (βm r)J0 (βm r ) b m=1 J1 (βm b) For homogeneous 1D transient conduction in an inﬁnite medium, the Green’s function is 1 (x − x )2 G(x, t|x , τ )|τ =0 = exp(− ) (4παt)1/2 4αt For homogeneous 1D transient conduction in a semiinﬁnite medium, the Green’s function is 1 (x − x )2 (x + x )2 G(x, t|x , τ )|τ =0 = [exp(− ) − exp(− )] (4παt)1/2 4αt 4αt For homogeneous 1D transient conduction in a slab of thickness L, ∞ 2 −αβm t 2 βm + H 2 G(x, t|x , τ )|τ =0 = 2 e cos(βm x) cos(βm x ) m=1 2 L(βm + H 2 ) + H 4.4 Applications The problem of ﬁnding ψ(x, t) for −∞ < x < ∞ satisfying ∂ 2ψ 1 ∂ψ 2 = ∂x α ∂t and subject to ψ(x, 0) = F (x) has the solution ∞ −1/2 (x − x )2 ψ(x, t) = [(4παt) exp(− ]F (x )dx = x =−∞ 4αt ∞ = G(x, t|x , τ )|τ =0 F (x )dx x =−∞ Therefore, the Green’s function of the problem of ﬁnding T (x, t) for −∞ < x < ∞ satisfying ∂ 2T 1 1 ∂T 2 + g(x, t) = ∂x k α ∂t and subject to T (x, 0) = F (x) 13 is simply −1/2 (x − x )2 G(x, t|x , τ ) = (4πα[t − τ ]) exp(− 4α[t − τ ] and the desired solution is (x − x )2 ∞ T (x, t) = (4παt)−1/2 F (x )dx exp(− x =−∞ 4αt α t ∞ (x − x )2 + dτ (4πα[t − τ ])−1/2 exp(− g(x , τ )dx k τ =0 x =−∞ 4α[t − τ ] The problem of ﬁnding ψ(r, t) for 0 ≤ r ≤ b satisfying ∂ 2 ψ 1 ∂ψ 1 ∂ψ 2 + = ∂r r ∂r α ∂t and subject to ψ(0, t) = 0 and ψ(r, 0) = F (r) has the solution b ∞ 2 2 J0 (βm r) ψ(r, t) = r[ e−αβm t 2 J0 (βm r )]F (r )dr r =0 b2 m=1 J1 (βm b) b = r G(r, t|r , τ )|τ =0 F (r )dr r =0 Therefore, the Green’s function of the problem of ﬁnding T (r, t) for 0 ≤ r ≤ b satisfying ∂ 2T 1 ∂T 1 1 ∂T + g(x, t) = ∂x2 r ∂r k α ∂t and subject to T (b, t) = f (t) and T (r, 0) = F (r) is simply ∞ 2 2 J0 (βm r) G(r, t|r , τ ) = e−αβm [t−τ ] 2 J0 (βm r ) b2 m=1 J1 (βm b) and the desired solution is b α t b T (r, t) = r G(r, t|r , τ )|τ =0 F (r )dr + dτ r G(r, t|r , τ )g(r , τ )dr r =0 k τ =0 r =0 t ∂G −α [r ]f (τ )dτ τ =0 ∂r 14 5 Goodman Approximate Integral Method The integral method provides approximate solutions to heat conduction problems in which the nergy balance is satisﬁed only in an average sense. 5.1 Linear Problems The integral method is usually implemented in four steps: • The heat equation is ﬁrst integrated over a distance δ(t) called the thermal layer to obtain the heat balance integral. The thermal layer is selected so that the eﬀect of boundary conditions is negligible for x > δ(t). • A low order polynomial is then selected to approximate the temperature distribution over the thermal layer. The coeﬃcients in the polynomial are in general functions of time and must be determined from the conditions of the problem. • The approximating polynomial is substituted in the heat balance integral to obtain a diﬀerential equation for δ(t). • The resulting δ(t) is substituted in the polynomial formula to obtain the required temperature distribution. Consider the case of a semiinﬁnite medium initially at Ti where at the boundary x = 0 the temperature is maintained at T0 > Ti . The formulation is, ﬁnd T (x, t) satisfying ∂2T 1 ∂T 2 = ∂x α ∂t subject to T (x, 0) = Ti and to T (0, t) = T0 First integrate the heat equation from x = 0 to x = δ(t) to obtain x=δ(t) ∂T ∂T ∂T )= d( |x=δ(t) − |x=0 = x=0 ∂x ∂x ∂x 1 x=δ(t) ∂T 1 d x=δ(t) dδ = dx = [ ( T dx) − T |x=δ(t) ] α x=0 ∂t α dt x=0 dt 15 where Leibnitz rule has be used on the right hand side. Since (dT /dx)|x=δ(t) = 0 and T |x=δ(t) = Ti ∂T d x=δ(t) −α |0 = ( T dx − Ti δ) ∂x dt x=0 this is the heat balance integral equation. Now assume that the temperature distribution in the thermal layer can be represented as a low order (fourth order) polynomial of x, i.e. T (x, t) = a + bx + cx2 + dx3 + ex4 The coeﬃcients a, b, c, d and e are determined by introducing the conditions T (0, t) = T0 T (δ, t) = Ti ∂T =0 ∂x x=δ ∂2T =0 ∂x2 x=0 and ∂ 2T =0 ∂x2 x=δ Substituting into the polynomial, the result is T (x, t) − Ti x x x = 1 − 2 + 2( )3 − ( )4 T0 − Ti δ δ δ Substituting this last expression for T (x, t) into the heat balance integral results in 20 dδ α=δ 3 dt which must be solved subject to δ(0) = 0 to give 40 δ(t) = αt 3 This expression for δ(t) can now be substituted into the expression for T (x, t) or in q(0, t) = −k(∂T /∂x)|x = 0 to give 2k q(0, t) = (T0 − Ti ) δ 16 5.2 Nonlinear Problems Consider now instead the case of a semiinﬁnite medium initially at Ti = 0 where the boundary x = 0 is subjected to a prescribed, time dependent heat ﬂux f (t) and in which the thermal properties are all functions of temperature. The formulation is, ﬁnd T (x, t) satisfying ∂ ∂T ∂T (k(T ) ) = ρ(T )Cp (T ) ∂x ∂x ∂t subject to T (x, 0) = 0 and to ∂T −k |x=0 = f (t) ∂x Introducing now a new variable U deﬁned as U = 0T ρCp dT the problem transforms into ∂ ∂U ∂U (α(U ) )= ∂x ∂x ∂t subject to U (x, 0) = 0 and to ∂U −(α(U ) )|x=0 = f (t) ∂x Integrating over the thermal layer δ(t) and since (∂U/∂x)|x=δ = U |x=δ = 0, leads to the heat balance integral equation d δ U dx = f (t) dt 0 Now, assuming that the distribution of U is cubic in x and using the conditions of the problem to deﬁne the coeﬃcients leads to δf (t) x U (x, t) = (1 − )3 3αx=0 δ Substituting this into the heat balance integral gives d δ 2 f (t) [ ] = f (t) dt 12αx=0 In this case, however, the solution depends on αx=0 which is unknown but can be determined from the value of U at x = 0. The ﬁnal result is √ 4 t (U α)|x=0 = [ f (t) f (t )dt ]1/2 3 0 This one can be used to ﬁnd αx=0 (t). From that, δ(t) can be calculated, which in turn gives U , which in turn gives T . 17