Lecture 8.  Diffusion - Part 2

 

 

 

by Reinaldo Baretti Machín

and Alfonso Baretti Huertas

 

www.geocities.com/serienumerica4

www.geocities.com/serienumerica3

www.geocities.com/serienumerica2

www.geocities.com/serienumerica

 

For questions and suggestions write to:

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.

 

 

Section 10. Spherical radial diffusion

 

To generalize the heat diffusion equation to spherical coordinates we make use of the Laplacian operator in spherical coordinates

 

 \Delta f 
= {1 \over r^2} {\partial \over \partial r}
  \left( r^2 {\partial f \over \partial r} \right) 
+ {1 \over r^2 \sin \theta} {\partial \over \partial \theta}
  \left( \sin \theta {\partial f \over \partial \theta} \right) 
+ {1 \over r^2 \sin^2 \theta} {\partial^2 f \over \partial \phi^2}.

                                       

                                                                                                    (55)

where our scalar function is T(r,θ,φ) instead of  f.

 

Suppose T is only a function of the radial coordinate r  , then all the previous DE where we had the second derivative with respect to x  become ,

 

K d2 C/dx2 → (1/r2) ∂ (r2  ∂C/∂r )/∂r = d2 C(r)/dr2 +(2/r)(dC/dr).

 

Instead of

∂C(x,t)/ ∂t  = K ∂2C(x,t) /∂x2   -R - λC      , we have for the steady state

 

ie. (∂C(x,t)/ ∂t  =0),

 

0 = d2 C(r)/dr2 +(2/r)(dC/dr)    -R - λC            or

 

 d2 C(r)/dr2 +(2/r)(dC/dr)    =  + R + λC                      .  (56)   

 

Suppose we know that C(r=0) ≠ ∞  then multilplying (56) by r and taking the limit as r goes to zero gives

 

lim r→ 0 ( r d2 C(r)/dr2 + 2(dC/dr)  ) =  lim r→ 0 ( r R+ r λC )  =0 ,

thus at the origin

 

                                  (dC/dr )r=0 = 0.                              (57)

 

The finite difference solution of (56) is

 

Cn = 2Cn-1 – Cn-2 + (∆r)2{ -(2/rn-1)( (Cn-1 –Cn-2)/ ∆r  ) + R/K+ (λ/K)Cn-1 } .                                                  

                                                                                                 (58)

 

Example : Let R=1, λ=1,K=1 , rsphere =1 with  initial conditions (IC)

 

 C(0)=1, dC(0)/dr=0.

Fig. 8.8 Spherical diffusion. IC C(0)=1, dC(0)/dr=0.

 

FORTRAN code

c spherical diffusion Simon page 184

      real L ,lscale,K

      data R,L,K,rsphere/1.,1.,1.,1./

      data nstep/1000/

      lscale=sqrt(K/L)

      dr=rsphere/float(nstep)

      c0=1.

      dcdr=0.

      c1=c0+dr*dcdr

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

      kount=kp

      print 100, 0.,c0

      do 10 i=2,nstep

      r=dr*float(i)

      c2=2.*c1-c0+dr**2*(-(2./(r-dr))*dcdr + (R/K) +(L/K)*c1 )

      if(i.eq.kount)then

      print 100, r, c2

      kount=kount+kp

      endif

      c0=c1

      c1=c2

      dcdr=(c1-c0)/dr

10    continue

      print 100, r, c2

100   format(1x,'r,C= ', 2(3x,e11.4))

      stop

      end

 

 

 

Section 10. Cylindrical diffusion and the Krogh equation

 

Laplace operator in cylindrical coordinates has the form

 

 \Delta f 
= {1 \over \rho} {\partial \over \partial \rho}
  \left( \rho {\partial f \over \partial \rho} \right) 
+ {1 \over \rho^2} {\partial^2 f \over \partial \theta^2}
+ {\partial^2 f \over \partial z^2 }. 
             ,    (58)

 

ρ = (x2 + y2 ) 1/2  ~  radial coordinate.

 

Assuming no dependence on either θ or z   and multiplying the laplacian by the parameter K we have the diffusion equation ,

 

      ∂C/ ∂t  = K [(1/ρ )∂ {ρ (∂C /∂ρ)}/ ∂ρ ] +   -R - λC                . (59)

 

