Dynamic Modelling of Batch Distillation Columns

Maria Nunes de Almeida Viseu

mariaviseu@ist.utl.pt

Instituto Superior Técnico, Lisbon, Portugal

November 2014

Abstract
Batch distillation is becoming increasingly important in specialty product industries where high added value
chemical compounds are produced with uncertain demands and lifetimes. Mathematical models that predict the
separation time and product purity thereby facilitating plant scheduling are required. The primary purpose of
this study is thus to develop a batch multi-staged distillation model in gPROMS ModelBuilder® based on mass
and energy balances, equilibrium stages and tray hydraulic relations. The multi-staged batch model was validated
for a methanol/water separation using literature data from an experimental pilot plant and from theoretical
results given by a model implemented in Fortran language and by commercial simulator Batchsim of Pro/II.
The results simulated in ModelBuilder for the batch multi-staged model reveal a 6.2% overestimation of the
experimental methanol recovery. A very good agreement is found between the ModelBuilder and Fortran
models: the methanol recovery predicted by ModelBuilder is only 2.3% lower. It is shown that the
ModelBuilder multi-staged batch model is robust with ±10% heat duty variations or ±0.5 reflux ratio
differences both affecting the total experiment time in approximately 12%. Differences of ±0.5 in the reflux
ratio are found to have a 2.2 to 6.4% absolute impact on the methanol recovery whilst this recovery is
practically not affected by 10% heat duty variations. This work offers a tool that may be applied to the
scheduling of batch chemical plants and aid industrial management at the planning level.
Keywords: Batch distillation; equilibrium stages; computational models; gPROMS.

1. Introduction negligible vapour and liquid tray holdups are useful if a

Batch distillation is becoming increasingly important in preliminary study of batch distillation is required.
industries such as fine chemical, food, biochemical and Luyben [1] studied a shortcut method to design batch
pharmaceutical companies. The flexibility of the process distillation columns, based on the assumption of
is a key issue in the separation performance of these constant relative volatility. A shortcut model with
industries where a wide range of chemical products with constant relative volatility was used by Barolo &
high added value are produced with uncertain demands Botteon [2] for a column at total reflux to separate pure
and lifetimes. Furthermore, the separation of multiple components from a binary mixture. Diweker &
components from a product mixture is one of the major Madhaven [3] developed a shortcut method for the two
difficulties in the production processes of these operating policies of constant distillate composition and
industries. By employing a single batch column a constant reflux ratio. The shortcut method of Sundaram
multicomponent mixture is separated into its pure & Evans [4] is also easy to employ; however, it applies
components whereas to serve the same purpose, the only to the case of constant reflux. In both latter cases
continuous process counterpart requires a train of the Fenske-Underwood-Gilliland (FUG) shortcut
distillation columns. method for continuous distillation is employed treating
Because the demand for efficiency and productivity in batch rectification as a sequence of continuous, steady-
the chemical industry is increasing, reliable mathematical state rectifications at successive time steps.
models that easily predict the separation degree of Stage-by-stage composition, temperature and flowrate
mixtures are needed. Hence, the main objective of this profiles are needed to perform a detailed simulation of a
work is to provide a rigorous dynamic distillation model batch multicomponent rectification. The first rigorous
in gPROMS ModelBuilder® based on mass and energy batch rectification model based on the assumptions of
balances, equilibrium equations and tray hydraulic negligible vapour holdup, constant liquid holdup and
correlations. perfect mixing of the phases on equilibrium trays was
Shortcut models for handling multicomponent mixtures developed by Meadows [5]. Distefano [6] extended this
under the assumptions of constant molar overflow and model by incorporating a computer-based-solution

1
procedure used to successfully simulate commercial 2. Dynamic multi-staged distillation model
batch-distillation columns. The multi-staged model implemented in gPROMS
If stage holdups are to be included in a dynamic column ModelBuilder® [10] consists of a countercurrent
model, the effective liquid holdup on the trays must be cascade of 𝑁 equilibrium stages (see Figure 1). Interstage
estimated. Bennett developed a clear liquid height liquid from the above stage 𝑗 − 1 with molar flowrate
calculation applicable to froth-type regimes based on 𝐿𝑗−1 and a mole fraction of component 𝑖 of 𝑥𝑖,𝑗−1 and
the weir height, the liquid flow and the froth density [7]. interstage vapour from the adjacent stage 𝑗 + 1 with
The Jeronimo & Sawistowski correlation was flowrate 𝐿𝑗+1 and a molar composition 𝑦𝑖,𝑗+1 of 𝑖 are
successfully used as a building block for spray regime entering stage 𝑗. Liquid and vapour feeds can also enter
entrainment, strictly predicting clear liquid heights at the stage 𝑗 at a total molar flowrate 𝐹𝑗 . The liquid and
froth to spray transition. However, it was shown that vapour phases are leaving each stage 𝑗 at temperature 𝑇𝑗
the clear liquid height in the spray regime is similar to and pressure 𝑃𝑗 .
the clear liquid height at this transition [8] and thus the
correlation is applicable to spray-type regimes.
The total pressure drop across a tray has to be taken
under consideration when implementing a pressure-
driven dynamic column system. To evaluate the
pressure drop of the vapour through the orifices of a
sieve tray, Fair and Van Winkle recommended the
correlation by Liebson et al. whereas Ludwig and Chase
suggested the Hughmark & O’Connell correlation [8].
For both the Liebson and the Hughmark & O’Connell
correlations, the orifice discharge coefficient is a
function of the ratio of tray thickness to hole diameter
and the fractional hole area. To determine the pressure
drop through the aerated liquid on the trays, Fair
followed a dimensionless tray aeration factor approach
correcting the clear liquid flow for aeration effects [8]. A
more recent and fundamental relationship to determine
the pressure loss through the aerated mass was
recommended by Bennett [9]. In this case, the pressure
drop through the aerated liquid is the sum of two
separate terms: the effective clear-liquid height used to
determine the liquid holdup and the residual pressure
drop due to surface tension. For 302 experimental data Figure 1 - Countercurrent cascade of N column stages.
points covering a wide range of systems the Bennett
pressure drop correlation gave an average relative error 2.1. Model assumptions
of ±0.35%. An error of approximately 5% was obtained The mathematical model is based on the following
for Fair’s correlation considering a similar data base [7]. assumptions:
In this work, tray hydraulic correlations comprising tray  Liquid and vapour feeds and products can
holdup and tray pressure drop estimations are first enter any tray;
implemented for a continuous methanol/water case  Neither of the existing phases entrains the
study. Having tested the hydraulic equations, a multi- other phase on each stage;
staged dynamic model is then validated against both  No heat is lost from any stage, i.e., the column
experimental data and theoretical models implemented is adiabatic;
in Fortran and Batchsim of Pro/II for a  The liquid and vapour phases are perfectly
methanol/water separation. The robustness of the mixed on each column stage;
multi-staged dynamic model is evaluated by performing  Phase equilibrium, that is, thermal, mechanical
a sensitivity analysis, analysing the effect of key variables and thermodynamic equilibrium occurs at each
such as the reflux ratio and the heat duty on separation stage (equilibrium stage model);
performance indicators, namely the total separation time  The downcomer holdup is negligible;
and the product recovery.  No chemical reactions occur.

