Lecture 8.  Diffusion

 

 

 

by Reinaldo Baretti Machín,

and Alfonso Baretti Huertas

 

 

 Free counter and web stats

http://www1.uprh.edu/rbaretti

http://www1.uprh.edu/rbaretti/methodsoftheoreticalphysics.htm

http://www1.uprh.edu/rbaretti/methodsoftheoreticalphysicsPart2.htm 

http://www1.uprh.edu/rbaretti/methodsoftheoreticalphysicsPart3.htm 

e-mail: reibaretti2004@yahoo.com

 

 Reference :

 

1. Mathematical Techniques for Biology and Medicine (Dover Publications- Paperback)

by William Simon ; http://www.amazon.com/Mathematical-Techniques-Biology-Medicine-William/dp/0486652475/ref=sr_1_1?ie=UTF8&s=books&qid=1223731765&sr=8-1

 

2. http://geocities.com/serienumerica4/Lecture1Biolologymedicine.doc

3. Mathematical Biology: I. An Introduction (Interdisciplinary Applied Mathematics) by James D. Murray

 

4. http://www.geocities.com/serienumerica4/Lecture3Biolologymedicine.doc

 

5. Schaum's Outline of Advanced Mathematics for Engineers and Scientists by Murray R Spiegel

6. Schaum's Outline of Laplace Transforms by Murray R Spiegel

 

7. http://www.mathworks.com/access/helpdesk/help/toolbox/symbolic/index.html?/access/helpdesk/help/toolbox/symbolic/laplace.html&http://search.yahoo.com/search?p=MATLAB+laplace+transfom&ei=utf-8&fr=b1ie7

 

8. Theoretical Physics (Paperback) by Georg Joos (Author), Ira M. Freeman , page 489.

 

9. Numerical Solution of Differential Equations by Ph. D., D. Sc. William Edmund Milne

 

10. Numerical Analysis by Richard L. Burden and J. Douglas Faires

 

11. Schaum's Outline of Numerical Analysis by Francis Scheid

 

 

Examples and problems are taken from:

Chapter 8 of  Reference 1. The problem numbers refer to the textbook listing.

 

 

 

The  heat flow equation is

 

            ∂2 u /∂ x2 = (1/α2) ∂ u /∂ t                                    (1)

 

where u can represent a  temperature  or a concentration of a diffusing substance like chemical or particles concentrations.  x ~ is position along the Xaxis , t time  and   α ~ length / time1/2  is the diffusivity constant which depends on specific properties of the material 

 

 

Consider an infinitesimal volume dV= dx dy dz =A dx or Ady or A dz from which heat is flowing out due to a temperature gradient at the walls.

 

For simplicity let the heat flow be only at the walls xith x= 0 and the the wall  x = dx. 

 

∂Q/∂t = k A ( (∂T /∂x)x=0 + dx - (∂T /∂x)x=0 )                (2)

 

where k~ energy/(time-length-kelvin) .

 

The substraction comes because the  area A and the the gradients are vectors.