Furthermore assuming , the steady stae has been reached (∂C/ ∂t  =0), and that there is no concentration dependendent loss (59) reduces to,

 

 K [(1/ρ )∂ {ρ (∂C /∂ρ)}/ ∂ρ ] = R                                            .  (60)

 

One can integrate (60 ) analytically in two steps. First

 

 

    d/dρ (ρ dC/dρ) = (R/K) ρ      , thus

 

 ρ dC/dρ = (1/2) (R/K) ρ2 + B1  ,                (61)

 

where  B1 is a constant of integration.  Solving for dC/d ρ ,

we get

                  dC/dρ =  (1/2) (R/K) ρ + B1/ ρ                      . (62)

 

Integrating a second time  gives

 

         C(ρ) =  (1/2) (R/K) ρ2 + B1 ln (ρ ) + B2                . (63)

 

If we demand a physical solution where the concentration at the origin , C(ρ=0) is finite the constant B1 has to be zero; recall that ln (0) = - ∞ .  

So we take

                          

               C(ρ) =  (1/2) (R/K) ρ2  + B2 

 

 and clearly B2 = C(0) . So ,

 

C(ρ) =  (1/2) (R/K) ρ2  + C(0).                                     (64)

 

 

Problem 8. The interior of a sphere of radius a=0.10 m is initially at 0.0 celsius while the surface is kept at 100 celsius. The diffusion constant is

 

 

          K = 2(69/π2) cm2/sec = 2(69/π2) E-4 m2/sec .,

                                            = 1.40E-3  m2/sec.

 

Obtain an estimate of T(r,t). Take the approximation

 

       T(r,t) = exp(-λ t) F( r)  + 100 , with only one effective  lambda , λ ~ K/a2 = 0.14 /sec . The boundary value T(a,t) =100 implies that F(a)=0.

As will be shown below λ is numerically adjusted such that the numerical solution satisfies F(r=a) =0.

 

The heat diffusion equation is then

 

∂T/∂t = -λ exp(-λ t) F( r) = K exp( -λ t ) {d2 F(r)/dr2 +(2/r)(dF/dr)}   ,

which reduces to the radial equation

 

                d2 F(r)/dr2 +(2/r)(dF/dr) = -( λ/K)F( r )

 

With regards to the  initial conditions we already have that F(a)=0 but dF(a)/dr is unknown. On the other hand from the last equation we find that

 

(dF/dr)r=0 =0 , while T(r=0,t=0)= 0 = F(0) exp(0) +100 , implies

 

F(0)= -100.

 

Therefore one can start the numerical solution of the DE at the origin with IC  F(0)=-100. ,  dF(0)/dr =0.

 

 

 

 

 

A fourth degree  polynomial fit is

F ( r ) = -1.84e5 r4 -7.32e4 r3 + 1.96e4 r2 -41.2 r -99.9

 

and the approximation for the temperature is

 

T(r,t) = exp(-1.40 t ) F ( r )    + 100.

 

=exp(-1.40 t ) { -1.84e5 r4 -7.32e4 r3 + 1.96e4 r2 -41.2 r - 99.9 }

  +100.

 

 

, where the effective lambda is  1.40 /sec , t ~ seconds .

 

 

 

c solution problem 8.8 Simon Plots temp as function of (r,t)

      real lambda

      data lambda,a,nstep/1.40, .1, 10/

      Temp(r,t)=exp(-lambda*t)*(-1.84e5*r**4 -7.32e4*r**3 + 1.96e4*r**2

     $  -41.2*r - 99.9 ) +100.

      dt=.5

      dr=a/float(nstep)

      do 10 it=0,4

      t=dt*float(it)

      do 20 ir=0,nstep

      r=dr*float(ir)

      print 100,t,r,temp(r,t)

20    continue

      print*,'  '

10    continue

100   format(1x,'t,r, temp=',3(3x,e10.3))

      stop

      end

     

     

 

 