2
2.2. Model equations correlation a dimensionless effective froth density ∅𝑒,𝑗 is
The tray model equations are contained in two subsets: calculated for each stage using a factor 𝐾𝑠,𝑗 :
the material, equilibrium, summation and energy
balances, i.e., the MESH equations; tray hydraulic ∅𝑒,𝑗 = exp[−12.55𝐾𝑠,𝑗 0.91 ] , 𝑗 = 1, … 𝐽. (7)
correlations covering tray holdup and tray pressure drop
predictions.
𝜌𝑗𝑉
𝐾𝑠,𝑗 = 𝑣𝑎,𝑗 ( )0.5 , 𝑗 = 1, … 𝐽. (8)
2.2.1. MESH equations 𝜌𝑗𝐿 − 𝜌𝑗𝑉
For 𝑐 components and 𝐽 stages the MESH equations are
as follows: where 𝑣𝑎,𝑗 refers to the vapour velocity (m/s) through
the active area 𝐴𝑎 . The froth height ℎ𝑓,𝑗 (m) is calculated
𝑑𝑀𝑖,𝑗
+ 𝐿𝑗 𝑥𝑖,𝑗 + 𝑉𝑗 𝑦𝑖,𝑗 = 𝐿𝑗−1 𝑥𝑖,𝑗−1 + 𝑉𝑗+1 𝑦𝑖,𝑗+1 + 𝐹𝑗𝐿 𝑥𝑖,𝑗
𝐹 and multiplied by the effective froth density giving the
𝑑𝑡 (1)
𝑉 𝐹
+ 𝐹𝑗 𝑦𝑖,𝑗 , 𝑖 = 1, … 𝑐; 𝑗 = 1, … 𝐽. effective clear liquid height ℎ𝑐,𝑗 𝐵 (m):
ℎ𝑐,𝑗 𝐵 = ℎ𝑓,𝑗 × ∅𝑒,𝑗 , 𝑗 = 1, … 𝐽. (9)
𝐿
𝑦𝑖,𝑗 ∅̂𝑖,𝑗
= , 𝑖 = 1, … 𝑐. ; 𝑗 = 1, … 𝐽. (2) 𝑞𝑗 2
𝑥𝑖,𝑗 ∅̂𝑖,𝑗 𝑉 ℎ𝑓,𝑗 = 10−3 {ℎ𝑤 + 15330𝐶( )3 } , 𝑗 = 1, … 𝐽. (10)
∅𝑒,𝑗

∑ 𝑥𝑖,𝑗 = 1, 𝑗 = 1, … 𝐽.
𝐶 = 0.0327 + 0.0286 exp(−0.1378 ℎ𝑤 ) (11)
(3)
𝑖
where ℎ𝑤 is the weir height (mm) and 𝑞𝑗 the liquid
volumetric flow in each stage (m3/s). In the Jeronimo &
∑ 𝑦𝑖,𝑗 = 1, 𝑗 = 1, … 𝐽. (4) Sawistowski correlation ℎ𝑐,𝑗 𝑆 (m) is calculated using the
𝑖 total hole area 𝐴ℎ (m2), the active area 𝐴𝑎 (m2) and the
hole diameter 𝑑ℎ (mm):
𝑑𝑈𝑗
+ 𝐿𝑗 ℎ𝑗 𝐿 + 𝑉𝑗 ℎ𝑗 𝑉 = 𝐿𝑗−1 ℎ𝑗−1 𝐿 + 𝑉𝑗+1 ℎ𝑗+1 𝑉 + ℎ𝑗 𝐹𝐿 𝐹𝑗𝐿
𝑑𝑡 (5) ℎ𝑐,𝑗 𝑆 = 10−3 {ℎ𝑐,𝐻2𝑂,𝑗 (
996 0.5(1−𝑛)
+ ℎ𝑗 𝐹𝑉 𝐹𝑗𝑉 , 𝑗 = 1, … 𝐽.
) }, 𝑗 = 1, … 𝐽. (12)
𝜌𝑗𝐿
𝐿 𝑉
The fugacity coefficients ∅̂𝑖,𝑗 and ∅̂𝑖,𝑗 from the
𝐴ℎ −0.791 0.833
equilibrium equations (2) are determined using the 0.497
𝐴𝑎
𝑑ℎ
ℎ𝑐,𝐻2𝑂,𝑗 = , 𝑗 = 1, … 𝐽. (13)
Infochem Multiflash™ software supplied by KBC 𝐴 −1.79
1 + 0.013𝐿𝑗 −0.59 ℎ
Advanced Technologies [11]. The NRTL model and the 𝐴𝑎

