Chemical Engineering Journal 76 (2000) 223–237

Modelling of gas–liquid reactors — implementation of the penetration model in

dynamic modelling of gas–liquid processes with the presence of a liquid bulk
E.P. van Elk a,∗ , P.C. Borman b , J.A.M. Kuipers c , G.F. Versteeg c
a Procede Twente BV, P.O. Box 217, 7500 AE Enschede, The Netherlands
b DSM Research, P.O. Box 18, 6160 MD Geleen, The Netherlands
c Department of Chemical Engineering, Twente University of Technology, P.O. Box 217, 7500 AE Enschede, The Netherlands
Received 14 April 1999; received in revised form 8 November 1999; accepted 29 November 1999

Abstract
Analytical solution of models for gas–liquid reactors is restricted to a few asymptotic cases, while most numerical models make use of
the physically less realistic stagnant film model. A model was developed that simulates the dynamic behaviour of gas–liquid tank reactors
by simultaneously solving the Higbie penetration model for the phenomenon of mass transfer accompanied by chemical reaction and the
dynamic gas and liquid phase component balances. The model makes it possible to implement an alternative for the well known Hinterland
concept, which is usually used together with the stagnant film model. In contrast to many other numerical and analytical models the present
model can be used for a wide range of conditions, the entire range of Hatta numbers, (semi-)batch reactors, multiple complex reactions
and equilibrium reactions, components with different diffusion coefficients and also for systems with more than one gas phase component.
By comparing the model results with analytical asymptotic solutions it was concluded that the model predicts the dynamic behaviour of
the reactor satisfactorily. It is shown that under some circumstances substantial differences exist between the exact numerical and existing
approximate results. It is also shown that for some special cases, differences can exist between the results obtained using the stagnant film
model with Hinterland concept and the implementation of the Higbie penetration model. ©2000 Elsevier Science S.A. All rights reserved.
Keywords: Gas–liquid reactors; Mass transfer; Dynamic penetration model

1. Introduction the overall gas and liquid phases, and the micro model, de-
scribing the gas–liquid interphase transport of mass and/or
Many chemical processes involve mass transfer of one heat in combination with the chemical reaction.
or more gaseous components to a liquid phase in which Often applied, idealised, macro models are the plug flow
a chemical reaction occurs. Typical examples include gas (PFR) model and the ideal stirred tank (CISTR) model. Other
purification, oxidation, chlorination, hydrogenation [1] and possibilities include the plug flow with axial dispersion and
hydroformylation processes. the tanks in series model respectively (see [3]).
For the selection of a reactor type it is important to clas- Frequently used micro models are the stagnant film model
sify gas–liquid transfer processes on the magnitude of the in which mass transfer is postulated to proceed via stationary
reaction-diffusion modulus (Hatta number [2]). The ‘fast re- molecular diffusion in a stagnant film of thickness δ [4],
actions’ (Ha>2) are considered to proceed predominantly the penetration model in which the residence time θ of a
near the gas–liquid interface, while the ‘slow reactions’ fluid element at the interface is the characteristic parameter
(Ha<0.2) are considered to occur mainly in the liquid bulk. [5], the surface renewal model in which a probability of
For reactions with 0.2<Ha<2, that are also frequently en- replacement is introduced [6] and the film-penetration model
countered in the chemical process industry, no distinct reac- which is a two-parameter model combining the stagnant film
tion region can be defined. and the penetration model [7,8].
In the general case the description of a gas–liquid reactor Solution of these micro models for mass transfer accom-
consists of the following parts: the macro model, describing panied by chemical reaction analytically is restricted to cases
in which many simplifying assumptions are made, e.g. reac-
∗ Corresponding author. Tel.: +31-53-4894480/+31-53-4894337; fax: tion kinetics are simple and the rate of the reaction is either
+31-53-4894774. very fast or very slow [9–12]. Other authors have devel-
E-mail address: [Link]@[Link] (E.P. van Elk). oped approximate analytical solutions [13–16]. For all other

1385-8947/00/$ – see front matter ©2000 Elsevier Science S.A. All rights reserved.
PII: S 1 3 8 5 - 8 9 4 7 ( 9 9 ) 0 0 1 6 2 - X
224 E.P. van Elk et al. / Chemical Engineering Journal 76 (2000) 223–237

situations numerical techniques are required for solving the zero. Other reactions, including multiple (in)dependent re-
coupled mass balances. actions, or other kinetics, like Langmuir–Hinshelwood, can
Overall dynamic gas–liquid reactor models can be solved easily be included in the present model. For simplicity rea-
in two fundamentally different ways: sequential and simul- sons only one reaction and only one gas phase component
taneous solution of the micro and macro model, respectively. is considered in the current work. Also for simplicity rea-
If the penetration model is used as micro model and the sons, isothermal models were assumed, the temperature can
gas and or liquid bulk concentrations may vary significantly however be implemented as additional component via im-
during the residence time of the liquid element at the mass plementation of the energy balance.
transfer interface (e.g. the liquid bulk is not at steady state The mass transfer in the gas phase is described with the
or at equilibrium), the simultaneous strategy constitutes the stagnant film model while for the liquid phase the penetra-
only correct approach. tion model is used. Further postulations in the actual dy-
Due to its simplicity and relative ease of implementa- namic reactor model are: (1) both the gas and the liquid bulk
tion in the overall gas–liquid reactor model (e.g. using a can be assumed to be a CISTR; (2) the reaction only takes
numerical analogy of the Hinterland concept [3]), the film place in the liquid phase; (3) the contact time is substantially
theory is often used in numerical studies [17,18]. Both these smaller than the liquid phase residence time.
articles deal with steady state models only. Romainen and
Salmi [19,20] have published a paper in which a dynamic 2.2. Micro model
film-penetration theory is adopted.
The film theory and Hinterland concept can, in many cases For the penetration model the balances for each species
be applied successful. However, the penetration model and for the phenomenon mass transfer followed by a chemical
the surface renewal model are to our opinion preferred since reaction yields the following set of equations:
they are physically more realistic models [3,21]. Since it is
not possible to implement the Hinterland analogy directly for ∂[A] ∂ 2 [A]
these models, a dynamic reactor model is required that uses = Da − Ra (3)
∂t ∂x 2
a somewhat different approach. Until now no such model
has been published in literature. ∂[B] ∂ 2 [B]
= Db − γ b Ra (4)
In this paper, a dynamic gas–liquid reactor model is ∂t ∂x 2
presented which simultaneously solves the Higbie penetra-
∂[C] ∂ 2 [C]
tion model and an instationary CISTR for the gas and the = Dc + γc Ra (5)
liquid phase, respectively. The model is suited for the ∂t ∂x 2
phenomenon mass transfer with complex (ir)reversible ∂[D] ∂ 2 [D]
chemical reaction(s) and is valid for the entire range of = Dd + γd Ra (6)
∂t ∂x 2
Hatta numbers. This paper concentrates on describing and
validating the model. A more practical application of the To be solved uniquely the four non-linear partial differen-
model for analysing the dynamics of gas–liquid processes tial equations (3)–(6) require one initial and two boundary
is given elsewhere [22]. conditions, respectively. The initial condition is given by

