Numerical Solution of the Diffusion Model

The solution of Eqs. (68) and (69) is calculated numerically by a finite difference

scheme. The origin of the coordinate system is located at the cylinder center. For

the radial concentration profile [Eq. (68)], the interval [–a, a] is divided into 2n

pieces of equal size Dh= a/n, for the axial concentration profile [Eq. (69)]the interval

[–Z/2, Z/2]is divided into 2n pieces of size Dh = Z/2n. Ci denotes the concentration

at iDh, thus Ci(t) = C(iDh,t).

The derivation of a smooth function can be approximated by a central difference

quotient:

df

dx _ x _≈ f _ x h _ _ f _ x _ h _

2 h

_ _70_

4.2 Kraft Pulping Processes 163

To obtain an approximation for a second-order derivation, the second-order derivation

is replaced by a central difference quotient of first-order derivations, after

which the first-order derivations are replaced by central difference quotients.

The resulting difference equations for Eq. (68) are

_ C

i _ t _≈ Dr

D h 2 __1

2 i _ Ci 1_ t _ _ 2 Ci _ t _ Ci _1_ t ___1 _

2 i __ i _ 1_ _____ n _ 1 _71_

and

_ C

i _ t _≈ Dz

D h 2 _ Ci 1_ t _ _ 2 Ci _ t _ Ci _1_ t __ i _ 1_ _____ n _ 1 _72_

for Eq. (67).

The condition that C is finite at r = 0 in Eq. (68) implies ∂ C

∂ r _0_ t _ _ 0 for t > 0

and symmetry of problem Eq. (69) implies ∂ C

∂ z _0_ t _ _ 0 for t > 0. After approximation

of C2, C1,C0 with a quadratic polynomial this conditions transform into:

C 0_ t _≈ 4

C 1_ t _ _

C 2_ t _ _73_

After inserting Eq. (73) into Eqs. (71) and (72) respectively, a system of ordinary

differential equations (ODE) is obtained which can be solved by any standard

numerical ODE solver with good stability properties.

Euler’s implicit method is used in the sample code. Only a set of linear equations

with a tridiagonal system matrix is solved each time step.

4.2.4