Numerical solution for diffusion equation
Consider the classical diffusion equation:
∂c ∂ 2c =D 2 ∂t ∂x
(1)
Let us normalise this equation: that is, let us divide the time parameter by a characteristic time T so that the normalised quantity t/T ranges from zero to unity; similarly, let us divide the distances by the characteristic distance L so that the normalised quantity x/L ranges from zero to unity . Let the normalised time be denoted by τ and the normalised distance by χ. Let us denote the normalised diffusivity by α = DT/L 2 . The following is the normalised diffusion equation: ∂c ∂2c =α 2 ∂τ ∂χ
(2)
For simplicity’s sake (without loss of generality) assume α = 1.0. Note that the Fourier law of heat conduction also leads to an identical equation (when the normalisation is carried out). The numerical solution of the diffusion equation above can be obtained using a finite difference schme, in which, the time and positions are discretized and the solution for each discrete position at different time steps are obtained by solving the corresponding difference equations (which are a set of algebraic equations as we show below). Let us consider the term on the left-hand side, namely, cretized using Forward-Euler scheme as follows:
∂c ∂τ
. It can be d is-
∂c cτ +δτ − cτ = ∂τ δτ Let us consider the term on the right-hand side, namely, discretized using a central difference formula as follows: ∂2c 1 = [cχ−δχ − 2cχ + cχ+δχ ] 2 ∂χ (δχ)2 1
(3) ∂ 2c ∂χ 2
. It can b e
(4)
Note that for partial derivative of time, we need to keep the position fixed, while, for that of position, we need to keep the tim e fixed. While fixing the position for partial derivative of time is rather straight-forward, for the partial derivative of position, the time can be fixed in three different ways: at τ (which is known as the explicit scheme), at τ + δτ (which is known as implicit scheme), and some at τ and some at τ + δτ (which is known as semi-implicit scheme). The explicit expressions for these three cases are as shown below:
Explicit scheme
cτχ+δτ − cτχ 1 τ = cχ−δχ − 2cτχ + cτχ+δχ 2 δτ (δχ) Implicit scheme cτχ+δτ − cτχ 1 τ +δτ τ +δτ τ +δτ = c − 2cχ + c χ − δχ χ + δχ δτ (δχ)2
(5)
(6)
Semi-implicit scheme
cτχ+δτ − cτχ 1 τ τ +δτ τ = c − 2c + c χ χ+δχ (δχ)2 χ−δχ δτ
(7)
Note that depending on the number of positional nodes, the number of equations in the above equations vary. From the expressions above, it is clear that knowing the composition profile at any time τ , one can obtain the same at the future time τ + δτ ; however, while in the case of implicit scheme, the composition profile at future time is obtained by solving a matrix equation, for explicit and semi-implicit schemes, the composition profile is obtained by solving an algebraic equation for each position independently. That is, to solve the above set of equations, the initial profile is essential. Further, note that at the two extremes, namely, at χ = 0 and χ = 1, the spatial discretization involves unknowns either to the left or to the right, which are not part of the calculation domain. These values are to be obtained
2
using the two boundary conditions at these two points. That is, to sovle the above set of equations, two boundary conditions are essential.
Problem 1 Consider the following boundary conditions: [c]χ=0(τ ) = 1.0; and,
∂c ∂χ
χ=1
(τ ) = 0.
The initial conditions is [ c]χ=0(τ = 0) = 1.0 and c(τ = 0) at all the other points is zero. (a) Write a script to solve the given partial differential equation with the given boundary conditions, using explicit method. The calculations should be done for a time of about 1.0 time unit and the solutions should be plotted for every 0.1 time units (on the same figure). Use a δχ value of 0.01. Check your implementation using δτ values of 1 × 10−6 , 1 × 10−4 and 1 × 10−2. For which value of δτ are your calculations stable? (b) Write a script to solve the given partial differential equation with the given boundary conditions, using implicit method. The calculations should be done for a time of about 1.0 time unit and the solutions should be plotted for every 0.1 time units (on the same figure). Use a δχ value of 0.01. Check your implementation using δτ values of 1 × 10−6 , 1 × 10−4 and 1 × 10−2. For which value of δτ are your calculations stable? (c) Write a script to solve the given partial differential equation with the given boundary conditions, using semi-implicit method. The calculations should be done timetime of about the solutions should be plotted for evfor erya 0.1 uni ts1.0 (ontime the unit sameand figure). Use a δχ value of −6 0.01. Check your implementation using δτ values of 1 × 10 , 1 × 10−4 and 1 × 10−2 . For which value of δτ are your calculations stable?
Problem 2 (Optional) You will be given bonus points if you solve this problem Consider the following boundary conditions: 3
[c]χ=0(τ ) = 1.0; and, [c]χ=1(τ ) = 0.0. The initial conditions is [ c]χ=0(τ = 0) = 1.0 and c(τ = 0) at all the other points is zero. (a) Write a script to solve the given partial differential equation with the given boundary conditions, using explicit method. The calculations should be done for a time of about 1.0 time unit and the solutions should be plotted for every 0.1 time units (on the same figure). Use a δχ value of 0.01. Check your implementation using δτ values of 1 × 10−6 , 1 × 10−4 and 1 × 10−2. For which value of δτ are your calculations stable? (b) Write a script to solve the given partial differential equation with the given boundary conditions, using implicit method. The calculations should be done for a time of about 1.0 time unit and the solutions should be plotted for every 0.1 time units (on the same figure). Use a δχ value of 0.01. Check −
−
−
your implementation using δτ values of 1 × 10 6 , 1 × 10 4 and 1 × 10 2. For which value of δτ are your calculations stable? (c) Write a script to solve the given partial differential equation with the given boundary conditions, using semi-implicit method. The calculations should be done for a time of about 1.0 time unit and the solutions should be plotted for ev ery 0.1 time uni ts (on the same figure). Use a δχ value of 0.01. Check your implementation using δτ values of 1 × 10−6 , 1 × 10−4 and 1 × 10−2 . For which value of δτ are your calculations stable?
Problem 3 (Home assignment) Consider the following boundary conditions: [c]χ=0(τ ) = 1.0; and,
∂c ∂χ
χ=1
(τ ) = 0.
The initial conditions is c(τ = 0) = 1.0 for all positions. (a) Write a script to solve the given partial differential equation with the given boundary conditions, using explicit method. The calculations should be done for a time of about 1.0 time unit and the solutions should be plotted for every 0.1 time units (on the same figure). Use a δχ value of 0.01. Check your implementation using δτ values of 1 × 10−6 , 1 × 10−4 and 1 × 10−2. For
4
which value of δτ are your calculations stable? (b) Write a script to solve the given partial differential equation with the given boundary conditions, using implicit method. The calculations should be done for a time of about 1.0 time unit and the solutions should be plotted for every 0.1 time units (on the same figure). Use a δχ value of 0.01. Check your implementation using δτ values of 1 × 10−6 , 1 × 10−4 and 1 × 10−2. For which value of δτ are your calculations stable? (c) Write a script to solve the given partial differential equation with the given boundary conditions, using semi-implicit method. The calculations should be done for a time of about 1.0 time unit and the solutions should be plotted for ev ery 0.1 time uni ts (on the same figure). Use a δχ value of 0.01. Check your implementation using δτ values of 1 × 10−6 , 1 × 10−4 and 1 × 10−2 . For which value of δτ are your calculations stable?
5
