Download Diffusion Equation: Derivation and Solutions in 1D and 2D and more Lecture notes Differential Equations in PDF only on Docsity!
Chapter 7
The Diffusion Equation
The 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 )
where u ( r , t ) is the density of the diffusing material at location r = ( x , y , z ) and time t. D ( u ( r , t ), r ) denotes the collective diffusion coefficient for density u at location r. If the diffusion coefficient doesn’t depend on the density, i.e., D is constant, then Eq. (7.1) reduces to the following linear equation:
∂ u ( r , t ) ∂ t
= D ∇^2 u ( r , t ). (7.2)
Equation (7.2) is also called the heat equation and also describes the distribution of a heat in a given region over time. Equation (7.2) can be derived in a straightforward way from the continuity equa- tion , which states that a change in density in any part of the system is due to inflow and outflow of material into and out of that part of the system. Effectively, no mate- rial is created or destroyed: ∂ u ∂ t
where Γ is the flux of the diffusing material. Equation (7.2) can be obtained easily from 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 is proportional to the local density gradient:
Γ = − D ∇ u ( r , t ).
59
7.1 The Diffusion Equation in 1D
Consider an IVP for the diffusion equation in one dimension:
∂ u ( x , t ) ∂ t
= D
∂ 2 u ( x , t ) ∂ x^2
on the interval x ∈ [ 0 , L ] with initial condition
u ( x , 0 ) = f ( x ), ∀ x ∈ [ 0 , L ] (7.4)
and Dirichlet boundary conditions
u ( 0 , t ) = u ( L , t ) = 0 ∀ t > 0. (7.5)
7.1.1 Analytical Solution
Let us attempt to find a nontrivial solution of (7.3) satisfying the boundary condi- tions (7.5) using separation of variables [4], i.e.,
u ( x , t ) = X ( x ) T ( t ).
Substituting u back into Eq. (7.3) one obtains:
1 D
T ′( t ) T ( t )
X ′′( x ) X ( x )
Since the right hand side depends only on x and the left hand side only on t , both sides 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 , (7.6) T ′( t ) + D λ T ( t ) = 0. (7.7)
Let us consider the first equation for X ( x ). Taking into account the boundary condi- tions (7.5) 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 ) = 0.
That is, the problem of finding of the solution of (7.3) reduces to the solving of linear ODE and consideration of three different cases with respect to the sign of λ :
- λ < 0: X ( x ) = C 1 e
√ − λ x (^) + C 2 e −
√ − λ x (^).
Fig. 7.1 Schematical repre- sentation of the FTCS finite difference scheme (7.9) for solving the 1-d diffusion equation (7.3).
u u u
Q t (^) j
Q
QkQ
^3
u t (^) j + 1
7.1.2 Numerical Treatment
The FTCS Explicit Method
Consider Eq. (7.3) with the initial condition (7.4). The first simple idea is an explicit forward in time, central in space (FTCS) method [28, 22] (see Fig. (7.1)):
u (^) ij +^1 − u (^) ij △ t
= D
u (^) ij + 1 − 2 u (^) ij + u (^) ij − 1 △ x^2
or, with α = D (^) △△ xt 2
u (^) ij +^1 = ( 1 − 2 α) u (^) ij + α ( u (^) ij + 1 + u (^) ij − 1 ). (7.9)
In order to check the stability of the schema (7.9) we apply again the ansatz (1.21) (see Sec. 1.3), considering a single Fourier mode in x space and obtain the following equation for the amplification factor g ( k ):
g^2 = ( 1 − 2 α) g + 2 g α cos( k △ x ) ,
from which
g ( k ) = 1 − 4 α sin^2
k △ x 2
The stability condition for the method (7.9) reads
| g ( k )| ≤ 1 ∀ k ⇔ α ≤
⇔ △ t ≤
△ x^2 D
Although the method (7.9) is conditionally stable, the derived stability condi- tion (7.10), however, hides an uncomfortable property: A doubling of the spatial resolution △ x requires a simultaneous reduction in the time-step △ t by a factor of four in order to maintain numerical stability. Certainly, the above constraint limits us to absurdly small time-steps in high resolution calculations. The next point to emphasize is the numerical dispersion. Indeed, let us compare the exact dispersion relation for Eq. (7.3) and relation, obtained by means of the schema (7.9). If we consider the perturbations in form exp( ikx − i ω t ) the dispersion
Fig. 7.2 Dispersion relation by means of the schema (7.9) for different values of α, compared with the exact dis- persion relation for Eq. (7.3).
0 0.2 0.4 0.6 0.8 1 0
5
10
15
20
k∆ x/π
iω∆
t/
α
exact α=0. α=0. α=0. α=0.
relation for Eq. (7.3) reads i ω = D k^2.
On the other hand, the FTCS schema (7.3) leads to the following relation
ei^ ω△ t^ = 1 − 4 α sin^2
k △ x 2
or, in other words
i ω△ t = − ln
1 − 4 α sin^2
k △ x 2
The comparison between exact and numerical disperion relations is shown on Fig. (7.2). One can see, that both relations are in good agreement only for k △ x ≪ 1. For α > 0 .25 the method is stable, but the values of ω can be complex, i.e., the Fourier modes drops off, performing damped oscillations (see Fig. (7.2) for α = 0. 3 and α = 0 .4). Now, if we try to make the time step smaler, in the limit △ t → 0 (or α → 0) we obtain
i ω△ t ≈ 4 α sin^2
k △ x 2
= k^2 D △ t
sin^2
k △ x 2
k △ x 2
i.e., we get the correct dispersion relation only if the space step △ x is small enougth too.
Fig. 7.3 Schematical rep- resentation of the DuFort- Frankel method (7.12). q t (^) j − 1 xi − 1
u
u u
xi
q xi + 1
Q t (^) j
Q
QkQ
^3
u t (^) j + 1
Fig. 7.4 Schematical rep- resentation of the implicit BTCS method (7.13).
u e u
u
Q t (^) j
Q
Q
Qk
t (^) j + 1
method. Indeed, consider the Taylor expansion for Eq. (7.3) by means of (7.12):
u (^) ij +^1 − u (^) ij −^1 2 △ t
= D
u (^) ij + 1 − u (^) ij +^1 − u (^) ij −^1 + u (^) ij − 1 △ x^2
ut −
△ x^2 3!
uttt +... =
D
△ x^2
△ x^2 uxx + 2 △ x^4 4!
uxxxx − △ t^2 utt − 2 △ t^4 4!
utttt +...
ut + O(△ t^2 ) = Duxx + O(△ x^2 ) − D
△ t^2 △ x^2
utt + O
△ t^4 △ x^2
In other words, the method (7.12) has order of accuracy
O
△ t^2 , △ x^2 ,
△ t^2 △ x^2
For cosistency, △ t /△ x → 0 as △ t → 0 and △ x → 0, so (7.12) is inconsistent. This constitutes an effective restriction on △ t. For large △ t , however, the scheme (7.12) is consistent with another equation of the form
D utt + ut = D uxx.
7.1.2.1 The BTCS Implicit Method
One can try to overcome problems, described above by introducing an implicit method. The simplest example is a BTCS (backward in time, central in space) method (see Fig. 7.4) [29]. The differential schema reads:
u (^) ij +^1 − u (^) ij △ t
= D
u (^) ij ++ 11 − 2 u (^) ij +^1 + u (^) ij −+ 11 △ x^2
Fig. 7.5 Schematical rep- resentation of the Crank- Nicolson method (7.14). q t (^) j − 1 xi − 1
u u u xi
q xi + 1
t (^) j + (^12)
u u u t (^) j + 1
or, with α = D △ t /△ x^2
− u (^) ij = α u (^) ij ++ 11 − ( 1 + 2 α) u (^) ij +^1 + α u (^) ij −+ 11. (7.13)
In this case the amplification factor g ( k ) is given by
g ( k ) =
1 + 4 α^ sin^2
k △ x^2 2
That is, the schema (7.13) is unconditionally stable. However, the method has order of accuracy O(△ t , △ x^2 ), i.e., first order in time, and second in space. Is it possible to improve it? The answer to is given below.
The Crank-Nicolson Method
An implicit scheme, introduced by J. Crank and P. Nicolson in 1947 [6] is based on the central approximation of Eq. (7.3) at the point ( xi , t (^) j + 12 △ t ) (see Fig. (7.5)):
u (^) ij +^1 − u (^) ij 2 △ 2 t
= D
u j + (^12) i + 1 −^2 u^
j + (^12) i +^ u^
j + (^12) i − 1 △ x^2
The approximation used for the space derivative is just an average of approxima- tions in points ( xi , t (^) j ) and ( xi , t (^) j + 1 ):
un i +^1 − uni △ t
= D
( un i ++ 11 − 2 un i +^1 + un i −+ 11 ) + ( uni + 1 − 2 uni + uni − 1 ) 2 △ x^2
Introducing α = D △ t /△ x^2 one can rewrite the last equation as
−^ α u^ ij ++ 11 + 2 ( 1 +^ α) u^ ij +^1 −^ α u^ ij −+ 11 =^ α unj + 1 +^2 (^1 −^ α) u^ ij +^ α unj − 1.^ (7.14)
All terms on the right-hand side of Eq. (7.14) are known. Hence, the equations in (7.14) form a tridiagonal linear system
A u = b.
Fig. 7.7 Contour plot of the diffusion of the initial Gauss pulse, calculated with the BTCS implicit method (7.13). x
t
−5^0 0
50
100
150
200
Example 2
Use the implicit BTCS method (7.13) to solve the one-dimensional diffusion equa- tion ut = uxx ,
on the interval x ∈ [− L , L ], if the initial distribution is a Gauss pulse of the form u ( x , 0 ) = exp(− x^2 ) and the density on both ends of the interval is given as ux (− L , t ) = ux ( L , t ) = 0. For the other parameters see the table below:
Space interval L = 5 Space discretization step △ x = 0. 1 Time discretization step △ t = 0. 05 Amount of time steps T = 200
Solution of the problem is shown on Fig. (7.7).
Example 3
Use the Crank-Nicolson method (7.14) to solve the one-dimensional heat equation
ut = 1. 44 uxx ,
on the interval x ∈ [ 0 , L ], if the initial heat distribution is u ( x , 0 ) = f ( x ) and again, the temperature on both ends of the interval is given as u ( 0 , t ) = Tl , u ( L , t ) = Tr. Other parameters are choosen as:
Space interval L = Space discretization step △ x = 0. 1 Time discretization step △ t = 0. 05 Amount of time steps T = 15 Boundary conditions Tl = 2 , Tr = 0. 5 Initial heat distribution f ( x ) = 2 − 1. 5 x + sin( π x )
Fig. 7.8 Contour plot of the heat distribution, calculated with the Crank-Nicolson method (7.14). x
t
(^00) 0.5 1
5
10
15
1
2
Numerical solution of the problem in question is shown on Fig. (7.8).
7.2 The Diffusion Equation in 2D
Let us consider the solution of the diffusion equation (7.2) in two dimensions
∂ u ∂ t
= D
∂ 2 u ∂ x^2
∂ 2 u ∂ y^2
where u = u ( x , y , t ), x ∈ [ ax , bx ], y ∈ [ ay , by ]. Suppose, that the inittial condition is given and function u satisfies boundary conditions in both x - and in y -directions. As before, we discretize in time on the uniform grid tn = t 0 + n △ t , n = 0 , 1 , 2 ,.. .. Futhermore, in the both x - and y -directions, we use the uniform grid
xi = x 0 + i △ x , i = 0 , · · · , M , △ x =
bx − ax M + 1
y (^) j = y 0 + j △ y , j = 0 , · · · , N , △ y =
by − ay N + 1
7.2.1 Numerical Treatment
The FTCS Method in 2D
In the case of two dimensions the explicit FTCS scheme reads
un i j +^1 − uni j △ t
= D
uni + 1 j − 2 uni j + uni − 1 j △ x^2
uni j + 1 − 2 uni j + uni j − 1 △ y^2
− α u^221 + γ u^211 − β u^212 = u^111 + α u^201 + β u^210 , − α u^222 + γ u^212 − β ( u^213 + u^211 ) = u^112 + α u^202 , − α u^223 + γ u^213 − β u^212 = u^113 + α u^203 + β u^214 , − α( u^231 + u^211 ) + γ u^221 − β u^222 = u^121 + β u^220 , − α( u^232 + u^212 ) + γ u^222 − β ( u^223 + u^221 ) = u^122 , − α u^221 + γ u^231 − β u^232 = u^131 + α u^241 + β u^230 , − α u^222 + γ u^232 − β ( u^233 + u^231 ) = u^132 + α u^142 , − α u^223 + γ u^233 − β u^232 = u^133 + α u^244 + β u^234.
Formally, one can rewrite the system above to the matrix form A u = b , i.e., γ − β 0 − α 0 0 0 0 0 − β γ − β 0 − α 0 0 0 0 0 − β γ 0 0 − α 0 0 0 − α 0 0 γ − β 0 − α 0 0 0 − α 0 − β γ − β 0 − α 0 0 0 − α 0 − β γ 0 0 − α 0 0 0 − α 0 0 γ − β 0 0 0 0 0 − α 0 − β γ − β 0 0 0 0 0 − α 0 − β γ
u^211 u^212 u^213 u^221 u^222 u^223 u^231 u^232 u^233
u^111 + α u^201 + β u^210 u^112 + α u^202 u^113 + α u^203 + β u^214 u^121 + β u^220 u^122 u^123 + β u^224 u^131 + α u^241 + β u^230 u^132 + α u^242 u^133 + α u^144 + β u^234
The matrix A is a five-band matrix. Nevertheless, despite of the fact that the schema is absolute stable, two of five bands are desposed so far apart from the main diagonal, that simple O( n ) algorithms like TDMA are difficult or even impossible to apply.
The ADI Method
The idea of the ADI-method ( alternating direction implicit ) is to alternate direction and thus solve two one-dimensional problem at each time step [17]. The first step keeps y -direction fixed:
un i j +^1 /^2 − uni j △ t / 2
= D
( (^) un + 1 / 2 i + 1 j −^2 u
n + 1 / 2 i j +^ u
n + 1 / 2 i − 1 j △ x^2
uni j + 1 − 2 uni j + uni j − 1 △ y^2
In the second step we keep x -direction fixed:
un i j +^1 − un i j +^1 /^2 △ t / 2
= D
u n + 1 / 2 i + 1 j −^2 u
n + 1 / 2 i j +^ u
n + 1 / 2 i − 1 j △ x^2
un i j ++^11 − 2 un i j +^1 + un i j +−^11 △ y^2
Both equations can be written in a triadiagonal form. Define
α = D △ t 2 △ x^2
, β = D △ t 2 △ y^2
Than we get:
− α un i ++ 11 j /^2 + ( 1 + 2 α) un i j +^1 /^2 − α un i −+ 11 / j^2 = β uni j + 1 + ( 1 − 2 β ) uni j + β uni j −(7.19) 1
− β un i j ++^11 + ( 1 + 2 β ) un i j +^1 − β un i j +−^11 = α un i ++ 11 j /^2 + ( 1 − 2 α) un i j +^1 /^2 + α un i ++ 11 / j^2.
Instead of five-band matrix in BTCS method (7.18), here each time step can be obtained in two sweeps. Each sweep can be done by solving a tridiagonal system of equations. The ADI-method is second order in time and space and is absolute stable [11] (however, the ADI method in 3D is conditional stable only).
7.2.2 Examples
Use the ADI method (7.19) to solve the two-dimensional diffusion equation
∂ t u ( r , t ) = △ u ( r , t ) ,
where u = u ( r , t ), r ⊆ R^2 on the interval r ∈ [ 0 , L ] × [ 0 , L ], if the initial distribution is a Gauss pulse of the form u ( x , 0 ) = exp(− 20 ( x − L / 2 )^2 − 20 ( y − L / 2 )^2 ) and the density on both ends of the interval is given as ur ( 0 , t ) = ur ( L , t ) = 0. Other parameters are choosen according to the table below.
Space interval L = 1 Amount of points M = 100, (△ x = △ y ) Time discretization step △ t = 0. 001 Amount of time steps T = 40
Solution of the problem is shown on Fig. (7.9).