t = 0 and x ≥ 0, [A] = [A]l,bulk , [B] = [B]l,bulk ,

2. Theory [C] = [C]l,bulk , [D] = [D]l,bulk (7)

2.1. Introduction The boundary condition for x=δ p is given by

The problem considered is a dynamic gas–liquid reactor t > 0 and x = δp , [A] = [A]l,bulk , [B] = [B]l,bulk ,
with mass transfer followed by an (ir)reversible chemical [C] = [C]l,bulk , [D] = [D]l,bulk (8)
reaction of general order with respect to both reactants and
products: where the thickness δ p of the liquid element is assumed to
be infinite with respect to the penetration depth of the gas
A(g) + γ b B(l)γc C(l) + γd D(l) (1) phase component. The concentrations of the liquid bulk used
with the following overall reaction rate equation: in Eqs. (7) and (8) follow from the macro model for the
liquid bulk.
Ra = kR,m,n,p,q [A]m [B]n [C]p [D]q The second associated boundary condition is obtained by
−kR,r,s,t,v [A]r [B]s [C]t [D]v (2) assuming that the species B, C and D are non-volatile and
that the flux of component A from the gas phase is equal
The reaction rate expression is based on order-kinetics. If to the flux of component A to the liquid phase. The use of
p, q, r and s are zero, a reversible 1,1-reaction is obtained. An the latter assumption instead of assuming that [A]=[A]l,i at
irreversible 1,1-reaction is obtained by putting kR ,r,s,t,v to x=0 is convenient in view of the applicability of the model
 E.P. van Elk et al. / Chemical Engineering Journal 76 (2000) 223–237 225

for cases where a part of the resistance against mass transfer up with the liquid bulk and replaced by a new fresh one.
is situated in the gas phase: The dimensions of the liquid element are assumed to be
  infinite compared to the penetration depth and therefore no
∂[A]
−Da = kg ([A]g,bulk − [A]g,i ), direct mass transport to the liquid bulk via the liquid element
∂x x=0 occurs.
     
∂[B] ∂[C] ∂[D] The mixing with the liquid bulk must be taken into account
= = = 0 (9) instantaneously after the contact time θ . It has been assumed
∂x x=0 ∂x x=0 ∂x x=0
that the contact time is much smaller than the liquid phase
2.3. Macro model residence time. Therefore the convection may also be taken
into account instantaneously. The liquid phase concentration
The generally applicable material balance of a component after these contributions is thus given by (i=A, B, C, D)
in a certain phase in a tank reactor is: Ni,elem + Ni,bulk + (8l,in [i]l,in − 8l,out [i]l,out )θ
[i]l = (18)
haccumulatedi = + hreactor feedi εl VR
− hreactor outleti where the first two terms of the numerator present the num-
(10)
+ hinterfacial mass transferi ber of moles present in the total liquid phase after the contact
+ hreaction in the bulk phasei time (see Fig. 1)
For the gas phase only component A has to be considered: Ni,elem + Ni,bulk = [i]bulk εl VR
Z δp
d[A]g
εg VR = 8g,in [A]g,in − 8g,out [A]g − Ja aVR (11) + ([i] − [i]bulk ) dx aVR (19)
dt 0

The liquid phase bulk component balances are given by and the last two terms of the numerator of Eq. (18) present
the convection in the liquid phase.
d[A]l The macro model Eqs. (11)–(15) to be used as boundary
ε l VR = 8l,in [A]l,in − 8l,out [A]l
dt conditions for the micro model need to be simplified to
+Ja aVR − Ra εl VR (12)
d[A]g,bulk [A]g,in [A]g,bulk Ja a
d[B]l = − − (20)
ε l VR = 8l,in [B]l,in − 8l,out [B]l − γ b Ra εl VR (13) dt τg,in τg,out εg
dt
d[A]l,bulk
d[C]l = −Ra,bulk (21)
ε l VR = 8l,in [C]l,in − 8l,out [C]l + γc Ra εl VR (14) dt
dt
d[B]l,bulk
d[D]l = −γ b Ra,bulk (22)
ε l VR = 8l,in [D]l,in − 8l,out [D]l + γd Ra εl VR (15) dt
dt d[C]l,bulk
= γc Ra,bulk (23)
where the mass transfer flux, Ja , is defined as: dt
Ja = (kov )kg ,kl ,Ea ,ma (ma [A]g,bulk − [A]l,bulk ) (16) d[D]l,bulk
= γd Ra,bulk (24)
dt
To solve the differential equations (11)–(15) uniquely, they
where the mass transfer flux, Ja , is being derived from
require an initial condition
boundary condition (9) of the micro model
   
