Chapter 7 The Diffusion Equation - uni-muenster.de

1 Chapter 7 The Diffusion EquationThe Diffusion Equation is a partial differential Equation which describes density fluc-tuations in a material undergoing Diffusion . The Equation can be written as: u(r,t) t= D(u(r,t),r) u(r,t) ,( )whereu(r,t)is the density of the diffusing material at locationr= (x,y,z)and (u(r,t),r)denotes the collectivediffusion coefficientfor densityuat the Diffusion coefficient doesn t depend on the density, ,Dis constant, thenEq. ( ) reduces to the following linear Equation : u(r,t) t=D 2u(r,t).( ) Equation ( ) is also calledthe heat equationand also describes the distribution ofa heat in a given region over ( ) can be derived in a straightforward way from thecontinuity equa -tion, which states that a change in density in any part of the system is due to inflowand outflow of material into and out of that part of the , no mate-rial is created or destroyed: u t+ =0,where is the flux of the diffusing material.

2 Equation ( ) can be obtained easilyfrom the last Equation when combined with the phenomenological Fick s first law,which assumes that the flux of the diffusing material in any part of the system isproportional to the local density gradient: = D u(r,t). The Diffusion Equation in 1 DConsider an IVP for the Diffusion Equation in one dimension: u(x,t) t=D 2u(x,t) x2( )on the intervalx [0,L]with initial conditionu(x,0) =f(x), x [0,L]( )and Dirichlet boundary conditionsu(0,t) =u(L,t) =0 t>0.( ) Analytical SolutionLet us attempt to find a nontrivial solution of ( ) satisfying the boundary condi-tions ( ) using separation of variables [4], ,u(x,t) =X(x)T(t).Substitutinguback into Eq.

3 ( ) one obtains:1DT (t)T(t)=X (x)X(x).Since the right hand side depends only onxand the left hand side only ont, bothsides are equal to some constant value ( sign is taken for convenience reasons).Hence, one can rewrite the last Equation as a system of two ODE s:X (x) + X(x) =0,( )T (t) +D T(t) =0.( )Let us consider the first Equation forX(x). Taking into account the boundary condi-tions ( ) one obtains (T(t)6=0 as we are loocking for nontrivial solutions)u(0,t) =X(0)T(t) =0 X(0) =0,u(L,t) =X(L)T(t) =0 X(L) = is, the problem of finding of the solution of ( ) reduces to the solving oflinear ODE and consideration of three different cases with respect to the sign of :1.

4 <0:X(x) =C1e x+C2e into account the boundary conditions one getsC1=C2=0, so for <0only the trivial solution =0:X(x) =C1x+C2 Again, due to the boundary conditions, one gets only trivialsolution of the prob-lem (C1=C2=0).3. >0:X(x) =C1cos( x) +C2sin( x).Substituting of the boundary conditions leads to the following equations for theconstantsC1andC2:X(0) =C1=0,X(L) =C2sin( L) =0 sin( L) =0 n= nL 2,n=1,2, ..Hence,X(t) =Cnsin nLx .That is, the second Equation for the functionT(t)takes the form:T (t) +D nL T(t) =0 T(t) =Bnexp D nL 2t ,whereBnis , the general solution of the problem ( ) can bewritten asu(x,t) = n=1 Ansin nLx exp D nL 2t ,An= order to findAnone can use the initial condition ( ).

5 Indeed, if we write thefunctionf(x)as a Fourier series, we obtain:f(x) = n=1 Fnsin nLx = n=1 Ansin nLx ,An=Fn=2 LZL0f( )sin nL d .Hence, the genetal solution of Eq. ( ) reads:u(x,t) = n=1 2 LZL0f( )sin nL d sin nLx exp D nL 2t .( )Fig. repre-sentation of the FTCS finitedifference scheme ( ) forsolving the 1-d diffusionequation ( ).uuu6tjQQQQk 3tj+ Numerical TreatmentThe FTCS Explicit MethodConsider Eq. ( ) with the initial condition ( ). The first simple idea is an explicitforward in time, central in space (FTCS) method [28, 22] (seeFig. ( )):uj+1i uji t=Duji+1 2uji+uji 1 x2,or, with =D t x2uj+1i= (1 2 )uji+ (uji+1+uji 1).( )In order to check the stability of the schema ( ) we apply again the ansatz ( )(see Sec.)

