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:
Reference :
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
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
(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
, (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