- Email: [email protected]
Numerical strategies in solving gas-liquid reactor models—1. Stagnant films and a steady state CSTR
Compufers dzem. Engng, Vol. 15, No. 11, pp. 769-781, 1991 Printed in Great Britain.All rightsreserved
009s1354/91
$3.00 + 0.00
CopyrightQ 1991 Pergamon Press plc
NUMERICAL STRATEGIES IN SOLVING GAS-LIQUID REACTOR MODELS-l. STAGNANT FILMS AND A STEADY STATE CSTR J. J.
ROMANAINEN
and T.
SALMI
Department of Chemical Engineering, Abe Akademi, SF-20500, (Received
7 November
199O;fiMI
revision received 2 July 1991;
receivedfor
Turku, Finland publication
11 Jzdy 1991)
Abstract-Two alternative numerical strategies are compared in solving the isothermal tank reactor model using the following test reactions: a first-order and a consecutive-competitive reaction scheme. Since the steady-state tank reactor model consistsof the componentconservationequationsof the stagnantfilms
and the bulk phase mass balances the numericalsolution of the equations can be performed either simultaneouslyor sequentially.The simultaneousstrategyimpliesthat the bulk phase mass balancesare includedin the boundaryconditionsof the film equations.In the sequentialstrategythe iXm and the bulk phase equationsare treatedseparately.The example simulationsindicatethat the simultaneousstrategy based on spline collocation is faster than the sequentialstrategybased on orthogonal collocation and Broyden’smethod for the filmequationsand the bulk phasemassbalances,respectively. The simultaneous strategysuffers,however,from stabilityproblems; the sequentialstrategyis thereforeconcludedto be a more robust way to solve the steady-statetank reactormodel.
1.
INTRODUff
ION
Chemical reactions between gas and liquid components are common in industry. Typical processes include absorption, oxidation, chlorination and hydrogenation (Doraiswamy and Sharma, 1984). These reactions are performed in either a tank or a column type of reactor. The selection of the reactor type depends on many factors, e.g. reaction kinetics, mass transfer requirements and thermal effects (Rase, 1977). Gas-liquid reactors have traditionally been modeled using analyticai solutions. This can, however, be done only in restricted cases, where reaction kinetics are simple enough. In the general case the gas-liquid reactor model consists of two basic parts: the reactor bulk phase balances, and the model for the gas-liquid contact. The resulting mathematical models must thus be solved numerically. This work concentrates on finding suitable numerical strategies to carry out the numerical solution of the forementioned models. This first part will focus on the model of an isothermal steady-state continuous stirred-tank reactor (SCSTR). The resulting mathematical model of a SCSTR consists of a two-point boundary value problem (BVP) and a nonlinear equation system (NLE) arising from the gas-liquid contact and the bulk phase balances, respectively. Basically three methods have been used to solve the BVP arising from simultaneous diffusion and reaction. The finite-difference (Lee, 1968; Deckwer, 1976) and the multiple shooting (Darde et al., 1983) methods were previously the most common algorithms. The latter approach is not a very recom769
mendable because it suffers from stability problems. Newer methods applied in this context are based on collocation (Ho and Lee, 1980; Huang et al., 1980; Shah and Paraskos, 1975; Nielsen and Villadsen, 1986), which has been found to be quite efficient in solving BVPs. In the current work both collocation and finite-difference methods are tested. The NLE system arising from thi bulk phase mass balances is a typical problem in chemical engineering computations. Thus numerous NLE solvers are available for this problem. In this work the Broyden algorithm including constraints for the variables is used. The benefits of the Broyden method are the speed and the low computer memory requirement.
2. STOICZ-IIOMETRY
AND
KINETICS
The general case of S reactions and N components in a chemical system can in most cases be described by the power law kinetics:
R,=k,f&F’, j=l,..., i= 1
s.
(1)
Using the stoichiometric coefficient matrix y, the rates of the components can be written in matrix notation: r=y-R.
(2)
The reaction rate constants are assumed to follow the
Arrhenius temperature dependence: k,=k,OexpG.
J. J.
770 3. SIMULTANEOUS REACTION AND IN STAGNANT FILMS
ROMANAINEN
DIFFUSION
Assuming that the diffusion through a volume element consists only of the molecular contribution, Fick’s law can be used. For a component diffusing through stationary boundaries the corresponding flux (No is given by equation (4):
Provided that the diffusion proceeds in one direction only in a fluid with a constant density, the wellknown balance for simultaneous diffusion and reaction in steady state is obtained: Diz+r,=o.
It is possible to get an analytical solution of equation (5) if the reaction rate ri is a linear function the corresponding concentration ci solely, or no reaction occurs in the stagnant element (ri = 0). This restricts the analytical solutions to first-order irreversible kinetics. Other, especially second-order bimolecular kinetics, have also been solved analytically, applying some additional assumptions or simplifications (e.g. van Krevelen and Hoftijzer, 1948). The film theory is commonly used to model the stagnant elements between the bulk phases and the gas-liquid interface. Film theory is used because of its simplicity and ease of implementation. Other theories, surface renewal (Dankwerts, 1970) and penetration theory (Higbie, 1935), have been developed, but they have not proved to be superior compared to the simple film theory. Therefore, the film theory is still used in most numerical studies concerning gas-liquid reactions. Applying the film theory for the stagnant elements, the following boundary conditions are valid: ci=c;
at
2 = 0,
ci=cF
at
2 -6,
(6)
where S is the thickness of the laminar layer. The gas-liquid interface is noted by z = 0 which sets z = 6 at the bulk phase boundary. This treatment is referred to as the film theory and it is commonly used to model the behavior of the stagnant elements in steady-state reactors. A schematic picture of the gas-liquid contact according to the film theory is given in Fig. 1. For numerical reasons the laminar layer spatial variable z is often scaled by the film thickness 6. The dimensionless variable x = z/6 is always within the interval [0, 11. Usually the chemical reaction takes place in only one of the phases. Here we assume that the reactions occur in the liquid phase, which is also the most common situation in industrial cases. Thus physical diffusion prevails in the gas film. Integrating equation (5) twice (ri= 0) and applying the boundary conditions (6) the flux can be obtained from equation (7): N&, = &(c!$;
- c;,,),
Since only physical diffusion is present in the gas phase, the fluxes equal on the film boundaries: N& = N$. Furthermore, the fluxes on gas-liquid interface are equal:
(9) both
N;, = -_k.
Llquld
Film
Bulk
b
Cl C
reaction
dlffuslon
6e
end
dfffusf on
0
sides
of
the
dCLi = N&, = kGi(c& - K&,), (10) L’( dx > x=0 where K, is the equilibrium constant between the component concentrations on both sides of the phase boundary and is related to the component vaporliquid equilibrium (VLE) constant or Henry constant if Henry’s law is used to describe the VLE. Equation (10) can be written in a more suitable form for
=a
\
(7)
where kGi is the gas-phase mass-transfer coefficient given by the film theory:
Liquid
Gas Film
Gas Bulk
and T. SALMI
%
Fig. 1. Schematicpicture of the gas-liquid contact according to the film theory.
b II
Solving gas-liquid reactormodels-1
f
T~_r
.
1 nPGi
771
A
.
Nb
n
a
L1
11 quid
*
4
c
n
FLi
I
film i
GAS BULK
6L
LIQUID BULK THE FILMS Fig. 2. Schematicpictureof the steady state CSTR.
numerical treatment, and the boundary conditions are thus given by:
- kLi dc,
--
k 01
d-x
-
c:, + K$,
= 0
CLi
czi
=
at
x = 0,
at
x=
1.
5. I. The sequential solution strategy (11)
As a measure of the accelerating effect of the chemical reaction on the absorption rate, the enhancement factor (E,) is defined as: Ei =
Ki k,(&
-
ct)
’
(12)
where the denominator gives the rate of the pure physical mass transfer. 4.
REACTOR
BULK
PHASE
BALANCES
Because the essential part of the current work concentrates on the numerical aspects of the problem, isothermal models were assumed to be sufficient in order to show the behavior of the numerical solutions. The stirred reactor consists of one gas and one liquid phase including one input and one output stream of both phases to and from the tank (Fig. 2). As mentioned previously, the reactions are restricted only to the liquid phase. Following these assumptions, the mass balances for the liquid and gas phases are given by equations (13) and (14), respectively: nFLi
-
. flPLi %GI
+Nb,fIt-(l-~~)V~r~=O, -
’ %x3r
5. NUMERICAL
-N&A=O.
the reactor bulk phase mass balances (13) and (14). The complete model is thus a combination of a two-point BVP and an NLE.
There are two possible strategies to solve this system of BVP and NLE. The traditionally used approach is to solve the film equations separately as a two-point BVP during the main iteration of the reactor mass balances. This leads to a sequential type of solution illustrated in Fig. 3. This traditional approach demands a very heavy computational effort, if the film equations must be solved every time the bulk phase model is called. A flag can be used based on the difference of the values of the variables since the previous call. This restricts the need to solve the BVP only, when the concentrations in the bulk phase have changed enough by the NLE solver in the main iteration, thus speeding the computation considerably. The changes resulting from this modification are shown in Fig. 4. The resulting main iteration loop (NLE) was solved by the Broyden algorithm (Broyden, 1965) and for the underlying flux update several BVP modifications were tested. These BVP algorithms included among others the COLSYS (Ascher et al., 1981) software and a new implementation of the
(13) (14)
APPROACH
The complete model of the steady-state CSTR consists of the combination of both the film equations (5) with the boundary conditions (11) together with
Fig. 3. The sequentialstrategyof solvingthe SCSTR model.
J. J. ROMANAINENand T.
772
Have variables changed enough from last call 7
YES
I Solve the Stagnant film model
( I
NO
SALMI
The boundary conditions are written partially inside the HSR solver. The nonvolatile components will be calculated using equation (19) for the second derivative at the gas-liquid interface. This method is called fake bounduries (Ames, 1977). The solution at x = 0 is approximated by a parabola having a maximum or a minimum at x = 0. Because of symmetry cLi(x_,) = +(x,), equation (18) reduces to: 2[cL&)
Fig. 4. The sequential strategy to solve the SCSTR model with a modification leading to less frequent stagnant film solutions.
finite-difference type of method HSR (Romanainen, 1989). Also the multiple shooting (DBVPMS) and finite difference (DBVPFD) algorithms provided by the IMSL (1987) for solution of BVPs were tested. Because of unreliability these codes were, however, discarded after preliminary tests. The COLSYS software applies orthogonal collocation on splines to approximate the solution of the BVP. The COLSYS requires the BVP model to be written in the following way for the second derivative:
f;(x)=$ -&
(1%
The boundary condition routine of COLSYS is divided into two logical sections according to the boundaries (x = 0, x = 1). At the gas-liquid interface the boundary conditions are given by equation (16): &(O)
=
_k,i dci k
dx
Or (1
_
+mLi-&i
(16)
Ax=-,
interface (x = 1) by equation n dp =
(17)
The initial guess for the component concentration profiles in the film was supplied by a linear function between the boundary concentrations. The Jacobian matrix was evaluated numerically for A(x) and analytically for both boundaries gi(0) and gj(l). In a finite-difference method using an even grid of discretization points, the second derivatives are normally evaluated with equation (18). The new implementation of a finite-difference type of solver (HSR) (Romanainen, 1989) is based on the idea of evaluating the second derivatives with an equation similar to (18), but at the same time using a refined grid. In addition, this was achieved without the need to add the number of discretization points in the grid: =CLi(XI--I)-2CLi(Xf)+CLI(Xl+1)~
X= X,
CLAXI)
-
Ax2
CLi(%)
9
(20)
where
x-0
&(l) = CL1- & = 0.
(19)
Equation (19) sets the gradient to zero at the gasliquid interface and thus no volatility can occur with the component in question. Now there are two possible approaches to the remaining boundary conditions. One is the dynamic approach in which the film concentrations are assumed time dependent, thus leading after discretization to a set of ordinary differential equations (ODES) instead of the original BVP [equation (S)]. This homotopy approach has very nice features including that the solution will always be within physically meaningful limits. The essential features of the method have been studied elsewhere (Haario and Seidman, 1988, 1989). The other possibility is to use the stationary approach in which the equations in the discretization points form an NLE system and solve this with an NLE solver. The problem with the HSR solver described above is the discretization. The even discretization inevitably implies that the gradient is very sensitive to the number of points used because the gradient is calculated by equation (20): Ax
= 0 and at the film-bulk (17):
- CLic%)l Ax= .
(18)
1 nLip - 1
number of discretization points,
(21)
In practical cases, ndpis between 11 and 41, giving Ax-values from 0.1 to 0.025. For very fast reactions, however, this will not be sufficient to compute the gradient accurately. To overcome this problem, there are basically two possibilities. First it is possible to use some kind of interpolation to calculate the gradient more accurately near the boundaries. Another choice is to make the discretization differently, e.g. by coordinate transformation (Versteeg et al., 1987). In this work the following discretization scheme was adapted. The scheme is based on dividing the remaining film into two equal parts as we approach each boundary (Fig. 5). The benefits of this discretization scheme are, on one hand, the accuracy near the boundaries and, on the other hand, the simplicity of the implementation. The accuracy is achieved because the discretization interval at the boundaries is now:
Solving gas-liquid reactormodels-l
I
III I 10 . . .
l/8
114
773
I
I
l/2
I III 7/a
3/b
. . .
1
X
Fig. 5. The half discretizationmethod. Ax=$
n,=-.
In these cases the integration of the film equations would have to be continued for a very long time causing the method to be impractical.
x’
% -1
5.2. The simultaneous solution strategy
”
L
Equation (22) will give the discretization interval Ax practical values from 2-5 to 2-20, which guarantees sufficient accuracy to the gradients at both boundaries (x = 0, x = 1). Usually the gas-liquid interface boundary is more important, since the profiles of the soluble, reactive components there are sharp (especially with fast reactions). This is also the case in the examples presented in this paper. On some other cases, however, the reactive liquid components will get sharp profiles near the bulk phase boundary because of fast and/or complex kinetics in the film. The half discretization scheme is thus using the most refined grid exactly where the greatest accuracy is required. The other benefit of the implementation is that this discretization allows us to use equation (18) for the second derivatives, since the point x,_, is always x0 and the point x,, , is the next point, when computing points x, from the center to the gas-liquid interface. The same applies for the other half of the film. The only drawback of this scheme is that the center of the film (x = 0.25, x = 0.5, x = 0.75) will be calculated always using Ax = 0.25, which may cause inaccuracy in the results. This can be avoided by combining the even and halving discretion schemes using the even scheme at the center of the film and the halving scheme near boundaries. Another method to improve the accuracy was also tested. It is a method commonly used to estimate the Auxes in catalyst particles (Froment and Bischoff, 1979) where simultaneous diffusion and reaction occur. This method is based on the total balances of the film elements:
Jo Estimating the integral on the right-hand side of equation (23) by the Simpson formula and using the other gradient that was near zero (Nkj in this case), corrected values for the gradients at the steeper end of the profile were obtained. Other methods such as splines on points near the boundaries were also tested, but they did not provide much improvement to the accuracy of the results. An additional problem with the dynamic approach of the HSR solver is the cases with slow reaction.
The alternative strategy to solve the film and bulk phase equations is to solve them simultaneously. Here the reactor model is converted to a boundary value problem. Equation (5) remains unchanged. The difference now is that the bulk phase mass balance equations (13) and (14) are included in the corresponding boundary conditions (16) and (17). The resulting BVP was solved by the collocation software COLSYS. The software requires the model equations to be written in the following way for the simultaneous strategy: (24) and for the film-bulk similarly:
g,(l) = -
2
boundary
x = 1 we
+ -cc;~~++~Li+(l.O-to)VRri
( >
&al/,
get
3 (25)
while equation (15) remains unchanged. The transformed concentrations CL, are obtained from: Cii(X) = & + c&r-(1 - x), where c,,, denotes concentration. 6.
the
maximum
(26) liquid
feed
EXAMPLES
All examples presented here were programmed in FORTRAN and executed on a DEC Micro VAX-II computer. 6. I. The sragnant jilm model solution Two different reaction schemes are used to test the different boundary value problem solvers. These examples deal only with the film equations since the reactor bulk balances are not needed in this section. 61.1. One first -or&r irreversible reaction. The reaction scheme is an irreversible first-order reaction where component A originally in the gas phase diffuses through the films to the liquid phase reacting simultaneously to a liquid component B: A(G) --+ B(L).
(27)
J. J. ROMANAINEN
774 Table
I.
Solution methods for the irrwrcrsible reaction scheme
first-order
NaaIC
Method
ANA CSN CST HSE HSH HDE HDH
The analytical solution COLSYS normal COLSYS with variable transformation HSR stationary with even discretization HSR stationary with half discretization HSR dynamic with even discretization HSR dynamic with half discretization
Table 2. Parameters used in simulation of the firstorder irreversible reaction scheme (example 6.1.1) Parameter K# k, kc, 6, ct cs*
Value 1.0 1.0 x lo-‘m s--l 1.O x 102’ m s- ’ ( = no gas film resistance) 1.0 x lo-‘m 0.0 mol dm-’ 0.01 mol dm-’
The solution methods studied for the film equations are listed in Table 1. The parameters used in simulations are presented in Table 2. To get a good coverage of the behavior of the solvers, a wide range of reaction rate constant values were used. Figure 6 shows the CPU-times demanded by the tested solvers relative to the analytical solution. The number of discretization points used in the HSR method was 21, which was found to be enough for reaction rate constant values less than 10s s-l.
and
T.
SALMI
Figure 6 shows that the stationary HSR method with either discretization scheme is the fastest overall solver in these cases. The COLSYS solver is also fast except for some individual values of the reaction rate constant (k, = lo’, l@, lo’, 10’s_‘). Thlb might indicate some problems connected to the implementation of the COLSYS software. Neither the discretization scheme in HSR methods nor the variable transformation in COLSYS seem to affect the computation time. The dynamic HSR methods required about four times the CPU-time compared with the stationary ones and they are thus discarded already at this point. The accuracy of the solution was examined by the enhancement factor, the value of which was compared to the analytical solution. The result of this comparison can be seen in Fig. 7. Despite the problems with COLSYS in computation time, the solution accuracy was found to follow the analytical solution extremely closely. The variable transformation did not effect the accuracy in any detail. The HSR methods, on the other hand, show a clear dependence on the discretization. The best version of the HSR solver-at least in this example--seems to be the stationary HSR solver with the half-discretization scheme. This scheme (HSH) is very accurate with 21 points up to reaction rate coefficient values k, = IO5s-l. The even discretization schemes with 21 points fail after values k, 3 10 s-l. There cannot be seen any differences in the accuracy of the stationary and the dynamic variations of the HSR solver.
Fig. 6. Relative computation time as a function of reactionrate coefficient.One first-orderirreversible reaction (example 4.1.1).
Solving gas-liquid reactormodels-l n.3
-A+ -t *ANA -A-
HDH HSH HSE HOE
Reaction
rate
coefficient
k,
Fig. 7. Enhancementfactor as a function of reactionrate coefficient.One first-order irreversiblereaction (example 6.1.1). The accuracy of the finite-difference methods (HSR) was further improved by the use of equation (23) and the Simpson formula. This approach was applied to the even discretization variations of the HSR solver (methods HSEm and HDEm). The enhancement factor values achieved with these modified routines are shown in Fig. 8. This method improved the enhancement factor values only up to rate coefficient values k, < lo3 s-‘. Therefore the even discretization versions were discarded at this point. The overshooting behavior of the enhancement factor in applying equation (23) is probably due to numerical problems. An examination of the concentration profiles in the film as function of the reaction rate constant values can be seen in Fig. 9. These profiles show the expected behavior of steepening curves as the reaction rate constant increases. The remaining solvers all follow the analytical solution quite nicely. 6.1.2. Two consecutive-competitive irreversible second-order reactions. This example was applied only to the remaining solvers which passed the tests in the first example. These solvers were the stationary HSR method with half discretization scheme (HSH) and the COLSYS method with the variable transformation (CST). The reaction scheme is presented in equation (28) and the corresponding parameters used in the simulation are listed in Table 3:
A(G) + B(L) -
C(L),
A(G) + C(L) +
D(L).
The simulations were performed with reaction rate coefficient values of lo-‘, 10’ and 103m3 kmol-‘s-l for both rate constants in every possible combination resulting in nine simulation runs. The resulting CPU-
775
time usages are depicted in Fig. 10, where it can be seen that COLSYS is again the fastest solver except for the last three cases where k, = lo3 m3 kmol-t s-l and k2 has all three values in turn. The accuracy could not be compared since an analytical solution is not obtainable. On the other hand, both methods seem to give similar results, so the accuracy is probably not a problem at least with these values of the rate constants. 6.1.3. Conclusions of the stagnant jilm model solu tion. COLSYS seems to be both the fastest and the most accurate method, except for some cases, where problems might occur. The variable transformation [equation (26)J did not seem to alter the behavior of COLSYS in any way. A few details of the COLSYS implementation were found to be worth noticing. The number of subintervals should always be set to one and let the software decide how to increase it. The problem should always be defined as nonregular although this might cause problems in some linear cases. The tolerances should be at least IO-*, otherwise the results will not be accurate. Finally it was noticed that giving the values from a previous failed calculation back to COLSYS as a new initial guess often lead to a visible solution. These observations lead us to believe that the problems encountered with COLSYS were mostly due to an inadequate initial guess of the film profiles. The problems arising from this can be minimized restricting the mesh size in COLSYS or using an estimate based on a previously successfully computed profile as an initial guess. These methods of improving the performance of COLSYS were found to be
lo:: Reaction
101
rate
105
106
coefficient
10'
ld
k,
Fig. 8. Enhancementfactor as a function of reaction rate coefficient. One first-order irreversible reaction. Finite-difference methods with even discretizstionscheme (HSE. HDE) and eauation(23) correction(HSEm. HDEml . corn&red td the an&ical solution (ANA). ’
J. J.
776
ROMANAINEN
and T. SALMI
(a) 0.010 0.009
-a-
CSN
0.008 0.W-f 0.006 0.005 0.004 0.003 0.002 0.001 0.000 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Dimensionless
film
thickness
,
x
(c)
(b)
Dimensionless
film
thickness
x
Dimensionless
film
thickness
x
Fig. 9. Component concentrationprofiles in the liquid film with differentvalues of the reaction rate coefficient(example 6.1.1): (a) k, = 1.0 x IO”s-‘; (b) k, = 1.0 x 102s-‘; and (c) k, = 1.0 x 104s-‘. successful; hence it was chosen as the stagnant film model solver. The HSR routines were still under development during this work, so they are not yet comparable with established software like COLSYS. However, Table 3. Parameters used in simulation of the two consecutive-competitive reaction schemes (example 6. I .2) PaIXmeter
4
b
k
V&E 1.0
= km OA
b,
=
km
1.0x
IO-‘m
s-’
1ti5 m s- ’ ( = no gas film resistance) I.0 x 10mSm 0.0mol dm-’ 0.01 mol dme3 1.O mol dm-’ 0.0 mol dm-’ 1 .O x
after some further development, they might prove to be competitive as stagnant film solvers. The half discretization scheme is obviously superior to the even one in terms of the result accuracy. The computation time is about the same with both methods. The accuracy is, however, not so high in the center of the film, and this might cause problems with some more complicated cases, where steep gradients are expected not only at the boundaries, but also at the center of the film. As discussed earlier, this can be helped combining both discretization schemes. Another way would be to develop some kind of adaptive scheme for the discretization during the computation, based, e.g. on the change of the derivatives. The dynamic HSR methods offer the benefit of approaching the solution always via a physically
777
Solving gas-liquid reactor models-1
k,=lO-’
k,=lO-’
k,=lO+
k,=lO
k+= IO-’
k2= IO’
k2= lo3
k*=lO-’
React
’
ion
k,=
10’
k2=70’ rate
k,=
IO’
kz= IO3
k,= lo=
k,=103
k,:
kglO-’
k,=lO’
k;lO”
IO3
coefficients
Fig. 10. Relative computation time as function of reaction rate coefficients. Two consecutive second-order reactions(example6.1.2).
meaningful path. The problem with these dynamic methods applied in stationary cases is the selection of the integration interval. When then reaction is a fast one, the choice of the integration interval is not critical, but with very slow reactions the integration would have to be continued for a very long time to obtain the steady state. This excludes this kind of methods from practical applications, although they might be appropriate for some special cases with difficult reaction-diffusion schemes. The stationary HSR method with the half discretization was the only one comparable to COLSYS. The benefit of this method over COLSYS is its easy implementation. The advantage of the physically viable path to the solution is regrettably lost in the stationary version. The final conclusion of the stagnant film solvers is that COLSYS seems to offer the best performance at the present moment. The stationary HSR model is the other possible choice, but it should be developed further.
6.2. The steady-state
CSTR
model solution
These latter examples 6.2.1 and 6.2.2 use variations of the two consecutive-competitive second-order reactions scheme [equation (2811. The examples are used to demonstrate the different characteristics of the sequential and simultaneous strategies introduced in connection with the steady state CSTR model. The example runs were obtained using COLSYS as the BVP solver.
6.2.1. Two consecutive-competitive second-order irreversible reactions. The reaction scheme used in these examples [equation (2911 was the same as applied by Darde and co-workers (Darde et al., 1983): A(G) + B(L) - ’
W-1 + G(G),
A(G) + C(L) - ’
W-1 + G(G),
(29)
where letters A-G correspond to the compounds listed in Table 4. The parameters used for the simulation were taken partly from Darde et al. (1983) and are listed in Table 5. Furthermore, all component properties were taken at T = 273 K and p = 0.1 MPa, for the vapor-liquid equilibrium Henry’s law was used and the gas phase was assumed to be ideal. The liquid phase mass transfer coefficients were set to 10m5m s-l for those components in the liquid phase whose real kLj values were not available. Tabk 4. Components in the example reaction [equation(29)) used in the simuktion of the steady stateCSTR model: chlorinationof ~aru-cresol
A
B C D E
ChlOCiIlC Pam-crrsol
Monochloro-para-c.resol Dichloro-porn-cresol Tetrachloromethane Nitrogen Hydrogenchloride
Gas Liquid Liquid Liquid Liquid GaS Gas
J. J. ROMANAINEN
778
and T. SALMI
Table 5. Parameters used in simulation of the steady state CSTR model. Taken partly from Darde ef al.
----
(1983)
of components Number of reactions KA kl k, k,, kLB kLc Number
k,, (i= A, B, C, D, E, F, G)
Simultaneous strata Seq”cntiol strategy
7 2 0.020 5.620 dm” mol-’ SK’* 0.107dm’mol-‘s-‘P 7.00 x 10mSm s-la 3.77 x lo-Srn S-In 3.50 x 10m5m s-la 1.00 x lo*sms-’ ( = no gas resistance) 1.486 x 10-5m’ 6.8714 mz rne3’ 0.2 1 kmol h-’ 40 kmol h-l 205 kmol h-’ 1700 kmol h-’
‘From Darde et al. (1983). Some values modified.
0
The rate of the reaction with respect to mass transfer can be described with the Hatta number. The Hatta number for the first reaction {equation (29)] corresponds to an intermediate or slow reaction (< 1.22) achieving the maximum value at the feed composition. The Hatta number for the second reaction is even smaller (c 0.3) indicating a slow reaction. The actual simulations were done in a sequence of an increasing liquid residence time in the reactor. This enabled the use of a good initial guess from the solution at the previously computed residence time. The results from the simulation sequence is
Residence
time
-r (s)
12. Component concentrations in the liquid bulk phase at different liquid residence times. Simulation of the steady state CSTR model in example 6.2.1. Irreversible reactions. Fig.
presented in Fig. 11, where the computation time is plotted against the liquid residence time. The obvious result from this is that the simultaneous strategy is considerably faster than the sequential one. The simultaneous strategy, however, displays some disturbing behavior at some residence times. This behavior can be seen in Fig. 12, where at z = 1260 s and
25C
B
Simultaneous
[77721
Sequential
strategy strategy
Residence
Fig.
11. Relative
computation
time
T
ts)
time of the steady state CSTR model with different liquid residence times (example 6.2.1). Irreversible reactions.
Solving gas-liquid reactor models-l z = 1440 s the correct solution can be achieved only after re-introducing values from an originally unsuccessful calculation and recalculating the simulation. Otherwise both strategies give comparable results. The physically interesting result given in Fig. 12 is the maximum in the concentration of the intermediate component C, which in this case happens to be the desired product. This simulation would indicate the suitable liquid residence time of about r = 1080 s for the operation of the actual reactor. 6.2.2. Two consecutive-competitive second-order reversible reactions. The steady state CSTR simulations were continued by modifying the original reactions [equation (29)J to include some reversibility. This resulted in schemes (30): A(G) + B(L) - ’ A(G)
+ C(L) +
C(L) + G(G), D(L) + G(G).
(30)
and (31): A(G) + B(L) --!A(G) + C(L) +
4
C(L) + G(G). D(L) + G(G).
(31)
These reaction sequences are hypothetical, but the physical properties from Table 5 were used for the sake of comparison. Here the simulations were performed using 0.1, 1 and 10 times the values of
the reaction rate coefficients of reactions 4 and 3 in respect to the values of reactions 1 and 2, respectively. The computation time with different reactor residence times is shown in Fig. 13 for reaction scheme (30). Since the computation time is approximately equal for both strategies, this indicates that the reversibility of reactions causes problems for the simultaneous solver. The computation time is especially high at the end of the simulation sequence which counts for long residence times in the reactor. Plotting the component concentrations at various residence times (Fig. 14) shows some instability of the simultaneous solver at higher residence times. In the last simulation sequence also the first reaction was made reversible [equation (3111 to verify the assumption of reversibility affecting the performance of the simultaneous solver. This resulted in a clear proof of the assumption, since the simultaneous solver was not able to solve more than few first points with low residence times, after which it failed. Recalculating the failed cases successively by returning the failed result to COLSYS did not provide the solution anymore as it did in the previous examples. The time used by the sequential solver was about the same as in the other reaction schemes tested. The resulting component concentration profiles as a function of the residence time are plotted in Fig. 15 for one example simulation, where k3 = 10 * k, and k4=
10.kl.
The difficulties with the simultaneous solver were first assumed to originate from a bad initial guess of
8OC Simultaneous
pf7
Sequential
strategy
strategy
6OC
;;
z i
779
4oc
E r 2 ” 2oc
0
Fig. 13. Relative computation time of the steady-state CSTR model with different liquid residence times (example 6.2.2). Irreversible and reversible reactions.
J. J. ROMANAINEN
780
and T. SALMI 0.25
0.25 n^
n^ E $
E $
0.20
.A.
.s s .-
E ._
0.15
CD
0.15
P E iti 5 0.10 0
5 E : 8
0.20
0.10
E E g E ;:
E fi n 0.05 E :: 500
1000
ResTdence
1500 time
2000 -r
2500
:
3
(s)
C*
CC 0.05
ca 0
500
1000
Residence
1500 time
2000 -r
2500
3
(s)
Fig. 14. Component concentrations in the liquid bulk phase at different liquid residence times. Simulation of the steady state CSTR model in example 6.2.2. Irreversible and reversible reactions.
Fig. 15. Component concentrations in the liquid bulk phase at different liquid residence times. Simulation of the steady state CSTR model in example 6.2.2. Reversible reactions with k3 = 10. k2 and k, = 10. k, in reaction (31).
the profiles. However, this is not the only source of problems, since the behavior of the simultaneous solver was predictible and got worse, when adding reversible reactions into the scheme. The simultaneous solver first started to deviate from the solution slightly and then failed. At first the failure was corrected by the successive guessing of the profile. When the residence time was still increased, this lead the whole method to fail. At first some of the concentrations became negative, the successive substitution found a physically nonfeasible solution and always converged to it, thus failing completely in computing the rest of the sequence. Because of this behavior, some modifications to the tolerances were tested with COLSYS. Making the tolerances lower resulted the method to be faster, and able to compute a few more simulations. However, this did not help very much, especially with the last scheme [equation (3111. In further test simulations, also different mesh size restrictions were tried, but without results. 6.2.3. Conclusion of the steady state CSTR model solution. The conclusion of the simultations were that the implementation of the simultaneous strategy was not a reliable one for an environment, where different types of reaction schemes are to be simulated, especially in the case of reversible reactions. In restricted use with a known and tested reaction scheme, however, it was sometimes superior to the sequential solution strategy. The problems with the simultaneous solution strategy were mostly connected with the problems of COLSYS software, and should not be taken as principal defects of the approach. The unstable behavior is probably connected with the boundary conditions,
and their very nonlinear nature, which is, furthermore, changing during the simulation. The sequential solution strategy implemented in this work using the Broyden method for the bulk phase mass balances and COLSYS for the stagnant film model equations performed well. The problems found with the COLSYS software were overcome by the NLE software, thus resulting in a reliable combination to solve steady-state CSTR models. Further studies should thus be directed into other boundary value problem solvers as well as to some other updating methods of the bulk concentrations and the general behavior of very nonlinear boundary conditions. The COLSYS software seems a very good BVP solver provided that the necessary modifications concerning very nonlinear boundary conditions could be added. The CSTR model can also be used to describe systems with more complex hydrodynamics (plug-flow, axial dispersion) by introducing the tanks-in-series concept. The numerical methods studied here can in principle be applied to these cases. Since the use of the tanks-in-series concept is restricted to cocurrent columns, the actual column simulation should be done with a dedicated model suitable for both cocurrent and countercurrent columns. NOMENCLATURE A = Interfacial area in reactor (m2)
0 = A/V, (m* m-‘) ci = Component concentration (km01 m-‘) c; = Transformed component concentration D, = Diffusion coefficient (m* s-l) Er I Enhancement factor
(kmol m-‘1
Solving gas-liquid reactor models-l E, = Activation energy of a reaction (kJ mol-‘)
f = Model function
g = Boundary condition K, = Gas-liquid equilbrium constant k, = Reaction rate coefficient; the dimension depends on kinetics, for n&order kinetics [(m’ kmol-I)(“- I) s-i] k,? = Arrhenius law frequency factor (same as k,) kLi = Liquid mass transfer coefficient (m s-‘) k, = Gas mass transfer coefficient (m s-‘) iV = Number of components N, = Flux (km01 rnT2 s-‘) n = Amount of substance (mol) ri = Flow of substance (mol s-‘) ” dp= Number of discretization points p = Pressure (MPa) R = Gas constant (8.314) (J k-l mall’) R = Vector of R,; the length of R is (S) (same as Rt) R, = Rate of reaction (km01 m-r s- ‘) r = Vector of r,; the length of r is (N) (same as r,) ri = Component reaction rate (kmol rno3 s-l) S = Number of reactions T = Temperature (K) V = Reactor volume (m3) v = Volumetric flowrate (m3 s-i) y = (N . S) Matrix of \I~ x = Film coordinate, drmensionless z = Film coordinate (m) Greek
symbols
a- = Component
$ = co = z= Y#=
order in a reaction Film thickness (m) Gas hold-up Residence time in a reactor (s) Stoichiometric coefficient
781
work was carried out as a part of the Computer Aided Reactor Design project financed by Teknologian kehittimiskeskus (TEKES), Kemira Oy, Neste Oy, University of Oulu and Abe Akademi.
Acknowledgement-This
REFERENCES Ames W.
F., Numerical Methods for Partial D#erentiaI 2nd Edn. Academic Press, New York
Equations,
(1977). Ascher U., J. Christiansen and R. D. Russell, ACM Trans. Math. Sofftware 7, 209 (1981). Broyden C. G., Math. Comput. 19, 577 (1965). Danckwerts P. V., Gas-Liquid Reactions. McGraw-Hill, New York (1970). Darde T., N. Midoux and J.-S. Charpentier, Entropic 19,92 (1983). Deckwer W. D., Chem. Engng Sci. 31, 309 (1976). Doraiswamy L. K. and M. M. Sharma, Heterogeneous Reactions, Vol. 2. Wiley, New York (1984). Froment G. F. and K. B. Bischoff, Chemical Reactor Analysis and Design. Wiley, New York (1979). Haario H. and T. I. Seidman, Proc. Minisymposium on Numerical
Methodr
for
Semiconductors
and
Magnets
(Neittaanmgki, Ed.). University of JyviiskylL, Dept of Mathematics, Report 42, JyvLskyll (1988). Haario H. and T. I. Seidman, Reaction and diffusion at a gas/liquid interface. Submitted for publication (1990). Higbie R., Trans. Am. Inst. Chem. Engnrs 35, 365 (1935).
Ho S.-P. and S.-T. Lee, Chem. Engng Sci. 35, 1139 (1980). Huang D. T.-J., J. J. Carberry and A. Varma, AIChE JI 26, 832 (1980).
Subscripts
F = G = i= i = I= L = P= R = A B
. i;-1
Reactor feed Gas phase Component Reaction Discretization point in the film Liquid phase Reactor outlet (product) Reactor
= Components A, B, . . .
Superscripts
b = Bulk phase s = Gas-liquid interface
IMSL, Znternational Mathematical and Statistical Libraries, Vol. 2, Houston, Texas (1987). van Krevelen D. W. and P. J. Hoftijzer, Rec. Trav. Chim. Pays-Bas 67, 563 (t948). Lee E. S., Quasilinearization
and
Invariant
ZmbeaWng.
Academic Press, New York (1968). Nielsen P. H. and J. Villadsen, Chem. Engng Sci. 41, 1655 (1986).
Rase H. F., Chemical Reacror Design for Process Plants, Vol. 1. Wiley, New York (1977). Romanainen J., Numerical strategies in solving gas-liquid tank reactor models, Lit. Technol. Thesis, Helsinki University of Technology, Helsinki (1989). Shah Y. T. and J. A. Paraskos, Chem. Engng Sci. 30, 465 (1975).
Versteeg G. F., J. A. M. Kuipers, F. P. H. van Beckum and W. P. M. van Swaaij, Chem. Engng Sci. 44, 2295 (1989).
009s1354/91
$3.00 + 0.00
CopyrightQ 1991 Pergamon Press plc
NUMERICAL STRATEGIES IN SOLVING GAS-LIQUID REACTOR MODELS-l. STAGNANT FILMS AND A STEADY STATE CSTR J. J.
ROMANAINEN
and T.
SALMI
Department of Chemical Engineering, Abe Akademi, SF-20500, (Received
7 November
199O;fiMI
revision received 2 July 1991;
receivedfor
Turku, Finland publication
11 Jzdy 1991)
Abstract-Two alternative numerical strategies are compared in solving the isothermal tank reactor model using the following test reactions: a first-order and a consecutive-competitive reaction scheme. Since the steady-state tank reactor model consistsof the componentconservationequationsof the stagnantfilms
and the bulk phase mass balances the numericalsolution of the equations can be performed either simultaneouslyor sequentially.The simultaneousstrategyimpliesthat the bulk phase mass balancesare includedin the boundaryconditionsof the film equations.In the sequentialstrategythe iXm and the bulk phase equationsare treatedseparately.The example simulationsindicatethat the simultaneousstrategy based on spline collocation is faster than the sequentialstrategybased on orthogonal collocation and Broyden’smethod for the filmequationsand the bulk phasemassbalances,respectively. The simultaneous strategysuffers,however,from stabilityproblems; the sequentialstrategyis thereforeconcludedto be a more robust way to solve the steady-statetank reactormodel.
1.
INTRODUff
ION
Chemical reactions between gas and liquid components are common in industry. Typical processes include absorption, oxidation, chlorination and hydrogenation (Doraiswamy and Sharma, 1984). These reactions are performed in either a tank or a column type of reactor. The selection of the reactor type depends on many factors, e.g. reaction kinetics, mass transfer requirements and thermal effects (Rase, 1977). Gas-liquid reactors have traditionally been modeled using analyticai solutions. This can, however, be done only in restricted cases, where reaction kinetics are simple enough. In the general case the gas-liquid reactor model consists of two basic parts: the reactor bulk phase balances, and the model for the gas-liquid contact. The resulting mathematical models must thus be solved numerically. This work concentrates on finding suitable numerical strategies to carry out the numerical solution of the forementioned models. This first part will focus on the model of an isothermal steady-state continuous stirred-tank reactor (SCSTR). The resulting mathematical model of a SCSTR consists of a two-point boundary value problem (BVP) and a nonlinear equation system (NLE) arising from the gas-liquid contact and the bulk phase balances, respectively. Basically three methods have been used to solve the BVP arising from simultaneous diffusion and reaction. The finite-difference (Lee, 1968; Deckwer, 1976) and the multiple shooting (Darde et al., 1983) methods were previously the most common algorithms. The latter approach is not a very recom769
mendable because it suffers from stability problems. Newer methods applied in this context are based on collocation (Ho and Lee, 1980; Huang et al., 1980; Shah and Paraskos, 1975; Nielsen and Villadsen, 1986), which has been found to be quite efficient in solving BVPs. In the current work both collocation and finite-difference methods are tested. The NLE system arising from thi bulk phase mass balances is a typical problem in chemical engineering computations. Thus numerous NLE solvers are available for this problem. In this work the Broyden algorithm including constraints for the variables is used. The benefits of the Broyden method are the speed and the low computer memory requirement.
2. STOICZ-IIOMETRY
AND
KINETICS
The general case of S reactions and N components in a chemical system can in most cases be described by the power law kinetics:
R,=k,f&F’, j=l,..., i= 1
s.
(1)
Using the stoichiometric coefficient matrix y, the rates of the components can be written in matrix notation: r=y-R.
(2)
The reaction rate constants are assumed to follow the
Arrhenius temperature dependence: k,=k,OexpG.
J. J.
770 3. SIMULTANEOUS REACTION AND IN STAGNANT FILMS
ROMANAINEN
DIFFUSION
Assuming that the diffusion through a volume element consists only of the molecular contribution, Fick’s law can be used. For a component diffusing through stationary boundaries the corresponding flux (No is given by equation (4):
Provided that the diffusion proceeds in one direction only in a fluid with a constant density, the wellknown balance for simultaneous diffusion and reaction in steady state is obtained: Diz+r,=o.
It is possible to get an analytical solution of equation (5) if the reaction rate ri is a linear function the corresponding concentration ci solely, or no reaction occurs in the stagnant element (ri = 0). This restricts the analytical solutions to first-order irreversible kinetics. Other, especially second-order bimolecular kinetics, have also been solved analytically, applying some additional assumptions or simplifications (e.g. van Krevelen and Hoftijzer, 1948). The film theory is commonly used to model the stagnant elements between the bulk phases and the gas-liquid interface. Film theory is used because of its simplicity and ease of implementation. Other theories, surface renewal (Dankwerts, 1970) and penetration theory (Higbie, 1935), have been developed, but they have not proved to be superior compared to the simple film theory. Therefore, the film theory is still used in most numerical studies concerning gas-liquid reactions. Applying the film theory for the stagnant elements, the following boundary conditions are valid: ci=c;
at
2 = 0,
ci=cF
at
2 -6,
(6)
where S is the thickness of the laminar layer. The gas-liquid interface is noted by z = 0 which sets z = 6 at the bulk phase boundary. This treatment is referred to as the film theory and it is commonly used to model the behavior of the stagnant elements in steady-state reactors. A schematic picture of the gas-liquid contact according to the film theory is given in Fig. 1. For numerical reasons the laminar layer spatial variable z is often scaled by the film thickness 6. The dimensionless variable x = z/6 is always within the interval [0, 11. Usually the chemical reaction takes place in only one of the phases. Here we assume that the reactions occur in the liquid phase, which is also the most common situation in industrial cases. Thus physical diffusion prevails in the gas film. Integrating equation (5) twice (ri= 0) and applying the boundary conditions (6) the flux can be obtained from equation (7): N&, = &(c!$;
- c;,,),
Since only physical diffusion is present in the gas phase, the fluxes equal on the film boundaries: N& = N$. Furthermore, the fluxes on gas-liquid interface are equal:
(9) both
N;, = -_k.
Llquld
Film
Bulk
b
Cl C
reaction
dlffuslon
6e
end
dfffusf on
0
sides
of
the
dCLi = N&, = kGi(c& - K&,), (10) L’( dx > x=0 where K, is the equilibrium constant between the component concentrations on both sides of the phase boundary and is related to the component vaporliquid equilibrium (VLE) constant or Henry constant if Henry’s law is used to describe the VLE. Equation (10) can be written in a more suitable form for
=a
\
(7)
where kGi is the gas-phase mass-transfer coefficient given by the film theory:
Liquid
Gas Film
Gas Bulk
and T. SALMI
%
Fig. 1. Schematicpicture of the gas-liquid contact according to the film theory.
b II
Solving gas-liquid reactormodels-1
f
T~_r
.
1 nPGi
771
A
.
Nb
n
a
L1
11 quid
*
4
c
n
FLi
I
film i
GAS BULK
6L
LIQUID BULK THE FILMS Fig. 2. Schematicpictureof the steady state CSTR.
numerical treatment, and the boundary conditions are thus given by:
- kLi dc,
--
k 01
d-x
-
c:, + K$,
= 0
CLi
czi
=
at
x = 0,
at
x=
1.
5. I. The sequential solution strategy (11)
As a measure of the accelerating effect of the chemical reaction on the absorption rate, the enhancement factor (E,) is defined as: Ei =
Ki k,(&
-
ct)
’
(12)
where the denominator gives the rate of the pure physical mass transfer. 4.
REACTOR
BULK
PHASE
BALANCES
Because the essential part of the current work concentrates on the numerical aspects of the problem, isothermal models were assumed to be sufficient in order to show the behavior of the numerical solutions. The stirred reactor consists of one gas and one liquid phase including one input and one output stream of both phases to and from the tank (Fig. 2). As mentioned previously, the reactions are restricted only to the liquid phase. Following these assumptions, the mass balances for the liquid and gas phases are given by equations (13) and (14), respectively: nFLi
-
. flPLi %GI
+Nb,fIt-(l-~~)V~r~=O, -
’ %x3r
5. NUMERICAL
-N&A=O.
the reactor bulk phase mass balances (13) and (14). The complete model is thus a combination of a two-point BVP and an NLE.
There are two possible strategies to solve this system of BVP and NLE. The traditionally used approach is to solve the film equations separately as a two-point BVP during the main iteration of the reactor mass balances. This leads to a sequential type of solution illustrated in Fig. 3. This traditional approach demands a very heavy computational effort, if the film equations must be solved every time the bulk phase model is called. A flag can be used based on the difference of the values of the variables since the previous call. This restricts the need to solve the BVP only, when the concentrations in the bulk phase have changed enough by the NLE solver in the main iteration, thus speeding the computation considerably. The changes resulting from this modification are shown in Fig. 4. The resulting main iteration loop (NLE) was solved by the Broyden algorithm (Broyden, 1965) and for the underlying flux update several BVP modifications were tested. These BVP algorithms included among others the COLSYS (Ascher et al., 1981) software and a new implementation of the
(13) (14)
APPROACH
The complete model of the steady-state CSTR consists of the combination of both the film equations (5) with the boundary conditions (11) together with
Fig. 3. The sequentialstrategyof solvingthe SCSTR model.
J. J. ROMANAINENand T.
772
Have variables changed enough from last call 7
YES
I Solve the Stagnant film model
( I
NO
SALMI
The boundary conditions are written partially inside the HSR solver. The nonvolatile components will be calculated using equation (19) for the second derivative at the gas-liquid interface. This method is called fake bounduries (Ames, 1977). The solution at x = 0 is approximated by a parabola having a maximum or a minimum at x = 0. Because of symmetry cLi(x_,) = +(x,), equation (18) reduces to: 2[cL&)
Fig. 4. The sequential strategy to solve the SCSTR model with a modification leading to less frequent stagnant film solutions.
finite-difference type of method HSR (Romanainen, 1989). Also the multiple shooting (DBVPMS) and finite difference (DBVPFD) algorithms provided by the IMSL (1987) for solution of BVPs were tested. Because of unreliability these codes were, however, discarded after preliminary tests. The COLSYS software applies orthogonal collocation on splines to approximate the solution of the BVP. The COLSYS requires the BVP model to be written in the following way for the second derivative:
f;(x)=$ -&
(1%
The boundary condition routine of COLSYS is divided into two logical sections according to the boundaries (x = 0, x = 1). At the gas-liquid interface the boundary conditions are given by equation (16): &(O)
=
_k,i dci k
dx
Or (1
_
+mLi-&i
(16)
Ax=-,
interface (x = 1) by equation n dp =
(17)
The initial guess for the component concentration profiles in the film was supplied by a linear function between the boundary concentrations. The Jacobian matrix was evaluated numerically for A(x) and analytically for both boundaries gi(0) and gj(l). In a finite-difference method using an even grid of discretization points, the second derivatives are normally evaluated with equation (18). The new implementation of a finite-difference type of solver (HSR) (Romanainen, 1989) is based on the idea of evaluating the second derivatives with an equation similar to (18), but at the same time using a refined grid. In addition, this was achieved without the need to add the number of discretization points in the grid: =CLi(XI--I)-2CLi(Xf)+CLI(Xl+1)~
X= X,
CLAXI)
-
Ax2
CLi(%)
9
(20)
where
x-0
&(l) = CL1- & = 0.
(19)
Equation (19) sets the gradient to zero at the gasliquid interface and thus no volatility can occur with the component in question. Now there are two possible approaches to the remaining boundary conditions. One is the dynamic approach in which the film concentrations are assumed time dependent, thus leading after discretization to a set of ordinary differential equations (ODES) instead of the original BVP [equation (S)]. This homotopy approach has very nice features including that the solution will always be within physically meaningful limits. The essential features of the method have been studied elsewhere (Haario and Seidman, 1988, 1989). The other possibility is to use the stationary approach in which the equations in the discretization points form an NLE system and solve this with an NLE solver. The problem with the HSR solver described above is the discretization. The even discretization inevitably implies that the gradient is very sensitive to the number of points used because the gradient is calculated by equation (20): Ax
= 0 and at the film-bulk (17):
- CLic%)l Ax= .
(18)
1 nLip - 1
number of discretization points,
(21)
In practical cases, ndpis between 11 and 41, giving Ax-values from 0.1 to 0.025. For very fast reactions, however, this will not be sufficient to compute the gradient accurately. To overcome this problem, there are basically two possibilities. First it is possible to use some kind of interpolation to calculate the gradient more accurately near the boundaries. Another choice is to make the discretization differently, e.g. by coordinate transformation (Versteeg et al., 1987). In this work the following discretization scheme was adapted. The scheme is based on dividing the remaining film into two equal parts as we approach each boundary (Fig. 5). The benefits of this discretization scheme are, on one hand, the accuracy near the boundaries and, on the other hand, the simplicity of the implementation. The accuracy is achieved because the discretization interval at the boundaries is now:
Solving gas-liquid reactormodels-l
I
III I 10 . . .
l/8
114
773
I
I
l/2
I III 7/a
3/b
. . .
1
X
Fig. 5. The half discretizationmethod. Ax=$
n,=-.
In these cases the integration of the film equations would have to be continued for a very long time causing the method to be impractical.
x’
% -1
5.2. The simultaneous solution strategy
”
L
Equation (22) will give the discretization interval Ax practical values from 2-5 to 2-20, which guarantees sufficient accuracy to the gradients at both boundaries (x = 0, x = 1). Usually the gas-liquid interface boundary is more important, since the profiles of the soluble, reactive components there are sharp (especially with fast reactions). This is also the case in the examples presented in this paper. On some other cases, however, the reactive liquid components will get sharp profiles near the bulk phase boundary because of fast and/or complex kinetics in the film. The half discretization scheme is thus using the most refined grid exactly where the greatest accuracy is required. The other benefit of the implementation is that this discretization allows us to use equation (18) for the second derivatives, since the point x,_, is always x0 and the point x,, , is the next point, when computing points x, from the center to the gas-liquid interface. The same applies for the other half of the film. The only drawback of this scheme is that the center of the film (x = 0.25, x = 0.5, x = 0.75) will be calculated always using Ax = 0.25, which may cause inaccuracy in the results. This can be avoided by combining the even and halving discretion schemes using the even scheme at the center of the film and the halving scheme near boundaries. Another method to improve the accuracy was also tested. It is a method commonly used to estimate the Auxes in catalyst particles (Froment and Bischoff, 1979) where simultaneous diffusion and reaction occur. This method is based on the total balances of the film elements:
Jo Estimating the integral on the right-hand side of equation (23) by the Simpson formula and using the other gradient that was near zero (Nkj in this case), corrected values for the gradients at the steeper end of the profile were obtained. Other methods such as splines on points near the boundaries were also tested, but they did not provide much improvement to the accuracy of the results. An additional problem with the dynamic approach of the HSR solver is the cases with slow reaction.
The alternative strategy to solve the film and bulk phase equations is to solve them simultaneously. Here the reactor model is converted to a boundary value problem. Equation (5) remains unchanged. The difference now is that the bulk phase mass balance equations (13) and (14) are included in the corresponding boundary conditions (16) and (17). The resulting BVP was solved by the collocation software COLSYS. The software requires the model equations to be written in the following way for the simultaneous strategy: (24) and for the film-bulk similarly:
g,(l) = -
2
boundary
x = 1 we
+ -cc;~~++~Li+(l.O-to)VRri
( >
&al/,
get
3 (25)
while equation (15) remains unchanged. The transformed concentrations CL, are obtained from: Cii(X) = & + c&r-(1 - x), where c,,, denotes concentration. 6.
the
maximum
(26) liquid
feed
EXAMPLES
All examples presented here were programmed in FORTRAN and executed on a DEC Micro VAX-II computer. 6. I. The sragnant jilm model solution Two different reaction schemes are used to test the different boundary value problem solvers. These examples deal only with the film equations since the reactor bulk balances are not needed in this section. 61.1. One first -or&r irreversible reaction. The reaction scheme is an irreversible first-order reaction where component A originally in the gas phase diffuses through the films to the liquid phase reacting simultaneously to a liquid component B: A(G) --+ B(L).
(27)
J. J. ROMANAINEN
774 Table
I.
Solution methods for the irrwrcrsible reaction scheme
first-order
NaaIC
Method
ANA CSN CST HSE HSH HDE HDH
The analytical solution COLSYS normal COLSYS with variable transformation HSR stationary with even discretization HSR stationary with half discretization HSR dynamic with even discretization HSR dynamic with half discretization
Table 2. Parameters used in simulation of the firstorder irreversible reaction scheme (example 6.1.1) Parameter K# k, kc, 6, ct cs*
Value 1.0 1.0 x lo-‘m s--l 1.O x 102’ m s- ’ ( = no gas film resistance) 1.0 x lo-‘m 0.0 mol dm-’ 0.01 mol dm-’
The solution methods studied for the film equations are listed in Table 1. The parameters used in simulations are presented in Table 2. To get a good coverage of the behavior of the solvers, a wide range of reaction rate constant values were used. Figure 6 shows the CPU-times demanded by the tested solvers relative to the analytical solution. The number of discretization points used in the HSR method was 21, which was found to be enough for reaction rate constant values less than 10s s-l.
and
T.
SALMI
Figure 6 shows that the stationary HSR method with either discretization scheme is the fastest overall solver in these cases. The COLSYS solver is also fast except for some individual values of the reaction rate constant (k, = lo’, l@, lo’, 10’s_‘). Thlb might indicate some problems connected to the implementation of the COLSYS software. Neither the discretization scheme in HSR methods nor the variable transformation in COLSYS seem to affect the computation time. The dynamic HSR methods required about four times the CPU-time compared with the stationary ones and they are thus discarded already at this point. The accuracy of the solution was examined by the enhancement factor, the value of which was compared to the analytical solution. The result of this comparison can be seen in Fig. 7. Despite the problems with COLSYS in computation time, the solution accuracy was found to follow the analytical solution extremely closely. The variable transformation did not effect the accuracy in any detail. The HSR methods, on the other hand, show a clear dependence on the discretization. The best version of the HSR solver-at least in this example--seems to be the stationary HSR solver with the half-discretization scheme. This scheme (HSH) is very accurate with 21 points up to reaction rate coefficient values k, = IO5s-l. The even discretization schemes with 21 points fail after values k, 3 10 s-l. There cannot be seen any differences in the accuracy of the stationary and the dynamic variations of the HSR solver.
Fig. 6. Relative computation time as a function of reactionrate coefficient.One first-orderirreversible reaction (example 4.1.1).
Solving gas-liquid reactormodels-l n.3
-A+ -t *ANA -A-
HDH HSH HSE HOE
Reaction
rate
coefficient
k,
Fig. 7. Enhancementfactor as a function of reactionrate coefficient.One first-order irreversiblereaction (example 6.1.1). The accuracy of the finite-difference methods (HSR) was further improved by the use of equation (23) and the Simpson formula. This approach was applied to the even discretization variations of the HSR solver (methods HSEm and HDEm). The enhancement factor values achieved with these modified routines are shown in Fig. 8. This method improved the enhancement factor values only up to rate coefficient values k, < lo3 s-‘. Therefore the even discretization versions were discarded at this point. The overshooting behavior of the enhancement factor in applying equation (23) is probably due to numerical problems. An examination of the concentration profiles in the film as function of the reaction rate constant values can be seen in Fig. 9. These profiles show the expected behavior of steepening curves as the reaction rate constant increases. The remaining solvers all follow the analytical solution quite nicely. 6.1.2. Two consecutive-competitive irreversible second-order reactions. This example was applied only to the remaining solvers which passed the tests in the first example. These solvers were the stationary HSR method with half discretization scheme (HSH) and the COLSYS method with the variable transformation (CST). The reaction scheme is presented in equation (28) and the corresponding parameters used in the simulation are listed in Table 3:
A(G) + B(L) -
C(L),
A(G) + C(L) +
D(L).
The simulations were performed with reaction rate coefficient values of lo-‘, 10’ and 103m3 kmol-‘s-l for both rate constants in every possible combination resulting in nine simulation runs. The resulting CPU-
775
time usages are depicted in Fig. 10, where it can be seen that COLSYS is again the fastest solver except for the last three cases where k, = lo3 m3 kmol-t s-l and k2 has all three values in turn. The accuracy could not be compared since an analytical solution is not obtainable. On the other hand, both methods seem to give similar results, so the accuracy is probably not a problem at least with these values of the rate constants. 6.1.3. Conclusions of the stagnant jilm model solu tion. COLSYS seems to be both the fastest and the most accurate method, except for some cases, where problems might occur. The variable transformation [equation (26)J did not seem to alter the behavior of COLSYS in any way. A few details of the COLSYS implementation were found to be worth noticing. The number of subintervals should always be set to one and let the software decide how to increase it. The problem should always be defined as nonregular although this might cause problems in some linear cases. The tolerances should be at least IO-*, otherwise the results will not be accurate. Finally it was noticed that giving the values from a previous failed calculation back to COLSYS as a new initial guess often lead to a visible solution. These observations lead us to believe that the problems encountered with COLSYS were mostly due to an inadequate initial guess of the film profiles. The problems arising from this can be minimized restricting the mesh size in COLSYS or using an estimate based on a previously successfully computed profile as an initial guess. These methods of improving the performance of COLSYS were found to be
lo:: Reaction
101
rate
105
106
coefficient
10'
ld
k,
Fig. 8. Enhancementfactor as a function of reaction rate coefficient. One first-order irreversible reaction. Finite-difference methods with even discretizstionscheme (HSE. HDE) and eauation(23) correction(HSEm. HDEml . corn&red td the an&ical solution (ANA). ’
J. J.
776
ROMANAINEN
and T. SALMI
(a) 0.010 0.009
-a-
CSN
0.008 0.W-f 0.006 0.005 0.004 0.003 0.002 0.001 0.000 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Dimensionless
film
thickness
,
x
(c)
(b)
Dimensionless
film
thickness
x
Dimensionless
film
thickness
x
Fig. 9. Component concentrationprofiles in the liquid film with differentvalues of the reaction rate coefficient(example 6.1.1): (a) k, = 1.0 x IO”s-‘; (b) k, = 1.0 x 102s-‘; and (c) k, = 1.0 x 104s-‘. successful; hence it was chosen as the stagnant film model solver. The HSR routines were still under development during this work, so they are not yet comparable with established software like COLSYS. However, Table 3. Parameters used in simulation of the two consecutive-competitive reaction schemes (example 6. I .2) PaIXmeter
4
b
k
V&E 1.0
= km OA
b,
=
km
1.0x
IO-‘m
s-’
1ti5 m s- ’ ( = no gas film resistance) I.0 x 10mSm 0.0mol dm-’ 0.01 mol dme3 1.O mol dm-’ 0.0 mol dm-’ 1 .O x
after some further development, they might prove to be competitive as stagnant film solvers. The half discretization scheme is obviously superior to the even one in terms of the result accuracy. The computation time is about the same with both methods. The accuracy is, however, not so high in the center of the film, and this might cause problems with some more complicated cases, where steep gradients are expected not only at the boundaries, but also at the center of the film. As discussed earlier, this can be helped combining both discretization schemes. Another way would be to develop some kind of adaptive scheme for the discretization during the computation, based, e.g. on the change of the derivatives. The dynamic HSR methods offer the benefit of approaching the solution always via a physically
777
Solving gas-liquid reactor models-1
k,=lO-’
k,=lO-’
k,=lO+
k,=lO
k+= IO-’
k2= IO’
k2= lo3
k*=lO-’
React
’
ion
k,=
10’
k2=70’ rate
k,=
IO’
kz= IO3
k,= lo=
k,=103
k,:
kglO-’
k,=lO’
k;lO”
IO3
coefficients
Fig. 10. Relative computation time as function of reaction rate coefficients. Two consecutive second-order reactions(example6.1.2).
meaningful path. The problem with these dynamic methods applied in stationary cases is the selection of the integration interval. When then reaction is a fast one, the choice of the integration interval is not critical, but with very slow reactions the integration would have to be continued for a very long time to obtain the steady state. This excludes this kind of methods from practical applications, although they might be appropriate for some special cases with difficult reaction-diffusion schemes. The stationary HSR method with the half discretization was the only one comparable to COLSYS. The benefit of this method over COLSYS is its easy implementation. The advantage of the physically viable path to the solution is regrettably lost in the stationary version. The final conclusion of the stagnant film solvers is that COLSYS seems to offer the best performance at the present moment. The stationary HSR model is the other possible choice, but it should be developed further.
6.2. The steady-state
CSTR
model solution
These latter examples 6.2.1 and 6.2.2 use variations of the two consecutive-competitive second-order reactions scheme [equation (2811. The examples are used to demonstrate the different characteristics of the sequential and simultaneous strategies introduced in connection with the steady state CSTR model. The example runs were obtained using COLSYS as the BVP solver.
6.2.1. Two consecutive-competitive second-order irreversible reactions. The reaction scheme used in these examples [equation (2911 was the same as applied by Darde and co-workers (Darde et al., 1983): A(G) + B(L) - ’
W-1 + G(G),
A(G) + C(L) - ’
W-1 + G(G),
(29)
where letters A-G correspond to the compounds listed in Table 4. The parameters used for the simulation were taken partly from Darde et al. (1983) and are listed in Table 5. Furthermore, all component properties were taken at T = 273 K and p = 0.1 MPa, for the vapor-liquid equilibrium Henry’s law was used and the gas phase was assumed to be ideal. The liquid phase mass transfer coefficients were set to 10m5m s-l for those components in the liquid phase whose real kLj values were not available. Tabk 4. Components in the example reaction [equation(29)) used in the simuktion of the steady stateCSTR model: chlorinationof ~aru-cresol
A
B C D E
ChlOCiIlC Pam-crrsol
Monochloro-para-c.resol Dichloro-porn-cresol Tetrachloromethane Nitrogen Hydrogenchloride
Gas Liquid Liquid Liquid Liquid GaS Gas
J. J. ROMANAINEN
778
and T. SALMI
Table 5. Parameters used in simulation of the steady state CSTR model. Taken partly from Darde ef al.
----
(1983)
of components Number of reactions KA kl k, k,, kLB kLc Number
k,, (i= A, B, C, D, E, F, G)
Simultaneous strata Seq”cntiol strategy
7 2 0.020 5.620 dm” mol-’ SK’* 0.107dm’mol-‘s-‘P 7.00 x 10mSm s-la 3.77 x lo-Srn S-In 3.50 x 10m5m s-la 1.00 x lo*sms-’ ( = no gas resistance) 1.486 x 10-5m’ 6.8714 mz rne3’ 0.2 1 kmol h-’ 40 kmol h-l 205 kmol h-’ 1700 kmol h-’
‘From Darde et al. (1983). Some values modified.
0
The rate of the reaction with respect to mass transfer can be described with the Hatta number. The Hatta number for the first reaction {equation (29)] corresponds to an intermediate or slow reaction (< 1.22) achieving the maximum value at the feed composition. The Hatta number for the second reaction is even smaller (c 0.3) indicating a slow reaction. The actual simulations were done in a sequence of an increasing liquid residence time in the reactor. This enabled the use of a good initial guess from the solution at the previously computed residence time. The results from the simulation sequence is
Residence
time
-r (s)
12. Component concentrations in the liquid bulk phase at different liquid residence times. Simulation of the steady state CSTR model in example 6.2.1. Irreversible reactions. Fig.
presented in Fig. 11, where the computation time is plotted against the liquid residence time. The obvious result from this is that the simultaneous strategy is considerably faster than the sequential one. The simultaneous strategy, however, displays some disturbing behavior at some residence times. This behavior can be seen in Fig. 12, where at z = 1260 s and
25C
B
Simultaneous
[77721
Sequential
strategy strategy
Residence
Fig.
11. Relative
computation
time
T
ts)
time of the steady state CSTR model with different liquid residence times (example 6.2.1). Irreversible reactions.
Solving gas-liquid reactor models-l z = 1440 s the correct solution can be achieved only after re-introducing values from an originally unsuccessful calculation and recalculating the simulation. Otherwise both strategies give comparable results. The physically interesting result given in Fig. 12 is the maximum in the concentration of the intermediate component C, which in this case happens to be the desired product. This simulation would indicate the suitable liquid residence time of about r = 1080 s for the operation of the actual reactor. 6.2.2. Two consecutive-competitive second-order reversible reactions. The steady state CSTR simulations were continued by modifying the original reactions [equation (29)J to include some reversibility. This resulted in schemes (30): A(G) + B(L) - ’ A(G)
+ C(L) +
C(L) + G(G), D(L) + G(G).
(30)
and (31): A(G) + B(L) --!A(G) + C(L) +
4
C(L) + G(G). D(L) + G(G).
(31)
These reaction sequences are hypothetical, but the physical properties from Table 5 were used for the sake of comparison. Here the simulations were performed using 0.1, 1 and 10 times the values of
the reaction rate coefficients of reactions 4 and 3 in respect to the values of reactions 1 and 2, respectively. The computation time with different reactor residence times is shown in Fig. 13 for reaction scheme (30). Since the computation time is approximately equal for both strategies, this indicates that the reversibility of reactions causes problems for the simultaneous solver. The computation time is especially high at the end of the simulation sequence which counts for long residence times in the reactor. Plotting the component concentrations at various residence times (Fig. 14) shows some instability of the simultaneous solver at higher residence times. In the last simulation sequence also the first reaction was made reversible [equation (3111 to verify the assumption of reversibility affecting the performance of the simultaneous solver. This resulted in a clear proof of the assumption, since the simultaneous solver was not able to solve more than few first points with low residence times, after which it failed. Recalculating the failed cases successively by returning the failed result to COLSYS did not provide the solution anymore as it did in the previous examples. The time used by the sequential solver was about the same as in the other reaction schemes tested. The resulting component concentration profiles as a function of the residence time are plotted in Fig. 15 for one example simulation, where k3 = 10 * k, and k4=
10.kl.
The difficulties with the simultaneous solver were first assumed to originate from a bad initial guess of
8OC Simultaneous
pf7
Sequential
strategy
strategy
6OC
;;
z i
779
4oc
E r 2 ” 2oc
0
Fig. 13. Relative computation time of the steady-state CSTR model with different liquid residence times (example 6.2.2). Irreversible and reversible reactions.
J. J. ROMANAINEN
780
and T. SALMI 0.25
0.25 n^
n^ E $
E $
0.20
.A.
.s s .-
E ._
0.15
CD
0.15
P E iti 5 0.10 0
5 E : 8
0.20
0.10
E E g E ;:
E fi n 0.05 E :: 500
1000
ResTdence
1500 time
2000 -r
2500
:
3
(s)
C*
CC 0.05
ca 0
500
1000
Residence
1500 time
2000 -r
2500
3
(s)
Fig. 14. Component concentrations in the liquid bulk phase at different liquid residence times. Simulation of the steady state CSTR model in example 6.2.2. Irreversible and reversible reactions.
Fig. 15. Component concentrations in the liquid bulk phase at different liquid residence times. Simulation of the steady state CSTR model in example 6.2.2. Reversible reactions with k3 = 10. k2 and k, = 10. k, in reaction (31).
the profiles. However, this is not the only source of problems, since the behavior of the simultaneous solver was predictible and got worse, when adding reversible reactions into the scheme. The simultaneous solver first started to deviate from the solution slightly and then failed. At first the failure was corrected by the successive guessing of the profile. When the residence time was still increased, this lead the whole method to fail. At first some of the concentrations became negative, the successive substitution found a physically nonfeasible solution and always converged to it, thus failing completely in computing the rest of the sequence. Because of this behavior, some modifications to the tolerances were tested with COLSYS. Making the tolerances lower resulted the method to be faster, and able to compute a few more simulations. However, this did not help very much, especially with the last scheme [equation (3111. In further test simulations, also different mesh size restrictions were tried, but without results. 6.2.3. Conclusion of the steady state CSTR model solution. The conclusion of the simultations were that the implementation of the simultaneous strategy was not a reliable one for an environment, where different types of reaction schemes are to be simulated, especially in the case of reversible reactions. In restricted use with a known and tested reaction scheme, however, it was sometimes superior to the sequential solution strategy. The problems with the simultaneous solution strategy were mostly connected with the problems of COLSYS software, and should not be taken as principal defects of the approach. The unstable behavior is probably connected with the boundary conditions,
and their very nonlinear nature, which is, furthermore, changing during the simulation. The sequential solution strategy implemented in this work using the Broyden method for the bulk phase mass balances and COLSYS for the stagnant film model equations performed well. The problems found with the COLSYS software were overcome by the NLE software, thus resulting in a reliable combination to solve steady-state CSTR models. Further studies should thus be directed into other boundary value problem solvers as well as to some other updating methods of the bulk concentrations and the general behavior of very nonlinear boundary conditions. The COLSYS software seems a very good BVP solver provided that the necessary modifications concerning very nonlinear boundary conditions could be added. The CSTR model can also be used to describe systems with more complex hydrodynamics (plug-flow, axial dispersion) by introducing the tanks-in-series concept. The numerical methods studied here can in principle be applied to these cases. Since the use of the tanks-in-series concept is restricted to cocurrent columns, the actual column simulation should be done with a dedicated model suitable for both cocurrent and countercurrent columns. NOMENCLATURE A = Interfacial area in reactor (m2)
0 = A/V, (m* m-‘) ci = Component concentration (km01 m-‘) c; = Transformed component concentration D, = Diffusion coefficient (m* s-l) Er I Enhancement factor
(kmol m-‘1
Solving gas-liquid reactor models-l E, = Activation energy of a reaction (kJ mol-‘)
f = Model function
g = Boundary condition K, = Gas-liquid equilbrium constant k, = Reaction rate coefficient; the dimension depends on kinetics, for n&order kinetics [(m’ kmol-I)(“- I) s-i] k,? = Arrhenius law frequency factor (same as k,) kLi = Liquid mass transfer coefficient (m s-‘) k, = Gas mass transfer coefficient (m s-‘) iV = Number of components N, = Flux (km01 rnT2 s-‘) n = Amount of substance (mol) ri = Flow of substance (mol s-‘) ” dp= Number of discretization points p = Pressure (MPa) R = Gas constant (8.314) (J k-l mall’) R = Vector of R,; the length of R is (S) (same as Rt) R, = Rate of reaction (km01 m-r s- ‘) r = Vector of r,; the length of r is (N) (same as r,) ri = Component reaction rate (kmol rno3 s-l) S = Number of reactions T = Temperature (K) V = Reactor volume (m3) v = Volumetric flowrate (m3 s-i) y = (N . S) Matrix of \I~ x = Film coordinate, drmensionless z = Film coordinate (m) Greek
symbols
a- = Component
$ = co = z= Y#=
order in a reaction Film thickness (m) Gas hold-up Residence time in a reactor (s) Stoichiometric coefficient
781
work was carried out as a part of the Computer Aided Reactor Design project financed by Teknologian kehittimiskeskus (TEKES), Kemira Oy, Neste Oy, University of Oulu and Abe Akademi.
Acknowledgement-This
REFERENCES Ames W.
F., Numerical Methods for Partial D#erentiaI 2nd Edn. Academic Press, New York
Equations,
(1977). Ascher U., J. Christiansen and R. D. Russell, ACM Trans. Math. Sofftware 7, 209 (1981). Broyden C. G., Math. Comput. 19, 577 (1965). Danckwerts P. V., Gas-Liquid Reactions. McGraw-Hill, New York (1970). Darde T., N. Midoux and J.-S. Charpentier, Entropic 19,92 (1983). Deckwer W. D., Chem. Engng Sci. 31, 309 (1976). Doraiswamy L. K. and M. M. Sharma, Heterogeneous Reactions, Vol. 2. Wiley, New York (1984). Froment G. F. and K. B. Bischoff, Chemical Reactor Analysis and Design. Wiley, New York (1979). Haario H. and T. I. Seidman, Proc. Minisymposium on Numerical
Methodr
for
Semiconductors
and
Magnets
(Neittaanmgki, Ed.). University of JyviiskylL, Dept of Mathematics, Report 42, JyvLskyll (1988). Haario H. and T. I. Seidman, Reaction and diffusion at a gas/liquid interface. Submitted for publication (1990). Higbie R., Trans. Am. Inst. Chem. Engnrs 35, 365 (1935).
Ho S.-P. and S.-T. Lee, Chem. Engng Sci. 35, 1139 (1980). Huang D. T.-J., J. J. Carberry and A. Varma, AIChE JI 26, 832 (1980).
Subscripts
F = G = i= i = I= L = P= R = A B
. i;-1
Reactor feed Gas phase Component Reaction Discretization point in the film Liquid phase Reactor outlet (product) Reactor
= Components A, B, . . .
Superscripts
b = Bulk phase s = Gas-liquid interface
IMSL, Znternational Mathematical and Statistical Libraries, Vol. 2, Houston, Texas (1987). van Krevelen D. W. and P. J. Hoftijzer, Rec. Trav. Chim. Pays-Bas 67, 563 (t948). Lee E. S., Quasilinearization
and
Invariant
ZmbeaWng.
Academic Press, New York (1968). Nielsen P. H. and J. Villadsen, Chem. Engng Sci. 41, 1655 (1986).
Rase H. F., Chemical Reacror Design for Process Plants, Vol. 1. Wiley, New York (1977). Romanainen J., Numerical strategies in solving gas-liquid tank reactor models, Lit. Technol. Thesis, Helsinki University of Technology, Helsinki (1989). Shah Y. T. and J. A. Paraskos, Chem. Engng Sci. 30, 465 (1975).
Versteeg G. F., J. A. M. Kuipers, F. P. H. van Beckum and W. P. M. van Swaaij, Chem. Engng Sci. 44, 2295 (1989).