Moving finite elements method applied to the solution of front reaction models: causticizing reaction

Pergamon

Computers them. Engng Vol. 19, Suppl., pp. S421S426.1995 Copyright 0 1995 Elsevier Science Ltd Printed in &eat Britain. All rights reserved 0098-1354(95)00052-6 0098-1354/95 $9.50 t 0.00

MOVING FINITE ELEMENTS METHOD APPLIED TO THE SOLUTION OF FRONT REACTION MODELS: CAUSTICIZING REACTION B.P.M. Duarte and A.A.T.G.

Portugal

Departamento de Engenharia Quimica Universidade de Coimbra 3000 Coimbra - Portugal

ABSTRACT In this paper a Moving Finite Elements Method based on cubic Hermite polynomials is presented. PDE’s solution is obtained by minimizing the square norm of discretized residuals all over the domain in order to the time derivatives of the solutions on the nodes and nodal velocities. Spatial discretization gives origin to a system of ODES, that can degenerate when singularities occur. These can be of two kinds: parallelism and nodal coalescence. Nodal coalescence is avoided by fixing the minimum internodal distance while parallelism is dealt with the addition of penalty functions to the objective function resulting from the application of Gale&n Method. Boundary conditions are previously time derived originating constraints to the minimization system arising from the basic formulation. The solution of this problem is obtained via the Lagrange Multipliers Method. In this paper this methodology is applied to the solution of a front reaction model - the causticizing reaction. The results obtained are in close agreement with those published in the literature where orthogonal collocation in Finite Elements was used. The adequacy of the Moving Finite Elements Method for the solution of front reaction models is demonstrated in this paper. Keywords:

Moving Finite Elements, Front reaction models, PDEs, Penalty functions

INTRODUCTION Front reaction models occur in chemical processes and have an undeniable practical interest. They appear in systems characterized by an initial spatial domain divided in two or more subdomains. Each of these subdomains have different dynamic equations or physical properties, being generally called phases. The point of contact between the phases, the interface, is time dependent. Its movement results from the direct application of conservation principles. There, the solution is discontinuous and the calculated position is obtained in accordance with the solution profiles in the subdomains. Murray (1959) solved a front reaction model that describes an ice block fusion subjected to a source of heat in a fixed boundary. Unequal conductivities in each phase were considered. Powers (1975) derived and solved a front reaction model for the reaction occuring in an hydride fixed bed. The numerical approach he used is based on Finite Differences and interpolation is used to locate the interface inside a constant subdomain scaling technique. Pais et0al. (1985) solved Powers’s model using a Fixed Finite Elements scheme. The position of the interface was defined by the solution of an additional equation resulting from the front spatial velocity discretization. The length of the subdomains was constantly updat.ed and the intcrfacc’s location was simultaneously computed. A very similar technique was applied by Portugal ct al. (1989) to the solution S421

European Symposium on Computer Aided Process Engineering-5

s422

of a model describing the causticizing reaction. The results obtained were validated by experimental data. The evaluation of the position of the interface was a central point to the work, implemented through an interpolation procedure included in the Fixed Finite Elements Method. FORMULATION

OF THE METHOD