Soave-Redlich-Kwong equation of state are used for the 𝑑ℎ

liquid and vapour phases, respectively. The energy 𝑛 = 0.00091 (14)
𝐴ℎ
holdup accumulation term is calculated using the energy 𝐴𝑎

balance (5) for each stage and thus the specific liquid 2.2.2.2. Tray pressure drop
and vapour molar enthalpies ℎ𝑗 𝐿 and ℎ𝑗 𝑉 must also be The total pressure drop across a sieve tray is the sum of
determined in Multiflash. the pressure drop due to vapour flow friction through
dry tray perforations, ℎ𝑑 , and the pressure drop through
2.2.2. Tray hydraulic equations the aerated liquid mass on a tray, ℎ𝑙 . In the present
2.2.2.1. Tray holdup work, the standard orifice equation, the Liebson
The effective clear liquid height ℎ𝑐,𝑗 gives a measure of correlation and the Hughmark & O’Connell correlation
the liquid level on the tray and can be used to calculate were used to estimate ℎ𝑑 . For sieve plates, the pressure
the liquid molar holdup for a given stage 𝑗, as follows: drop through the dispersers ℎ𝑑,𝑗 (m) is given by:
𝐿
ℎ𝑐,𝑗 × 𝐴𝑎 × 𝜌𝑚,𝑗 = 𝑀𝑗𝐿 , 𝑗 = 1, … 𝐽. (6) 50.8 𝜌𝑗𝑉
ℎ𝑑,𝑗 = 𝑣ℎ,𝑗 2 , 𝑗 = 1, … 𝐽. (15)
𝐶𝑣 2 𝜌𝑗𝐿

where 𝐴𝑎 refers to the active area and 𝜌𝑚,𝑗

𝐿
and 𝑀𝑗𝐿 are the where 𝑣ℎ,𝑗 refers to the vapour velocity through the
liquid molar density and liquid molar holdup of stage 𝑗, holes (m/s). Note that SI units should be used when
respectively. In this study the Bennett clear liquid height estimating ℎ𝑑,𝑗 using the above equation. If the simple
correlation and the Jeronimo & Sawistowski correlation standard orifice equation is required, a constant value of
𝐶𝑣𝑆 = 0.75 should be given. The more accurate Liebson
have been used. In the Bennett clear liquid height
correlation uses directly the fractional hole area 𝐴𝐴ℎ to
𝑎

3
estimate the orifice coefficient 𝐶𝑣𝐿 shown in the above
472000𝜎𝑗 𝑔(𝜌𝑗𝐿 − 𝜌𝑗𝑉 ) 1⁄
equation. 𝐶𝑣𝐿 is given by: ℎ𝑅,𝑗 = [ ] 3, 𝑗 = 1, … 𝐽. (24)
𝑔𝜌𝑗𝐿 𝑑ℎ × 𝜎𝑗 × 106
𝐴ℎ 𝑡
0.29 𝑡⁄𝑑 −0.56 (16)
𝐶𝑣𝐿 = 0.74 +𝑒 ℎ
(25)
𝐴𝑎 𝜎𝑗 = 𝑓(𝑇𝑗 , 𝑃𝑗 , ⃗⃗⃗
𝑥𝑗 , ⃗⃗⃗
𝑦𝑗 ), 𝑗 = 1, … 𝐽.

where 𝑑ℎ is the hole diameter and 𝑡𝑡 the tray thickness. 3. Results and discussion
The Hughmark & O’Connell correlation determines an 3.1. Column hydraulics
orifice coefficient 𝐶𝑣𝐻 fit by the following equation: Table 1 shows the operating and geometrical data for
the continuous methanol/water case study.
𝑑ℎ 𝑑ℎ
𝐶𝑣𝐻 = 0.85032 − 0.04231 + 0.0017954( )2 (17)
𝑡𝑡 𝑡𝑡
Table 1 - Case study methanol/water data.

ℎ𝑑,𝑗 𝐻 (m), is then given by: Feed flowrate (kg/s) 6.05

