# Solution of the Diffusion Equation

Document Sample

```					                                   College of Engineering and Computer Science
Mechanical Engineering Department
Notes on Engineering Analysis
Larry Caretto November 17, 2011

Solution of the Diffusion Equation
Introduction and problem definition

The diffusion equation describes the diffusion of species or energy starting at an initial
time, with an initial spatial distribution and progressing over time. The simplest example
has one space dimension in addition to time. In this example, time, t, and distance, x,
are the independent variables. The solution is obtained for t  0 and 0  x  xmax. The
dependent variable, species or temperature, is identified as u(x,t) in the equations
below. The basic diffusion equation is written as follows.

u    2u
 2                                             [1]
t   x

Here, , is a species or thermal diffusion coefficient with dimensions of length squared
over time. The initial condition specifies the value of u at all values of x at t = 0. This
initial condition is usually written as follows:

u(x,0) = u0(x)                                       [2]

The boundary conditions at x = 0 and x = xmax, can be functions of time in general.
These are written as follows (with subscripts L and R denoting left and right).

u(0,t) = uL(t)        u(xmax,t) = uR(t)                         [3]

There is actually a discontinuity in the initial conditions and the boundary conditions.
The initial conditions are assumed to hold for all times less than or equal to zero. The
boundary conditions hold for t > 0. This discontinuity only occurs at the boundaries. For
example, the initial condition may be that u = 0 for all x and the boundary conditions
may be that u = 100 at both boundaries for t > 0. This creates a discontinuity at t = 0 for
x = 0 and x = xmax. Fortunately, this discontinuity does not affect our ability to solve the
problem.

For the simplest cases we will consider both the initial and boundary conditions to be
constant. However, in the most general case the boundary conditions may be given by
any functional forms for u0(x). uL(t), and uR(t)

Basic solution process
The solution of the diffusion proceeds by a method known as the separation of
variables. In this method we postulate a solution that is the product of two functions,

Jacaranda (Engineering) 3333               Mail Code                Phone: 818.677.6448
E-mail: lcaretto@csun.edu                    8348                     Fax: 818.677.7062
Solution of the diffusion equation                L. S. Caretto, November 17, 2011                       Page 2

T(t) a function of time only and X(x) a function of the distance x only. With this
assumption, our solution becomes.

u(x,t) = X(x)T(t)                                             [4]

We do not know, in advance, if this solution will work. However, we assume that it will
and we substitute it for u in equation [1]. Since X(x) is a function of x only and T(t) is a
function of t only, we obtain the following result when we substitute equation [4] into
equation [1].

u X ( x)T (t )             T (t )
                  X ( x)          
t         t                   t
[5]
 2u     2  X ( x)T (t )            2 X ( x)
 2                         T (t )
x             x 2                      x 2

If we divide the final equation through by the product X(x)T(t), we obtain the following
result.

1 1 T (t )     1  2 X ( x)
                                                       [6]
 T (t ) t   X ( x) x 2

The left hand side of equation [6] is a function of t only; the right hand side is a function
of x only. The only way that this can be correct is if both sides equal a constant. This
also shows that the separation of variables solution works. In order to simplify the
solution, we choose the constant1 to be equal to2. This gives us two ordinary
differential equations to solve.

1 1 dT (t )                        1 d 2 X ( x)
  2                                2                           [7]
 T (t ) dt                      X ( x) dx 2

dT (t )
The first equation becomes               2T (t ) . The general solution to this equation is
dt
 2t
known to be T (t )  Ae
2
. The second differential equation in [7] can be written as
d 2 X ( x)
2
 2 X ( x)  0 . The general solution to this equation is X(x) = Bsin(x) +
dx
Ccos(x). Thus, our general solution for u(x,t) = X(x)T(t) becomes

u( x, t )  T (t ) X ( x)  Ae t B sin(x)  C cos(x)  e  t C1 sin(x)  C2 cos(x)
2                                     2
[8]

The choice of – for the constant as opposed to just plain  comes from experience. Choosing the
1                  2

constant to have this form now gives a more convenient result later. If we chose the constant to be
simply , we would obtain the same result, but the expression of the constant would be awkward.
2
As usual, you can confirm that this solution satisfies the differential equation by substituting the solution
into the differential equation.
Solution of the diffusion equation                       L. S. Caretto, November 17, 2011                Page 3

We first consider a simple set of boundary conditions where u(0,t) = 0 and u(xmax,t) = 0.
In order to satisfy these boundary conditions on u for all times, we have to have
boundary conditions on X(x) that X(0) = X(xmax) = 0. With these homogenous boundary
conditions on X(x) and the form of the original differential equation for X(x), the solution
for X(x) is the solution of a Strum-Liouville problem. Accordingly, expect that we will be
able to obtain an infinite set of solutions for X(x) that will give a complete set of
orthogonal eigenfunctions.

If we substitute the boundary condition at x = 0 into equation [8], get the following result
because sin(0) = 0 and cos(0) = 1.

u(0, t )  0  e  t C1 sin(0)  C2 cos(0)  C2 e  t
2                                     2
[9]

Equation [9] will be satisfied for all t if C2 = 0. With C2 = 0, we can apply the solution in
equation [8] to the boundary condition at x = xmax.

u( xmax, t )  0  C1e  t sin(xmax )
2
[10]

Equation [10] can only be satisfied if the sine term is zero.3 This will be true only if xmax
is an integral times . If n denotes an integer, we must have

n
xm ax  n       or                                         [11]
xm ax

We have an infinite number of solutions for X(x), as the integer n goes from 1 to infinity.
We can write an individual solution as Xn(x) = sin(nx/xmax). Because the differential
equation and boundary conditions for X(x), form a Sturm-Liouville problem, we know
that the solutions to this problem, Xn(x), are an infinite set of orthogonal eigenfunctions.
Since any integral value of n gives a solution to the original differential equations, with
the boundary conditions that u = 0 at both boundaries, the most general solution is one
that is a sum of all possible solutions, each multiplied by a different constant, C n. This
general solution is written as follows

n
u ( x, t )   Cn e nt sin(n x)                        n 
2
with                                [12]
n 1                                                  xm ax

We still have to satisfy the initial condition that u(x,0) = u0(x). Substituting this condition
for t = 0 into equation [12] gives


 nx                     nx         
u0 ( x)  u ( x,0)      Cn e sin 
0
x
 max

    
 n 1
Cn sin 
x
 max



[13]
n 1

3
Of course, the boundary condition is also satisfied if C1 = 0, but this gives u(x,t) = 0, not a very
interesting solution.
Solution of the diffusion equation                             L. S. Caretto, November 17, 2011                                 Page 4

Since the eigenfunctions sin(nx/xmax) are a complete set of orthogonal eigenfunctions
in the interval x ≤ x ≤ xmax, we can expand any initial condition function, u0(x) in terms of
these eigenfunctions. We can obtain the usual equation for eigenfunction expansion
coefficients of Sturm-Liouville solutions by using the orthogonality relationships for
integrals of the sine. If we multiply both sides of equation [13] by sin(mx/xmax), where
m is another integer, and integrate both sides of the resulting equation from a lower limit
of zero to an upper limit of xmax – the upper and lower limits of the Sturm-Liouville
problem for X(x) – we get the following result.

 mx                               mx   nx 
x m ax                            x m ax 

0
 x dx 
u0 ( x) sin      
 max 
 C
0     n 1
n        x  sin  x dx
sin            
 max   max 
[14]
 x m ax
 mx   nx                                     mx 
x m ax

                      x  sin  x dx  Cm
Cn sin                                               x dx
sin 2      
n 1 0             max   max                       0             max 

In the second row of equation [14] we reverse the order of summation and integration
because these operations commute. We then recognize that the integrals in the
summation all vanish except for the one where n = m, because the eigenfunctions,
sin(nx/xmax), are orthogonal. This leaves only one integral to evaluate. Solving for Cm
and evaluating the last integral4 in equation [14] gives the following result.

 mx 
x max

 u ( x) sin x
0
dx

max 


x
2 max               mx 
xmax 
Cm  xmax
0
                      x dx
u 0 ( x) sin                                               [15]
 mx                              max 
  xmax 
sin 2      dx            0

0             

For any initial condition, then, we can perform the integral on the right hand side of
equation [15] to compute the values of Cm and substitute the result into equation [12] to
compute u(x,t). We recognize that equation [15] is the usual equation for eigenfunction
expansions in Sturm-Liouville problems.

As an example, consider the simplest case where u0(x) = U0, a constant. In this case
we find Cm from the usual equation for the integral of sin(ax).

x max                                        x max                                                        x
 mx                               mx       2U 0  xmax     mx 
max
2                                     2U 0
Cm 
xmax          u 0 ( x) sin
x    dx 

 max      xmax                   
sin      dx 

 xmax 

xmax  m

cos      

 xmax  0

0                                             0

Using a standard integral table, and the fact that the sine of zero and the sine of an integral multiple of 
4

are zero, we find the following result:
x max                                                           x max
2 mx      x x        2mx                                    xmax xmax     2mxmax                   x
   sin 
x    dx    max sin

 max 

 2 4m  xmax 

                                
2
     sin
4m  xmax

  0  max

        2
0                                                           0
Solution of the diffusion equation                         L. S. Caretto, November 17, 2011                                Page 5

2U 0  xmax     mxmax                   xmax     m0  2U 0
Cm              cos
 x                       
 m   cos
x     
       1  cosm                           [16]
xmax  m       max                               max  m

Since cos(m) is 1 when m is even and –1 when m is odd, the factor of 1 – cos(m) is
zero when m is even and 2 when m is odd. This gives the final result for Cm shown
below.

 4U 0
 m                                 m odd

Cm                                                                            [17]
0                                 m even


Substituting this result into equation [12] gives the following solution to the diffusion
equation when u0(x) = U0, a constant, and the boundary temperatures are zero.

ent sin(n x)                                                      n
2

4U 0
u ( x, t )       
 n 1,3,5         n
with         n 
xmax
[18]

If we substitute the equation for n into the summation terms we get the following result.
n 2  2 t

 nx 
x 
sin 
2
x m ax

e                               
4U 0                                             m ax 
u ( x, t ) 


n 1, 3, 5                         n
[19]

Equation [19] shows that the ratio u(x,t)/U0 depends on the dimensionless parameters
t/(xmax)2 and x/xmax. Thus we can compute a universal plot of the solution as shown in
Figure 1. Because the solution is symmetric about the midpoint x/L = 0.5, only the left
half is plotted here.

Although we can obtain reasonable solutions using equation [19] for “large” times, we
will not get a correct solution, within the limits of practical computers, if we try to
evaluate equation [19] for times near zero. Figure 2 shows the results of using equation
[19] to compute u(x,0)/U0 using various numbers of terms for t = 0. In this case we are
trying to use the sine series to evaluate X(x) = 1. We see that taking more terms makes
the solution more accurate away from x = 0, but no matter how many terms we take,
there are still strong oscillations near x = 0. This result is well known and is called the
Gibbs phenomena.

Nonzero boundary conditions

The use of zero boundary conditions to give us a Sturm-Liouville problem appears to be
a major limitation on this approach. We can solve a broader class of problems, still
using the same approach, by appropriate definitions of variables.
Solution of the diffusion equation                                                               L. S. Caretto, November 17, 2011                                      Page 6

Figure 1 -- Diffusion Equation Solutions
t = t/xmax2
1
t = 0.01                      t = 0.02
t = 0.0001
0.9                                                 0   05                                                         t = 0.03
0.             0   08                                                                   t = 0.04
t=               0.
0.8                                                     t=
t = 0.05
0.7

0.6
t = 0.08
u/U

0.5
t = 0.10
0.4
t = 0.125
0.3                                                                                                                                 t = 0.15
t = 0.175
0.2
t = 0.2

0.1
t = 0.30
0
0        0.05          0.1          0.15                     0.2          0.25      0.3     0.35         0.4      0.45                 0.5
x/L

Figure 2 -- Sine expansion for f(x) = 1
1.3

1.2

1.1

1

0.9

0.8
1 terms
f(x) = 1

0.7                                                                                                                                          2 terms
10 terms
0.6
25 terms
0.5                                                                                                                                          100 terms

0.4

0.3

0.2

0.1

0
0        0.02       0.04         0.06           0.08                   0.1       0.12    0.14    0.16        0.18         0.2

x

First we consider the case were the boundary potentials, u(0,t) and u(xmax,t) are both set
to a fixed value uB. In this case we can make the following simple variable
transformation from u(x,t) to v(x,t).
Solution of the diffusion equation                       L. S. Caretto, November 17, 2011                        Page 7

v(x,t) = u(x,t) – uB                                          [20]

This transformation gives v(0) = u(0,t) – uB = uB – uB = 0 and v(xmax) = u(xmax,t) – uB = uB
– uB = 0. Thus, v has homogenous boundary conditions.

Since uB is a constant, when we substitute equation [20] into equation [1], we get the
same equation in terms of v instead of u. Since the boundary conditions on v are that v
= 0 at x = 0 and x = xmax, the differential equation for u is that of equation [1], and
boundary conditions for v are the same as those used to derive the solution for u in
equation [12]. Accordingly, equation [12] is the solution for v(x,t). When we apply the
initial condition that u(x,0) = u0(x) to equation [12] for this revised problem, we get the
following equation in place of equation [13].


 nx 
u0 ( x)  u ( x,0)  v( x,0)  uB   Cn sin
x                                [21]
n 1     max 

Starting with this equation, we can repeat the derivation of equation [15] for Cm, which
gives the following result for this case where both x = 0 and x = xmax boundaries have
the same common boundary value, uB.

                                      mx 
u 0 ( x)  uB sin mx dx  2                                       cosm  1
x max                                         x max
2
Cm 
xmax    
0
x 
 max     xmax     
0
 x dx  uB
u 0 ( x) sin     
 max          m
[22]

Thus, the solution for u(x,t) for the general initial condition, u(x,0) = u 0(x) and the
boundary conditions that u(0,t) = u(xmax,t) = uB is given by the following equation.

     x m ax                       2
   u 0 ( x)  uB sin n x dxent sin(n x)
2
u ( x, t )  u B                                                                            [23]
xm ax n 1  0
                             


We can extend the procedure used above to consider the case where the two
boundaries are set at two different constant temperatures. That problem is specified by
the following equations, where uL and uR are constants.

u    2u
 2            with u ( x,0)  u0 ( x), u (0, t )  uL , u ( xm ax, t )  uR               [24]
t   x

The solution to this problem proceeds by defining u(x,t) as the sum of two functions
v(x,t) and w(x), both of which satisfy the differential equation in [24].

u(x,t) = v(x,t) + w(x)                                           [25]

Again, the function v(x,t) is defined to have a zero value at x = 0 and x = xmax, to give us
the Sturm-Liouville problem, which provides us with a complete set of eigenfunctions.
This zero boundary condition on v(x,t) tells us what the boundary conditions on w(x)
must be set to satisfy the boundary conditions in equation [24].
Solution of the diffusion equation                             L. S. Caretto, November 17, 2011           Page 8

w(0) = u(0,t) – v(0,t) = uL – 0 = uL                  w(xmax) = u(xmax,t) – v(xmax,t) = uR – 0 = uR   [25a]

Setting w(0) = uL and w(xmax) = uR means that the solution for u(x,t) proposed in
equation [25] satisfies the boundary conditions of equation [24]. Of course the
proposed solution must also satisfy the differential equation. Since both v(x,t) and w(x)
satisfy the differential equation in [24], their sum, u(x,t) = v(x,t) + w(x), which is the
solution that we seek, will also satisfy the differential equation. We then examine the
solution to the differential equations for both v(x,t) and w(x).

Since w(x) does not depend on t, we have the following result when we substitute w(x)
into the differential equation.

w       2w                             d 2w
 0  2                                  0                      [26]
t       x                              dx 2

That is, we have a simple, second-order, ordinary differential equation to solve for w.5
Integrating this equation twice gives the following general solution for w(x).

w(x) = C1x + C2                                     [27]

Since this equation has to satisfy the boundary conditions that w(0) = uL and w(xmax) =
uR, we have the following solution for w(x) that satisfies the differential equation and
boundary conditions.

uR  uL
w( x )  u L               x                             [28]
xm ax

Since v is defined to be zero at both boundaries, the solution for v(x,t) is the same as
the solution originally found for u(x,t) in equation [12]. Thus the solution for u(x,t) =
v(x,t) + w(x) is found by adding equation [12] for v(x,t) and equation [28] for w(x) to give.


 nx           uR  uL
 Cne t sin xmax   uL 
2
u ( x, t )  v( x, t )  w( x)                  
n
                          x        [29a]
n 1                                   xmax

Setting t = 0 in this equation gives the initial condition for u(x,t).


 nx         uR  uL
u ( x,0)  v( x,0)  w(0)             Cn sin xmax   uL 



             xmax
x             [29b]
n 1

Applying the general initial condition that u(x,0) = u0(x) gives the following result.

5
The w function does not depend on time. If the solution for v(x,t) is such that this function approaches
zero for long times, then u(x,t) = v(x,t) + w(x) becomes equal to w(x). Thus, w(x) represents the steady
state solution and we can regard the definition of u(x,t) = v(x,t) + x(x,t) as dividing the solution into a
transient part that decays to zero and a steady-state part.
Solution of the diffusion equation                                L. S. Caretto, November 17, 2011                 Page 9


 nx            uR  uL
u0 ( x )     Cn sin xmax   uL 



                 xmax
x                          [29c]
n 1

Repeating the derivation of equation [15] for Cm, in the same manner used to find Cm for
the boundary condition of two equal, non-zero boundary values in equation [22], gives
the following result in this case.

x m ax
2                             uR  uL       mx 
Cm 
xmax            u 0 ( x)  u L 
                  xmax
x  sin 
  dx 

  xmax 
0                                                                 [30]
2u 1  cos m 2u R  u L  cos m
x m ax
2                               mx 
xmax               u 0 ( x) sin 
x    dx  L

 max           m

m
0

Thus, the solution to the differential equation and boundary conditions in equation [24] is
given by the following equation, where equation [30] is used to find the values of C n.

uR  uL      
n
x   Cn e   nt sin(n x)                       n 
2
u ( x, t )  u L                                                           with                      [31]
xm ax      n 1                                                  xm ax

We see that the e  nt terms in the summation approach zero for large times leaving the
2

solution as uL + (uR – uL)x/xmax which is the steady state solution, defined in equation
[25] as w(x).

As in the case of zero boundaries, the easiest problem to examine is the one where
u0(x) is a constant, U0. In this case equation [30] becomes

2u 1  cos m 2u R  u L  cos m
x m ax
2                              mx 
Cm 
xmax               U 0 ( x) sin 
x    dx  L

 max            m

m                         [32]
0
2U 0  u L 1  cos m 2u R  u L  cos m
                          
m                   m

Since cos mπ = 1 for even m and -1 for odd m, we have the following result for Cm.

 2uR  uL 
 m                       m even

Cm                                                                    [33]
 4U 0  uL  2uR  uL 
 m                      m odd
                 m

We can substitute this result into equation [31] and divide by U0 – uL to get the following
dimensionless result.
Solution of the diffusion equation                                                                                   L. S. Caretto, November 17, 2011                                            Page 10

                                                                                                
2
 n 
u ( x, t )  u L u R  u L  x        2  1  xm ax 
                                                     t        nx 
                         6, n e  
U 0  u: L  xm ax  n  2, 4,

sin        

U 0  u: L                                                                                                          xm ax 
                                                                                                          [34]
2
 n 
2    u  uL   1                                            
x      t
              nx 
2  R                                                                      sin 
                                                                                       x 
 m ax 
e                                                             
  U 0  u: L n 1,3,5, n
                                                                                 m ax 

In this equation the dimensionless ratio, [u(x,t) – uL] / (U0 – uL) depends on the
dimensionless distance x/xmax, dimensionless time,  = αt/(xmax)2, and the ratio of initial
and boundary conditions, (uR – uL) / (U0 – uL). A plot of this result for (uR – uL) / (U0 – uL)
= 0.5 is shown in Figure 3. This plot shows that the steady-state solution is a linear
profile with the expected boundary values for [u(x,t) – uL] / (U0 – uL); these are zero at x
= 0 and (uR – uL) / (U0 – uL) = 0.5 at x = xmax. If we changed the value of the parameter
(uR – uL) / (U0 – uL) from 0.5 to some other value, the plot would be similar, but the
location of the dimensionless parameter at x = xmax would shift.

Figure 3 – Solution of diffusion equation for nonzero boundaries
1
02

0   5
t = 0.0001

8

0.0              1
0.0
.000

t=                0.0               2
0.9                                             t=                 0.0
t=

t=                3
t=0

0.0
t=                          5
0.0
0.8                                                                          .04        t=
t    =0

0.7                                                                                          8
0.0
(u(x,t) - uL)/(U0 - uL)

t=
0
0.1
0.6
t=
0.5
5
0.1
t=                             0.2
0.4                                                                                                             t=
0
0.3
0.3
t=
0
0.8
t=
0.2
(uR - uL) / (U0 - uL) = 0.5
t = (time)/(xmax)2
0.1

0

0                  0.1             0.2                 0.3              0.4                0.5               0.6         0.7               0.8       0.9   1
x

Cylindrical geometry

In a cylindrical coordinate system, 0 ≤ r ≤ R, the diffusion equation has the following
form.
Solution of the diffusion equation               L. S. Caretto, November 17, 2011        Page 11

u    1  u
     r                                   [35]
t    r r r

The most general initial and boundary conditions for the radial problem are

u(r,0) = u0(r)     u/r|r=0,t = 0      u(R,t) = uR(t)          [36]

We will first consider the conditions where uR is a constant. We will later see that we
will have to use an expansion of Sturm-Liouville eigenfunctions to match the initial
condition. Anticipating this, we will obtain a problem with zero boundary conditions by
splitting u(r,t) into two parts, as we did in equation [20].

u(r,t) = v(r,t) + uR                           [37]

The new function, v(r,t) satisfies the differential equation in [35] and has v/r = 0 at r =
0 and v = 0 at r = R. Since uR is a constant, the proposed solution in [37] will satisfy the
original differential equation and the boundary conditions from equations [35] and [36].
As before, we seek a solution for v(r,t) using separation of variables. We postulate a
solution that is the product of two functions, T(t) a function of time only and P(r) a
function of the radial coordinate, r, only. With this assumption, our solution becomes.

v(r,t) = P(r)T(t)                             [38]

We substitute equation [38] for u in equation [35]. Since P(r) is a function of r only and
T(t) is a function of t only, we obtain the following result.

u P(r )T (t )          T (t )
               P( r )         
t        t                t
1  P(r )T (t )
[39]
1  u                                     1  P(r )
     r          r               T (t )       r
r r r     r r       r                  r r   r

If we divide the final equation through by the product P(r)T(t), we obtain the following
result.

1 1 T (t )    1 1  P(r )
            r     2                         [40]
 T (t ) t   P(r ) r r   r

The left hand side of equation [40] is a function of t only; the right hand side is a function
of r only. The only way that this can be correct is if both sides equal a constant. As
before, we choose the constant to be equal to2. This gives us two ordinary
differential equations to solve.

dT (t )
 2T (t ) has the general solution T (t )  Ae t . The
2
The first equation becomes
dt
 P(r )
second differential equation in [40] can be written as       r      2 rP (r )  0 . The
r   r
Solution of the diffusion equation                            L. S. Caretto, November 17, 2011              Page 12

general solution to this equation is P(r) = BJ0(r) + CY0(r), where J0 and Y0, are the
Bessel functions of the first and second kind with zero order.6 Thus, our general
solution for u(x,t) = X(x)T(t) becomes

v(r, t )  T (t ) P(r )  Ae t BJ0 (r )  CY0 (r )  e t C1J 0 (r )  C2Y0 (r )
2                                       2
[41]

As r → 0, Y0(r) → -∞; to keep the solution finite, we require that C2 = 0. (We can also
show that the first derivative of J0 is zero at r = 0 so that we satisfy the zero gradient
condition in equation [36]. To satisfy the condition that v(R,t) = 0, we must satisfy the
following equation.

v( R, t )  0  C1e t J 0 (R)
2
[42]

Equation [42] defines a transcendental equation similar to that in equation [10] which
required sin(λxmax) = 0. That equation has a simple solution since we know that the
sin(nπ) = 0 when n is an integer. We do not have such a simple result for the equation
that J0(λR) = 0. However, the points at which J0 = 0 can be determined and are
available in tables. We use the symbol αmn to indicate the mth point where Jn is zero.
That is we define αmn by the following equation.

J n ( mn )  0             for        m  1,,                         [43]

For example, the first five points at which J0 is zero are for α10 = 2.40483, α20 = 5.52008,
α30 = 8.65373, α40 = 11.79153, and α50 = 14.93092. There are an infinite number of
such values such that J0(λmr), with λmR = αm0 provides a complete set of orthogonal
eigenfunctions that can be used to represent any other function over 0 ≤ r ≤ R. These
mn values are called the zeros of Jn.

We can now write the solution to our radial diffusion, from equation [38], as the sum of
v(r,t) and uR.


    Cm e   m t J 0 ( m r )  u R
2
u (r , t )                                                  m R   m0            [44]
m 1

6
Bessel functions, like sines and cosines, are solutions to differential equations. Although Bessel
functions occur less frequently than sines and cosines, they are functions that can be found in tables or in
computer function calculations. The general form of Bessel’s equation used for obtaining series solutions
d 2 y ( x) 1 dy( x) x 2  2
is                              y ( x)  0 . This equation can be transformed into
dx 2      x dx         x2
d  dy    2        2 
 z dz     z  k z  y  0 whose solution is AJ(kz) + BY(kz). We see that this second equation
dz                      
             
 P(r )
has the same form as        r      2 rP (r )  0 , provided that we set  = 0. This gives the result above
r   r
2
that the solutions are J0 and Y0. The factor of r in the λ rP(r) term is a weighting factor that must be
included in the definition for the inner product of solutions to the radial portion of the diffusion equation.
Solution of the diffusion equation                            L. S. Caretto, November 17, 2011                          Page 13

We still have to satisfy the initial condition that u(r,0) = u0(r). We can do this by using an
eigenfunction expansion.


u (r ,0)  u R  u0 (r )  u R             C m J 0 ( m r )                          [45]
m 1

The values of Cm are found from the general equation for orthogonal eigenfunction
expansions which includes a weighting function. Here the weighting function p(r) is
equal to r. For any initial condition, u0(r), we find the values of Cm from the following
equation. (In the second equation below, the integral in the denominator has been
evaluated using integral tables and the fact that J0(mR) = 0.)
R                                          R

 rJ 0 (m r )u0 (r )  uR dr             rJ     0   (m r )u0 (r )  u R dr
Cm      0
R
   0
[46]
R2
 rJ   0   (m r ) dr
2
J1 (m R)2
0
2

If we consider the simple case where u0(r) is a constant equal to U0, we have the
following result for Cm.
R                                                                               rR
 rJ1 ( m r ) 
 rJ 0 ( mr )U 0  uR dr                             
2U 0  u R    m  r  0

Cm     0
                                                 
R2
J1 ( m R)2                           R2               J1 ( m R)2
2

RJ1 ( m R)
2U 0  u R 
m               2U 0  u R    2U 0  u R 
Cm                                                                                                 [47]
R J1 ( m R)
2                 2             m RJ1 ( m R)  m0 J1 ( m0 )

With this result for Cm our solution for u(r,t) with the constant initial condition is written
as shown below after some algebra.

2U 0  uR   m 0 R 2 
                               t
2
r
u (r , t )                      e    J 0  m 0   uR                                   [48]
m 1  m 0 J1 ( m 0 )               R

This equation shows the important variables are r/R and αt/R2. We can rearrange this
equation to get a dimensionless form, plotted in Figure 4.
t
u (r , t )  uR   
2             m 0 2
2
    r
                      e       R
J 0  m0                              [49]
U 0  uR       m 1  m 0 J1 ( m 0 )                    R
Solution of the diffusion equation                                                      L. S. Caretto, November 17, 2011                                                             Page 14

Figure 4 – Solutions to radial diffusion with constant initial condition
1.0

tau
ta

tau
tau
u

=0
=
0.9                                                                                                           0.

.00
ta

=0
tau

=0
tau                                  u                    01
tau                                                        =

01
=0

.00

.00
=0
=0                    .06                               0.
.08                                                  02

2
04
0.8

.00
ta
tau                                  u

6
=
=0                                     0.
.1                                  04
[u(r,t) - uR] / (U0 - uR)

0.7

0.6                               tau
= 0.
15

0.5             tau = 0.2

0.4

0.3             tau = 0.3

0.2
tau = 0.4

0.1              tau = 0.5
tau = 0.6
tau = 0.8
0.0
0   0.1       0.2       0.3              0.4           0.5       0.6              0.7                   0.8                0.9             1
r/R

In Figure 4, the notation “tau” denotes the dimensionless time αt/R2. The results in
Figure 4 are similar to those in Figure 1 for the rectangular geometry. Here the plot
from 0 to R covers half the cylinder and is equivalent to the plot in Figure 1 that covered
only half of the one-dimensional region.

Mixed boundary condition

Boundary conditions on differential equations can be expressed in terms of the value of
u at the boundary, the gradient of u at the boundary, or a combination of the two. The
boundary condition which is an equation relating the value of u and its gradient at the
boundary is called a mixed boundary condition.

In this example we consider a problem of one dimensional heat conduction. Here we
will explicitly use the temperature, T, as the dependent variable (instead of u) and note
that α in this case is called the thermal diffusivity. We consider the case where a slab,
initially at a constant temperature, T0, is placed in a medium at a different temperature,
T∞. Convective heat transfer occurs on each side such that the results in the slab are
symmetric about the center. In such a case, we can solve for only half the geometry. If
we assume that the slab has a thickness of 2L, where –L ≤ x ≤ L, we solve for a region
between 0 and L using a symmetry boundary condition that the gradient is zero at the
center line. The convection boundary condition at x = L simply states that the
conduction heat flux in the solid at the boundary is the same as the convection heat flux
leaving the boundary. This is written as follows.
Solution of the diffusion equation                     L. S. Caretto, November 17, 2011              Page 15

dT
qx  L  k               h(TL  T )                         [50]
dx   xL

For this problem, we will define the following dimensionless variables that will give us a
Sturm-Liouville problem for the x solution in the separation of variables result.

T  T                     x              t
                                                                 [51]
T0  T                    L              L2

With these definitions, the diffusion equation and the initial and boundary conditions
may be written in the following dimensionless form.

  2                                                         hL
                ( ,0)  1                  0                     1  0      [52]
  2                                   0              1 k

We can perform our usual separation of variables solution to obtain the following
general solution.

( , )  e  C1 sin( )  C2 cos( )
2
[53]

The derivative of this solution is shown below.


 e    C1 cos( )  C2 sin( )
2
[54]


The gradient Θ/ξ will be zero at ξ = 0 only if C1 is zero. With C1 = 0, the boundary
condition at ξ = 1 becomes

       hL                               hL  2
     1  e   C2 sin( )     e C2 cos( )  0
2
[55]
  1 k                                  k

We can rearrange this equation to obtain a transcendental equation for λ.

k    cos( )
          cot( )                                     [56]
hL    sin( )

The roots of this equation can be visualized as the points were the plot of a straight line
through the origin with a slope of k/hL intersects with the plot of the cotangent. Such a
plot is shown in Figure 5 for hL/k = 1.
Solution of the diffusion equation                        L. S. Caretto, November 17, 2011                      Page 16

Figure 5 – Eigenvalues for convection problem
25

20

15

10

5
Functions

Line
0                                                                             Cotangent
Intersections
-5

-10

-15

-20

-25
0   2   4       6       8     10      12       14   16      18     20

x

In this problem, the eigenfunctions are cos(λmx), but the values of λm do not have a
simple equation such as they did in previous problems with simpler boundary
conditions. Instead the values of λm for this problem, similar to those for the
eigenvalues for J0, are not evenly spaced and have to be found by numerical solutions.
The first few values of λm for hL/k = 1 are7 λ1 = 0.860333589, λ2 = 3.425618459, λ3 =
6.437298179, λ4 = 9.529334405, λ5 = 12.64528722, λ6 = 15.77128487, and λ7 =
18.90240996. Changing the value of hL/k would change the slope of the line in Figure 5
and consequently change the intersection points. As m increases, the values of λm
increase and the intersections come closer to an integer value of π where the cotangent
goes to infinity. For example, λ7 = 18.90240996 is close to 6 π = 18.85.

Even though the eigenvalues, λm, are not evenly spaced, the part of the separation-of-
variables problem that was solved for the ξ dependence is a Sturm-Liouville problem, so
we know that we can expand any function, representing the initial condition, in terms of
these eigenfunctions.

As usual, the most general solution for Θ is a sum of all the eigenfunctions that satisfy
the differential equations and the boundary conditions.

hL
( , )   Cme  m cos(m )        m       cot(m )
2
[57]
m 1                                  k

7
You can observe that these values of l are the intersections indicated by the red dots in Figure 5.
Solution of the diffusion equation                         L. S. Caretto, November 17, 2011                           Page 17

If, for the moment, we regard T0 as some arbitrary measure of a variable initial
temperature profile, we can use an eigenfunction expansion to fit any initial condition,
Θ(ξ,0) = Θ0(ξ) = (T0(L) – T)/(T0 - T∞).

hL
( ,0)  0 ( )   Cm cos(m )                               m        cot(m )              [58]
m 1                                             k

We have the usual equation for the eigenfunction expansion, where we can find the
integral in the denominator from integral tables.
1                                              1

  ( ) cos(  )d
0              m                  2m  0 ( ) cos(m )d
Cm    0
          0
[59]
1
cos(m ) sin(m )  m
 cos (  )d
2
m
0

In particular for a constant initial temperature, T0, Θ0 = 1, and Cm is found as follows.
1
1
 sin(m ) 
2m  cos(m )d       2m            
Cm       0
       m  0             2 sin(m )
[60]
cos(m ) sin(m )  m cos(m ) sin(m )  m cos(m ) sin(m )  m

Substituting this expression for Cm into equation [57] gives our solution for the
dimensionless temperature, Θ, assuming a constant initial temperature.

2 sin(m )em cos(m )
2

hL
( , )                                                              m       cot(m )          [61]

m 1 cos( m ) sin(m )  m                                          k

The solution is plotted for hL/k = 1 in Figure 6 and for hL/k = 10 in Figure 7.
Both solutions are different from the earlier ones in which the boundary temperature
was set to a new value at zero time. Here, the boundary temperature decreases over
time because of convection heat transfer with the environment. With hL/k = 10, the
convection is ten times as strong as it is with hL/k = 1. This leads to a more rapid
change in temperature and a greater temperature gradient (i.e., greater heat transfer) at
the x = L boundary.

Consider the diffusion equation with zero gradient boundary conditions shown below.

u    2u                                u                     u
 2           0 x L                                                      0 u  x,0   f ( x)      [62]
t   x                                  x      x  0, t       x   x  L,t
Solution of the diffusion equation                                                  L. S. Caretto, November 17, 2011                                                     Page 18

Figure 6. Convection Boundary Condition, hL/k = 1
1

tau = 0.02
0.9                                            tau = 0.2            tau = 0.1              tau = 0.05

0.8
tau = 0.5
Dimensionless Temperature

0.7

tau = 0.75
0.6

tau = 1
0.5

0.4
tau = 1.5

0.3
tau = 2
0.2

tau = 3
0.1
tau = 4
0
0   0.1     0.2     0.3           0.4         0.5     0.6           0.7                  0.8            0.9                1
x/L

Figure 7. Convection Boundary Condition, hL/k = 10
1

tau = 0.1
ta

tau
u
ta

0.9                                                                                  ta
=
u

u
=0
0.
=

=
00
0.
0.

ta
.00
02
01

5

0.8                                                                   u
tau = 0.2                                      =                                                      3
0.
05
Dimensionless Temperature

0.7

0.6
tau = 0.35

0.5

tau = 0.5
0.4

0.3
tau = 0.75

0.2
tau = 1
0.1
tau = 1.5
0
0   0.1     0.2      0.3           0.4          0.5     0.6              0.7                  0.8             0.9               1

x/L

The same process that was used to solve equation [1] can be applied to equation [62] to
obtain the following solution.
Solution of the diffusion equation                 L. S. Caretto, November 17, 2011                              Page 19

2
              n 
   t
 nx 
u ( x , t )  C0   C n e        L 
cos                                      [63]
n 1                        L 
The constants in equation [63] are given by the following equation.

 mx 
L                                             L
1                           2
C0   f ( x)dx For m  1 : Cm   f ( x) cos     dx                                             [64]
L0                          L0            L 

As t → , the exponential terms will vanish and the only term that will be left in the
L
solution is the C0 term. This gives the result that ux, t    C0 
1
L    f ( x)dx .
0
Physically, this result tells us what happens when there is no flux into the region.
(Recall that the flux of physical quantities described by the diffusion equation is
the region, the final potential at all points in the region is simply the average of the initial
condition.

Nonzero gradients on two sides of a one-dimensional geometry

This section should be omitted at first reading. Go directly to the heading “Mixed
boundary condition in a cylinder” on page 25. If we want to consider two nonzero
gradients we have a problem that does not have a steady-state solution unless both
gradients are the same. Even if both gradients are the same the usual notion of a
solution in terms of two variables, u(x,t) = v(x,t) + w(x) will not have a unique solution.
To see this, consider the following problem.

u    2u                         u                     u
 2           0 x L                          g0                    gL   u  x ,0   f ( x )      [65]
t   x                           x   x  0, t          x    x  L,t

A solution of the form u(x,t) = v(x,t) + w(x) where v(x,t) satisfies the diffusion equation
with zero gradient boundary conditions and w(x) satisfies the equation d2w/dx2 = 0 with
the boundary conditions that dw/dx = g0 at x = 0 and dw/dx = gL at x = L will satisfy the
differential equation. However, the equation for w cannot satisfy its boundary
conditions. To see this, we integrate the equation d2w/dx2 = 0 two times to obtain the
general solution w = A + Bx.

Taking the first derivative of this general solution gives du/dx = B. This shows that
du/dx is a constant which must be the same at both boundaries so that B = g 0 = gL.
However, after we find the value for B, we have no other boundary condition for the
constant A. Thus, any solution u = A + g0x = A + gLx will satisfy the differential equation
and the boundary conditions only when g0 = gL.
Solution of the diffusion equation                             L. S. Caretto, November 17, 2011       Page 20

This mathematical result parallels physical reasoning. In the diffusion equation
derivation the gradient of u represents a flux. If we have a one dimensional region in
which we have a different flux on two sides, there cannot be a steady solution. How can
we find a solution to the problem posed in equation [65]?

Since this problem does not admit a steady state solution we can augment our usual
division of u(x,t) into a function v(x,t) that satisfies the diffusion equation with zero
boundary conditions and a function, w(x), that depends only on the space coordinate.
We have to include another function, (t) that provides the long term solution. We write
our proposed solution as follows.

u(x,t)=v(x,t)+ w(x) + (t)                        [66]

Here v satisfies the diffusion equation with zero gradient boundary conditions, and the
dw
function w(x) has the following gradient boundary conditions:            g 0 and
dx x  0
dw
 g L . We can show that the solution in equation [66] satisfies the original
dx x  L

u                   v                  (t )     w
                                         0  0  g0  g0          [67]
x    x 0,t         x   x 0,t          x x0,t x x0,t

u                    v                  (t )       w
                                             0  0  gL  gL      [68]
x    x  L ,t        x   x  L ,t        x x L ,t x x L ,t

The initial condition for v can de derived from the initial condition for u:

v( x,0)  u( x,0)   (0)  w( x)  f ( x)   (0)  w( x)             [69

If we substitute equation [66] into the diffusion equation and note that w(x) is a function
of x only and (t) is a function of time only, we obtain the following result.

u    2 u v  w       2 v  w   v   2v  2w
 2                                    2   2  0 [70]
t   x         t               x 2      t t x    x
Since v satisfies the diffusion equation, the v terms in the last expression cancel leaving
the following relationship between  and w.

(t )     2 w( x)
                                     [71]
t         x 2
Solution of the diffusion equation                L. S. Caretto, November 17, 2011              Page 21

We see that equation [71] has a function of t on one side and a function of x on the
other side. Just as in the solution of an original partial differential equation by
separation of variables, the only way that a function of time can equal a function of x is if
both functions equal a constant. Thus, we can rewrite this equation in the same way as
we usually apply separation of variables.

(t )     2 w( x)
           A                                  [72]
t         x 2
This gives two ordinary differential equations whose solutions are known.

d
A               At  B
dt                                                           [73]
d 2w A           Ax 2
         w       Cx  D
dx 2            2
We can apply the boundary conditions on w to evaluate A and C.

dw       Ax    
g0              C                             C  g0
dx x0         x 0                                                      [74]
dw       Ax                         g C   g  g0
gL              C                  A L    L
dx x L        x L                   L      L

With the constants just found we can write the solutions for (t) and w(x) as follows.

g L  g0                        g L  g0 2
(t )  At  B                  tB          w( x)            x  g0 x  D          [75]
L                              2L
The solution for v(x,t) is the solution to the diffusion equation with zero gradient
boundary conditions. This solution is an infinite series in the cosine of nx/L, which was
given in equation [63].
2                                      2
           n                                 n 
  t
 nx                     t
 nx 
v ( x, t )   C n e     L 
cos       C0   C n e     L 
cos           [76]
n 0                      L           n 1                    L 
The solution for u(x,t) = v(x,t) + w(x) + (t) is then found by combining equations [73]
and [76].
Solution of the diffusion equation              L. S. Caretto, November 17, 2011                Page 22

2
           n 
   t
 nx 
u ( x, t )  v( x, t )  w( x)  (t )  C0   C n e     L 
cos     
n 0                      L         [77]
A 2
    x  Cx  D  At  B
2
Before substituting the results of equation [74] for A and C, we will consider the initial
condition on v: v(x,0) = f(x) – (0) – w(x) = f(x) – B – (A/)(x2/2) – Cx – D. Satisfying this
initial condition requires


Ax 2                          nx 
v( x,0)  f ( x)  B            Cx  D  C0   Cn cos                        [78]
2                   n 1     L 
We can determine the coefficients Cn by an eigenfunction expansion. We use a
separate integration for the n = 0 eigenfunction, cos(0) = 1, because it has a different
normalization constant.

L 
L
              Ax2                  L
 nx 
  f ( x)  B 

              2                     
 Cx  D  (1)dx  C0 (1)dx 



n 1
C n cos
 L 
dx
0                                      0           0
L                                             [79]
L
                                                  L
A x3    x2                     L   nx 
   f ( x)dx   Bx 
       2 3
C
2
 Dx  C0 L  C n
0                

sin     C0 L
n   L  0

0                                              n 1

We thus find the following result for C0

L
1                 A L2   L
C0 
L 
f ( x)dx  B 
2 3
C  D
2
[80]
0

To get the remaining values of Cn, we multiply both sides of the initial condition equation
by cos(mx/L) and integrate from 0 to L. Applying the orthogonality relationship for the
cosine eigenfunctions gives the following result for m > 0.

L
              Ax 2            mx 
   f ( x)  B 
2
 Cx  D  cos      dx 
0                               L                                 [81]
L 
 nx   mx                  2  mx 
L
Cm L

 n0 Cn cos L  cos L dx   Cm cos  L   2
                                   
0                                    0

The integral for the terms B and D vanishes as shown below.
Solution of the diffusion equation                      L. S. Caretto, November 17, 2011                   Page 23

L                                                              L

 B  Dcos mx dx  B  Dsin mx   0
                   L                                             [82]
0
 L                        0

Combining equations [81] and [82] gives the following expression for Cm when m ≥ 1.

L
2             A 2        mx 
Cm 
L    
 f ( x) 
2
x  Cx  cos
  L 
dx m  1                                      [83]
0

If we substitute equations [80] and [83] into equation [77] we obtain the following
equation for u(x,t).

L
1                         A L2   L
u ( x, t ) 
L      f ( x)dx  B 
2 3
C  D 
2
0
2
 n 
 2 L                     mx     L                      t
 nx 
 
A 2
  f ( x)      x  Cx  cos   dx e                               cos           [84]
L 0
n0
2           L                                             L 

A 2
    x  Cx  D  At  B
2

The B and D terms in this equation cancel so these terms do not enter the solution.8
The following steps can be used to modify equation [84]: (1) use the equations in [74] to
replace A and C; (2) use the definition of the average initial condition,
1L
U 0   f ( x)dx , to replace the integral in equation [84]; and (3) use the symbol, Cm,
L0
in place of the full expression used in [84]. With these steps the following result obtains
after some algebra.
2
 n 
 g  g0 g0              t
 nx 
u ( x , t )  U 0  L L            C n e  L  cos     
    6     2  n1                 L                                   [85]
g  g0 2           g  g0
 L      x  g0 x   L        t
2L                   L
If we consider a problem with such a large value of t that the terms in the summation
are negligible, we can integrate equation [85] to find the (spatial) average value of u(x,t)
at any time. In this integration there is a large amount of cancellation leading to the
result shown below.

In many approaches to this problem the functions  and w are defined so that (0) = w(0) = 0. This
8

drops the constants B and D early in the derivation. However, there is no physical reason for setting
these conditions; they are only a convenience to rewrite the derivation in a simpler manner once we know
the solution.
Solution of the diffusion equation                  L. S. Caretto, November 17, 2011                  Page 24

1L                     g  g0
u (t )   u ( x, t )  U 0   L     t                                 [86]
L0                        L

In this case of large , the average value of u is increasing or decreasing with time
depending on the sign of gL – g0.

For large times, the exponential term in the summation of equation [85] will vanish;
removing this summation from equation [85] gives the large time limit for the solution
shown in equation [87].

 g  g0 g0  g L  g0 2          g  g0
Lim u ( x, t )  U 0  L L                x  g0 x   L     t                             [87]
t                        6    2      2L                  L

u   g  g0
At this point the change in u with time is                L     . This derivative is not a
t      L
function of x or t.9 At any value of x then, the change in u(x,t) at any x location, u(x,t2) –
u(x,t1), will be the same for any given time difference, t2 – t1.

g L  g0
Lim u ( x, t2 )  u ( x, t1 )              (t2  t1 )                    [88]
t1 , t 2                                  L
At large times then, the spatial profile of u(x,t) will not change with x. At any value of x
the change in u from one time point to another will be the same. This change in u will
be an increase if gL > g0; it will be a decrease if gL < g0, and there will be no change in u
if gL = g0. In the case where g0 = gL, equation [87] gives

    L
Lim u ( x, t )  U 0  g 0  x                    g L  g0                   [89]
t                            2

In this case where the flux is the same constant at both boundaries, the long-time
potential is the average of the initial condition plus a linear gradient whose slope is the
constant flux term, g0 = gL.

9
A comparison of equations [86] and [87] shows that two related quantities have the same value:
u (t )   g  g0        u ( x, t )
 L      Lim
t          L     t     t
The derivative on the left is the average (over x) at any time; the derivative on the right is the value at any
x for large times. The two equations are consistent. At large times, the derivative at any x has a constant
value. Thus the average (over x) will have the same value. Equation [86] tells us that the average (over
x) will have this same value at any (large) time.
Solution of the diffusion equation                L. S. Caretto, November 17, 2011                    Page 25

Mixed boundary condition in a cylinder

The problem of diffusion in a cylindrical coordinate system, 0 ≤ r ≤ R, for a fixed
boundary condition at the outer radius was treated above, starting with equations [35]
and [36]. If there is a mixed boundary condition at the outer radius of the cylinder, the
initial and boundary conditions for this problem become.

u(r,0) = u0(r)        u/r|r=0,t = 0             –ku/r|r=R,t = h(u – u)             [90]

In order to solve this problem, we will first define a new variable, v(r,t) as we did in
equation [37], to give a zero boundary condition at r = R.

v(r,t) = u(r,t) – u                                          [91]

With this definition, the boundary condition at r = R becomes

–kv/r|r=R,t = hv(R,t)        or         (k/h)v/r|r=R,t + v(R,t) = 0               [92]

The solution of the differential equation [35] proceeds in the same manner as done
previously, leading to the following equation, similar to equation [42].

v(r , t )  Ce     t
2
J 0 (r )                              [93]

We apply this solution to the boundary condition at r = R, using the result that dJ0(x)/dx
= –J1(x) to give the following result.

k v                          k
 v( R, t )  0   Ce   t J1 (R )  Ce   t J 0 (R)
2                     2
[94]
h r r  R                    h

Dividing by the constant, C, and the exponential term and defining  = R gives the
following results.

k                                     hR
      J1 ()  J 0 ()  0                 J 0 ()  J1 ()  0               [95]
hR                                      k

The roots of this equation give the eigenvalues  = /R. These eigenvalues will depend
on the value of hR/k. The first ten eigenvalues for various values of hR/k are shown in
the table below.

Table of Eigenvalues,  n = nR, for Various Values of hR/k
Index       hR/k =       hR/k =    hR/k =        hR/k =     hR/k =     hR/k =              hR/k =
,n         0.001        0.01       0.1            1         10        100                 1000
1      0.04471577 0.141245 0.441682           1.25578    2.17950    2.38090             2.40242
2          3.83197    3.83431   3.85771       4.07948    5.03321    5.46521             5.51456
3          7.01573    7.01701   7.02983       7.15580    7.95688    8.56783             8.64508
4          10.1736    10.1745   10.1833       10.2710    10.9363    11.6747             11.7797
Solution of the diffusion equation                                     L. S. Caretto, November 17, 2011                                  Page 26

5               13.3238           13.3244          13.3312           13.3984         13.9580           14.7834       14.9160
6               16.4707           16.4712          16.4767           16.5312         17.0099           17.8931       18.0530
7               19.6159           19.6164          19.6210           19.6667         20.0829           21.0036       21.1904
8               22.7601           22.7605          22.7645           22.8040         23.1710           24.1147       24.3281
9               25.9037           25.9041          25.9075           25.9422         26.2698           27.2264       27.4660
10               29.0469           29.0472          29.0503           29.0812         29.3767           30.3387       30.6040

With the eigenvalues, the general solution to the differential equation for v(r,t) is the sum
of all the individual eigenvalue solutions. The solution for u(r,t) is the sum of v(r,t) and
u .


 Cm e        m t                               m R   m hR / k 
2
u (r , t )                           J 0 ( m r )  u                                                          [96]
m 1

The constants in this equation are found by satisfying the initial condition that u(r,0) =
u0(r).


u (r ,0)  u  u0 (r )  u                 C m J 0 ( m r )            m R   m hR / k                        [97]
m 1

We use the orthogonality condition for the Bessel functions to compute the constants.
Multiplying equation [97] by rJ0(nr) and integrating from r = 0 to r = R gives the
following result, using Matlab to get the integral in the denominator. Note that this
integral has an additional term, [J0(mR)]2, that is not present in equation [46]. In that
equation, the eigenvalue equation, J0(mR) = 0, sets this term to zero.
R                                                   R

 rJ 0 ( m r )u0 (r )  u dr                     rJ 0 ( m r )u0 (r )  u dr
Cm      0
     0
[98]
J (          2  J1 ( m R)2 
R                                         2
R
 rJ 0 ( m r ) dr
2
0    m R)
2
0

As before, we consider the simplest case where u0(r) is a constant equal to U0; this
gives the following result for Cm.

R                                                                                           rR
 rJ1 ( m r ) 
   rJ 0 ( m r )U 0  u dr                                                  
2U 0  u           m  r 0
Cm              0

J (            2  J1 ( m R)2                       J 0 ( m R)2  J1 ( m R)2
2
R                                                       R2
0    m R)                                                                                                                         [99]
2
RJ1 ( m R)
2U 0  u 
m                                    2U 0  u J1 ( m R)                             2U 0  u J1 ( m )


R 2 J 0 ( m R)2  J1 ( m R)2                

 m R J 0 ( m R)2  J1 ( m R)2             

 m J 0 ( m )2  J1 ( m )2   
Solution of the diffusion equation                 L. S. Caretto, November 17, 2011              Page 27

Substituting this result for Cm into equation [96] for u(r,t) gives the solution for the
constant initial condition shown below after some algebra.

t

2U 0  u J1 (m )       2
   r
  J ( )2  J ( )2 e
m
u (r , t )                                            R2   J 0  m   u    [100]
m 1 m 0 m         1 m                              R

This equation can be rearranged into the following form.

                                            t
u ( r , t )  u                   2 J1 (m )              m
2
R2 J  r 
                                   e          0 m       [101]
m 1  m J 0 ( m )  J1 ( m )
U 0  u                               2             2
   R

Equation [101] shows the explicit dependence of a dimensionless potential difference
u (r , t )  u
on a dimensionless radial distance, r/R, and a dimensionless time, t/R2. In
U 0  u
addition, this solution depends on the ratio hR/k since the eigenvalues, m depend on
this parameter. The results of calculations for hR/k = 1 and hR/k = 10 are shown in the
Figures 8 and 9 as plots of dimensionless potential versus dimensionless radius with
dimensionless time, tau = t/R2 as a parameter. The comparison between these two
figures is similar to the comparison between Figures 6 and 7. Comparing Figures 6 and
7 showed that a higher values of the parameter hL/k give higher gradients at x = L and
shorter times for the values of u at any x location to get closer to u . In Figure 9, with a
higher value of hR/k than Figure 8, the gradients at r = R are steeper than in Figure 8
and the times required to reach a given dimensionless u value are shorter.

Spherical geometry

In a spherical coordinate system, 0 ≤ r ≤ the one-dimensional diffusion equation (no
angular dependence) has the following form.

u     1  2 u
        r                                      [102]
t    r 2 r   r

The most general initial and boundary conditions for the radial problem are

u(r,0) = u0(r)     u/r|r=0,t = 0      u(R,t) = uR(t)              [103]

We will first consider the conditions where uR is a constant. We will later see that we
will have to use an expansion of Sturm-Liouville eigenfunctions to match the initial
condition. Anticipating this, we will obtain a problem with zero boundary conditions by
splitting u(r,t) into two parts, as we did in equation [37] for cylindrical geometry.

u(r,t) = v(r,t) + uR                               [104]
Solution of the diffusion equation                                                                                  L. S. Caretto, November 17, 2011                  Page 28

Figure 8. Dimensionless Potential Difference in a Cylinder hR/k = 1.0

1

0.9

0.8
Dimensionless Potential Difference

tau = 0.01
0.7                                                                                        tau = 0.05
tau = 0.1
0.6                                                                                        tau = 0.2
tau = 0.3
0.5                                                                                        tau = 0.4
tau = 0.5
0.4                                                                                        tau = 0.6
tau = 0.8
0.3                                                                                        tau = 1
tau = 1.5
0.2

0.1

0
0       0.1    0.2   0.3      0.4      0.5     0.6      0.7   0.8   0.9    1

Figure 9. Dimensionless Potential Difference in a Cylinder hR/k = 10.0

1

0.9

0.8
Dimensionless Potential Difference

tau = 0.01
0.7                                                                                                         tau = 0.05
tau = 0.1
0.6                                                                                                         tau = 0.2
tau = 0.3
0.5                                                                                                         tau = 0.4
tau = 0.5
0.4                                                                                                         tau = 0.6
tau = 0.8
0.3                                                                                                         tau = 1
tau = 1.5
0.2

0.1

0
0          0.1    0.2    0.3      0.4     0.5      0.6      0.7   0.8    0.9       1

Solution of the diffusion equation            L. S. Caretto, November 17, 2011            Page 29

The new function, v(r,t) satisfies the differential equation in [102] and has v/r = 0 at r =
0 and v = 0 at r = R. Since uR is a constant, the proposed solution in [104] will satisfy
the original differential equation and the boundary conditions from equations [102] and
[103]. As before, we seek a solution for v(r,t) using separation of variables. We
postulate a solution that is the product of two functions, T(t) a function of time only and
P(r) a function of the radial coordinate, r, only. With this assumption, our solution
becomes.

v(r,t) = P(r)T(t)                                [105]

We substitute equation [105] for u in equation [102]. Since P(r) is a function of r only
and T(t) is a function of t only, we obtain the following result.

u P(r )T (t )         T (t )
              P(r )         
t     t                  t
[106]
1   u         1       P(r )T (t )           1    P(r )
 2 r2         2 r2                       T (t ) 2 r 2
r r r         r r             r               r r    r

If we divide the final equation through by the product P(r)T(t), we obtain the following
result.

1 1 T (t )    1 1  2 P(r )
              r     2                        [107]
 T (t ) t   P(r ) r 2 r   r

The left hand side of equation [107] is a function of t only; the right hand side is a
function of r only. The only way that this can be correct is if both sides equal a constant.
As before, we choose the constant to be equal to2. This gives us two ordinary
differential equations to solve.

dT (t )
 2T (t ) has the general solution T (t )  Ce   t .
2
The first equation becomes
dt
The second ordinary differential equation in [107] can be written as follows.

d 2 dP(r ) 2 2
r        r P( r )  0 .                           [108]
dr   dr

We can show that the following proposed solution satisfies this differential equation.

sin r    cosr
P( r )  A          B       .                         [109]
r        r

Differentiating this proposed solution gives.

dP(r )      cosr      sin r   sin r cosr
 A        B        A 2 B 2 .                         [110]
dr           r           r       r     r
Solution of the diffusion equation                      L. S. Caretto, November 17, 2011                 Page 30

Substituting this first derivative of the proposed solution into the first term in the
differential equation in [108] gives.

d 2 dP(r ) d 2  cosr                sin r     sin r   cosr 
r          r  A              B        A 2 B 2 
dr      dr      dr          r            r         r        r 
 Ar cosr  Br sin r  A sin r  B cosr  
d
.                 [111]
dr
A cosr  2 Ar sin r  B sin r  2 Br cosr  A cosr  B sin r
 2 Ar sin r  2 Br cosr

Substituting this result and the proposed solution into equation [108] shows that the two
terms in the differential equation sum to zero; thus, the proposed solution satisfies the
differential equation.

d 2 dP(r )
r          2 rA sin r  2 rB cos r
dr        dr                                               .             [112]
 sin r     cos r 
2 r 2 P(r )  2 r 2  A      B            2 rA sin r  2 rB cos r
     r         r 

At r = 0, the cos(r)/r term in the solution becomes infinite. In order to retain a finite
solution at the center of the sphere we must set B = 0. The sin(r)/r term becomes zero
at the center of the sphere, providing a finite term in the solution. With B = 0, the
product solution, v(r,t) = P(r)T(t) becomes.

sin r            sin r
v(r , t )  T (t ) P(r )  ACe        t
 De  t
2                      2
.          [113]
r                 r

We now have to satisfy the condition that v(r,t) = 0 at r = R,

sin R
v( R, t )  T (t ) P( R)  0  De  t
2
.                [114]
R

This gives the expected eigenvalue solution that R must be an integral multiple of .

n
R  n         or           n                              [115]
R

The general solution for v(r,t) is a sum of all eigenvalue solutions, each multiplied by a
different constant, Dn. From equation [104], we know that the solution we want for u(r,t)
= v(r,t) + uR. Thus, our general solution for u(r,t) is written as follows


sin( n r )                                      n
 Dne  t
2
u (r , t )                n                 uR              with        n          [116]
n 1
r                                            R
Solution of the diffusion equation                      L. S. Caretto, November 17, 2011                 Page 31

We still have to satisfy the initial condition that u(r,0) = u0(r). Substituting this condition
for t = 0 into equation [116] gives


Dn  nr 
u0 (r )  u (r ,0)         r
sin      uR
 R 
[117]
n 1

As usual, the eigenfunctions (1/r)sin(nr/R) form a complete set of orthogonal
eigenfunctions in the interval 0 ≤ r ≤ R; however, the Sturm-Liouville problem defined in
equation [108] has a weighting factor of r2 that must be used in the orthogonality
integral. If we multiply both sides of equation [117] by r2(1/r)sin(mr/R)dr, where m is
another integer, and integrate from zero to R, we get the following result. The steps
used here are similar to those in equation [14].

 mr                       nr   mr 
R                       sin          R           sin       sin    
 R  dr                    R   R dr
   u0 (r )  u R r                        
2                           2
Dn r
r         n 1
r           r
0                                            0                                          [118]
 R                                              R
 mr   nr            2  mr 
     Dn sin

 sin 
R   R 
dr  Dm sin 
 R dr
n 1 0                                 0

Solving for Dm and evaluating the final integral gives the following result, which is similar
to equation [15].
R
 mr 
 ru0 (r )  u R sin
    R 
dr           R
 mr 
r u0 (r )  u R sin 
2
Dm      0
R
 mr 

R                       R 
dr    [119]

 sin 2 
 R 
dr                   0

0

As an example, consider the simplest case where u0(r) = U0, a constant. In this case
we find Dm as follows.

2U 0  u R 
R                                                    R
 mr                           mr 
r u0 (r )  u R sin 
2
Dm 
R                       R 
dr 
R
r sin 
 R dr
0                                                    0
2U 0  u R   R    mr   mr   mr 
2                                         R
                      sin             cos                              [120]
R       m    R   R   R  0
2 RU 0  u R 
sin m  sin(0)  m cosm  0   2 RU 0  u R cosm
m2                                                           m

Since cos(m) is 1 when m is even and –1 when m is odd, we can write the final result
for Dm as shown below.
Solution of the diffusion equation                      L. S. Caretto, November 17, 2011                              Page 32

 2 RU0
 m                m odd
2 RU 0  u R  1m 1       
Dm                                                                                   [121]
m                     2 RU0
  m             m even


Substituting this result into equation [116] gives the following solution to the diffusion
equation when u0(x) = U0, a constant, and the boundary temperatures are zero.

2 RU 0  u R    1m 1  2 t sin( n r )                                         n
u (r , t ) 
              m
e n
r
 uR             with           n 
R
[122]
n 1

We can rearrange this equation to introduce nr in the summation and create the
following dimensionless form.

n 2 t
u (r , t )  u R     
 1n 1 e  2n t sin( r )  2   1n 1 R e                       nr
U 0  uR
2       nr
n            nr
R2      sin(
R
)    [123]
n 1                                   n 1

The dimensionless ratio on the left of this equation depends on the dimensionless
radius r/R and the dimensionless time, expressed as a Fourier number, t/R2.

Spherical geometry with mixed boundary condition

We still have the differential equation in [102], but we now have the following boundary
conditions:

u(r,0) = u0(r)      u/r|r=0,t = 0           –ku/r|r=R,t = h(ur=R – u)                     [124]

In order to solve this problem, we will first define a new variable, v(r,t) as we did in
equation [91], for the mixed boundary condition for the cylinder, to give a zero boundary
condition at r = R.

v(r,t) = u(r,t) – u                                                [125]

The boundary conditions for this new variable, v(r,t) are shown below.

v(r,0) = u0(r) – u v/r|r=0,t = 0               –kv/r|r=R,t = hvr=R                      [126]

The general solution for v(r,t) is will be the same as the solution to the problem with
fixed temperature given in equation [113]. This equation is copied below.

sin r            sin r
v(r , t )  T (t ) P(r )  ACe    t
 De  t
2                      2
.                       [113]
r                 r
Solution of the diffusion equation                     L. S. Caretto, November 17, 2011                     Page 33

In order to satisfy the boundary condition at r = R, –kv/r|r=R,t = hvr=R, with equation
[113], we have to satisfy the following equation.

v                                  cos R sin R     sin R 
 hvr  R  0  De  t k  
2
kv                                                   2 
h          .                  [127]
r r  R                               R     R           R  

Dividing the final equation by (De-2tk/R2)sin R gives.

cos R     hR     R       hR              hR 
R           1            1     0   R  1     tan R                             [128]
sin R      k   tan R      k                 k 

The roots of this equation are the eigenvalues of the solution to the present problem.
These eigenvalues will depend on the dimensionless group hR/k. We can compare this
equation with equation [56], which gives the eigenvalues for the rectangular slab with a
convection boundary condition. The roots of equation [56] are the values of  that
satisfy  = hL/k cot(). If we define  = R in equation [128], the roots of that equation
are the values of  that satisfy  = (1 – hR/k) tan(). For hR/k = 1, this equation can be
written as  cot() = 0. Thus, the eigenvalues for hR/k = 1 are simply the zeros of the
cotangent, (2n – 1)/2, where n is an integer greater than or equal to one.10 For other
values of hR/k, we can write the eigenvalue equation as /(1 – hR/k) = tan(). For hR/k
< 1, the roots of the equation [128] will be the intersection of the tangent curve with a
line, starting at zero and having a positive slope. This will yield eigenvalues as shown in
Figure 10. The straight line for values of hR/k greater than one will have a negative
slope as shown in Figure 11.

Regardless of the behavior of the eigenvalues, the resulting solution will still have the
same form as equation [116]; however, this solution has u in place of uR, and the
eigenvalues are different.


sin( n r )                                nR
 Dne t
2
u (r , t )           n                  u           with                  tan  n R    [129]
r                                        hR
n 1                                                   1
k

We still have to satisfy the initial condition that u(r,0) = u0(r). Substituting this condition
for t = 0 into equation [129] gives



Dn
u0 (r )  u (r ,0)            sin  n r  u                        [130]
n 1
r

As in equation [117], the eigenfunctions (1/r)sin(nr) form a complete set of orthogonal
eigenfunctions in the interval 0 ≤ r ≤ R, with a weighting factor of r2 that must be used in
the orthogonality integral.

Note that  = 0 is not a root of cot() = 0. As  → 0, cot() → 1.
10
Solution of the diffusion equation                                L. S. Caretto, November 17, 2011              Page 34

Figure 10. Sphere Convection Eigenvalues for hR/k = 0.1

20

15

10

5
Functions

Line
0                                                                                Tangent
Intersections

-5

-10

-15

-20
0       2          4   6     8    10      12     14     16      18    20



Figure 11. Sphere Convection Eigenvalues for hR/k = 2.0

20

15

10

5
Functions

Line
0                                                                                      Tangent
Intersections

-5

-10

-15

-20
0       2          4       6    8     10      12     14      16     18     20


Solution of the diffusion equation                           L. S. Caretto, November 17, 2011                      Page 35

Repeating the process used in obtaining equation [118] – multiplying both sides of
equation [130] by r2(1/r)sin(mr)dr, where m is another eigenvalue, and integrating from
zero to R – gives the following result.

R                                              R 

    u0 (r )  u r sin  m r dr 
2
  Dn r
2    sin  m r sin  n r
dr
r                                      r        r
0                                             0 n 1
[131]
 R                                            R
      Dn sin  mr sin  n rdr  Dm  sin 2  mrdr
n 1 0                                          0

Solving for Dm and evaluating the final integral gives the following result, which is similar
to equation [119].

R                                             R

   u0 (r )  u r sin  m r dr
2
 u0 (r )  u r
2   sin  m r
dr
r                                                 r
Dm      0                                            0                                        [132]
R                                       R sin  m R cos m R

 sin        m rdr
2
2         2 m
0

The integral of sin2(mr) was found by symbolic integration using Matlab. This integral
has an additional term that was not present in equation [119], where mR = m which
gives sin(mR) = 0.

As usual we will apply this equation to the simplest case where u0(r) = U0, a constant.
In this case we find the numerator in the integral for Dm as follows.

R                                                     R
 R sin  m R cos m R 
 
2                      Dm  r u0 (r )  u sin  m r dr  U 0  u  r sin  m r dr
                                                                         
           2 m            0                                           0
2
 1 
   sin  m r   m r cos m r 0 
 U 0  u     
R                                     [133]
 m
U 0  u
sin  m R  sin(0)   m R cos m R  0  U 0  u sin  m R   m R cos m R
m
2
2m

Solving for Dm and performing some algebraic rearrangement gives the final result for
Dm.

sin  m R   m R cos  m R
Dm  2U 0  u                                                                   [134]
R2   m sin  m R cos  m R
m

Substituting this expression for Dm into equation [129] gives the following solution for
u(r,t) in this case.
Solution of the diffusion equation                             L. S. Caretto, November 17, 2011                                Page 36


sin  R   n R cos n R                            sin( n r )
 R2  
u (r , t )  2U 0  u                                                   e   n t
2
n                                                                      u       [135]
n 1          n       n sin  n R cos n R                           r

We can obtain a dimensionless form by introducing m = nR and performing some
algebraic rearrangement.

 r 
                                  2
t   sin  n 
u (r , t )  u        sin  n   n cos n                               R 

n
2                          e                   R2                                    [136]
U 0  u                sin  n cos n                               n r
n 1 n
R

u ( r , t )  u
This equation shows that the dimensionless ratio,                                             , depends explicitly on the
U 0  u
dimensionless radius, r/R and the dimensionless time (Fourier number) t/R2. In
addition, there is an implicit dependence on the dimensionless group hR/k, known as
the Biot number, through the eigenvalue n. A plot of equation [136] for hR/k = 1 is
shown in Figure 12.

Figure 12. Dimensionless Potential Difference in a Sphere with hR/k = 1.0
1

0.9

0.8
tau = 0.01
tau = 0.1
0.7
tau = 0.2
tau = 0.3
0.6                                                                                                                tau = 0.4
tau = 0.5
0.5                                                                                                                tau = 0.6
tau = 0.8
0.4                                                                                                                tau = 1
tau = 2
0.3                                                                                                                tau = 3
tau = 4
tau = 5
0.2

0.1

0
0        0.1     0.2        0.3         0.4      0.5      0.6        0.7            0.8        0.9      1

Dimensionless time, t/R  2

```
DOCUMENT INFO
Shared By:
Categories:
Tags:
Stats:
 views: 12 posted: 11/17/2011 language: English pages: 36