Lecture 8. Diffusion
by Reinaldo Baretti Machín,
and Alfonso Baretti Huertas
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 :
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.
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