𝐻 𝜌𝑤𝑎𝑡𝑒𝑟 Methanol–20
ℎ𝑑,𝑗 = 5.1204 × 10 −5
𝑣ℎ,𝑗 2 𝜌𝑗𝑉 [1 Feed composition (mass %)
𝜌𝑗𝐿 (18) Water–80
2
𝐴ℎ Feed temperature (°C) 95
− ( ) ] 1⁄ 𝐻 2 , 𝑗 = 1, … 𝐽.
𝐴𝑎 𝐶𝑣 Feed pressure (bar) 1
The pressure drop through the aerated mass ℎ𝑙 was Trays 13
Feed tray 6
calculated using the Fair correlation and the Bennett Fractional active area 0.82
aerated liquid pressure drop correlation. In the Fair Fractional hole area 0.128
correlation the pressure drop through the aerated mass Hole diameter (mm) 4.5
on stage 𝑗, ℎ𝑙,𝑗 𝐹 , is calculated using an aeration factor 𝛽: Hole diameter/pitch length 0.376
Hole diameter/tray thickness 2.3
ℎ𝑙,𝑗 𝐹 = 𝛽𝑗 ℎ𝑑𝑠,𝑗 , 𝑗 = 1, … 𝐽. (19) Plate spacing (m) 0.6
Weir type Segmental
𝑞𝑗 0.5
Weir length/column diameter 0.7
𝛽𝑗 = 0.0825 𝑙𝑛 − 0.269 ln (𝑣ℎ,𝑗 𝜌𝑗𝑉 ) + 1.679, Weir height (m) 0.04
𝑙𝑤 (20)
𝑗 = 1, … 𝐽.

In the above equation, ℎ𝑑𝑠,𝑗 is the sum of the weir height Figure 2 shows the liquid and vapour flowrate profiles
ℎ𝑤 and the height of the crest over the weir ℎ𝑜𝑤,𝑗 . ℎ𝑜𝑤,𝑗 is
for the present case study. It is interesting to note that
calculated both for segmental (section 3.1) and circular the liquid and vapour flowrates remain relatively
(section 3.2) weirs using the Francis weir equation and constant above and below feed tray 6, where a regime
its modification for different weir types. For a segmental transition is observed. Figure 2 also shows the froth-
weir and for height ℎ𝑜𝑤,𝑗 𝑆 in meters of liquid: spray correlation presented by Porter & Jenkins, where
the flow parameter 𝐹𝑙𝑣 , the active area 𝐴𝑎 (m2) and the
𝑞𝑗 2 weir length 𝑙𝑤 (m) are used to determine the froth-spray
ℎ𝑜𝑤,𝑗 𝑆 = 0.664( )3 , 𝑗 = 1, … 𝐽. (21)
𝑙𝑤 regime limit, as follows:

The modified Francis equation for circular weirs was 𝐿𝑗 𝜌𝑗𝑉 0.5
𝐹𝑙𝑣,𝑗 = ( ) , 𝑗 = 1, … 𝐽. (26)
used to determine ℎ𝑜𝑤,𝑗 (m), as follows: 𝑉𝑗 𝜌𝑗𝐿

𝑞𝑗
ℎ𝑜𝑤,𝑗 = 44.3(
𝑑𝑤 × 1000
)0.704 , 𝑗 = 1, … 𝐽. (22) 𝐹𝑙𝑣,𝑗 𝐴𝑎
= 0.07, 𝑗 = 1, … 𝐽. (27)
𝑙𝑤

where 𝑑𝑤 (m) refers to the weir diameter and 𝑞𝑗 (m3/s) 𝐹 𝐴

In the rectifying section 𝑙𝑣,𝑗 𝑎
< 0.07 and therefore the
is the volumetric liquid flow. In the Bennett aerated 𝑙𝑤

liquid pressure drop correlation the pressure drop in column is operating under a spray regime. Indeed, a
spray-type regime is attained under high gas rates and
stage 𝑗 through the aerated mass ℎ𝑙,𝑗 𝐵 is the sum of the
low liquid rates, where jets are formed at each
effective clear liquid height ℎ𝑐,𝑗 𝐵 (calculated according to perforation, with the fraction of the holes that are
the Bennett effective clear liquid height correlation) and jetting increasing with vapour velocity. On the contrary,
the residual pressure drop ℎ𝑅,𝑗 due to surface tension 𝜎𝑗 in the stripping section the inequality 𝑙𝑣,𝑗
𝐹 𝐴𝑎
> 0.07 is
𝑙𝑤
(N/m): observed and consequently the column is operating
under a froth-type regime. Here, the dispersion on the
ℎ𝑙,𝑗 𝐵 = ℎ𝑐,𝑗 𝐵 + ℎ𝑅,𝑗 , 𝑗 = 1, … 𝐽. (23) tray is in the form of a bubbly, or aerated liquid.

4
 8 0.1
7
0.08
Flowrate (kg/s)

Flv.Aa/Lw (m)
6
5 0.06
4
3 0.04
2 0.02
1 L
0 V 0
1 2 3 4 5 6 7 8 9 10 11 12 13 1 2 3 4 5 6 7 8 9 10 11 12 13
Trays Trays
Figure 2 - Liquid and vapour flowrate profiles and Porter & Jenkins froth-spray prediction.