t = 0, [A]g = [A]0g , [A]l = [A]0l , ∂ [A] [A]x=0
Ja = −Da = kg [A]g,bulk − (25)
[B]l = [B]0l , [C]l = [C]0l , [D]l = [D]0l (17) ∂ x x=0 ma
The Eqs. (20)–(24) are only valid provided that the mass
2.4. Overall reactor model transfer and the convection terms that have been left out are
processed instantaneously according to Eq. (18) after each
In literature [23,24] it is sometimes asserted that the macro contact time θ .
model can be implemented as boundary conditions for the The solution of the overall reactor model thus proceeds
micro model. For the penetration model, however, this is not in three steps, where N=1 at the start (see Fig. 1): (1) the
straightforward in case that the liquid phase bulk concen- micro model (Eqs. (3)–(9)) and the simplified macro model
trations may vary significantly during the contact time. This (Eqs. (20)–(25)) are solved simultaneously from t=(N−1)θ
situation occurs when the liquid phase is not at equilibrium to t=Nθ ; (2) at t=Nθ the liquid element is mixed up with
and/or in case of (ir)reversible reactions when Ha<2. the liquid bulk instantaneously using Eqs. (18) and (19); (3)
According to the penetration model a liquid element is N is increased by one and the process is repeated for the
exposed at the gas–liquid interface for a period θ during next period θ . Therefore the overall reactor model produces
which mass transfer takes place. Next, the element is mixed a result after each contact period.
226 E.P. van Elk et al. / Chemical Engineering Journal 76 (2000) 223–237

Fig. 1. The overall reactor model.

2.5. Enhancement factor and utilisation factor of the liquid bulk. The degree of saturation of the liquid
bulk is the ratio of the liquid phase concentration of A to
The enhancement factor, Ea , defined as the ratio of the liquid phase concentration of A when the liquid phase is
the mass flux of component A through the interface with saturated
chemical reaction and driving force ([A]l,i −[A]l,bulk ) to the [A]l,bulk [A]l,bulk
mass flux through the interface without chemical reaction, ηa ≡ ≈ (27)
[A]l,sat ma [A]g
but with the same driving force, can be obtained from the
calculated concentration profiles. If gas phase mass trans- For increasing reaction rate or decreasing mass transfer rate
fer resistance can be neglected (i.e. [A]g,bulk =[A]g,i ) the the degree of saturation will approach to a certain minimum
enhancement factor, Ea , is estimated by (often zero) and for decreasing reaction rate or increasing
mass transfer rate the degree of saturation will approach a
Ja,with reaction certain maximum.
Ea ≡
Ja,without reaction The degree of utilisation of the liquid phase is the ratio
R Nθ of the actual conversion rate of A to the conversion rate
(1/θ ) (N−1)θ − Da (d[A]/dx)x=0 dt
≈ R Nθ (26) of A that would occur if the entire reaction phase were in
(1/θ ) (N −1)θ kl ([A]l,i − [A]l,bulk ) dt equilibrium with the interface
R Nθ
For slow reactions (Ha<0.2), the influence of the reaction (1/θ ) (N −1)θ ξa dt
η ≈ R Nθ (28)
in the liquid bulk has a significant effect, but since the en-
(N −1)θ Ra,sat dt εl VR
hancement factor is always 1.0 for slow reactions it supplies
no additional information. For increasing reaction rate or decreasing mass transfer rate
Therefore two additional dimensionless parameters are the degree of utilisation will approach to a certain minimum
used: the degree of saturation and the degree of utilisation (often zero) and for decreasing reaction rate or increasing
 E.P. van Elk et al. / Chemical Engineering Journal 76 (2000) 223–237 227

mass transfer rate the degree of utilisation will approach to

a certain maximum.

3. Numerical treatment

In the penetration model, the concentration profiles

are time-dependent: they develop a solution of a sys-
tem of coupled non-linear parabolic partial differential
equations subject to specified initial and two point
boundary conditions. The approach used to solve these
models is based on the method presented by Versteeg
et al. [23], however, in the discretisation scheme two ad-
ditional gridlines are introduced for the ordinary differen- Fig. 3. Development of the component concentrations in the reactor during
Case 1.
tial equations of the macro model that have to be solved
simultaneously.
The implicit discretisation method used is known as the 4. Validation results
Baker and Oliphant [25] discretisation. For the time deriva-
tive a three-point backward discretisation is used (with su- 4.1. Introduction
perscript j indicating the time level) leading to (i=A, B,
C, D)
In order to validate the model, runs of the numerical solu-
tion method have been carried out on three fictitious cases.
∂ [i] 3[i]j +1 − 4[i]j + [i]j −1 The numerical solutions are compared with analytical solu-
→ (29)
∂t 21t tions to check the validity of the numerical method used.
The cases are fictitious, and chosen so that, using physically
The finite difference form of the reactor model thus leads to realistic data, all relevant limit cases are tested. Emphasis is
relations between concentrations in five grid points for the put on validation of the micro model, since this is the most
micro model, four grid points for the gas phase macro model complex part of the overall reactor model. Validation of the
and three grid points for the liquid phase macro model, clus- macro model is shown in Figs. 3, 5 and 8. The cases are:
tered as ‘molecules’ as shown in Fig. 2. Only for j=0 these 1. physical absorption: A(g) → A(l) with a dynamic gas
molecules are not possible, because no grid points with time phase and a static liquid phase,
index −1 exist. Therefore in the first step, a two-point back- 2. absorption and equilibrium reaction: A(g) → A(l), A(l)+
ward discretisation (Euler) is used at the cost of lower order B(l)C(l) + D(l) and Ra =k1,1 [A][B]−k−1,−1 [C] [D] in
truncation error. a batch reactor,
The integrals of Eqs. (19), (26) and (28) are calculated 3. absorption accompanied by first-order irreversible chem-
using the Simpson’s 1/3 rule. ical reaction: A(g) → A(l), A(l) → P(l) with Ra =k1 [A].
This case is compared with the analytical solution of the
film model.

4.2. Physical absorption

The first simulation was carried out to check that the

dynamic gas phase macro model (Eq. (20)) was correctly
implemented in the previously tested [23] micro model. The
simulation is described in Table 1.
The results of Case 1 are presented in Fig. 3. It can be
seen that the numerical solution coincides with the analyt-
ical solution. Studying the profiles in the liquid element in
Fig. 4 it can be concluded that the decrease of the concen-
tration at x=0 corresponds to the decrease of the gas phase
concentrations. This is achieved by the implementation of
the gas phase macro model as boundary condition of the
Fig. 2. Discretisation scheme. micro model.
228 E.P. van Elk et al. / Chemical Engineering Journal 76 (2000) 223–237

