Heat Conduction by Laplace Transform, Duhamel, Green Function and by klz10308

VIEWS: 0 PAGES: 17

									                    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 briefly 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 differential 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    Definitions
Consider a real function F (t). The Laplace transform of F (t), L[F (t)] = F (s) is defined 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 defined 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 finding 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) finally 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 find 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 final 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    Definitions
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, find Φ(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. Specifically, 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    Definitions
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 diffuses throughout the body. The mathematical formulation
of this problem consists of finding 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 infinite). 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 defined 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 infinite 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 semiinfinite 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 finding ψ(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 finding 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 finding ψ(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 finding 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 satisfied only in an average sense.

5.1      Linear Problems
The integral method is usually implemented in four steps:

    • The heat equation is first integrated over a distance δ(t) called the thermal layer to
      obtain the heat balance integral. The thermal layer is selected so that the effect 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 coefficients 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
      differential equation for δ(t).

    • The resulting δ(t) is substituted in the polynomial formula to obtain the required
      temperature distribution.

   Consider the case of a semiinfinite medium initially at Ti where at the boundary x = 0
the temperature is maintained at T0 > Ti . The formulation is, find 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 coefficients 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 semiinfinite medium initially at Ti = 0 where the boundary
x = 0 is subjected to a prescribed, time dependent heat flux f (t) and in which the thermal
properties are all functions of temperature.
   The formulation is, find 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 defined 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 define the coefficients 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 final result is
                               √            4          t
                            (U α)|x=0 = [ f (t) f (t )dt ]1/2
                                            3        0
This one can be used to find αx=0 (t). From that, δ(t) can be calculated, which in turn gives
U , which in turn gives T .

                                              17

								
To top