The area at the right extreme (x=0) points left and the area at x=dx points to the right.

 

 Expanding (∂T /∂x)x=0 + dx ≈ (∂T /∂x)x=0  + (dx)( (∂ 2T /∂x2)x=0 converts (2) in

 

∂Q/∂t = k A (dx)( (∂ 2T /∂x2)x=0 = k  (∂ 2T /∂x2)x=0 dV    .   (3)

 

On the other hand   the calorimetry formula relates heat with mass, specficic heat c and change in temperature

 

                         ∆ Q = m c ∆ T = ρ ∆ V c ∆ T                     (4)

 

where ρ ~mass/volume , ∆ V element of volume , c ~energy/(mass-kelvin).

 

Passing to infinitesimals (4) becomes

 

 

                                  ∂Q/∂t = ρ c dV ∂T/∂t                      (5)

Equating (3) with (5) gives

 

                          k  (∂ 2T /∂x2) = ρ c ∂T/∂t              ,

or

                            (∂ 2T /∂x2) = (ρ c/k) ∂T/∂t           ,           

 

                                             ≡   (1/α2) ∂ T/∂ t              (6)

 

where     α2 = k/(ρ c)  ~ length2/time.

 

Section 12. Time dependent diffusion.

 

The analytical solution of (6) will be obtained for a uniform bar ( or infinite slab) that extends from x=0 to x= L . The initial temperature of the bar is 100 degrees , T(x,t=0)=100. The end point points ,(“ boundary conditions” ) are kept  at  0 degrees,  i.e.  T(x=0,t) =T(x=L,t) =0.  Heat diffuses out of the   bar from both ends and eventually will reach the equilibrium temperature of 0 degrees.

 

Fig 8.1 In this example an insulated rod or an infinite slab  is initially, at its interior points,  at a temperature T(x,0)=100 degrees while its extremes are kept fixed at T(x=0,t)=T(x=L,t) = 0. degrees at all times. Heat flow occurs only at the end  points.

 

For  convenience we rewrite now  the diffusion equation as

 

                                    ∂ T/∂ t     = K  (∂ 2T /∂x2)    ,     (7)

 

with K ~ length2/time .

 

 

The standard analytical method for the solution of (6) consists in assuming separation of variables. One takes

 

     T(x,t) =  g (t) h(x)                                                .        (8) 

 

Then eq(7) leads to

 

           h(x) d g(t)/dt = K g(t) d2 h(x)/dx2                      .     (9)

 

Divide (9) by g(t) h(x) to get   

 

(1/g) dg/dt   = K (1/h) d2 h(x)/dx2                   .   (10)

 

The left side of (10) is a function of t only and the right side is a function of the position x.

 

A  separation parameter  λ is introduced to satisfy

 

(1/g) dg/dt   = - β   =  K (1/h) d2 h(x)/dx2               .    (11)    

 

As will be seen shortly there is a whole set of parameters β ~1/time  , that are necessary when the boundary conditions are applied. For this reason  we append a subscript (n)  to read

 

             (1/g) dg/dt   = - β n  = K (1/h) d2 h(x)/dx2       (12)

 

 

The solution will in the from of a series     

 

                  T(x,t)  = ∑ n gn(t) hn(x)                                          . (13) 

 

Let         

                             gn (t)=exp(-β n  t)                               .     (14)

The DE for h(x) is

 

                 d2 h(x)/dx2  =    - (β n /K ) h(x) 

                                    ≡   - λ2 n   h(x)                             ,    (15)

where        λ2 n   =   β n /K ~ 1/length2  .

Equation (15) is the DE of the “harmonic oscillator” , whose  general solutions are the sum of sine and cosine functions.

 

Thus

hn(x) = An cos(λ n x) + Bn sin(λ n x )                          .     (16)

 

At this point we will justify the inclusion of the subscript n.

 

Applying the boundary conditions

T(x=0,t) = 0 = g (t) {An cos(λ n x) + Bn sin(λ n x )}x=0 ,   (17)

gives

                 0    = An = 0  .

At x=L  we have that ,

T(x=L,t) = 0 =g(t) { Bn sin(λ n x )  }x=L  .

 

Bn can not be zero because we are left without a solution.

So the trigonometric function has to satisfy

                  sin(λ n L) = 0    ,which implies

 

                  λ n L = n π     ( n=1,2,3,4…..)

 

and therefore    

                                   λ n = (n π /L )    .                              (18)

 

So the expression (13) is

 

          T(x,t) = ∑ n=1 Bn exp(-β n t) sin(λ n x)          ,             (19-a)

 

or

                 T(x,t) = ∑ n=1 Bn exp(-K(λ n )2 t) sin(λ n x)   .     (19-b)

       

with Bn to be determined from the initial temperature distribution T(x,t=0).                                        

                   

   Let T(x,t=0) = f(x) .

Then from (19)

               f(x) =  ∑ n=1 Bn  sin(λ n x)                       .     (20)

 

Multiply by a specific   sine term, say  sin (λ m x)  and integrate from from x=0 to x = L. We get

 

             L 0 f(x) sin (λ m x)  dx = (L/2) Bm  ,

or

               B m = (2/L) ∫ L 0 f(x) sin (λ m x)  dx             .   (21)

 

In our example f(x) =100.  The coefficients turn out to be

 

Bn = (4/(nπ)) (100)     for  odd n

 

Bn = 0         for   even n  .

 

Inserting    Bn  in (19) gives the final expression

 

T(x,t) = ∑ n=0 {400/((2n+1)π)} exp(-K (λ 2n+1 )2 t) sin(λ 2n+1 x)    .       (22)  

 

Example :  In eq(22) let  L= 1 meter , K=1.0E-3 m2/s , find T(x,t) various time instants. Tscale = L2/K = 1E3 seconds.

 

 

   

Fig 8.2 Plots of T(x,t) –eq (22) with L=1 meter , K= 1.E-3 m2/s

 

 

Fig 8.3 shows the result of applying a parabolic approximation to the diffusion problem given here. It is based on the forward difference method. From eq (7) we have  ( see references 9,10)

 

T(x,t+∆t) = T(x,t) +

                              (∆t K/ (∆x)2 ){ T(x+∆x,t) -2T(x,t)+T(x-∆x,t) }    . (23)

 

In eq(23) T(x,t) is known.

 

Let  x= L/2 and  ∆x=L/2  , taking into account the BC, T(x+∆x,t)= T(L,t)=0 and T(x-∆x,t)=T(0,t)=0. We have to calculate only one value of T at a time

t + ∆t , which is (since T(L/2,t) is known ) ,

 

 T(L/2,t+∆t) = T(L/2,t) + (∆t K/ (∆x)2 ){  -2T(L/2,t) }         .    (24)

 

The intervals ∆t and ∆x have to satisfy the inequality

                               (∆t K/ (∆x)2 << ½                                  .   (25)

 

 

          

Fig 8.3 Forward difference method.

 

 

At a given instant we have three values of the temperature

g0 = T(0,t)   , g1 = T(L/2,t)    , g2 = T(L,t)    at fixed points

x0 = 0 , x1=L/2 and x2=L  .

 

A Lagrange interpolation polynomial of the second order (see Ref 11 ,chapter 12) would be

 

T(x,t) ≈ g0*(x-x1)*(x-x2)/((x0-x1)*(x0-x2)

 

            + g1*(x-x0)*(x-x2)/((x1-x0)*(x1-x2))

      

            + g2*(x-x0)*(x-x1)/((x2-x0)*(x2-x1))                 (26)

 

 

In the actual example g0 = g2 =0 .

 

The parabolic approximation to T(x,t) is

 

T(x,t) = g1 ( x-x0) (x-x2)/ {(x1-x0)(x1-x2)}              ,

 

           =  g1 ( x) (x-L)/ {(L/2)(L/2-L)}           ,

 

           =    g1 ( x) (x-L)/ {(L/2)(L/2-L)}           ,

   

             =  -(4/L2) g1x(x-L)                 .                          (27)

 

Fig 8.2 Plot with a parabolic approximation. Compare with the exact solution in fig 8.1 .

 

 

 

 

FORTRAN code for eq (22).

c heat diffusion in a bar  or slab

      real L  ,k

      data L,k/1.,1.e-3/

      data nt,nx,nterms/4,10, 10/

      pi=2.*asin(1.)

      tscale=L**2/k

      deltat=.2*tscale/float(nt)

      deltax=L/float(nx)

      do 10 it=1,nt

      t=deltat*float(it)

      do 20 ix=0,nx

      x=deltax*float(ix)

      call temp(pi,nterms,L,K,x,t,tempx)

      print 100 , t,x,tempx

20    continue

      print*,'  '

10    continue

100   format('t,x,Temp=',3(3x,e10.3))

      stop

      end

     

      subroutine temp(pi,nterms,L,K,x,t,tempx)

      real L , K

      sum=0.

      do 10 n=0,nterms

      an=2.*float(n)+1.

      bn=(400./(an*pi))

      al=an*pi/L

      arge=k*al**2

      sum=sum+bn*exp(-arge*t)*sin(al*x)

10    continue

      tempx=sum

      return

      end

 

 

 

    FORTRAN code for a parabolic approximation

c parabola approximation to heat diffusion equation

      real L  ,k

      dimension  told(0:200),tnew(0:200)

      data L,k/1.,1.e-3/

      data nx,nx2,x0,x1,x2/2,10,0.,.5,1./

      pi=2.*asin(1.)

      deltax=L/float(nx)

      dx=L/float(nx2)

      tscale=deltax**2/(2.*k)

      dt=tscale/50.

      tfinal=200.

      nt=int(tfinal/dt)

      told(0)=0.

      told(1)=100.

      told(2)=0.

      do 10 it=1,nt

      t=dt*float(it)

      tnew(1)=told(1)+(dt*K/(deltax)**2)*(-2.*told(1))

      told(1)=tnew(1)

      tnew(0)=told(0)

      tnew(2)=told(2)

10    continue

      g1=tnew(1)

      g0=told(0)

      g2=told(2)

      do 20 ix=0,nx2

      x=dx*float(ix)

      call tlagr(g0,g1,g2,x0,x1,x2,x,tempx)

      print 100 , t,x,tempx

20    continue

      print*,'  '

100   format('t,x,Temp=',3(3x,e10.3))

      stop

      end

 

      subroutine tlagr(g0,g1,g2,x0,x1,x2,x,tempx)

      real Lg0,LG1,Lg2

      Lg0(x)=g0*(x-x1)*(x-x2)/((x0-x1)*(x0-x2))

      Lg1(x)=g1*(x-x0)*(x-x2)/((x1-x0)*(x1-x2))

      Lg2(x)=g2*(x-x0)*(x-x1)/((x2-x0)*(x2-x1))

      tempx=lg0(x)+lg1(x)+lg2(x)

      return

      end

 

 

 

Section 4. Chemical diffusion

 

Chemical diffusion through a thin layer, with no losses, obeys eq (1)

 

 ∂C(x,t)/ ∂t  = K ∂2C(x,t) /∂x2                                      (28)

 

where C~ moles/volume   , K ~ length2/time.

 

Experimentally it is found that there may be two kinds of losses of the substance C. The substance may decay and be lost at a rate – λC  , where

λ ~ 1/time. This is a a concentration dependent loss.

The rod or slice can be porous and a concentration independent loss exists at a rate  -R  ( R~ moles/(volume-time) )  may occur.

These two terms must be included in the right hand side of (28),

 

 

∂C(x,t)/ ∂t  = K ∂2C(x,t) /∂x2   -R - λC                        .  (29)                                  

 

Generalizing eq. (23) the forward difference solution to (29) is

 

C(x,t+∆t) = C(x,t)  

                            +  (∆t K/ (∆x)2 ){ C(x+∆x,t) -2C(x,t)+C(x-∆x,t) } 

                          

                            + ∆t ( -R – λC(x,t) )                           .     (30)

 

The inequality involving  ∆t and ∆x is now ,

 

              ∆t { 2K/  (∆x)2  +  λ } << 1                       .          (31)

 

Examples of steady state one dimensional diffusion:

In the steady state  ∂C(x,t)/ ∂t  =0 and (29 ) becomes an ordinary DE

 

             d2 C(x)/dx2  = + (λ/K)C (x) + (R/K)                     (32)

 

We may visualize the flow along the X axis in a slab like that of figure 8.1.

The flow , F(x1), at a plane x= x1 is defined by

 

F(x1) = - A K dC(x1)/dx  ~ length2 (length2/time)(moles/volume-length)

         

         ~ moles/time

 

 

a) Let the IC be C(0) = C0 , dC(0)/dx = C0prime , the solution of (32) are the complementary plus the particular solution.

 

C(x) = A exp(r x) +B exp(-rx)  -R/λ                           ,     (33)

 

where   r = (λ/K)1/2 .

 

Applying the IC gives

C0= A + B    -R/λ                                           ,        (34)

C0prime = r  A – r B   , or

C0prime/r = A – B                                         .         (35)

Adding or substracting (34)-(35) leads to

 

A = (1/2) { C0 + C0prime/r + R/λ }                ,       (36)

 

B = (1/2) { C0 – C0prime/r + R/λ }                .      (37)

 

If one fixes the IC to be C0=0 and C0prime =0 , then

 

A=B= R/(2 λ)     and

 

C(x)= R/(2 λ) { exp((λ/K)1/2 x)  +  exp(-(λ/K)1/2 x) } -R/λ       .    (38)            

 

Matlab code for plot of (38)

K=1; L=1; R=1;r=(L/K)^(1/2);

x=[0:.05:2.5];

 

C= -R/L + (R/(2*L))*(exp(r*x)+exp(-r*x));

plot(x,C),xlabel('x '),ylabel('Concentration')

 

                   

 

Fig 8.4 Concentration in a semi infinite medium with arbitrary values of λ=1, K=1 , R=1 and IC C(0)=0 , dC(0)/dx=0.

                    

   

 

Section 7 : Threshold dependent loss

 

 

Suppose we have two regions where eq (32) is modified such that

          

  d2 C(x)/dt2  = + (λ/K)C (x)             ,          x ≤ 0       (39)

 

and   

             d2 C(x)/dt2  = + (λ/K)C (x) + (R/K)   , x ≥ 0 .    (40)

 

 We call the C(x=0), the concentration  threshold

 

               C(x=0) = CT      (R=0) for  C ≤ CT    .

Let C1(x) be the solutions of (39) and (40 ) respectively.

The solutions must satisfy

 

 C1(x=0) = C2(0) = CT                                       (41)

 

and 

 

dC1(x=0) /dx = dC2(x=0)/dx                             (42)

 

 

The general solution of (39) is 

 

 C1(x) = A exp ((λ/K)1/2 x)  + B exp ( -(λ/K)1/2 x )     . (41)

 

But by requiring that C1(x = -∞ ) =0 and C1(x = 0) = CT , (41) reduces to

 

C1(x) = CT exp ((λ/K)1/2 x)         , x ≤ 0                    .          (42)

 

The derivative at the origin is

 

 dC1(0)/dx = ( λ/K)1/2 CT                                .      (43)

 

 

C2(x)  is the sum of two exponentials and the particular solution,

 

C2(x)  = D exp ((λ/K)1/2 x)  + E exp ( -(λ/K)1/2 x )   -R/λ  .  (43)

 

Applying the IC we get

 

C2(0) = CT = D + E  - R/λ                                                (44)

 

dC1(0)/dx = ( λ/K)1/2 CT =  ( λ/K)1/2 (D – E )     , or

          

                              CT = D – E .                               .        (45)  

Thus

 D= (  CT + R/(2λ) )   ; E =   R/(2λ)                           .      (46)

 

Inserting in (43)   gives

                                                                

 

C2(x)  = ( CT + R/(2λ) )exp ((λ/K)1/2 x)  

                                        +  ( R/(2λ)) exp ( -(λ/K)1/2 x )   -R/λ  ,

 

                                                            x ≥ 0 .  (47)

 

                   

   Fig 8.5 The arbitary values chosen are CT =1moles/volume ,

R=1 mole/(volume-time), L =1/time  , K= 1 length2 .

 

Fortran code for figure 8.5

c Simon plot fig 8.8

      real k , L ,lscale

      data r,L,K,ct/1.,1.,1.,1./

      lscale=sqrt(k/l)

      xi=-2.*Lscale

      xf= Lscale

      nstep=60

       dx=(xf-xi)/nstep;

      do 10 i=0,nstep

      x=xi+dx*float(i)

      if(x.le.0.)c=CT*exp(sqrt(L/K)*x)

      if(x.gt.0.)c=((CT+R/(2.*L))*exp(sqrt(L/K)*x)

     $+(R/(2.*L))*exp(-sqrt(L/K)*x)-(R/L) )

      print 100,x,c

10    continue

100   format(1x,'x,C(x)=',2(3x,e10.3))

      stop

      end

                         

 

 

Section 8. Semi-Infinite Section with Internal Generation

 

Recall eq (29)

 

∂C(x,t)/ ∂t  = K ∂2C(x,t) /∂x2   -R - λC   , where R represents a concentration independent loss. Suppose that instead of a loss, we have a concentration independent gain   I ~ moles/(volume-time) .

 

Then (29 would be

  ∂C(x,t)/ ∂t  = K ∂2C(x,t) /∂x2   + I  - λC      .          (47)  

 

In the steady state ∂C(x,t)/ ∂t  has become zero and

 

                      d2C(x) /dx2  -( λ/K) C      =  -(I/K)      .  (48) 

 

If we do not apply the criteria of a threshold ,the general solution is the particular solution  plus the complementary.

 

C(x) = (I/λ) + B1 exp( ( λ/K)1/2 x ) + B2 exp( -( λ/K)1/2 x )   . (49)

 

To find the coefficients  B1 and B2 one must know two conditions on C(x) , very much like in eqs (44) and (45) above.

 

Given the the two IC

C(0) =   C0 =      (I/λ) + B1  + B2                          (50)

 

 

           dC(0)/dx = ( λ/K)1/2 (B1  - B2 )                    (51)

 

allows us  to solve for B1 , B2 .            

In the next example  we solve the finite difference equation

 

C(xn) = 2 C(xn-1) – C(xn-2) + (∆x)2 {( λ/K) C(xn-1)    -(I/K) }      (52)

 

choosing    (∆x) < < (K/λ)1/2  . 

 

 

 

 

Example : IC , C(0) = 1.0E-2 moles/volume , dC(0)/dx =1 moles/length4 .         

Fig 8.6 Semi infinite section with internal generation I .

 

FORTRAN code

c Simon page 180

      real L , K ,lscale, Intern

      data L,k,Intern/1.,1.,1./

      data nstep/2000/

      lscale=sqrt(k/L)

      xi=0

      xf=2.*lscale

      dx=(xf-xi)/float(nstep)

      dcdx=1.

      c0=1.e-2

      c1=c0+dx*dcdx

      kp=int(float(nstep)/60.)

      kount=kp

      print 100, xi,c0

      do 10 i=1,nstep

      x=xi+dx*float(i)

      c2=2.*c1-c0+dx**2*((L/K)*c1-intern/K)

      if(i.eq.kount)then

      print 100, x, c2

      kount=kount+kp

      endif

      c0=c1

      c1=c2

10    continue

      print 100, x, c2

100   format(1x,'x,c(x)=',2(3x,e10.3))

      stop

      end

 

Section 9. Finite section with internal generation

 

Let the size of the slab section or rod be 2D. Assume again there is internal generation I and the DE is

                    d2C(x) /dx2  -( λ/K) C      =  -(I/K)   ,  -D ≤ x ≤ D.

 

The choice for the range of x is a matter of conevenience.

 

The general solution is ( see eq(48) )

C(x) = (I/λ) + B1 exp( ( λ/K)1/2 x ) + B2 exp( -( λ/K)1/2 x )   .  (53)

 

The IC may consist in specifying arbitrary values of  C(-D) and C(D) .

Or alternatively one could specify at any point the value of C(x) and its derivative at such point  e.g.  C(-D) , dC(-D)/dx.

 

One possible choice is the symmetric case C(-D) = C(D) ≡ A.

Then from (53) we get

 

      A  = (I/λ) + B1 exp( -( λ/K)1/2 D ) + B2 exp( ( λ/K)1/2 D ) , x=-D ,

      A  = (I/λ) + B1 exp( ( λ/K)1/2 D ) + B2 exp( -( λ/K)1/2 D ) , x = +D ,

implying B1 = B2 = B  .

 

Then

              A = (I/λ) + 2 B cosh(( λ/K)1/2 D)                          ,         (53)

where the hyperbolic cosine is  cosh(u) = (exp(u)+exp(-u))/2.

 

Then B ={ A- (I/λ) }/(2 cosh(( λ/K)1/2 D)       and

 

C(x) = (I/λ) + [A- (I/λ) }/(2 cosh(( λ/K)1/2 D) ] 2cosh( ( λ/K)1/2 x ) . (54)

 

Example: Plot (54) with the values A=1.e-2 , I=1 , λ=1 , K=1 , D=1 .

 

 

 

Fig 8.7 Finite section with internal generation.

FORTRAN code for the example

c Simon page 182

      real inter,K,L

      data inter,A,K,D,L/1.,1.e-2,1.,1.,1./

      data nstep/50/

      C(x) = (Inter/L) + ((A- (inter/L))/(2.*cosh(sqrt(L/K)*D) ))*

     $  2.*cosh(sqrt(L/K)*x )

      xi=-d

      xf=d

      dx=(xf-xi)/float(nstep)

      do 10 i=0,nstep

      x=xi+dx*float(i)

      print 100,x,c(x)

10    continue

100   format(1x,'x,C(x)=',2(3x,e10.3))

      stop

      end

 

 

 

 

Section 10. Spherical radial diffusion

 

(January 20 , 2009)

 

UNDER CONSTRUCTION