Table 1
Parameters used in Case 1
Case A(g) → A(l)
[A]g,in 1 mol m−3
[A]0g (initial state) 1 mol m−3
[A]l (fixed) 0 mol m−3
kl 5×10−5 m s−1
kg 100 m s−1 (no gas resistance)
Da 10−9 m2 s−1
ma 0.5
VR 10 m3
εg 0.5
a 100 m2 m−3
8g,in/out 0.01 m3 s−1
θ 0.51 s

Fig. 5. Development of the component concentrations in the reactor during

Case 2.

tion will be at equilibrium (k1,1 /k−1,−1 =[C]l [D]l /[A]l [B]l )

and the gas and liquid phases will also be at equilibrium
([A]l =ma [A]g ).
The results of Case 2 are presented in Fig. 5. The numer-
ical results approach and achieve the analytical steady state
solution.
From the profiles in the liquid element (Fig. 6) it can be
seen that at the beginning component B is consumed within
the liquid element and products C and D are produced. Near
the gas–liquid interface, where the concentration of A has the
highest value the reaction rate is fastest and the consumption
of component B is maximum. Once the stationary state is
Fig. 4. Concentration profiles in the liquid element at the end of a contact
period at different times during Case 1. achieved (t>500) the concentration of all components remain
constant over the entire liquid element. The reaction is at
equilibrium and there is no net gas–liquid mass transfer.
4.3. Equilibrium reaction in a batch reactor
4.4. Absorption and irreversible 1,0-reaction
The second case (Table 2) was carried out to check that the
chemical reaction and the dynamic liquid phase macro model For the third case there is no real analytical solution of the
(Eqs. (18)–(24)) were correctly implemented in the model. A penetration model available. Therefore, a comparison has
reversible reaction in a batch reactor was chosen so that the been made between the numerical model and the analytical
analytical steady state solution is known: the chemical reac-

Table 2
Parameters used in Case 2
Case (batch reactor) A(g) → A(l), A(l) + B(l)C(l) + D(l);
Ra =kR,1,1 [A][B]−kR,−1,−1 [C][D]
kR,1,1 10 m6 mol−2 s−1
kR,−1,−1 1 m6 mol−2 s−1
[A]0g (initial state) 1 mol m−3
[B]0l (initial state) 1 mol m−3
[A, C, D]0l (initial state) 0 mol m−3
kl 5×10−5 m s−1
kg 100 m s−1 (no gas resistance)
Di 10−9 m2 s−1
ma 0.5
VR 10 m3
εg 0.5
a 100 m2 m−3
Fig. 6. Concentration profiles in the liquid element at the end of a contact
θ 0.51 s
period at different times during Case 2.
 E.P. van Elk et al. / Chemical Engineering Journal 76 (2000) 223–237 229

Fig. 7. Comparing analytical (Hinterland concept) and numerical (80 gridpoints) results of saturation and enhancement factors under steady state conditions
of Case 3.

Table 4
solution of the film model using the Hinterland concept [3].
Comparing analytical and numerical results (80 gridpoints) of saturation
Both solutions must result in similar, but not necessarily and enhancement factors under steady state conditions of Case 3
equal results. The conditions for the simulation are described
Ha Al Ha2 Ea,an,pen Ea,num ηa,an,film ηa,num
in Table 3. Fig. 7 and Table 4 give the steady state results of
the saturation (left axis) and the enhancement factor (right 0.001 0.0001 1.000 1.001 0.926 0.925
axis) as a function of both the Hatta number and the value 0.003 0.0009 1.000 1.001 0.925 0.925
0.01 0.01 1.000 1.002 0.917 0.917
of (Al−1)Ha2 , where the Hinterland ratio Al is defined as
0.03 0.09 1.000 1.003 0.855 0.855
the ratio between the total reaction phase volume and the 0.06 0.36 1.002 1.005 0.695 0.695
reaction phase film volume 0.1 1 1.004 1.008 0.482 0.481
0.2 4 1.017 1.020 0.197 0.198
0.5 25 1.103 1.106 0.037 0.040
Table 3 1 100 1.379 1.381 0.008 0.012
Parameters used in Case 3 2 400 2.196 2.199 0.001 0.005
Case A(g) → A(l), A(l) → P(l); Ra =kR,1,0 [A] 4 1600 4.098 4.102 0.000 0.002
kR,1 2.5×10−6 ≤kR,1 ≤250 10 10000 10.04 10.05 0.000 0.001
0.001≤Ha≤10
[A]g,in 1 mol m−3
[i]l,in 0 mol m−3
[A]0g (initial state) 1 mol m−3 εl εl kl
Al = = (30)
[P]0l (initial state) 0 mol m−3 aδ aDa
[A]0l (initial state) 0 mol m−3
kl 5×10−5 m s−1 The physical meaning of (Al−1)Ha2 is the ratio of the
kg 100 m s−1 (no gas resistance) maximum conversion in the liquid bulk to the maximum
Di 10−9 m2 s−1
transport through the film.
ma 0.5
VR 10 m3 The analytical solutions shown in Fig. 7 can be found in
εg 0.5 the book of Westerterp et al. [3]. Several conclusions can be
a 250 m2 m−3 drawn from the figure:
8g,in/out 0.01 m3 s−1 1. For this case the numerical calculated saturation co-
8l,in/out 0.01 m3 s−1
incides the analytical solution of the film model with
θ 0.51 s
Hinterland concept.
230 E.P. van Elk et al. / Chemical Engineering Journal 76 (2000) 223–237

pronounced differences can occur between the stagnant film

model with Hinterland concept [3] and application of the
Higbie penetration model.
The applications discussed here concentrate on micro
scale parameters (concentration profiles) and steady state
results (degree of utilisation and saturation). Application of
the model on macro scale parameters (dynamics of concen-
tration and temperature at the reactor outlet) are discussed
elsewhere [22].

5.2. Absorption and irreversible 1,0-reaction

Additional simulations were performed for Case 3, vary-