3.1.1. Tray holdup the fractional hole area, the tray thickness and the sieve
Figure 3 shows the liquid molar holdup profiles tray hole diameter. Nonetheless, for this particular case
calculated using Bennett’s and Jeronimo & study using the simpler standard orifice equation would
Sawistowski’s clear liquid height correlations.
suffice. Indeed, from Figure 4 it can be noted that a
In the rectifying section the column is operating under a
similar and constant dry vapour head loss profile is
spray regime and therefore the Jeronimo & Sawistowski
predicted by the standard orifice equation, the
should be used if precise liquid holdups are to be
Hughmark & O’Connell correlation and the Liebson
determined. Nevertheless, in this particular case the two
correlation for the two distinct column sections.
correlations predict a similar and approximately
Because in the rectifying section the column is operating
constant liquid holdup with a maximum relative
under a spray-type regime with high vapour velocities,
difference of 6.9% for tray 1. In the stripping section a
the predicted vapour head loss is significantly higher
froth type regime is observed and thus the Jeronimo &
than in the stripping section, where a froth-type regime
Sawistowski correlation is not applicable. In this case,
with lower vapour rates is established. For example, an
the Bennett correlation should be used. Note that both
average dry vapour head loss of 65.7 mm is given by the
correlations estimate an approximately constant liquid
standard orifice equation for the rectifying section,
molar holdup profile for this lower column section.
whereas this value drops to 10.2 mm in the bottom
However, using the correct correlation is crucial since
column stages, i.e., an 84.5% reduction in the dry
the Bennett holdup correlation generates a significantly
vapour head loss is observed.
high average liquid molar holdup of 2517 moles
An approximately constant aerated liquid head loss
whereas the Jeronimo & Sawistowski correlation
profile is also predicted for the two column sections
estimates an average liquid holdup of only 1658 moles.
according to the Fair and the Bennett correlations.
Here, the effect of the liquid flowrate on the aerated
4000
liquid head loss is noticeable: for the rectifying section
Liquid holdup (mol)

Bennett
3000 Jeronimo & Sawistowski the liquid flowrates are significantly lower and therefore
for both correlations the liquid head loss is almost half
2000 of the value calculated for the stripping section.
1000 Theoretically the Fair correlation should estimate liquid
head losses more accurately in the upper column section
0 since using the Bennett correlation is not recommended
1 2 3 4 5 6 7 8 9 10 11 12 13
Trays for spray-type regimes. However, for the rectifying
Figure 3 – Liquid molar holdup profile. section, an average aerated liquid head loss of 17.0 mm
is given by the Bennett correlation resulting in only a
7.0% relative difference from the value predicted by the
3.1.2. Tray pressure drop
Fair correlation. Thus, for this particular case study
The Hughmark & O’Connell and the Liebson
using either correlation would not compromise the
correlations should theoretically predict a more accurate
liquid head loss profile prediction in the upper column
dry vapour head loss than the standard orifice equation,
section. A larger discrepancy can be noted for the
since the latter assumes a constant orifice coefficient of
bottom column stages where the Fair correlation
0.75 while in the first two cases this same orifice
estimates an average liquid head loss of 34.1 mm, with a
coefficient is estimated using geometrical inputs such as

5
 80 40
70 35
60 30

h aerated (mm)
h vapour (mm)

50 25
40 Standard orifice equation 20 Bennett
30 Hughmark & O'Connell 15 Fair
20 Liebson 10
10 5
0 0
1 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7 8 9 10 11 12
Trays Trays

Figure 4 – Dry vapour and aerated liquid head loss profiles.

16.5% relative difference from the 29.2 mm value relative volatility of 4.02 and in Batchsim using the
predicted by the Bennett correlation. In this case, NRTL/Virial thermodynamic models.
although both estimations apply to froth-type regimes,
the Bennett correlation should predict liquid head losses Table 2 - Fortran and Batchsim model specifications.
more accurately. Indeed, unlike the Fair correlation the Specification Fortran/Batchsim
Bennett aerated liquid loss correlation does not correct Initial reboiler charge (mol) 184.32
the clear liquid flow for aeration effects by including a Initial reboiler charge composition Methanol – 20
tray aeration factor. (mol%) Water - 80
Column diameter (mm) 50
Number of column stages 15
3.2. Batch distillation Tray 1 pressure (bar) 1.013
This section shows the implementation of a multi- Fortran: α = 4.02
Vapour-liquid equilibrium
staged batch distillation system in ModelBuilder based Batchsim:NRTL/Virial
on a pilot plant experiment lead by Bonsfills [12]. Stage modelling Equilibrium stages
Column liquid flowrate Constant
Bonsfills carried out a pilot plant experiment for a Column vapour flowrate Constant
methanol/water batch separation with a 15 trays Tray liquid holdup Constant
Oldershaw column of 3.75 m height and 50 mm inner Tray vapour holdup Negligible
diameter. The column worked with an average absolute Average holdup per stage (mol) 0.175
Column heat loss Adiabatic
pressure of 1 atm. To reduce heat losses, an insulation Reflux ratio 3
jacket was added. A 4 L 20 mol% methanol/80 mol% Condenser type Total
water mixture was charged to a 6 L reboiler with a Condenser holdup (mol) 0.175
maximum heat duty of 1400 W. The distillation still was
heated up with a constant heat duty of 933.3 W under 3.2.1. ModelBuilder flowsheet
total reflux until steady-state conditions were attained. Figure 5 shows the flowsheet used to simulate the
Then, the column operated under a finite constant multi-staged batch system in ModelBuilder. A
reflux ratio regime with product being cooled in a total methanol/water mixture contained in reboiler E-114 is
condenser and withdrawn at the top of the column. heated using a constant heat duty of 736.4 W, which is
Bonsfills compared the experimental results with the theoretical duty required to generate an average
simulations obtained from a batch distillation model distillate flowrate of 0.28 mol/min. The fact that the
implemented in Fortran language and with results given experimental Oldershaw column is not adiabatic may
by commercial simulator Batchsim of Pro/II. In this explain this discrepancy between the ModelBuilder heat
work the Fortran and Batchsim models will also be used duty of 736.4 W and the experimental value of 933.3 W.
to validate the ModelBuilder column batch model. Indeed, if a 21.1% heat loss equivalent to a 196.9 W
Table 2 shows the specifications used in the Fortran and heat loss is assumed for the Bonsfills experiment, the
Batchsim models. Both models assume 15 theoretical distillate flowrate profiles become identical and the
trays with constant liquid and negligible vapour validity of the ModelBuilder column model may be
holdups. However, the vapour-liquid equilibrium is assessed.
determined in the Fortran model using a constant