6 , considering a single Fourier mode inxspace and obtain the followingequation for the amplification factorg(k):g2= (1 2 )g+2g cos(k x),from whichg(k) =1 4 sin2k stability condition for the method ( ) reads|g(k)| 1 k 12 t 12 x2D.( )Although the method ( ) is conditionally stable, the derived stability condi-tion ( ), however, hides an uncomfortable property: A doubling of the spatialresolution xrequires a simultaneous reduction in the time-step tby a factor offour in order to maintain numerical stability. Certainly, the above constraint limitsus to absurdly small time-steps in high resolution next point to emphasize is the numerical dispersion. Indeed, let us comparethe exact dispersion relation for Eq.

7 ( ) and relation, obtained by means of theschema ( ). If we consider the perturbations in form exp(ikx i t)the dispersionFig. relationby means of the schema ( )for different values of ,compared with the exact dis-persion relation for Eq. ( ). x/ i t/ exact = = = = for Eq. ( ) readsi =D the other hand, the FTCS schema ( ) leads to the following relationei t=1 4 sin2 k x2 ,or, in other wordsi t= ln 1 4 sin2 k x2 .The comparison between exact and numerical disperion relations is shown onFig. ( ). One can see, that both relations are in good agreement only fork x > the method is stable, but the values of can be complex, , theFourier modes drops off, performing damped oscillations (see Fig.)

8 ( ) for = = ). Now, if we try to make the time step smaler, in the limit t 0 (or 0) we obtaini t 4 sin2 k x2 =k2D tsin2 k x2 k x2 2, , we get the correct dispersion relation only if the space step xis small Richardson MethodThe first idea to improve the approximation order of the schema is to use the centraldiferences for the time derivative of Eq. ( ), namely [28]uj+1i uj 1i2 t=Duji+1 2uji+uji 1 x2,or, with =D t/ x2uj+1i=uj 1i+2 uji+1 2uji+uji 1 .( )Unfortunately, one can show that the schema ( ) is unconditional unstable. In-deed, amplification factorg(k)in this case fulfilles the following Equation :g2+2 g 1=0, =4 sin2k x2,givingg1,2= p 2+ |g2(k)|>1 for all values ofk, the schema ( ) is absolut DuFort-Frankel MethodLet us consider one of many alternative algorithms which have been designed toovercome the stability problems of the simple FTCS and Richardson methods.

9 Wemodify Eq. ( ) as (see Fig. ( )) [28]uj+1i uj 1i2 t=Duji+1 2uj+1i+uj 1i2+uji 1 x2,which can be solved explicitly foruj+1i:uj+1i=1 1+ uj 1i+ 1+ uji+1+uji 1 ,( )where =2D t/ x. When the usual von Neumann stability analysis is appliedto the method ( ), the amplification factorg(k)can be found from(1+ )g2 2g cos(k x) + ( 1) = can be easily shown, that stability condition is fulfilledfor all values of , so themethod ( ) is unconditionally stable. However, this does not imply that xand tcan be made indefinitely large; we must still worry about the accuracy of theFig. rep-resentation of the DuFort-Frankel method ( ).tj 1qxi 1uuu6xiqxi+1tjQQQQk 3tj+1uFig. rep-resentation of the implicitBTCS method ( ).

10 Euuu6tjQQQQk 3tj+1method. Indeed, consider the Taylor expansion for Eq. ( )by means of ( ):uj+1i uj 1i2 t=Duji+1 uj+1i uj 1i+uji 1 x2 ut x23!uttt+..=D x2 x2uxx+2 x44!uxxxx t2utt 2 t44!utttt+.. ut+O( t2) =Duxx+O( x2) D t2 x2 utt+O t4 x2 .In other words, the method ( ) has order of accuracyO t2, x2, t2 x2 .For cosistency, t/ x 0 as t 0 and x 0, so ( ) is inconsistent. Thisconstitutes an effective restriction on t. For large t, however, the scheme ( )is consistent withanotherequation of the formD utt+ut=D The BTCS Implicit MethodOne can try to overcome problems, described above by introducing an implicitmethod. The simplest example is a BTCS (backward in time, central in space)method (see Fig.)