ing the specific contact area a from 2.5 to 2500 m−1 . The
Fig. 8. Development of the component concentrations in the reactor during steady state results are shown in Fig. 9 together with the an-
Case 3. alytical solution of the stagnant film model with Hinterland
concept.
For a specific contact area of 2.5, 25 and 250 m−1 the
2. The numerical calculated enhancement factors coincide solution of the stagnant film model seems to coincide the
with the analytical solution of the penetration model. solution of the penetration theory. However, at a specific
3. For Ha<0.2 (slow reactions) the enhancement factor higher contact areas, e.g. 2500 m−1 , there is a substantial
equals unity for all values of (Al−1)Ha2 . The mass discrepancy between the solutions of both the models. The
transfer is not enhanced. In that case the saturation is an penetration theory predicts substantial larger liquid phase
important parameter. concentrations for Ha > 0.2 than predicted by the film model
4. For Ha > 2.0 (fast reaction) the enhancement factor equals with Hinterland concept. For Ha<0.2 no discrepancy is
the Hatta number. The mass transfer is enhanced by the found.
reaction. The reason for this discrepancy is that the film model with
5. For (Al−1)Ha2 1 the saturation reaches an upper limit. Hinterland concept assumes a stagnant film, not being part
6. For (Al−1)Ha2 1 the saturation goes to zero. In that of the liquid bulk, whereas the penetration model assumes
case the enhancement factor is an important parameter. liquid elements that are mixed up with the liquid bulk after
From Table 4 it can be concluded that the analytical and each contact time. Therefore, the species present inside the
the numerical solutions are similar, although a small dis- liquid element are considered as part of the liquid phase and
crepancy is found between the saturation factors at large for specific conditions this results in a larger liquid phase
Hatta numbers (Ha > 0.2). This is discussed in more detail in concentration.
Section 5.2. For processes controlled by the rate of the kinetics
For Ha > 2.0, a dynamic analytical solution is also avail- (Ha<0.2), these two different methods give equal results
able. The development of the concentrations in time for Case because most species are present inside the liquid bulk and
3 is shown in Fig. 8 for Ha=2.0. Fig. 8 demonstrates that the excess of species present in the liquid element (the
the dynamic solution obtained analytically gives exactly the last term of Eq. (19)) can be neglected. For processes con-
same results as the numerical solution obtained using the trolled completely or partly by the rate of the mass transfer
presented model. (Ha > 0.2), these two different methods give different results
because the excess of species present in the liquid element
(the last term of Eq. (19)) can no longer be neglected com-
5. Applications pared to the total amount of species present in the liquid
phase (Eq. (19)).
5.1. Introduction The above mentioned explanation means that the discrep-
ancy must exist for all situations where Ha > 0.2, indepen-
In the previous section the numerical results of the model dent of the specific contact area a. How to explain that in
were validated by comparison with exact, analytical solu- Fig. 9 a discrepancy is only visible for a=2500 m−1 ? This
tions. This section discusses cases for which only approx- is because of the scale of the y-axis in Fig. 9. The cases
imate analytical or no solutions at all are available. Again, for specific contact areas of 2.5, 25 and 250 m−1 have the
all data used is fictitious and chosen so that some interesting same discrepancies between the film model with Hinterland
phenomena can be stressed. concept and the presented implementation of the penetra-
It is shown that under certain defined conditions substan- tion theory, but the differences are not visible because of the
tial deviations exist between the exact numerical and the scale of the y-axis. From Table 5 it is clear that there is a
approximate analytical results. It is also demonstrated that discrepancy for all cases once Ha > 0.2.
 E.P. van Elk et al. / Chemical Engineering Journal 76 (2000) 223–237 231

Fig. 9. Comparing analytical (Hinterland concept) and numerical results of saturation under steady state conditions of Case 3, varying the specific contact
area from 2.5 to 2500 m−1 .

Table 5  s
∼ Db [B]bulk Da
Saturation under steady state conditions of Case 3 obtained with the Ea∞ = 1+ (31)
presented model divided by the analytical solution of film model with γb Da [A]i Db
Hinterland concept. A value of 1.0 means that both models give equal
solutions For situations where the diffusion coefficients Da and Db are
Ha a=2.5 m−1 a=25 m−1 a=250 m−1 a=2500 m−1 the same the van Krevelen and Hoftijzer method gives very
0.001 1.0 1.0 1.0 1.0
good results (Table 7). The conditions for this case were
0.003 1.0 1.0 1.0 1.0 chosen so that the enhancement factor of the van Krevelen
0.01 1.0 1.0 1.0 1.0 and Hoftijzer method equals Ea∞ so that they are easy to
0.03 1.0 1.0 1.0 1.0 verify.
0.06 1.0 1.0 1.0 1.0 However, if the diffusion coefficients of components A
0.1 1.0 1.0 1.0 1.0
0.2 1.0 1.0 1.0 1.0
and B are different, the deviations increase significantly
0.5 1.1 1.1 1.1 1.1
1 1.4 1.4 1.4 1.4
Table 6
2 3.6 3.6 3.6 3.6
Parameters used in Cases 4 and 5
4 27 27 27 27
10 ∼11000 ∼11000 ∼11000 ∼11000 Case A(g) → A(l), A(l)+B(l) → P(l); Ra =kR,1,1 [A][B]
kR,1,1 250 000 (Ha2) (Case 4)
2.5×10−6 ≤kR ,1,1 ≤2.5×1012 (Case 5)
[A]g,in 1 mol m−3
[B]l,in 1 mol m−3
5.3. Absorption and irreversible 1,1-reaction kl 5×10−5 m s−1
kg 100 m s−1 (no gas resistance)
For this case (see Table 6) an exact analytical solution Da 10−9 m2 s−1
is not available. van Krevelen and Hoftijzer [13] gave an Db,p 10−9 or 5×10−9 m2 s−1
ma 1≤ m≤3 (Case 4)
approximate analytical solution method. The essence of their
m=0.5 (Case 5)
method is the approximation of the concentration profile of VR 10 m3
component B by a constant [B]1 all over the reaction zone, εg 0.5
so that the reaction becomes essentially pseudo first-order. a 100 m2 m−3
For calculation of the enhancement factor their method 8g,in/out 0.01 m3 s−1
8l,in/out 0.01 m3 s−1
requires the maximum enhancement factor, which for the
θ 0.51 s
penetration model is at least in good approximation given by
232 E.P. van Elk et al. / Chemical Engineering Journal 76 (2000) 223–237

