Numerical Solution of the Kinetic Model

The numerical approximation for the solution of Eqs. (119–124) is calculated by a

finite difference scheme. After replacing the spatial derivations with difference

quotients, a system of ordinary differential equations for the concentration C at

discrete points is obtained.

The origin of the coordinate system at the chip center is located and the onedimensional

wood chip is divided into 2n slices with the width Dh = s/2n. Ci

denotes the concentration at height 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

_ _126_

The difference quotient is applied consecutively in Eq. (119), with h= Dh/2 obtaining

the following difference equations

_ C

i _ t _ ≈ D

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

To simplify expressions, it was assumed that D does not depend on the spatial

direction; the general case, however, can be solved using the same principle.

After approximating C2, C1, C0 with a quadratic polynomial and minding

[Eq. (121)], we obtain

C 0_ t _ ≈ 4

C 1_ t _ _

C 2_ t _ _128_

The same approximation for Cn, Cn – 1, Cn – 2 results in

∂ C _ s _2_ t _ ∂ z ≈ 1

2D h _3 Cn _ t _ _ 4 Cn _1_ t _ Cn _2_ t __ which, after combining with

Eqs. (122) and (124), yields

Cn _ t _ ≈ CBulk _ t _ _

D

k 2D h _3 Cn _ t _ _ 4 Cn _1_ t _ Cn _2_ t __ _129_

and

_ C

Bulk _ t _ ≈

VChip D

s VBulk D h _3 Cn _ t _ _ 4 Cn _1_ t _ Cn _2_ t __ _130_

Equations (127–130) define a system of differential algebraic equations (DAEs).

After elimination of C0(t) and Cn(t) by inserting Eq. (128) and Eq. (129) into Eqs.

228 4 Chemical Pulping Processes

(127) and (130), the DAEs simplify to a system of ordinary differential equations

(ODE) which can be solved by any standard numerical ODE solver that has good

stability properties, for example, an implicit Runge Kutta method. Euler’s – which

has excellent stability properties – is used in the sample code, and although a set

of linear equations must be solved for every time step, the method is very fast

because the system matrix is almost trigonal.

4.2.6