Moving Finite Elements Method (MFEM) was originally presented by K. Miller and R.N. Miller (1981) with the purpose of solving PDEs that develope typical steep moving fronts. Initially K. Miller derived an approach based on piecewise linear polynomials, evaluating later the use of higher order spaces of approximations. As conclusions he stated that higher order approximations were particularly indicated for higher order PDEs. The methodology is based on the minimization of square norm of residuals of the discrete equation all over the domain in order to time derivatives of the solution on the nodes and node velocities. The nodes are allowed to move in order to decrease spatial gradients. Since these structures occur in connection with fronts (moving or not) the nodes tend to migrate to its vicinity. It is the mobility of the nodes the main characteristic of the Moving Finite Elements Methods. The discretization originates an ODES system solved by an implicit integrator. Sometimes the system becomes singular due to nodal coalescence or parallelism. Parallelism occurs when the solution in two neighbour elements is a straight line originating linear dependent equations. Nodal coalescence occurs due to shocks between neighbour nodes. While nodal coalescence is overcome by fixing the minimum internodal distance, parallelism is dealt by the addition of a penalty function. Later, Wathen and Baines (1988) derived a global formulation for MFEM and presented the techniques adopted to overcome singularities. Jimack (1991) derived a new approach based on Wathen and Baines’s method using cubic Hermite and Lagrange polynomials. The results obtained were promising. Following K. Miller’s methodology K. Pipilis (1991) implemented a general approach based on cubic Hermite polynomials whose results were very robust and stable when tested with a wide range of problems. Let us consider the MFE method applied to the general equations a = &(u, u,, grr, Z, t) subject to FL&u,, z, t) = 0, ER(g,v,, z, t) = 0 and ~(0, z) = go(z). Approximating the dependent variable m in a finite element k by cubic Hermite polynomials, we may write

up = 2 HI(z

(1)

I=1

Being Rk the scaled residual in the finite element k, u the scaled spatial variable, SI the spatial position of node I, aIJ the 1a coefficient of Hermit,e polynomials in element J and N the number of elements, the PDE solution is obtained by the minimization of

that yields I = 1,..,4

bR Rk- “du=O &k

I = 1,..,4

(3)

(4)

The implicit ODES system formed by spatial discretization of equations (3) c (4) can be writ.ten in the form &y)jr g(y), and then solved by a stiff integrator. In the present, case DASOLV -_ = --

s423

European Symposium on Computer Aided ProcessEngineering-5

(ml (ml o(l) (Jarvis and Pantelides, 1992) was used. The solution of the system y = [o(l) 0,ll 0,27 **,ao,1 9a0,2 7WY (1) (1) w+llT give us the node location as well the solution on the ‘-3%+1,17%+1,2~ **y4$1,1r 4$,2r nodes. The occuring singularities are dealt using Miller’s approaches. Nodal coalescence is eliminated by the definition of a minimum inter-nodal distance (Amin). If nodes become closer than A min Ythey start moving at equal velocity. When their separation increases and gets larger than A min they become independent again. Parallelism is avoided by the addition of a penalty term (4ir+1 - i~)~) to equation (2). The tunning parameters involved in the method formulation are Amin and e. Unfortunately there is no general way to choose the value of these parameters. Their selection is based on the user’s experience and their value may require an initial trial of the program. As an example, for solutions with deep fronts or schocks E may take values from low5 to 10-l to prevent parallelism. Fortunately, the solution uirnl does not depend on the value of the parameters used. Only nodal positions are affected by Amin and E. The treatment of boundary conditions needs a previous time derivation originating constraints to the minimization problem (2). The classical Lagrange Multipliers method is used to solve this problem. CAUSTICIZING

REACTION

MODEL

Causticizing reaction is an important step of the kraft process of pulp production being described by a front reaction model. The reaction occurs when lime mud (Coo) particles react with sodium carbonate (Na2C03) producing sodium hydroxide (NaOH) and calcium carbonate (CaCO3). The reaction is heterogeneous and occurs in a solid/liquid interface. The limiting step is the inward diffusion of CO;and the difusion of OH- outward the particles, as fig. 1 shows. For the modelling of the process it was assumed Cd,’ that diffusion controls the global process forcing the reaction to take place in a moving interface of infinitesimal thickness. The interface moves towards the centre of the particles and the carbonate and hydroxide ions are considFigure 1: Causticizing reaction scheme ered in chemical equilibrium. The following equations describe the process:

62COp- + at

aT2

e

E

zacOH-

th

CT

Z’(t)

&2,,2-

3

bt

=

a2Cco*-

1 Ve E

82

3

De 2accoa+TF-

dr

and

o
<

T

<

and

0 < T < d(t) z’(t)

< T <

(5)

1

1

(6)

subject to the following boundary conditions acOH-(o~

t,

=

o

(7)

a ac,,:-

&

acO,-

(9, t) (1, t)

a

=o =

(8) $

(@OH-

-

COH-

(0, t))

European Symposium on Computer Aided Process Engineering-5

s424

acco,l- 074

a

ac,,- -ar

=

acOH-

a

zI+

2

Kco2-

%oz& 3

-2

~

.I-

,I+

*I-

and initial conditions COHcCO,2-(? %(OH)2 %(OH)&,

=

cOH,T

(13)

0)

=

cco:,

(14)

CT, O)

=

CCa(OH)2i,

O)

=

o

(rye)

~~(0) =

0

2

z’(t)

T

<

z’(t)

<

T

5

R

d(O) = 0.8R

(15) (16) (17)

the reaction front velocity is given by

%o:-

dx’ dt=

z1+

&

1

I1 II-

(18)

CCa(OH)ain

0.228 mol/l 0.223 mol/l

600 s < t 5 2400 s 2400 s < t < 4200 s

0.211 mol/l 2.144 mol/l 2.254 mol/l

4200 s < t 5 6000 s 600 s < t 5 2400 s 2400 s < t 5 4200 s

2.279 mol/l 0 mol/l

4200 s < P < r 2

0.337 moljl 0 moljl

I‘ <

2.025 molfl

r 2

t 5

6000 s

z’(t) r’(t) z’(t) z’(t)

27 mol jl 7.2 x lO_”

ml/s

50 molfl 5 X IO-’ m/s 0.1 JI

X lO-g m

Table 1: Physical data for the causticizing reaction model MODEL IMPLEMENTATION Since the reaction front models are a particular case of distributed parameter models, the solution of this class of problems by the MFE will be described. OH- and CO:- concentrations are discountinuous variables at the interface. Contrarely, cubic Hermite approximations sssure the continuity of the solution on the nodes. Therefore, The the spatial derivatives at z’+ and z’- are different for both variables (COH- and Cco:-). interface is described by a node numbered as NFR. If ahFR,l is the solution on zr+ for variable is its spatial i and ahFR+l 2 its spatial derivative and ahFR,! is the solution on .zI- and ahFR,* derivate on zi-,

then ahFR+l,l = akFR,s and

a),FR+1,2