Fig. 10. Numerical results of saturation under steady state conditions of Case 5, varying the specific contact area from 2.5 to 2500 m−1 .

Table 7 Table 9
Comparing approximate analytical solution of Eq. (31) with numerical Saturation under steady state conditions of Case 5 obtained with the
results (160 gridpoints) for Case 4 with Da =Db =10−9 m2 s−1 presented model for large Hatta numbers (Ha>2.0)
ma Eq. (31) Ea,num a=2.5 m−1 a=25 m−1 a=250 m−1 a=2500 m−1
1.0 2.000 2.002 0.00004 0.0004 0.004 0.04
2.0 1.500 1.502
3.0 1.333 1.335

contact area a from 2.5 to 2500 m−1 . The simulation con-

(Table 8). This is probably caused by a wrong approxima- ditions are given in Table 6 and the steady state results are
tion of the infinite enhancement factor (Eq. (31)). For a case shown in Fig. 10.
with [B]bulk /γ b [A]i not much larger than one, Lightfoot [26] For a specific contact area of 2.5, 25 and 250 m−1 the so-
suggests the following equation: lution is at first glance as expected: at low Hatta numbers the
s saturation approaches a maximum and at high Hatta numbers
∼1+ D b [B]bulk Da the saturation approaches zero. This was explained in Sec-
Ea∞ = (32)
γb Da [A]i Db tion 4.4, conclusions 5 and 6. However, at a specific contact
area of 2500 m−1 , the saturation clearly does not approach
which is also given in Table 8.
zero (Fig. 10), even at very high Hatta numbers (Ha=1000).
From Table 8 it can be concluded that large discrepancies
The penetration theory predicts still substantial concentra-
exist between Eq. (31) and (32) and the exact numerical
tions of species A present in the liquid phase for Ha > 2.0. A
solution under certain conditions.
saturation was found of 0.04 (4%) for Ha > 7.0, independent
Additional simulations were performed with Case 5,
of the Hatta number (7.0<Ha<70 000). The Hinterland ratio
which differs slightly from Case 4, varying the specific
Al equals 10 for Case 5 at a=2500 m−1 , which means that
the process is still characterised by a significant liquid bulk.
Table 8 The ratio of the maximum conversion in the liquid bulk to the
Comparing approximate analytical solutions of Eqs. (31) and (32) with
numerical results (160 gridpoints) for Case 4 with Da =10−9 m2 s−1 and
maximum transport through the film, (Al−1)Ha2 , is much
Db =5×10−9 m2 s−1 larger than 1 for Ha > 7.0 ((Al−1)Ha2 =441 at Ha=7.0).
More detailed analysis of the cases for a specific contact
ma Eq. (31) Ea,num Eq. (32)
area of 2.5, 25 and 250 m−1 shows that here also the sat-
1.0 2.681 2.892 3.232 uration does not converge to zero at high Hatta numbers.
2.0 1.565 1.871 2.117
This is not visible from Fig. 10 because of the scale of the
3.0 1.193 1.554 1.745
y-axis, but from Table 9 it is clear that there is a limit slightly
 E.P. van Elk et al. / Chemical Engineering Journal 76 (2000) 223–237 233

above zero for all cases. The reason for the saturation not The reason that for a 1,0-reaction the area under the con-
approaching zero is that the penetration model assumes liq- centration profile converges to zero, while for a 1,1-reaction
uid elements that are mixed up with the liquid bulk after this is not the case, is explained below:
each contact period. Therefore, the species present inside 1. For a 1,0-reaction, species A is the only required reactant
the liquid element are considered part of the liquid phase. and increasing the Hatta number by increasing the reac-
We are aware that Ha=70 000 are not usually found in tion rate constant will under all circumstances increase
process industry. We only used these extreme numbers to the conversion of species A inside the liquid element and
show that even at these extremely high Hatta numbers, the lower the concentration of species A. At large enough
saturation does not converge to zero. In fact we show that, rate constants the concentration of species A will con-
at a realistic Ha=7, the same saturation is found as for verge to zero.
Ha=70 000 (Fig. 10). This is true, no matter what the value 2. For a 1,1-reaction, species A reacts with species B. This
of the specific contact area is (Table 9). means that component B is required to convert compo-
Why does for a 1,0-reaction (Section 5.2, Fig. 9) the sat- nent A into the products. Increasing the Hatta number by
uration of the liquid phase with species A converge to zero increasing the reaction rate constant will initially increase
with increasing Hatta number, while for a 1,1-reaction (Fig. the conversion of species A. However, above a certain
10) the saturation converges to a certain constant slightly reaction rate the supply of component B coming from
above zero, even if the Hatta number is increased up to the liquid bulk becomes limiting — a so-called instan-
70 000? This can be explained by studying the concentra- taneous reaction — and further increase of the reaction
tion profiles of the liquid element at the end of the contact rate does not increase the conversion. At large enough
period. rate constants the concentration of species A will con-
The concentration profiles of the 1,0-reaction at various verge to a constant value larger than zero. From the con-
Hatta numbers are shown in Fig. 11. From this figure it is centration profiles shown in Fig. 12 it is obvious that for
clear that the amount of species present in the liquid element Case 5 at Ha>7 the supply of component B for the liquid
(the area under the graph) converges to zero with increasing phase becomes limiting (for x<0.4 component B is not
Hatta number and therefore the liquid phase saturation will available and the reaction is instantaneous).
also converge to zero.
The concentration profiles of the 1,1-reaction at various
Hatta numbers are shown in Fig. 12. From this figure it is 5.4. Instantaneous reaction without enhancement
clear that the amount of species A present in the liquid ele-
ment (the area under the graph of species A) does not con- A special case occurs for instantaneous reactions (Ha2)
verge to zero with increasing Hatta number. For this specific without chemical enhancement of mass transfer (Ea =1).
case, the area under the graph does not change for Hatta Case 4 (see Table 6) is an example of such a case if
numbers larger than about 7. Therefore the liquid phase sat- a=2500 m−1 , Db =5×10−9 m2 s−1 and m=3.0. For these
uration does converge to a constant value, corresponding to conditions the enhancement factor equals the maximum
the amount of species A marked by the grey triangle in Fig. enhancement factor Ea∞ =1.09≈1.
12. The size and shape of this triangle is the same for all The concentration profile in the liquid element for this
Ha>7 for this case. case is shown in Fig. 13. The penetration theory does not
have any difficulties with calculating and processing this pro-
file. The corresponding profile for the stagnant film model is
shown schematically in Fig. 14. From this figure it is clear
that the concentration profile of the film model has a dis-
continuity at x=δ for this specific case. This discontinuity
may lead to serious difficulties when applying the stagnant
film model for such a case.