6
Vapour is formed in the reboiler and fed to the bottom a negligible pressure drop are maintained. For the
of column D-110. A vapour-liquid equilibrium process ModelBuilder simulations the dynamic mode of the
occurs on all column stages, with liquid flowing from condenser was turned off which in practice means that a
the top of the column to the bottom and vapour negligible holdup was assumed in E-121.
ascending throughout the column section. The top For a given reboiler heat duty, two distinct
section of the column is enriched with the lighter ModelBuilder simulations are scheduled. First, a
component methanol as the bottom section simulation is designed to achieve steady-state conditions
composition in water increases. The liquid leaving the by operating the column under total reflux. In this case,
bottom of the column enters the reboiler, closing the a reflux ratio of 𝑅𝑅 = 𝐷𝐿 = ∞ is specified via splitter M-
bottom loop. 123 and no distillate is withdrawn to product receiver F-
140. Note that the reflux ratio is perfectly controlled by
specifying the required split fraction and un-assigning
the pump pressure increase. For a particular reflux
flowrate, reflux valve V-124 generates a given pressure
drop which is perfectly compensated by pump L-126,
where the required pressure increase is established. All
the initial conditions are given at time 𝑡 = 0 in the first
total reflux simulation. The second simulation uses a
particular split fraction based on the desired finite reflux
ratio. This simulation starts operating in a steady-state
mode, using a saved variable set from the first
simulation, that is, all the values assigned to the
variables at the end of the first simulation (when steady-
state is achieved) are saved and used to initialise the
second simulation.
Globally the system is pressure-driven, that is, the mass
flowrates are calculated in the valve model based on
differential pressures. For this reason, valves were
introduced before and after each equipment.
The input specifications for the ModelBuilder dynamic
model are listed in Table 3 and Table 4.
Figure 5 – Batch distillation ModelBuilder flowsheet.
Table 3 – Column specification data.
H-119 represents the column sump, i.e., the liquid
Number of stages : 15
holdup that is held on the bottom of the distillation Column
Dummy feed stage: 5
column, below the column stages. A cascade control Tray type: sieve
scheme was implemented using the level and the flow Column diameter: 50 mm
Active area fraction: 0.8
controllers to maintain a constant ratio of liquid level to Design Hole area fraction: 0.1
sump height of 0.6. In this arrangement, the constant Weir fraction: 0.7
liquid level of 3.1 cm is maintained in the sump by Weir height: 4.46 mm
varying the sump outlet flowrate, instead of directly Plate spacing: 0.25 m
Dry vapour pressure drop correlation:
changing the sump outlet valve V-117 stem position, standard orifice equation
which would be the case in a simple (non-cascade) Aerated liquid pressure drop correlation:
Pressure
control scheme. In the primary loop, the level controller hydrostatic
Liquid height correlation:
monitors and compares the liquid level fraction to the
Francis formula for circular weirs
specified set-point of 0.6, changing the outlet sump Mode: pressure-driven
flowrate set-point accordingly. In the secondary loop, Initial liquid level fraction trays 1-15: 0.02185
the flow controller reads and compares the flowrate Dynamics Initial pressure trays 1-15: 1.113 to 1.2205
Initial holdup composition trays 1-15: 80.8 to
input value to the flowrate set-point, manipulating the 64.3 methanol mol%
sump outlet valve stem position as necessary. Efficiency Tray modelling: equilibrium
The vapour distillate leaving the upper section of the
column is condensed in total condenser E-121 where a
constant sub-cooling temperature difference of 5ºC and

7
 Table 4 – Sump, reboiler and product receiver specifications. plateau ends when the methanol purity falls below 0.99.