#

ahFR,4. This assumption implies the

EuropeanSymposiumon ComputerAided ProcessEngineering-5

s425

For the system addition of two variables (c~j,,~~+~,~= !*l,,+ and o$FR+1,2 = ylz,+). to be determined two additional equations from the boundary conditrons on the interface are required. The first, results from the spatial discretization and time derivation of (11) that yields

(19) The second is defined after spatial discretization relation

and time derivation of the equilibrium

2.4 2.3

COH

0.05

Figure 2: Solution profiles of Co,-

Figure 3: Solution profiles of C,,;-

1

0 ~~~~'/~~~'~~,~'~~~~'~~~~'~ 800 1,600 2,600 3,600

t

4.600

5.600

Figure 5: Front reaction evolution

Figure 4: Grid evolution

The velocity of the front is obtained by spatial discretization of the interface velocity -De BNFR

=

(&FR+1,2

-

a&FR,4)

(21)

CC~Kmi,

The grid comprises 13 nodes. In each subdomain 6 nodes are equidistributed. The subdomains are limited by the interface and fixed bounds, being the interface initially placed at r/R= 0.8 since an important part of the particles react outside the reactor. A node is placed coincident with the initial position of the front reaction. The penalty parameter used c = 10-l and Amin = 10A5. CONCLUSIONS The results obtained are in good agreement with t*he experimental data as well as the results

European Symposium on Computer Aided Process Engineering-5

S426 published

by Portugal

an

Pais (1989) and Levenspiel (1972) for a similar system.

Cc.:-

in