5.5. Equilibrium reaction

Case 6 in Table 10 defines a case with an equilibrium

reaction. The equilibrium constant is 10 and varying the
reaction rate constant of the forward reaction varies the Hatta
number (based on the forward reaction) from 0.001 to 100.
The steady state saturation is shown in Fig. 15 as a function
Fig. 11. Concentration profiles in the liquid element at the end of a
contact period for Case 3 at steady state conditions and a=2500 m−1 .
of Hatta number and specific contact area. The model has
The reaction rate constant was varied to obtain Hatta numbers of 1, 2, 4 no difficulties handling the area where the liquid phase is
and 10. not at equilibrium (Ha<1.5).
234 E.P. van Elk et al. / Chemical Engineering Journal 76 (2000) 223–237

Fig. 12. Concentration profiles in the liquid element at the end of a contact period for Case 5 at steady state conditions and a=2500 m−1 . The reaction
rate constant was varied to obtain Hatta numbers of 1, 3, 8, 685 and 70 494.

Fig. 13. Concentration profiles in the liquid element at the end of a contact period for Case 4 at steady state conditions and a=2500 m−1 , Db =5×10−9 m2 s−1
and m=3.0 (instantaneous reaction without enhancement).

Comparing Fig. 15 with Fig. 10 shows the influence of increases the saturation. For very low Hatta numbers the
the backward reaction. The backward reaction decreases the saturation converges to the same limits in both Figs. 10 and
net conversion of species A in the liquid phase and thereby 15. At large Hatta numbers the influence of the backward
reaction forces the saturation of species A to increase.
Table 10
Parameters used in Case 6
Case A(g) →A(l), A(l) + B(l)C(l) + D(l),
Ra =kR,1,1 [A][B]−kR ,−1,−1 [C][D]
kR ,1,1 2.5×10−6 ≤kR ,1 ≤250 000 m6 mol−2 s−1
Keq =kR ,1,1 /kR ,−1,−1 10
[A]0g (initial state) 1 mol m−3
[B]0l initial state) 1 mol m−3
[A, C, D]0l (initial state) 0 mol m−3
kl 5×10−5 m s−1
kg 100 m s−1 (no gas resistance)
Di 10−9 m2 s−1
ma 0.5
VR 10 m3
εg 0.5
a 2.5, 25, 250 and 2500 m2 m−3
8g,in/out 0.01 m3 s−1
Fig. 14. Schematic concentration profile for the stagnant film model,
8l,in/out 0.01 m3 s−1
corresponding to the profile for the penetration model presented in Fig.
θ 0.51 s
13 (instantaneous reaction without enhancement).
 E.P. van Elk et al. / Chemical Engineering Journal 76 (2000) 223–237 235

Fig. 15. Numerical results of saturation under steady state conditions of Case 6, varying the specific contact area from 2.5 to 2500 m−1 .

6. Conclusions The numerical model presented in this study is especially

useful for complex systems for which no analytical solu-
The main conclusion is that the dynamic behaviour of ide- tions exist and for which the stagnant film model in not
ally stirred gas–liquid reactors can be simulated successfully accurate.
over a wide range of conditions with the numerical solution
method presented in this study. From comparison of the nu-
merical results with analytical solutions it was concluded
that the errors can be neglected. 7. Nomenclature
Contrary to many other numerical and analytical models
the present model can be used over a wide range of con- a specific surface area (m2 m−3 )
ditions. The model is suitable for the entire range of Hatta A component A
numbers for both reversible and irreversible reactions. Al Hinterland ratio (defined by εl /δa)
Multiple (in)dependent reactions, even with more than one B component B
gas phase component, can be implemented in the model as C component C
well, by extending the number of components and the kinetic D component D
rate expression. The model can also be used for simulation Dsubscript diffusivity (m2 s−1 )
of (semi-)batch reactors. Esubscript enhancement factor
Since the Hinterland concept [3] requires a constant mass Esubscript∞ enhancement factor for instantaneous reaction
transfer flux from the liquid film or element to the liquid Ha Hatta number defined by (kR,m,n,p,q [A]m−1
bulk phase, the Hinterland concept is not directly applica- [B]n [C]p [D]q Da )0.5 /k1
ble to the penetration models. The model presented in this Jsubscript molar flux (mol m−2 s−1 )
study successfully uses a different approach to implement kg gas phase mass transfer coefficient (m s−1 )
the simultaneous solution of the Higbie penetration and the kl liquid phase mass transfer coefficient (m s−1 )
dynamic gas–liquid macro model. kov overall mass transfer coefficient (m s−1 )
Under some circumstances substantial differences exist kR ,subscript reaction rate constant (m3(m +n +p +q −1)
between the exact numerical and existing approximate re- mol−(m +n +p +q −1) s−1 or m3(r +s +t +v −1)
sults. For specific cases substantial differences can exist be- mol−(r +s +t +v −1) s−1 )
tween the results obtained using the stagnant film model m reaction order
with Hinterland concept and the presented implementation msubscript gas–liquid partition coefficient
of the Higbie penetration model. n reaction order
236 E.P. van Elk et al. / Chemical Engineering Journal 76 (2000) 223–237

N counter with start value 1 at t=0 Superscripts