Because Bonsfills obtained an experimental plateau time
Sump Reboiler Product receiver of 80 min, ModelBuilder gives the closest estimation to
this value, with an error of 17.2%. The methanol plateau
Cylindrical vessel Cylindrical vessel Cylindrical vessel duration is directly related to the distillate flowrate
with flat heads with flat heads with flat heads
D = 5 cm D = 19.5 cm D = 21.7 cm profile since higher distillate flowrates in this period
V = 0.1 L V=6L V=8L mean that more pure methanol is being withdrawn per
unit of time and consequently shorter first plateau
Initial holdup Initial holdup Initial holdup durations are obtained. This fact explains why the
-Methanol=0.312 -Methanol=19.9 -Methanol≈0 mol
mol mol% -Water≈0 mol
ModelBuilder estimation gives the closest
-Water=4.440 mol -Water=80.1 -T=66.3ºC approximation to the experimental methanol plateau
-T=76.9ºC mol% duration: both the experimental setting and the
Control -Liquid Liquid
level fraction=0.6 ModelBuilder model use a constant heat duty; because
-PI controllers surface pressure
-Cascade control -P=1.221 bar -P=0.9999 bar methanol has a lower heat of vaporization than water,
-Set-point LC- to maintain an average distillate flowrate of 0.28
119=0.6 liquid mol/min, the experimental and ModelBuilder distillate
level fraction
flowrate profiles must reveal higher flowrates in the
pure methanol withdrawal period (>0.28 mol/min)
3.2.2. Model validation
when compared to the pure water withdrawal period
The left-hand side of Figure 6 shows the molar
(<0.28 mol/min) and thus shorter methanol periods are
methanol fraction profile in the distillate for the
to be expected in these cases. Both the Fortran and
ModelBuilder simulation with a constant reflux ratio
Batchsim models assume a constant distillate flowrate
𝑅𝑅 = 3 and a reboiler heat duty of 736.4 W. The same
of 0.28 mol/min and thus the pure methanol
profile is also illustrated for the Bonsfills experiment,
withdrawal period lasts longer.
Fortran model and the Batchsim model for a reflux
The second phase of the distillation process generates
ratio of 3 and for the heat duty required to generate an
off-cuts. The shifted transformations show that the
average distillate flowrate of 0.28 mol/min. Globally
ModelBuilder and the Batchsim slopes are almost
three distinct time periods can be noted: a first period
parallel in this second time period, predicting very sharp
where a plateau occurs for a molar methanol fraction of
decreases in the methanol composition, i.e., efficient
1; a second phase where the methanol molar fraction
separations linked to short off-cut periods. The models
decreases and a final period where pure water is
used to determine the VLE data might explain the off-
withdrawn from the column.
cut slopes: both the ModelBuilder and the Batchsim
During the first period pure methanol is withdrawn
models use the NRTL equations to determine the
from the column and collected in the product receiver.
activity coefficients of the liquid phase. The Bonsfills
The estimated duration of this first plateau is 94 min for
experiment gives a significantly longer off-cut period
the ModelBuilder model, 105 min for the Fortran model
(>174 min).
and 112 min for the Batchsim model. The first plateau
duration is defined by the methanol purity, i.e., the first
1 1
0.9 0.9
Fortran Fortran
0.8 0.8 ModelBuilder
y methanol (mol frac)
y methanol (mol frac)

ModelBuilder
0.7 Batchsim 0.7 Batchsim
Experimental
0.6 Experimental 0.6
0.5 0.5
0.4 0.4
0.3 0.3
0.2 0.2
0.1 0.1
0 0
0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5
Time (hr) Time (hr)

Figure 6 - Distillate methanol composition profiles with RR=3 and 𝐃 ̅ = 0.28 mol/min. The right-hand side figure
represents the same profiles with an xx axis shift for each of the 3 functions to G(X) = F(X+T1stplateau).

8
The assumption of equilibrium stages for all the models 3.2.3. Sensitivity analysis
may have affected the distillate profiles generating The relative importance of key operating conditions was
optimistic (short) off-cut periods and efficient assessed by performing a sensitivity analysis. The effect
separations. The three models predict a higher methanol of the reflux ratio and the heat duty on the total
recovery when compared to the experimental value of separation time and on the methanol recovery was
79.7%: 85.4% for ModelBuilder, 86.7% for Fortran and analysed in a total of 9 simulations with ±0.5 reflux
92.7% for Batchsim. ratio differences and ±10% heat duty relative variations
Now that the behavior of the ModelBuilder simulations starting from the base case scenario with a reflux ratio
has been analysed for the distillate stream, let us of 3 and 736 W heat duty.
consider an intermediate column plate 5. Figure 7 shows Figure 9 shows the total separation time (TST) for all
that the three time periods analysed for the distillate simulations. The total separation time starts in the
profile are still very clearly defined for plate 5 with pure steady-state regime and ends when the reboiler is run
methanol being obtained during 81 min for the dry, i.e., it is the period of time where a constant finite
ModelBuilder simulation, 86 min for Fortran and 103 reflux separation is established. Higher distillate
min according to Batchsim. As expected, the period flowrates generate a shorter total separation time since
where pure methanol is obtained in the distillate is more product is withdrawn per unit of time. The
longer than for plate 5. However, for plate 5 the highest average distillate flowrate of 0.35 mol/min is
duration of the phase where pure water is obtained obtained combining the lowest reflux ratio of 2.5 with
increases. As the stage column number increases the the highest heat duty of 810 W. On the contrary, the
methanol is held during shorter times while the water is lowest average distillate flowrate of 0.22 mol/min is
retained for longer periods. generated for the highest reflux ratio and the lowest
heat duty combination (3.5 and 663 W, respectively).
1
Consequently, the shortest total separation time of 7.8
x methanol (mol frac)

0.8 hours is given by the 2.5 RR – 810 W simulation

0.6 Fortran
whereas the longest total separation time of 12.1 hours
ModelBuilder is predicted by the 3.5 RR – 663 W simulation. Another
0.4 Batchsim conclusion that may be drawn from Figure 9 is that
0.2 ±10% heat duty variations with fixed reflux ratios and
0 ±0.5 reflux ratio differences with fixed heat duties
0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 affect the total experiment time approximately 12%.
Time (hr)

Figure 7 - Methanol composition for plate 5 with RR=3 and

̅ = 0.28 mol/min. Q=810 W
𝐃
Figure 8 compares the temperature profiles obtained for
plate 5. The same three time periods can be noted: in
Q=736 W
the first plateau pure methanol is distilled and the
RR=3.5
temperature given by all the models is 64.7 °C, which is
RR=3
equivalent to methanol’s boiling point; the temperature Q=663 W RR=2.5
increases and the second phase ends when the distillate
temperature reaches 100 °C, when the process of 6 8 10 12 14
withdrawing pure water from the column starts. TST (hr)

110
100 Figure 9 – Total separation time (TST) for Q = 663 W, 736 W
Temperature (°C)