L,L/k=  1.3982321  999.999939

 

 r,f=    0.00000E+00   -0.10000E+03

 r,f=    0.16000E-02   -0.99955E+02

 r,f=    0.32000E-02   -0.99824E+02

 r,f=    0.48000E-02   -0.99608E+02

 r,f=    0.64000E-02   -0.99309E+02

 r,f=    0.80000E-02   -0.98925E+02

 r,f=    0.96000E-02   -0.98457E+02

 r,f=    0.11200E-01   -0.97906E+02

 r,f=    0.12800E-01   -0.97274E+02

 r,f=    0.14400E-01   -0.96560E+02

 r,f=    0.16000E-01   -0.95765E+02

 r,f=    0.17600E-01   -0.94892E+02

 r,f=    0.19200E-01   -0.93942E+02

 r,f=    0.20800E-01   -0.92915E+02

 r,f=    0.22400E-01   -0.91814E+02

 r,f=    0.24000E-01   -0.90640E+02

 r,f=    0.25600E-01   -0.89396E+02

 r,f=    0.27200E-01   -0.88082E+02

 r,f=    0.28800E-01   -0.86700E+02

 r,f=    0.30400E-01   -0.85254E+02

 r,f=    0.32000E-01   -0.83745E+02

 r,f=    0.33600E-01   -0.82175E+02

 r,f=    0.35200E-01   -0.80547E+02

 r,f=    0.36800E-01   -0.78864E+02

 r,f=    0.38400E-01   -0.77127E+02

 r,f=    0.40000E-01   -0.75339E+02

 r,f=    0.41600E-01   -0.73503E+02

 r,f=    0.43200E-01   -0.71623E+02

 r,f=    0.44800E-01   -0.69699E+02

 r,f=    0.46400E-01   -0.67736E+02

 r,f=    0.48000E-01   -0.65737E+02

 r,f=    0.49600E-01   -0.63703E+02

 r,f=    0.51200E-01   -0.61639E+02

 r,f=    0.52800E-01   -0.59546E+02

 r,f=    0.54400E-01   -0.57429E+02

 r,f=    0.56000E-01   -0.55290E+02

 r,f=    0.57600E-01   -0.53132E+02

 r,f=    0.59200E-01   -0.50959E+02

 r,f=    0.60800E-01   -0.48773E+02

 r,f=    0.62400E-01   -0.46577E+02

 r,f=    0.64000E-01   -0.44375E+02

 r,f=    0.65600E-01   -0.42169E+02

 r,f=    0.67200E-01   -0.39963E+02

 r,f=    0.68800E-01   -0.37760E+02

 r,f=    0.70400E-01   -0.35563E+02

 r,f=    0.72000E-01   -0.33374E+02

 r,f=    0.73600E-01   -0.31197E+02

 r,f=    0.75200E-01   -0.29034E+02

 r,f=    0.76800E-01   -0.26889E+02

 r,f=    0.78400E-01   -0.24764E+02

 r,f=    0.80000E-01   -0.22662E+02

 r,f=    0.81600E-01   -0.20585E+02

 r,f=    0.83200E-01   -0.18537E+02

 r,f=    0.84800E-01   -0.16519E+02

 r,f=    0.86400E-01   -0.14534E+02

 r,f=    0.88000E-01   -0.12586E+02

 r,f=    0.89600E-01   -0.10675E+02

 r,f=    0.91200E-01   -0.88045E+01

 r,f=    0.92800E-01   -0.69765E+01

 r,f=    0.94400E-01   -0.51928E+01

 r,f=    0.96000E-01   -0.34556E+01

 r,f=    0.97600E-01   -0.17666E+01

 r,f=    0.99200E-01   -0.12777E+00

 r,f=    0.10000E+00    0.67235E+00

 

c Simon problem 8.8  d^2F(r)/dr^2 +(2/r)(dF/dr) = -( L/K)*F( r )

c a factor is necessary to multiply -( L/K)*F( r )

      real L, k

      data a/.1 /

      data nstep, factor/1000, 10./

      g(r)=-(2./r)*(f1-f0)/dr -(L/k)*f1

      pi=2.*asin(1.)

      k=2.*(69./pi**2)*1.e-4

      L=factor*k/a**2

      dr=a/float(nstep)

      f0=-100.

      f1= f0

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

      kount=kp

      print*,'L,L/k=',L, L/k

      print*, '  '

      print 100 ,0.,f0

      do 10 i=2,nstep

      r=dr*float(i)

      f2=2.*f1-f0+dr**2*g(r-dr)

      if(i.eq.kount)then

c      print*,'f2,dr,g(r-dr)=',f2,dr,g(r-dr)

      print 100 ,r ,f2

      kount=kount+kp

      endif

      f0=f1

      f1=f2

10    continue

      print 100 ,r ,f2

100   format(1x,'r,f=',2(3x,e12.5))

      stop

      end

     

     

     

 

 

 

 

 

(February  4 , 2009)

 

UNDER CONSTRUCTION