Nsubscript number of moles (mol)
0 initial value
p reaction order
j time level
q reaction order
r reaction order
Rsubscript reaction rate (mol m−3 s−1 ) Acknowledgements
s probability of replacement according to
surface renewal model
These investigations were supported by DSM Research
s reaction order
Geleen.
t reaction order
t simulation time variable (s)
v reaction order References
VR reactor volume (m3 )
x place variable (m) [1] L.K. Doraiswamy, M.M. Sharma, Heterogeneous Reactions, Vol. 2,
√ Wiley, New York, 1984.
xdimensionless place variable defined as x/ 4Da θ [2] S. Hatta, Tech. Rep. Tohoku Imp. Univ. 10 (1932) 119.
[] concentration (mol m−3 ) [3] K.R. Westerterp, W.P.M. van Swaaij, A.A.C.M. Beenackers,
[]subscript concentration (mol m−3 ) Chemical Reactor Design and Operation, Wiley, New York, 1990.
[4] W.G. Whitman, Preliminary experimental confirmation of the
two-film theory of gas absorption, Chem. Metall. Eng. 29 (1923)
146–148.
Greek letters [5] R. Higbie, The rate of absorption of a pure gas into a still liquid
δ film thickness according to film model (m) during short periods of exposure, Trans. Am. Inst. Chem. Eng. 35
(1935) 36–60.
δp assumed thickness of liquid element (m) [6] P.V. Danckwerts, Significance of liquid-film coefficients in gas
εg gas phase hold-up absorption, Ind. Eng. Chem. 43 (1951) 1460–1467.
εl liquid phase hold-up [7] W.E. Dobbins, in: M.L. McCable, W.W. Eckenfelder (Eds.),
γ subscript stoichiometry number Biological Treatment of Sewage and Industrial Wastes, Part 2-1,
Reinhold, New York, 1956.
η degree of utilisation [8] H.L. Toor, J.M. Marchello, Film-penetration model for mass transfer
ηsubscript degree of saturation and heat transfer, AIChE J. 4 (1958) 97–101.
θ contact time according to penetration [9] T.K. Sherwood, R.L. Pigford, Absorption and Extraction,
model (defined by 4Da /πkl2 ) (s) McGraw-Hill, New York, 1952, p. 332.
[10] P.V. Danckwerts, A.M. Kennedy, Kinetics of liquid-film process in
τg gas phase residence time (s) gas absorption. Part 1: Models of the absorption process, Trans. Inst.
τl liquid phase residence time (s) Chem. Eng. 32 (1954) 49–53.
ξ subscript absolute conversion (mol) [11] D.R. Olander, Simultaneous mass transfer and equilibrium chemical
8g gas phase volume flow rate (m3 s−1 ) reaction, AIChE J. 6 (1960) 233–239.
[12] C.-J. Huang, C.-H. Kuo, Mathematical models for mass transfer
8l liquid phase volume flow rate (m3 s−1 ) accompanied by reversible chemical reaction, AIChE J. 11 (1965)
901–910.
[13] D.W. van Krevelen, P.J. Hoftijzer, Kinetics of gas–liquid reactions.
Subscripts Part 1: General theory, Rec. Trav. Chim. 67 (1948) 563–568.
[14] K. Onda, E. Sada, T. Kobayashi, M. Fujine, Gas absorption
a component A accompanied with complex chemical reactions — I. Reversible
an analytical solution chemical reactions, Chem. Eng. Sci. 25 (1970) 753–760.
[15] W.J. DeCoursey, Enhancement factors for gas absorption with
b component B reversible reaction, Chem. Eng. Sci. 37 (1982) 1483–1489.
bulk at bulk conditions [16] H. Hikita, S. Asai, A. Yano, H. Nose, Kinetics of absorption of
c component C carbon dioxide into aqueous sodium sulfite solutions, AIChE J. 28
d component D (1982) 1009–1015.
[17] L. Landau, Desorption with a chemical reaction, Chem. Eng. Sci.
elem at liquid element conditions 47 (1992) 1601–1606.
film film model [18] S.J. Parulekar, N.A. Saidina Amin, Complex gas–liquid reactions:
g gas phase feedback from bulk liquid to liquid-side film, Chem. Eng. Sci. 51
i interface (1996) 2079–2088.
[19] J.J. Romanainen, T. Salmi, Numerical strategies in solving gas–liquid
i species i
reactor models — 1. Stagnant films and a steady state CSTR, Comput.
in at inlet conditions Chem. Eng. 15 (1991) 769–781.
l liquid phase [20] J.J. Romanainen, T. Salmi, Numerical strategies in solving gas–liquid
num numerical solution reactor models — 2. Transient films and dynamic tank reactors,
out at outlet conditions Comput. Chem. Eng. 15 (1991) 783–795.
[21] G.F. Versteeg, P.M.M. Blauwhoff, W.P.M. van Swaaij, The effect of
p penetration element diffusivity on gas–liquid mass transfer in stirred vessels. Experiments
pen penetration model at atmospheric and elevated pressures, Chem. Eng. Sci. 42 (1987)
sat at saturation conditions 1103–1119.
 E.P. van Elk et al. / Chemical Engineering Journal 76 (2000) 223–237 237

[22] E.P. van Elk, P.C. Borman, J.A.M. Kuipers, G.F. Versteeg, Modelling complex reversible chemical reactions in gas–liquid systems: an
of gas–liquid reactors — stability and dynamic behaviour of overview, Chem. Eng. Sci. 47 (1992) 3181–3195.
gas–liquid mass transfer accompanied by irreversible reaction, Chem. [25] G.A. Baker, T.A. Oliphant, An implicit, numerical method for solving
Eng. Sci. 54 (1999) 4869–4879. the two-dimensional heat equation, Quart. Appl. Math. 17 (1960)
[23] G.F. Versteeg, J.A.M. Kuipers, F.P.H. van Beckum, W.P.M. van 361–373.
Swaaij, Mass transfer with complex reversible chemical reactions — [26] E.N. Leightfoot, Approximate expressions for predicting the effect
I. Single reversible chemical reaction, Chem. Eng. Sci. 44 (1989) of fast second-order chemical reaction on interphase mass-transfer
2295–2310. rates, Chem. Eng. Sci. 17 (1962) 1007–1011.
[24] W.P.M. van Swaaij, G.F. Versteeg, Mass transfer accompanied with