the unreacted lime zone since even there the equilibrium must be compulsorily satisfied. The simulation starts at 600 s because the preceeding step of causticizing, the slaking, takes approximately 10 minutes. Notice that different concentrations in the boundaries correspond to different reactors. The jump of particles along the reactors is compared to different step perturbations suffered by the concentration in the boundaries. The retention time in each reactor is sufficient for reaching a new steady state characterized by stable solution profiles and constant movement of front reaction (Fig. 2 and 3). The front reaction evolution is almost linear, except in the vicinity of T N 0 where it becomes faster due to the increasing role of diffusion phenomena (Fig. 5). Figure 4 shows a very smooth grid evolution due to pure diffusion profiles in both phases. The retention time in the three causticizers of a real pulp mill is 5400 s. The numerical solution indicates that for an equivalent reaction time the particles react to almost completeness. This fact proves the validity of the mathematical model and also the agreement between numerical and industrial data. This work shows that MFEM is particularly suited for this kind of problems. The nodes displacement give a good description of the interface’s location. This procedure is more acurate and “natural” than interpolation strategies coupled with scaling domain techniques used in fixed grid strategies. A CPU time of 224.64 s in a SUN SPARClO computer was obtained for this problem. _ NOMENCLATURE ar,k SI ,“I HI,k c C.(OQin

Ica coefficient of Hermite polynomials on 141 element position of node I time derivative of s interface’s position l& Hermite polynomial for element I

diffusion coefficient mass transfer coefficient equilibrium constant concentration of ions OHconcentration of CO:-

initial concentration

concentration

Acknowledgment: Investiga+io

of Ca(OH)*

This work was supported

by grant BD 1094/90

- RM from JNICT

of i in the bulk

(Junta

National

de

I. Theory”,

J.

Cientffica e Tecnolbgica).

REFERENCES [l] Baines, M.J., A.J. Wathen, “Moving finite element methods for evolutionary Comput. Phys., 79, 245, (1988)

problems.

[2] Jarvis, R.B., C.C. Pantelides, “DASOLV - a differential algebraic equation solver - version 1.2”, Centre for Process Systems Engineering, Imperial College of Science, Technology and Medicine, London, (1992) [3] Jima&, P.K., A.J. Wathen, “Temporal derivatives in the finite-element method on continuosly deforming grids”, SIAM J. Numer. Anal., 28, 990, (1991) [4] Levenspiel, O., “Chemical Reaction Engineering” - 2”6 edt., John W&y [5] Miller, K., R.N. Miller, “Moving finite elements. I”, SIAM J. Numer. [S] Miller, K., “Moving finite elements. II”, SIAM J. Numer.

Anal.,

d Sons, Inc., New York (1972)

Anal.,

18, 1019, (1981)

18, 1033, (1981)

[7] Murray, W.D., F. Landis, “Numerical and machine solutions of transient heat-conduction problems involving melting and freezing, Part I - Method of analysis and sample solutions”. Transactions of the ASME, (1959) [S] Pais, F.I.C.C. and Portugal, A.A.T.G., “A mathematical reactions”, Comp. Chem. Eng. (submitted) [9] Portugal, A.A.T.G, and Pais, F.I.C.C., “A mathematical CHEMPOR’89 - * Int. Chem. Eng. Conf., Lisbon (1989)

model for non-catalytic

liquid solid reversible

model for the causticizing

reaction”,

Proc.

[lo] Pipilis, K.G., “Higher order moving finite element methods for systems described by partial differentialalgebraic equations”, Ph.D. Thesis, Dept. of Chemical Engineering and Chemical Technology, Imperial College of Science, Technology and Medicine, (1991) [ll]

Powers, C.J., D.L. Cummings, “The dynamics of energy and mass transfer with reaction in hydride storage systems”, Hydrogen Energy, eds. T.N. Veziroglo em Proc. Hydrogen Econ. Miami, Energy Conf. Miami, Plinium Press, New York, (1975)

[12] Wathen, A.J., M.J. Baines, “On the structure Numerical Analysis, 5, 161, (1985)

of the moving finite element. equations”,

IMA Journal

of