and 810 W and RR = 2.5, 3 and 3.5.

90
80 Table 5 shows the methanol recoveries for all the
Fortran
70 ModelBuilder simulations under study. The methanol recoveries are
Batchsim affected by ±0.5 differences in the reflux ratio. For
60 Experimental
instance, a 6.8% relative decrease and a 2.6% relative
50 increase in the methanol recovery is obtained using the
0 1 2 3 4 5 6 7 8 9 10 11
Time (hr) 2.5 RR – 736 W and the 3.5 RR – 736 W simulations,
̅ = 0.28
Figure 8 - Plate 5 temperature with RR=3 and 𝐃 respectively, when compared to the base case 3.0 RR –
mol/min. 736 W simulation. On the contrary, a ±10% heat duty

9
variation practically does not modify the methanol average distillate flowrate of 0.28 mol/min. It was
recovery since the off-cut slopes were found to not be found that the assumption of equilibrium stages
affected and hence the separation efficiency was generated for all models a non-negligible overestimation
maintained. This fact is demonstrated, for example, in of the separation efficiency and thus the methanol
the fixed RR = 3 simulations where a negligible recovery. The experimental methanol recovery of 80%
methanol recovery relative variation of -0.7% is given was thereby lower than the ModelBuilder prediction of
for the 3.0 RR – 663 W simulation and of -0.1% for the 85%, the 87% Fortran estimation and the Batchsim
3.0 RR – 810 W simulation, when compared to the 3.0 93% methanol recovery.
RR – 736 W base case. In conclusion, the methanol The robustness of the dynamic ModelBuilder model
recovery varies between 2.2 and 6.4% in absolute values was attested by performing a sensitivity analysis with
when fixing the heat duties and modifying the reflux 10% heat duty variations for the reflux ratios of 2.5, 3
ratio in ±0.5 increments; on the contrary, a maximum and 3.5. It was shown that ±10% heat duty variations
absolute change of only 1.2% is observed in the with a fixed reflux ratio and ±0.5 reflux ratio
methanol recovery when the reflux ratio is maintained differences with a fixed heat duty affect the total
and the heat duties are varied ±10%. experiment time in approximately 12%. Furthermore,
±0.5 variations in the reflux ratio were found to have an
Table 5 – Methanol recovery for Q = 663 W, 736 W and 810 absolute impact of 2.2 to 6.4% on the methanol
W and RR = 2.5, 3 and 3.5.
recovery. On the contrary, this recovery was practically
Methanol not affected by the 10% heat duty variation.
recovery (%) 663 W 736 W 810 W
Works cited
RR=2.5 80.4 79.6 78.9 [1] W. Luyben, “Multicomponent batch distillation,”
RR=3 84.8 85.4 85.4 vol. 27, 1988.
RR=3.5 88.7 87.6 88.8 [2] M. Barolo and F. Botteon, “Simple method of
obtaining pure products by batch distillation,” vol.
4. Conclusions 43, 1997.
The demand for productivity in the fine chemical [3] U. Diwekar and K. Madhavan, “Multicomponent
industry is increasing and therefore robust models that batch distillation column design,” vol. 30, 1991.
can predict separation performance indicators such as [4] S. Sundaram and L. Evans, “Shortcut procedure for
operation time and component purity are required. The simulating batch distillation operations,” vol. 32,
1993.
main contribution of this study has been the
[5] E. Meadows, “Chem. Eng. Progr. Symp. Ser. No.
implementation of a dynamic multi-staged distillation
46,” vol. 59, 1963.
model in gPROMS ModelBuilder® based on
[6] G. Distefano, “Mathematical modelling and
equilibrium stages and including the effects of non-
numerical integration of multicomponent batch
negligible tray holdup and of tray pressure drop. distillation equations,” vol. 14, 1968.
Tray holdup and pressure drop equations were [7] R. H. Perry and D. W. Green, Perry's Chemical
incorporated in the rigorous dynamic multi-staged Engineers' Handbook, McGraw-Hill, 1999.
model to ensure a pressure-driven system with non- [8] H. Kister, Distillation Design, New York: McGraw-
negligible tray holdups. A methanol/water continuous Hill, 1992.
distillation case study with suitable geometrical [9] J. Seader, E. Henley and K. Roper, Separation
specifications was used to preliminarily test the Process Principles, New York: John Wiley & Sons,
implemented hydraulic correlations. Here, the type of 2011.
flow regime was found to significantly affect the liquid [10] “Process System Enterprise Limited,” 1997.
holdup, dry vapour head loss and aerated liquid head [Online]. Available:
loss profiles. After testing the hydraulic relations for the http://www.psenterprise.com/gproms/platform.
methanol/water continuous distillation case study, a html. [Accessed 8 September 2014].
multi-staged dynamic distillation model based on an [11] “KBC Advanced Technologies,” 2014. [Online].
experimental pilot plant assembled by Bonsfills [12] was Available: http://www.kbcat.com/infochem-
software/flow-assurance-software-multiflash.
implemented in ModelBuilder for a methanol/water
[Accessed 08 September 2014].
mixture. Model validation was successfully carried out
[12] A. Bonsfills and L. Puigjaner, “Batch distillation:
by comparing the ModelBuilder simulations with
simulation and experimental validation,” Chemical
experimental data and theoretical results given by Engineering and Processing, vol. 43, pp. 1239-1252,
Bonsfills for a Fortran model and for simulator 2004.
Batchsim of Pro/II for a reflux ratio of 3 and an

10