Received: 15 May 2020 Revised: 19 June 2020 Accepted: 22 June 2020

DOI: 10.1002/cjce.23857

MINI-REVIEW

Experimental methods in chemical engineering: Process

simulation

Jacopo De Tommaso1 | Francesco Rossi2 | Nooshin Moradi3 | Carlo Pirola4 |

Gregory S. Patience1 | Federico Galli1

1
Chemical Engineering, Polytechnique
Montréal, Montréal, Québec, Canada
Abstract
2
School of Chemical Engineering, Purdue Process simulation software designs equipment, simulates operations, opti-
University, Davidson School of Chemical mizes a plant’s configuration (heat exchangers network, for example), esti-
Engineering, Forney Hall of Chemical
mates operating and capital expenses, and serves as an educational tool.
Engineering, West Lafayette, Indiana
3 However, mastering the theoretical background minimizes common mistakes
School of Chemical and Petroleum
Engineering, Shiraz University, such as applying an incorrect thermodynamic method, selecting improper
Shiraz, Iran algorithms in the case of tear systems, and setting irrational system specifica-
4
Dipartimento di Chimica, Universitá tions. Engineers and researchers will exploit this tool more often in the future
degli Studi di Milano, Milan, Italy
as constant advancements in simulation science as well as new models are
Correspondence released continually. Process simulators make it easier to build digital twins
Federico Galli, Chemical Engineering,
and thus will facilitate the implementation of the industry 4.0 guidelines. We
Polytechnique Montréal, C.P. 6079, Succ.
“CV”, Montréal, H3C 3A7 Québec, highlight the mathematical and technical features of process simulators, as
Canada. well as the capabilities and the fields of application. A bibliometric map of key-
Email: federico.galli@polymtl.ca
words from articles citing Aspen+, Aspen plus, Hysys, and Pro/II indexed by
Web of Science between 2017 and 2020 identified the main research clusters,
such as design, optimization, energy or exergy, biomass; H2 and CO2 capture,
thermodynamics; and separations and techno-economic analysis.

KEYWORDS
cost estimation, design, optimization, process simulation, thermodynamics

1 | INTRODUCTION University of Houston developed CHESS in collaboration

with industry in 1968. It simulated the phase equilibria of
Process simulators solve mathematical equations not only more than 70 hydrocarbons.[1] PROCESS, by Simulation
to characterize reactors, distillation columns, heat Science, was distributed in 1966 and simulated distilla-
exchangers, pumps, compressors, and other unit opera- tion columns. It then became PRO/II, first distributed by
tions but also for process design and optimization. Fur- Scheider-Electrics and then by AVEVA. Fortran
thermore, they offer engineers a tool to estimate plant expanded the capabilities of these programs in the 1970s.
costs, both operating expenses (OPEX) and capital In the late 1970s, MIT developed Advanced System for
expense (CAPEX). Users input operating pressure, tem- Process ENgineering (ASPEN) that AspenTech commer-
perature, flow rates, composition, and thermodynamic cialized in the 1980s. Concurrently Profs. Bishnoi and
packages and the software computes mass and energy Svercek at the University of Calgary developed HYSYS.
balances for each stream and module. Examples of the In the late 1980s and early 1990s, PC-based simulators
first software date back to 1960s with PACER developed were developed.[2] In 2002, AspenTech acquired
by McMaster University for educational purposes. The Hyprotech; however, the Federal Trade Commission

Can J Chem Eng. 2020;98:2301–2320. wileyonlinelibrary.com/journal/cjce © 2020 Canadian Society for Chemical Engineering 2301
2302 DE TOMMASO ET AL.

0
required AspenTech to divest HYSYS because they g1 ðx 1 , x 2 ,x 3 , …,x N Þ = 0
deemed it contravened anti-trust laws. Honeywell B g ðx , x ,x , …,x Þ = 0
B 2 1 2 3 N
became the owner of HYSYS and created UniSim. Even- B , g ðx Þ = 0 ð1Þ
@ 
tually, AspenTech re-acquired HYSYS in 2016 and now
gN ðx 1 , x 2 , x 3 , …, x N Þ = 0
91 simulators are on the market, including open source
software.
This mini-review on process simulation is part of a The objective is to solve this nonlinear system (xs),
series of articles dedicated to experimental methods in which corresponds to a root of all the functions gi
chemical engineering.[3] Although practicing engineers (gi(xs) = 0). There exist no analytical method to accom-
and students use these simulators for plant design and plish this task, so we only estimate an approximate
optimization, researchers also apply them to calculate value of xs through numerical algorithms, Newton
thermodynamics, identify equilibrium compositions of type and fixed-point methods. Both these iterative pro-
reactions, verify the feasibility of separation operations, cedures rely on a user-supplied initial value, x0, and
and thus optimize experimental designs. Here we discuss iteration (xn)n ≥ 1 until they approach the solution vec-
these features but address the mathematical approaches tor xs. These two families of algorithms differ with
to solve these large systems of equations, resolution algo- respect to the procedure they apply to compute
rithms, and applications. We describe the main sources (xn)n ≥ 1, their efficiency and robustness, and their
of error and how to avoid common mistakes that arise at application domain.
the set up of the flowsheet and include a bibliometric sur- Starting from the last known iteration, xn, Newton-
vey that highlights the major applications. type methods calculate xn + 1 with the following
procedure:

2 | THEORY 1. Exact or approximate linearization of the nonlinear

systems at xn. This step requires calculation of the
Steady-state models of most unit operations are Jacobian matrix at xn or of an estimate thereof (the
nonlinear systems of either algebraic or differential- matrices J(xn) and Bn, respectively).
algebraic equations (eg, plug flow reactors and pipes). 2. Solution of the linearized system of equations and cal-
Discretization techniques convert differential-algebraic culation of the search direction dn, according to the
equations into a system of algebraic equations; there- first line of Equation (2) and Equation (3).
fore, any steady-state process flowsheet is equivalent to 3. Computation of the optimal step length an through a
a large-scale system of nonlinear equations. Therefore, line search strategy, performed along dn.[6]
in this section of the article, we first review the ratio- 4. Calculation of xn + 1 from xn, an, and dn, according to
nale of methods that solve these systems of nonlinear the second line of Equation (2) and Equation (3).
equations, and then explain how simulation packages
use them to solve process flowsheets. We discuss both
0 1
modular and equation-oriented strategies for flowsheet ∂g1 ∂g1 ∂g1
B ∂x 1  
simulation, comment on their advantages and disad- B ∂x 2 ∂x N C
C
B ∂g ∂g2 ∂g2 C
vantages, and provide guidelines on when/how to  B 2 C
J ðx n Þdn = − gðx n Þ B   C
apply each approach.[4,5] ^J ðx Þ = B ∂x ∂x 2 ∂x N C ð2Þ
B 1 C
x n + 1 = x n + a n dn B .. .. .. .. C
B . . . . C
B C
@ ∂gN ∂gN ∂gN A
2.1 | Numerical methods for systems of  
∂x 1 ∂x 2 ∂x N
nonlinear equations

Bn dn = −gðx n Þ
The most general mathematical formulation of a system ^BðnÞ≈J ðx n Þ ð3Þ
of nonlinear equations, Equation (1), where every func- x n + 1 = x n + a n dn
tion gi represents a nonlinear expression of the unknown
variables xi (in a process flowsheet, the functions gi repre- These steps are repeated until convergence criteria
sent steady-state mass, energy, and momentum balances, are met, for example, the magnitude of the search direc-
and discretization thereof or design specifications, tion (kdnk) becomes sufficiently small, or the sequence
whereas the quantities xi are process variables, for exam- (xn)n ≥ 1 diverges (Figure 1). The sequence (xn)n ≥ 1 fol-
ple, temperatures, pressures, molar/mass fraction, s and lows the same type of pattern (Figure 1) for systems of
flowrates) is as follows: any number of equations.
DE TOMMASO ET AL. 2303

(A) (A)

(B) (B)

FIGURE 1 Typical sequence of iterations generated by a F I G U R E 2 Typical sequence of iterations generated by direct
Newton-type method: A, convergent sequence; and B, divergent substitution: A, convergent sequence; and B, divergent sequence.
sequence. From the initial guess x0, the algorithm calculates the From the initial guess x0, the algorithm calculates xi and f(xi),
direction vector di and the step length ai, moving towards xS or moving towards xS or without converging
without converging

Newton-type methods do not guarantee conver- After this preliminary step, starting from the last
gence to a solution of the nonlinear system, starting known iteration xn, the next iteration xn + 1 is calculated
from any random initial guess as this family of algo- through the following procedure:
rithms is only locally convergent. However, their rate
of convergence is superlinear. Therefore, they are effec- 1. Calculation of the functions fi at xn.
tive when an adequate initial guess is available. 2. Selection of values for all the relaxation factors ai,n.
Another feature of Newton-type methods is their capa- 3. Application of Equation (5).
bility of solving nonlinear systems of arbitrarily large
scale, provided that adequate computing power is x i,n + 1 = ai,n f i ðx n Þ + ð1 −ai,n Þx i,n 8i½1, N  ð5Þ
available.
Fixed-point methods require a reformulation of the These steps are repeated until convergence criteria
original system of nonlinear equations, which consists are met, for example, the norm of the difference between
of re-writing every equation gi(x) = 0 in the same form two consecutive iterations (kxn + 1 − xnk) becomes suffi-
as Equation (4) (this reformulation is always possible, ciently small, or the sequence ((xn)n ≥ 1 diverges (the
although not necessarily unique, and may affect the direct substitution recursion matches Equation (5) with
convergence properties of the algorithm): all the relaxation factors ai,n = 1). The sequence (xn)n ≥ 1
follows a similar pattern for systems of two or more equa-
x i = f i ðx Þ 8i½1, N  ð4Þ tions (Figure 2).
2304 DE TOMMASO ET AL.

As with Newton-type methods, fixed-point schemes 2.3 | Modular solution (SM) strategies
do not guarantee convergence to a solution of the for process flowsheets
nonlinear system, starting from any random initial
guess (these algorithms are only locally convergent), Any chemical process is a set of unit operations, primarily
but they are usually less sensitive to the initial value. connected in series as the number of recycles is usually
On the other hand, their rate of convergence is slower significantly smaller than the total number of connections.
than Newton-type methods. Thus, they are particularly Modular solution strategies for process simulation leverage
suitable for situations in which no adequate initial this sequential nature of chemical processes to first sim-
guess can be computed. Note that, as opposed to plify and then perform flowsheet calculations. The ratio-
Newton-type algorithms, fixed-point methods can only nale of these solution strategies is as follows:
solve small-scale/medium-scale nonlinear systems,
regardless of the amount of computing power avail- 1. Given a certain process flow diagram (PFD), we first
able. This limitation considerably restricts their appli- convert it into a directed graph (the simulation flow
cation domain. diagram [SFD]), whose vertices and edges represent
units and streams, respectively (Figure 3).
2. Then, well-established algorithms, for example, the
2.2 | Computational strategy connection matrix method[13] or the depth-first search
and backtrack method,[14] detect all the simple cycles
The barebones of a flowsheeting software differs from in the flowsheet graph (Figure 3), which corresponds
those of others based on the computational strategies. to identifying all the recycles in the original process
The main techniques are: flowsheet.
3. Subsequently, we break all the simple cycles by tear-
• Sequential modular approach (SM). ing a number of edges in the flowsheet graph
• Equation oriented approach (EO). (Figure 4), thus transforming the latter into a
• Simultaneous modular approach. directed tree (the flowsheet tree). This graphical oper-
ation corresponds to breaking all the recycles in the
The SM solves the process units in sequence, starting original process flowsheet, by replacing specific pro-
from the feed and tear the common streams in case of cess streams (the tear streams) with pairs of new
recycle. Where an inlet stream is given, each block com- streams (the artificial streams), in order to generate a
putes its outlet stream. However, in the presence of a simplified process flowsheet, in which all the unit
recycle, an iteration is required. Sturdy and reliable, this operations can be solved in series once half the artifi-
is the default method of commercial and general soft- cial streams have been assigned (the remaining artifi-
ware, especially for steady-state calculations. On the cial streams become simulation outputs). Note that
other hand, the SM is weak and time-consuming for the every recycle stream is literally torn at least once,
following[7]: thus the name tear stream.
4. Finally, we exploit the connectivity information,
• Highly recycled processes. stored in the flowsheet tree, to solve all the unit opera-
• Highly-heat integrated processes. tions and compute all the process streams in the sim-
• Optimization. plified process flowsheet, under the additional
• Simultaneous flowsheet and design specification loops. requirement that every pair of artificial streams be
• Dynamic simulations. identical (these artificial constraints ensure that all
the recycles are enforced). This is mathematically
In the EO approach the whole process is solved equivalent to solving Equation (6) (the recycle prob-
simultaneously as a system of nonlinear algebraic equa- lem), in which the vector x represents the process var-
tion. EO requires a more experienced simulator, is diffi- iables of the artificial streams that must be assigned,
cult to debug, and works well only when the initial and the function f(x) denotes the process variables of
value is close to the solution. Nevertheless, it works bet- the artificial streams that are simulation outputs. f(x)
ter where SM is weak. The simultaneous modular is an implicit function of x, as we can only evaluate it
approach is a combination of both SM and EO.[8] The by solving some/all of the unit operations in the sim-
flowsheet is solved in an EO fashion, while the singular plified process flowsheet.
unit are solved sequentially.[9] So far, this approach
finds application in academia but not in commercial
x − f ðx Þ = 0 ð6Þ
software.[10–12]
DE TOMMASO ET AL. 2305

FIGURE 3 Simple process flowsheet complemented with its graph representation (S-N: Nth process stream; C-N: Nth simple cycle or
th
N recycle)

(A) (B)

F I G U R E 4 Reduction of the flowsheet graph to a directed tree by elimination of all the simple cycles: A, optimal tear streams
according to Barkley and Motard[15] and Christensen and Rudd[16]; and B, optimal tear streams according to Upadhye and Grens[17] and
Motard and Westerberg[18]

Equation (6) contains only a small fraction of the influence the size and numerical properties of the recycle
nonlinear equations, which describe the original process problem (Equation (6)). Therefore, many different tearing
flowsheet of interest, provided that the latter contains a algorithms have been proposed over the years. However,
small number of recycles (the scale of Equation (6) is all of these algorithms can be grouped in three different
indeed proportional to the number of recycles). Therefore, families, which rely on as many different underlying
modular solution strategies are suitable for simulation of ideas:
process flowsheets of arbitrary scale, which contain few
recycles. On the other hand, they should never be used to • According to Barkley and Motard,[15] the tear streams
solve process flowsheets with many recycles, as the should be selected as the smallest set of process
numerical solution of Equation (6) becomes impractical. streams that break all the recycles in the process flo-
wsheet (Figure 4).
• According to Christensen and Rudd,[16] the tear
2.3.1 | Selection of the optimal tear streams should be selected as the set of process streams
streams that break all the recycles in the process flowsheet and
minimize the size of the recycle problem (Figure 4).
The identification of optimal tear streams is an element • According to Upadhye and Grens[17] and Motard and
of any modular solution strategy for process flowsheets Westerberg,[18] the tear streams should be selected as
because the number and features of the tear streams the set of process streams that break all the recycles in
2306 DE TOMMASO ET AL.

the process flowsheet, under the additional require- simultaneous alternatives for specific types of unit opera-
ment that the number of times every recycle is torn be tions, for example, distillation, absorption and stripping
minimum (Figure 4). towers (for completeness, the most common partitioning
solution scheme for separation towers is the inside-out
Although none of these three alternatives is fully sat- method by Boston and Sullivan Jr.,[19] Saeger and
isfactory, the best option depends on the numerical Bishnoi,[20] and Jelínek[21]). As a rule of thumb, we rec-
method, used to solve the recycle problem. More specifi- ommend partitioning methods for solving mildly non-
cally, if we solve the recycle problem with Newton-type ideal separations and simultaneous methods for solving
methods, then the second option represents the most highly non-ideal or reactive separations. In addition,
appropriate choice from a computational perspective. On simultaneous methods may be more appropriate for spe-
the other hand, if we solve the recycle problem with cific tasks, such as running case studies and/or per-
fixed-point algorithms, then the third option is the most forming flowsheet optimization, because of their higher
suitable from a computational point of view (this tearing computational efficiency. Additional recommendations
scheme maximizes the rate of convergence of fixed-point on which numerical solution scheme to choose for spe-
methods, when they are applied to Equation (6)). For cific simulation problems are usually available in the user
computational (and historical) reasons, all the main soft- guides of process simulation packages.
ware were originally designed as SM simulators. Hence,
the first developers focused on the implementation of a
single tearing algorithm that worked best in a SM envi- 2.3.3 | Solution of the recycle problem
ronment. Therefore, most process simulation packages
implement a single tearing algorithm, which often resem- The last key component of any modular solution strategy
bles the third, so the user’s choices are limited to what is for process flowsheets is the solution of the recycle prob-
available in existing process simulation software. Due to lem (Equation (6)), which apply either fixed-point or
the mathematical complexity of the tearing algorithm, we Newton-type methods. In the 1960s and 1970s, we
recommend a non-expert user to operate the default witnessed the birth of several new numerical algorithms,
method. specifically designed to solve this problem. Nowadays,
those that are still implemented in process simulation
packages are the direct substitution method, the
2.3.2 | Solution of the single unit Wegstein method,[22] the dominant eigenvalue
operation models method,[23] the Newton-Raphson method, and the
Broyden method.[24] The rationale of each are as follows:
A key step of all modular solution strategies for process
simulation is the solution of the recycle problem • The direct substitution method is the simplest and
(Equation (6)), which in turn requires repeated calcula- oldest fixed-point algorithm, whose recursion matches
tions of single unit operations. Since models of different Equation (5) with all the relaxation factors ai,n = 1.
types of unit operations exhibit unique features, dedi- • The Wegstein method is an accelerated version of
cated numerical methods to solve each and every one of direct substitution, in which the relaxation factors ai,n
them have been developed over the years. However, all are computed by Equation (7). We emphasize that
these numerical algorithms belong to either partitioning Equation (7) relies on the assumption that every vari-
methods or simultaneous methods. Partitioning methods able of the recycle problem is independent of the
decompose unit operation models into blocks of equa- others; therefore, this numerical method may not per-
tions, which are then solved iteratively, in a predefined form satisfactorily when the Jacobian matrix of f(x)
sequence, until convergence criteria are met or the itera- (Equation (6)) is not diagonally dominant.
tions diverge. Internally, they rely on fixed-point • The dominant eigenvalue method is an improved ver-
methods, and thus exhibit the same numerical properties sion of Wegstein, in which the relaxation factors ai,n
as this family of algorithms. Simultaneous methods do are computed by Equation (8). In this case,
not decompose unit operation models but rather solve Equation (8) relies on the assumption that the Jaco-
them as nonlinear systems. Since they apply Newton-type bian matrix of f(x) (Equation (6)) has a dominant
methods (possibly complemented with some parametric eigenvalue; therefore this numerical method may not
continuation approach), they exhibit similar numerical perform satisfactorily when several variables of the
properties as these types of algorithms. Most process sim- recycle problem are equally important.
ulation packages implement simultaneous methods for • The Newton-Raphson method is the simplest and
unit operations but offer both partitioning and oldest Newton-type scheme, whose recursion matches
DE TOMMASO ET AL. 2307

Equation (2). It relies on the actual Jacobian of Equa- attempt to simplify flowsheet calculations by exploiting
tion (6), so it may be more computationally efficient the sequential nature of chemical processes. They con-
than Broyden when the recycle problem is a small- struct the large-scale sparse system of nonlinear equa-
scale/medium-scale nonlinear system. tions, which describes the process flowsheet of interest
• The Broyden method is a Newton-type scheme, whose and includes unit operation models, process streams, and
recursion matches Equation (3). It relies on an esti- design specifications, and solve it using dedicated numer-
mate of the Jacobian of Equation (6), computed by ical methods. The numerical schemes, suitable for flo-
Equation (9); therefore, it may be more computation- wsheet simulation, are Newton-type methods, combined
ally efficient than Newton-Raphson when the recycle with some merit function minimization strategy and,
problem is a large-scale nonlinear system. often, some parametric continuation approach, which
helps the sequence of iterations converge whenever equa-
tions are highly nonlinear, their Jacobian matrix is locally
1 fi ðx n Þ −fi ðx n− 1 Þ rank deficient, and/or an adequate initial guess cannot
ai , n = ^si, n = 8i½1, N  ð7Þ
1 −si, n x i, n −x i, n− 1 be computed.[6]
Regardless of any specific numerical details,
1 k x n −x n − 1 k equation-oriented solution strategies solve any type of
ai,n = ^λMAX = 8i½1, N  ð8Þ
1 −λn
MAX n
k x n− 1 −x n− 2 k process flowsheets, independently of their scale and
topology (note that the performance of these solution
gðx n Þ −gðx n − 1 Þ −Bn − 1 dn− 1 T schemes is unaffected by the number of recycles in the
Bn = Bn − 1 + dn − 1 ^B0 = J ðx 0 Þ process flowsheet, provided that adequate computing
k dn− 1 k2
power is available). However, they should be primarily
ð9Þ
used to solve process flowsheets, which contain many
recycles, as modular solution approaches are more com-
All of these numerical approaches exhibit advantages putationally efficient for simulation of process flowsheets
and limitations, which make them more adequate for with few recycles.
specific types of process flowsheets. Table 1 collects a
general set of recommendations to help select a numeri-
cal alternative for most real-world simulation problems. 3 | A P P L I C A TI O N S

The program VoSViewer generated a bibliometric map of

2.4 | Equation-oriented (EO) solution keywords generated from articles indexed by Web of Sci-
strategies ence Core Collection (WoS) from 2017 to April
2020.[25,26] We queried the database four times with
As opposed to modular solution approaches, equation- ASPEN +, ASPEN Plus, HYSYS, and PROII as topics,
oriented solution strategies for process simulation do not while simultaneously excluding forest (to avoid work

TABLE 1 Numerical methods for flowsheet simulation

Numerical method Classification Robustness Efficiency Recommendations

Direct substitution Fixed-point Low Low Neither reliable nor efficient. Avoid its use.
Wegstein Fixed-point Medium Low Suitable for flowsheets that contain weakly interacting
recycles.a
Dominant eigenvalue Fixed-point Medium Low Applicable to process flowsheets with weakly and/or strongly
interacting recycles.a,b
Use when Wegstein fails.
Newton–Raphson Newton-like Low High Fast but sensitive to initial guesses.
Useful for flowsheet optimization and/or running case
studies.
Broyden Newton-like Low Medium Equivalent to Newton-Raphson; more efficient for large-scale
process flowsheets.
a
A set of recycles exhibits weak interactions when changes in the properties of any recycle stream induce small variations in the properties of
all the other recycle streams.
b
A set of recycles exhibits strong interactions when changes in the properties of any recycle stream cause consistent variations in the proper-
ties of all the other recycle streams.
2308 DE TOMMASO ET AL.

TABLE 2 Most employed thermodynamic models to calculate components’ fugacities or activity coefficients

Model Advantages Limitations

Equation of State (EoS)
IDEAL Good for mixtures of similar fluids at P < 0.3MPa
low P and T
Chao-Seader[107] Generalized correlation for P < 14MPa, T < 260∘C
hydrocarbon mixture Dissolved gas <20% by mole
Lee-Kesler[108]  S
Correlation for H,  and ρ Not good for highly polar mixtures
Redlich-Kwong[109] Calculates hydrocarbons VLE Not good for polar compounds
accurately
Peng-Robinson[110] Better than SRK with polar components Worse for hydrocarbons mixtures
[111]
UNIWAALS Predictive method Components’ liquid volumes required
Activity coefficients
Wilson[112] Polar or associating components Does not calculate LLE
[113]
NRTL Strongly nonideal mixtures Very different molecules dimensions
[114]
UNIQUAC Very good with most of systems No electrolytes
[115]
UNIFAC Predictive method T < < TC, P < 8.5MPa

related to Aspen trees and forests rather than process recovery. Xiao et al[29] dimensioned a three-column pro-
simulation). We combined the four databases into one, cess to transesterify ethylene glycol with methyl acetate,
eliminated all duplicates, and were left with 2582 articles. regressing reaction kinetics parameters, and optimized
The WoS category energy & fuels had 930 articles, while the design minimizing the TAC. Wu et al[30] designed a
the chemical engineering category was a close second dividing-wall column, together with its control strategy,
with 824 articles, followed by thermodynamics (314 arti- for reformed gasoline. Biodiesel and esterification are
cles), environmental engineering (224), and green & sus- concepts less related to process simulation as they are far
tainable science and technology (215). from the map centre. This is because first- and
We retained about 90 of the keywords that were cited second-generation biofuels interest has been dropping in
most frequently in these articles but excluded simulation, recent years in favour of other green technologies for the
performance, column, behaviour, and growth. The conversion of green house gases, such as CO2 into
VOSViewer grouped the keywords into five clusters with fuels.[31–33] In fact, the magenta cluster centered around
Aspen+ at the centre of the red cluster with the most H2 (227 articles) concentrates on gasification (biomass
articles (433) and the cluster with the most keywords at +steam with 239 articles), syngas (141), and reactors and
29 (Figure 5). This cluster includes topics related to bio- fluidized beds (84).
mass (193 articles): ethanol and bioethanol (132), ligno- CO2 capture (193 articles), with CO2 (162), and CH4
cellulosics, and cellulose (64). The two other major topics dominate the yellow cluster that includes combustion
relate to catalysis/kinetics (162) and techno-economic (87), absorption (82), and HYSYS (66). This cluster spans
analysis/economics (227), in addition to other subjects a large area that covers technology related to carbon cap-
related to biomass. ture and sequestration (CCS). Duhoux et al,[34] for exam-
Process simulation discriminates between process ple, calculated the optimal pressure and flowrate of a
alternatives (design: 246 articles), which comprises pressurized fluidized bed combustion that sequestrates
modelling (187 articles) and determines the dimensions carbon dioxide via Calcium looping. Joule has published
of unit operations like separation (126), including distilla- the most cited article since 2019 (101 citations) that
tion (137 articles) (reactive and extractive), heat integra- describes a plant to capture 1 Mt CO2 per year with aque-
tion (33), and dehydration (26). Hung et al,[27] for ous KOH sorbent and a calcium caustic recovery loop.[35]
example, demonstrated that recovering dilute acetic acid In fact, the article spans many of the keywords in the
water solutions (between 30-70% by mass of acid) with map: unit operations, capital cost estimates, and heat and
amyl alcohols among C1-C5 minimizes the total annual mass balances with ASPEN+.
cost (TAC), without any stream pretreatment. Pirola The optimization (266 articles) cluster is connected to
et al[28] simulated an extractive distillation column both design and CO2 capture (Figure 5). Many of the sub-
employing p-xylene as entrainer to maximize acetic acid jects relate to process intensification, which includes
DE TOMMASO ET AL. 2309

F I G U R E 5 VOSViewer keyword bibliometric map based on keywords (Aspen+, Aspen plus, HYSYS, and Pro/II) derived from articles
that WoSindexed from 2017 to April 2020.[25,26] The text font size and diameter of the circles are proportional to the number of occurrences
in the 2583 articles indexed during this time. Aspen+ (largest circle) appeared in 433 articles. The smallest circles appear in 26 articles

energy (228 articles), systems (215), energy/efficiency creative environment that promotes team work is
(127), and exergy (163).[36,37] Dynamic simulation ana- generated.[46]
lyzes process start-up and shut downs and closely con- Process simulation requires a solid theoretical back-
nects to process design.[38] ground as well as a basic knowledge of the computational
Process simulation is an educational tool to demon- methods the software are based on. For this reason, lec-
strate many facets of design.[39,40] Steady-state simulation tures on the simulation suite and capabilities are propae-
represents the plant at stationary conditions, while deutic to understand and operate correctly.
dynamic simulation is required for the plant control strat- Dynamic simulation is required for control strategies,
egies, for start up and showdown operations, and for for start up and showdown operations, and for operator
operator training by plant virtualization.[41] The training by plant virtualization.[41] The detailed represen-
suggested approach for undergraduate students is steady tation of commercial control algorithms, the high-
state simulation,[42] while dynamic simulation is appro- realistic models of several unit operations and the repro-
priate for graduate-level courses.[43] For instance, the duction of the real behaviour of pipes, vessels, and valves
Universidad Complutense de Madrid, introduced steady are the basis of the most advanced dynamic simulation
state simulation software in a fourth-year course of the tools.[43] These options allow students to interact with a
chemical engineering degree,[44] proposing a case study virtual plant and to operate instruments. There, the stu-
on the reactive distillation for the production of methyl- dents engage at a practical level (starting up of a pump,
tertbutyl ether (MTBE). Simulation in chemical plant sampling, reading gauges, etc.) and in plant management
courses guides students through the principles of unit and control. This immersive simulation software com-
operations and is a support in the development of a bines augmented reality, computational dynamics, and
whole process.[45] For example, simulations reveal theo- computer graphics.[47]
retical trends of processes, verify the kinetic and thermo- Augmented reality simulation examples concerning
dynamics of reactions, and calculate fluid phase mixing simulated chemical plants are rare, while robotics and
and equilibria.[42] When proposing experimental activi- physical sciences already developed commercial applica-
ties in parallel with theoretical lectures and simulation, a tions and university courses.[48] Nevertheless, new
2310 DE TOMMASO ET AL.

3.1 | Cost estimation

The log of the investment (logI) of a process unit and

even an entire plant is linearly correlated with the log of
its capacity (logQ). The slope of the curve depends on the
equipment type, specifically on its characteristic
dimension,[63] and it varies from 0.33 to 1.0 for equip-
ment and from 0.38 to 0.90 for plants, averaging to a
value of 0.6[64] (Equation (10)):

 0:6
FIGURE 6 Operator training simulators (OTS) are an Q2
I2 = I1 ð10Þ
educative tool for master’s and undergraduate students. A dynamic Q1
simulation of the plant calculates the consequences of each
operator’s actions and the student learns through experience
For modular units, the learning elasticity is proportional
to this exponent and the cost of these units decreases with
educational proposals are on-going, such as for example multiple units (learning/experience).[65] Lang[66–68] elabo-
the Eye4Edu project at the University of Milan (Italy). rated a simplified method to calculate the total installed cost
Eye4Edu applies EYESIM software, from AVEVA, that of equipment based on factors that vary according to the
proposes the immersive virtualization of a crude distilla- type of plant (solid handling, solid–fluid, or fluids).
tion unit (CDU) plant with the representation of the A more detailed approach, proposed by Guthrie,[69,70]
detailed realistic 3D graphics of the whole plant estimates the bare module cost of a unit, CBM, based on a
(Figure 6).[49] DYNSIM software simulates the plant and characteristic variable, A (volume for a reactor or a col-
supplies the physico-chemical properties. The main prob- umn, surface for an heat exchanger, etc.), corrected with
lems and limitations of the software in education are: pressure coefficients, Ci, and material coefficients, Bi
(a) the imperfect reproduction of the reality, that is, there (Equation (11)):
is always some discrepancy compared to the real world. 8
It is therefore impossible to recreate all the layouts of a >
> C BM = C 0p ðB1 + B2 F M F P Þ
<
real plant; (b) simulators are expensive and require main-
log10 C0p = K 1 + K 2 log10 A + K 3 ðlog10 AÞ2 ð11Þ
tenance and constant updates; and (c) operators and >
>
: 2
teachers need training. log10 F P = C1 + C 2 log10 P + C 3 ðlog10 PÞ
Process simulators also train senior engineers and
field operators.[50–52] Commercial software simulate FM is the material factor, that is unity for carbon steel
mostly reactions, separations, and design heat integra- and >1 for other materials, and FP is a pressure factor. In
tion, and offer a variety of ideal reactors.[53,54] Depending the case of vessels, a different equation calculates FP,vessel
on the field of application, the software is integrated in (Equation (12)):
the design of a new process, in the simulation of an exis-
8 1
ting plant, or both.[55] For instance, if we simulate and > ðP + 1ÞD
>
> + 0:00315 C
benchmark the energy consumption of an existing plant, >
<F 2850 −0:6P + 1
P,vessel = C
we increase its overall efficiency and save money. In 2017 0:0063 C ð12Þ
> C
>
> F = 1 if t < 6:3mm A
and 2018, researchers demonstrated the efficacy of the >
:
P,vessel vassel

process integration in the cement[56] and the petrochemi- F P,vessel = 1:25 if P < −0:05MPag
cal industry,[57] with a total duty saving of 24% and 14%,
respectively. Cold and hot utilities consumption, energy which is valid for a corrosion allowance of 3.15 mm.
requirements and emissions of existing oil and gas,[58,59] Ulrich and Vasudevan[71,72] improved the method of
electric power,[60] biofuels,[61] chemical, pharmaceutical, Guthrie, extended his correlation to 58 unit operations of
and urban systems[62] plants can all be reduced with a the chemical industry, and updated the correlations with a
retrofit based on a simulated process. Optimizing a pro- larger database. All the estimations are then actualized by
cess, instead, requires setting economic, process, or envi- an inflation index like the Marshall and Swift equipment
ronmental objective functions to satisfy. Commercial cost index, the chemical engineering plant cost index
software, such as AspenPlus, HYSYS, PRO II, or Super- (CEPCI), the engineering news-record (ENR) construction
Pro Designer, are designed to work with a fully defined index, or the Nelson-Farrar refinery construction cost
problem (sum of degrees of freedom equal to zero). index (NFRCI), where CEPCI and NFRCI refer to the
DE TOMMASO ET AL. 2311

entire plant and the others to either plant or equip- environmental and/or social impact, and the error over
ment.[73] Simulation packages have tools to estimate a time of a control model. When two or more conflicting
units’ costs, which compute a combination of the methods objective functions must be maximized/minimized simul-
reported with vendors’ quotes/estimates.[74] They prefer to taneously, a dedicated optimization method is required to
implement algorithm-based methods (Guthrie, Ulrich, solve this multi-objective optimization problem. For
Seider[75]) rather then cost-to-capacity or the Lang-based example, Patel and Padhiyar[77] solved a bio-reactor
methods. For instance, the Aspen Capital Cost Estimator design problem in which they simultaneously minimized
(ACCE) (formerly APEA), the in-house Aspen estimation the batch time and maximized the process yield.
software, couples design and equipment cost modules with Additionally, Equation (13) may have multiple solu-
proprietary time-proven industry-based direct and indirect tions, called local optima, and finding the best of these,
factors based on company project history and vendor that is, the global optimum, is essential. In this case, the
quotes. These factors are updated regularly.[75] Industry user selects special types of optimization algorithms that
and academia developed dozens of integrated and stand- offer global optimality guarantees:
alone software to estimate cost. Cleopatra and ACCE are
the built-in software of PRO/II and AspenOne, respec- • Linear problems (LP): linear objective function and
tively. Another example is the built-in Capital Cost Esti- constraints, and continuous decision variables.
mator of the CAPE-OPEN project DWSIM. EstPro is a • Mixed-integer linear problem (MILP): linear objective
stand-alone process plant cost estimation package from function and constraints, and both discrete and contin-
Gulf Production while CapCost, CCEP, and DFP are uous decision variables.
included with books.[73,75,76] EconExpert is web-based soft- • Non-linear problems (NLP): non-linear objective func-
ware. Even software heavily based on updated vendor tion and/or constraints, and continuous decision
quotes have errors of at least 25%.[74] We recommend variables.
operating with one model (or software) to evaluate differ- • Mixed-integer non-linear problems (MINLP): non-
ent plant designs to generate comparable results. linear objective function and/or constraints, and both
discrete and continuous decision variables.

3.2 | Optimization Scheduling, supply-chain, and planning problems are

usually LPs or MILPs. NLPs and MINLPs characterize
Most process synthesis problems have multiple feasible dynamic optimization, design, and process synthesis prob-
solutions. Quantitative optimization techniques select the lems. Gradient-based and derivative-free algorithms solve
best of these solutions, which corresponds to a process con- NLPs and MINLPs. Both the algorithms are iterative; they
figuration that maximizes a problem-specific measure of start from an user-supplied initial guess and progressively
process performance. Any process synthesis problem can be move towards a solution of the optimization problem
reformulated as an equivalent optimization problem (a set of optimal values of the decision variables).
(Equation (13)) composed of three principal components: Gradient-based methods rely on the values and the deriva-
an objective function, f(x, y) (a problem-specific measure of tives of constraints and objective functions, while
process performance), some decision variables, continuous derivative-free methods rely on the value of the functions.
X (process variables), or binary y (they specify the presence Gradient-based algorithms are suitable for mildly non-
or absence of certain process units), and a set of constraints, linear optimization problems with smooth objective func-
inequality g(x, y) (eg, minimum product purity, maximum tions and constraints, given that the initial estimates are
allowable equipment cost, environmental emission limits) coherent. Derivative-free methods are appropriate for
and equality h(x, y) (eg, mass, momentum, and energy bal- highly non-linear small-scale optimization problems with
ances, mass and heat transfer correlations, and phase equi- both smooth and non-smooth objective function and con-
libria). The numerical solution of this optimization problem straints, even if an adequate initial guess is unknown.
requires discrete or continuous solvers: These optimization algorithms are partially complemen-
8 tary, which is the reason why process simulation packages
< Minimize f ðx,yÞ
> implement a few gradient-based and derivative-free
Subject to gðx, yÞ ≤ 0, hðx, yÞ = 0 ð13Þ methods, and let the user excogitate on the one to select.
>
: The LP simplex method for LPs,[78] sequential qua-
xX, yf0, 1g
dratic programming (SQP) for NLPs,[79] and the branch
and bound (BB) method for MINLPs are examples of
Examples of common objective functions include gradient-based methods.[80] SQP is a quasi-Newton algo-
total cost, profit, energy expenditure, exergy loss, rithm that is computationally efficient but it relies on first
2312 DE TOMMASO ET AL.

and second order derivatives of the objective function

and constraints of the optimization problem. For this rea-
son, SQP solves smooth problems best such as LPs and
NLPs (Aspen Plus and PRO/II feature SQP as default
optimization method). Additionally, SQP does not guar-
antee convergence to the global optimum.
Derivative-free optimization strategies include BOX,
which is a basic implementation of the Nelder-mead sim-
plex algorithm,[81] and COMPLEX, whose implementa-
tion details are property of Aspen HYSYS and Aspen
Plus. These built-in algorithms handle simple optimiza-
tion problems with a single objective function, provided
that a decent initial guess is available.[82,83] However, F I G U R E 7 Example of hot (red) and cold (blue) composite
complex or multi-objective optimization problems may curves. The blue area is the minimum cooling duty required Qc and
require external optimization packages (eg, CPLEX and the red area is the minimum heating duty required Qh
XPRESS for solving LPs and MILPs, CONOPT, IPOPT
and KNITRO for solving NLPs, and DICOPT, BARON
and ANTIGONE for solving MINLPs), interfaced with in terms of the heat content, we have, for each stream, a
the simulation flowsheet. For instance, Galli et al[84] curve whose slope is 1/CP. For a single couple of hot-cold
applied Monte Carlo optimization of operating costs and streams this approach is of little use because to a k hot
environmental impact of a plant for producing oxygen- load increase, corresponds a k cold required duty
enriched air, simulated with PRO/II. Quiroz-Ramírez increase.[90] However, the graphical representation comes
et al[85,86] optimized a bio-butanol plant with AspenPlus, in handy when there are several streams involved. In
MATLAB, and a home-made VBA script. Eslick and this case:
Miller[87] optimized an amine absorption process, with
Excel, Aspen HYSYS, and the solver NSGAII. Ponce- • First, plot all the streams separately on a T/H graph.
Ortega and Hernandez-Perez[83] developed a procedure Since we are interested in the H, we only have to place
to tune any commercial simulator to any multi- the stream on the right interval of the y axes.
optimization algorithm with MS Excel. In fact, they • Then, add up all the cold and heat duties separately.
implemented a hybrid stochastic algorithm called I- Since the slope of each stream represents its CP,
MODE to maximize the gross annual profit and to mini- whenever two or more streams share a temperature
mize the CO2 emissions of a power and a biodiesel interval T1 − T2, the heat available in this interval
plant.[88] Interfacing commercial simulation packages will be (CP,i + CP,ii + …CP,n)(T1 − T2). The outcome of
with external optimization tools is a considerable task this procedure are the hot and cold composite curves.
with respect to computational power, and time; therefore,
we recommend this technique only as a last resort. When we plot these together, we identify (Figure 7):

• The minimum cooling duty required QC.

3.3 | Energy integration intensification • The minimum heating duty required QH.
• The recovered duty where the curves overlap QREC.
Since 1970,[89] pinch analysis has guided the placement of • The point of closest approach, the pinch point.[89]
heat exchanges and other unit operations (distillation, for
example) to minimize energy consumption.[9] We visualize This means that as long as the heat exchange across
each non-reactive stream by a heat content H (kW), tem- the pinch operates at the corresponding minimum ΔT, it
perature couple. When a differential heat flow dQ is added is possible to design an heat exchange network (HEN)
to a stream, it increases its enthalpy H by CPdT: who will recover QREC. The only external heat for the
process will be QH and QC. Lindhoff and Flower[91]
ðT2
expanded the graphical concept of the composite curve,
Q= C P dT = CP ðT 2 −T 1 Þ = H ð14Þ
T1 to the problem table method. This latter method, algebra-
ically splits the process in a cascade of temperature inter-
vals. To do so, we adjust the hot or cold stream
where T2 and T1 are the target and the supply tempera- temperatures of a negative or positive fraction of ΔTpinch.
ture, respectively. In this way, plotting the temperature For instance, we either shift the hot streams by -0.5Tpinch
DE TOMMASO ET AL. 2313

and the cold streams by +0.5Tpinch, or the hot streams by variations, for example).[96–98] For this reason, we recom-
−ΔTpinch and leave the cold streams untouched. Setting mend understanding the dynamics of all the operations
the temperatures this way guarantees that each interval in the simulation before integrating energy.
will either have a positive or a negative duty balance.
Foreach interval i we will, therefore, have a balance:
X 
3.4 | Process simulation as research tool
ΔH i = ðT i − T i + 1 Þ CP,h −CP,c ð15Þ
i Process simulators started as academic tools but soon
thereafter industry adopted this software for design. Now
When we sum up the ΔH for all the intervals, we this tool is returning to academia and research. Saidi and
obtain a positive and a negative heat residual. The only Kadkhodayan[99] integrated an experimental Taguchi opti-
way to close the energy balance is to introduce a cold and mization method to a process simulation to bridge indus-
a hot utility. We have to supply a QC and QH at a thermal trial operation with laboratory scale experiments. Jafari
level lower and higher of the lowest and highest interval et al[100] simulated a circulating fluidized bed with a
of the cascade. That is, if the highest interval is at sequential modular approach and showed that it could be
T = 160∘C, we include a hot utility above 160∘C. added to a commercial simulation. Commercial software
Brown[92] demonstrated how pinch technology achieved developers are expanding their offering with non-
an average 25% cut on the site energy bills, equivalent to conventional components, such as electrolytes, solids, or
30 million pounds in 1989, for all the main sectors in plastic.[55] Process simulators support the experiments and
chemical engineering. Although other methods like the approximate thermodynamic parameters with equations
composite curve method are equivalent,[75,93] the inher- of state for those parameters difficult to measure. Shen
ent modularity of the problem table method makes it per- et al[101] studied a primary cooler for a coke oven gas LiBr
fect to implement it in process simulators. heat pump. They relied on simulation results to estimate
Process simulators perform pinch point analysis operating parameters they were unable to measure experi-
(PPA) to design an optimal heat exchange network mentally to build a refrigerator pilot plant. They retroac-
(HEN) and minimize the energy requirement of the pro- tively validated the simulation with pilot plant data and
cess and the utilities consumption.[53] Several software reported an error of 0% to 10% depending on the parame-
packages on the open market incorporate a PPA ter. Process simulators are also tools to calculate thermo-
section such as HYSIM, HEXTRAN, ADVENT by Aspen dynamic properties (enthalpy, entropy, Gibb’s free
Tech, TARGET, or, more recently, FI2EPI.[94] Simulators energies, densities, viscosity) as input to estimate reaction
couple PPA analysis with a UA heat exchange model, kinetics, heat duties, and phase changes. For example, to
and in this way they work on the design of the network define an experimental plan to study methane partial oxi-
outside the simulation environment. They share the pos- dation, we estimate the theoretical bound on conversion
sibility of importing the H&M balance results from out- and selectivity as a function of temperature and pressure,
side, and the presence of at least a manual and an and then determine the relationship between coking and
automatic design. However, we recommend completing oxygen partial pressures. These data define the boundaries
individual PPA before integrating them all into the entire of the plan. PRO/II simulated a Gibbs reactor (RGibbs)
network. In fact, the automatic mode will achieve the tar- with the Peng-Robinson equation of state. A case study
get (minimum number of heat exchangers or maximum varied the O2/CH4 ratio in the range 0.2 to 1.0 and the
energy saving), regardless the process engineering of the temperature of the reactor from 1000 to 1200 K (Figure 8).
system. For instance, it might suggest you to exchange
between stream 1 and stream 1000 of the plant, not con-
sidering the physical distance. Alternatively, it might sug- 4 | UNCER TAINTY
gest to split one stream into five sub-streams, to
maximize its heat exchange. Based on what you are 4.1 | Limitations
looking for, it might also propose uneconomical heat
exchangers. On the other hand, this software will not Commercial software licenses are expensive and to simu-
suggest you where to locate an equipment with respect to late non-standard operation often additional code is
the pinch.[9,75] Furthermore, energy intensification works required. The commercial software models and algo-
as long as the plant operates at steady-state. However, rithms are robust, and most of the errors come from user-
while the steady-state defines the baseline operation,[95] added inputs. Integrating custom code into process simu-
in real life a plant alternates between steady-state and lators is cumbersome and takes longer than writing a
unsteady state (start-up, shut-down, feed composition program from scratch.
2314 DE TOMMASO ET AL.

inputs the main kinetic parameters (Arrhenius constants,

Apre, activation energies, Ea, and reaction orders, ε) for a
kinetic model or by defining a pseudo-homogeneous
kinetic scheme (Equation (16)):
   Y
Ea 1 1
r = Apre exp − − T n  aεi ð16Þ
R T To

where a the activity of reagents and products (calculated

as the product concentration, molar fraction, or according
to a thermodynamic model).
PRO/II permits the user to supply the values of these
F I G U R E 8 Thermodynamic (maximum) methane conversion parameters and the volume and the pressure UOMs, as
(dotted lines) and obtainable H2/CO ratio (full lines) for the the r is always defined as (mol  volume−1  time−1). Also,
catalytic partial oxidation of methane at 1000 K (black lines) and the value of the gas constant R is 8.314 J  mol−1  K−1.
1200 K (red lines). The simulator also identifies at which ratio coke
Therefore, the UOMs of pressure and volume in the reac-
starts to form (points). We calculated the H2/CO ratio at which
tion window must be chosen to give J when multiplied,
coke forms every 20 K. The simulations took less than 8 minutes to
that is, or kPa and L, or Pa and m3, respectively.
complete
On the other hand, when defining a kinetic proce-
dure, the user writes a FORTRAN or pseudo-FORTRAN
To set up a simulation file, the following procedure code and calculates the system of differential equation
applies: that the software will integrate using a numerical method
(typically Runge-Kutta). In this case, the user should pro-
• Define the units of measure, according to a predefined vide the proper rate expression UOM to the solver, in
set (metric, English, or SI), or to a customized set. volume−1  time−1.
• Select the components (from a databank).
• Choose a thermodynamic model.
• Check the simulation tolerances. 4.2.2 | Components selection
• Build the flowsheet. In case of recycles, see section 2.
The selection of components means, for the software, to
Errors compound as the simulations proceeds from import all the pure component properties, the fixed
one step of the process to the next. Simulators have no (molecular weight, acentric factor, normal boiling point,
error assessment capability. critical point, Van der Waals area and volume, etc.) and
the temperature dependent properties (density, vapour
pressure, viscosity, heat capacity, etc.).
4.2 | Sources of error We always recommend employing components pre-
sent in the software database, as these are usually
4.2.1 | Units of measure (UOM) updated with one or more methods to estimate the com-
ponent’s properties. Besides, a good practice is to operate
One of the most frequent error source is incorrect data with less than 40 pure compounds per simulation.
entry, that is, the human error. Even though most of the In the case of new components (those not yet in any
commercial simulation suites have a window dedicated database), all the physical properties should be provided
to the selection of UOMs, it is possible to change these as well; we recommend to retrieve the physical and ther-
inside all unit operations. For example, if the ENGLISH modynamic properties on databases such as NIST chem-
system is set as default, but for a specific distillation istry WebBook[102] or Knovel.[103]
application the user knows that the distillate flowrate is There are dozens of equation to estimate thermody-
100 kmol  h−1, inserting a value of 100 without cor- namic properties.[104] This reference covers also electro-
recting the UOM will result in setting a distillate flow of lytes systems and solids solubilities. We recommend
100 klbmol  h−1, resulting in a error of 1/2.2. The most paying particular attention while operating with new
likely error is the simulation will not converge to a materials and processes such as:
solution.
Another common mistake is when users define units • Green and biological processes.
of measure in the reaction kinetics window. The user • Molten organic and inorganic salts.
DE TOMMASO ET AL. 2315

FIGURE 9 Simplified decisional flow diagram for the selection of a thermodynamic method (Table 2).[116,117] Reprinted from A. C.
Dimian, C. S. Bildea, A. A. Kiss, in Integrated Design and Simulation of Chemical Processes, 2nd ed. (Eds: A. C. Dimian, C. S. Bildea, A. A.
Kiss), Elsevier, Amsterdam, 2014, p. 248 (Copyright 2020) with the permission of Elsevier

• Nonconventional solids (eg, solid with unknown for the gas phase association of acetic acid (Hayden-
chemical formula). O’Connell equation[118]) or regressing with a robust method
• Electrolytes. experimental data led to a sound simulation.[28] Therefore,
• Novel materials (nanomaterials, composite materials, we always recommend a literature analysis of the phase equi-
copolymers, etc.). libria of the systems involved. Gmehling et al[119] published a
book regarding chemical thermodynamics and process simu-
Commercial software as ASPEN, PRO/II, and ProSim lation. Gani and O’Connell[116] created a decision matrix to
estimate new component’s properties from UNIFAC select proper physical methods and thermodynamic models,
group contribution methods.[105] In this case, the flaws of depending on the unit operation simulated and the operative
the simulators rely in the estimation of the physical prop- parameters.
erties, which are either imprecise,[106] or incomplete.

4.2.4 | Tolerances
4.2.3 | Thermodynamics
Similar to any other iterative calculation algorithm, the
One of the most critical points is the selection of a proper solver of a process simulator needs specific tolerances to
thermodynamic method to calculate mixture properties terminate the calculations. Generally the default parame-
(excess Gibbs free energy or fugacity). Analyzing the system ters (1 × 10−3 for temperature and pressure (relative),
(Figure 9) leads to the selection of a γ − ϕ or a ϕ − ϕ 1 × 10−6 for compositions (absolute), and 1 × 10−4 for
approach. However, this is a simplified scheme as each com- duties (relative)) are sufficient to guarantee the robust-
ponent mixture may show different non-idealities. For exam- ness of the termination criteria.
ple, the system water and acetic acid (fully miscible,
condensed phase and no supercritical gases dissolved)
should be treated with an activity coefficient model. How- 4.2.5 | Flowsheet
ever, with UNIQUAC equations and the binary parameters
available in PRO/II database only, the software calculates an Any unit operation requires specifications to close the
non-existent azeotrope (Figure 10). Only adding a correction mass and energy balances. Generally these specifications
2316 DE TOMMASO ET AL.

vary from one unit to another. For an exhaustive descrip-

tion of each method consult the software manual.
For complex simulations, involving several units and
recycling streams, we recommend choosing initial condi-
tions judiciously, or estimate them based on shortcut cal-
culations or literature data, to facilitate the algorithm
convergence to the solution. Poor initial guesses will
increase the number of iterations unnecessarily (if the
simulation converges at all).
To initiate the initial guess for tear streams, an ade-
quate strategy is to simulate the flowsheet with no recycle,
and then connect the recycles, so that the software stores
F I G U R E 1 0 The correct selection of a thermodynamic after the first simulation the initial values of flow rate,
method (blue line) avoids the calculation of non-existent azeotropes temperature, pressure, and compositions of these streams.
(red line) for the binary water-acetic acid system, P = 101.325kPa.
UNIQUAC parameters: aij = − 118.309, bij = − 0.124345K−1,
aji = 402.126, and bji = − 0.355725K−1 4.3 | Sources of error in the EO
Approach
are temperatures, pressures, or material flows (like the
flow of distillate and residue for a distillation column) While thermodynamic, component selection, and
but any unit may have peculiar specifications (reactor flowsheets are common for every computational architec-
dimensions for a plug flow or time for a batch reactor, ture, the EO is intrinsically free of any sources of error
reflux ratio for a distillation column, etc.). Converging to related to iteration. However, since the EO solves all the
a viable solution depends on setting these specifications units at the same time, the initial values are essential. As
correctly. Most of the error in this case are related to the a rule of thumb, we recommend first running the simula-
unit of measures of the specifications or to gross errors tion in a SM environment, and then implement the SM
such as setting a too high reflux ratio or a distillate flow outcome as initial value for the EO. This does not mean
higher than the feed flow rate of the column. that the SM needs to converge, but at least each block
Some specifications make the solver’s life more diffi- has to be solved once. The closer the SM solution is to
cult. If the simulation never converges to a solution, the the real solution, the easier it will be for the EO to con-
user should replace some of them with equivalent specifi- verge without error. While this initialization technique is
cations (eg, instead of decreasing the molar fraction of the one currently present in commercial simulators, new
impurities in the distillate, increase the reflux ratio until approaches such as the graphic approach,[120] or the
the composition reaches the desired value). Another com- pseudo-transient,[121] have been recently emerging.
mon mistake is to set the composition of distillate or resi-
due to a value impossible to obtain with the number of
theoretical trays specified (such as imposing a distillate 5 | CONCLUSIONS
or residue purity higher than the minima/maxima
azeotropic compositions) for distillations, or imposing by Engineers design, control, optimize, retrofit, and estimate
mistake to concentrate the heavy component in the distil- the cost of a plant and equipment with process simula-
late rather than in the residue stream. tors, which calculate thermodynamic properties, model
If the user imposes extreme constraints (ie, a separa- vapour- and/or liquid-liquid equilibria of mixtures and
tion in which the molar fraction of impurities in the prod- close mass and energy balances. Process simulation will
uct is less than 1 ppm), then they should start with a more become ever more present in a chemical engineer’s arse-
relaxed constraint (±5%) and tighten it up gradually. nal of tools to address society’s desire for sustainable
products that minimize waste and maximize recyclability
as the mathematical methods and computational power
4.2.6 | Initialization improve such that any PC can run these multi-objective
optimization problems. The feedstocks of the future will
All numerical methods for the resolution of unit opera- become more heterogeneous (bio-based, waste-based, for
tions require initial values to start iterating. However, all example) and because of this complexity, new models,
commercial packages possess several initialization unit operations, and thermodynamic correlations will be
methods for each unit operation, and these algorithms needed. Quantification of the uncertainty on model
DE TOMMASO ET AL. 2317

predictions is another key aspect that should be further ACCE Aspen capital cost estimator
improved in future simulators as well as adding stochas- Apre Arrhenius pre-exponential factor
tic optimization capabilities to these tools. ai activity of component i
Prospectives for future innovations include: B() Jacobian matrix estimation
Bi material coefficients
• generating kinetics databases (including packages like CBM cost of bare module
Chemkin;[122]) CCS carbon capture and sequestration
• handling solids and other unconventional materials CDU crude distillation unit
(plastic mixtures, urban waste, sludge, and poorly CEPCI chemical engineering plant cost index
characterized complex mixtures); Ci pressure coefficients
• better dynamic simulators that consider the time scale CP specific heat
of reaction vs hydrodynamics contribution; D diameter of a vessel
• new operating functions for process intensification of dn search direction vector
current processes;[123,124] Ea reaction activation energy
• more reliable safety assessments that automatically ENR engineering news-record
identify hazardous operating conditions (HAZOP auto- EO equation oriented approach
matic calculation of lower and upper explosion EoS equation of state
limits,[125] identification of risk zones, improving the FM material factor
hazard identification when streams mix); FP pressure factor
• integrating computational fluid dynamics; i
H enthalpy of component i in mixture
• dynamically query other databases on the inter- HOC Hayden-O’Connell
net (NIST); I investment
• integrating artificial intelligence and artificial neural J() Jacobian matrix
network algorithms for control and optimization; Ki cost coefficient of an equipment
• improving the thermodynamic packages for new mole- LLE liquid-liquid equilibrium
cules (pharma and green chemistry); LP linear problem
• improving the predictive control of unit operations; MILP mixed integer liner problem
• expanding virtual reality that reproduces realistic plant MINLP mixed integer non-liner problem
dimensions; and MTBE methyl-tertbutil ether
• improving the integration between home-made code NFRCI Nelson-Farrar refinery construction cost index
and process simultators. NLP non-linear problem
OTS operator training simulator
Also, the integration of multiobject optimization on P pressure
performance key indicators not only based on economic PPA pinch point analysis
parameters, such as the CO2 equivalent or the water foot- PR Peng-Robinson
print of a process (environmental impact) or the hazard Q plant capacity
and operability analysis, which at the moment is feasible Qc minimum cooling duty
via open platform communication.[126] Qh minimum heating duty
As stated by Coon et al[127] in 1998, “it is the responsi- Qrec recovered duty
bility of the engineer to validate the quality of process R gas constant
design, troubleshoot, and optimization results from simu- r reaction rate
lation tools” and that wrong answers are “the results of i
S entropy of component i in mixture
getting the correct answer to a poorly chosen question.” SM sequential modular approach
SQP successive quadratic programming
ACK NO WLE DGE MEN TS SRK Soave-Redlich-Kwong
The authors acknowledge professor Rahmat Sotudeh T temperature
Gharebagh for his valuable comments and suggestions, TC critical temperature
and the fruitful discussions. tvessel vessel wall thickness
TAC total annual costs
UOMs unit of measures
NO MEN CLATU RE VLE vapour-liquid equilibrium
A characteristic variable of an equipment xn solution vector at the iteration n
a optimal step length xs solution vector
2318 DE TOMMASO ET AL.

G R EE K L E TT ER S [27] W.-J. Hung, I.-K. Lai, Y.-W. Chen, S.-B. Hung, H.-P. Huang,
γ activity coefficient M.-J. Lee, C.-C. Yu, Ind. Eng. Chem. Res. 2006, 45, 1722.
ε reaction order [28] C. Pirola, F. Galli, F. Manenti, M. Corbetta, C. L. Bianchi,
Ind. Eng. Chem. Res. 2014, 53, 18063.
ϕ fugacity
[29] Y. Xiao, W. Cai, H. Sun, F. Shi, G. Li, Can. J. Chem. Eng.

ρi density of component i in mixture 2018, 96, 722.
[30] S. Wu, K. Guo, C. Liu, W. Qi, T. Zhang, H. Liu, Can. J. Chem.
ORCID Eng. 2018, 96, 2431.
Carlo Pirola https://orcid.org/0000-0002-7013-5424 [31] C.-K. Kuei, M.-D. Lee, Can. J. Chem. Eng. 1991, 69, 347.
Gregory S. Patience https://orcid.org/0000-0001-6593- [32] N. A. Bahari, W. N. R. Wan Isahak, M. S. Masdar, Z. Yaakob,
7986 Int. J. Energy Res. 2019, 43, 5128.
[33] D. Previtali, M. Longhi, F. Galli, A. D. Michele], F. Manenti,
M. Signoretto, F. Menegazzo, C. Pirola, Fuel 2020, 274, 117804.
R EF E RE N C E S
[34] B. Duhoux, R. T. Symonds, R. Hughes, P. Mehrani, E. J.
[1] J. J. McKetta, Encyclopedia of Chemical Processing and Anthony, A. Macchi, Can. J. Chem. Eng. 2020, 98, 75.
Design, CRC Press, New York 1980, p. 150. [35] D. W. Keith, G. Holmes, D. S. Angelo, K. Heidel, Joule 2018,
[2] D. C. Y. Foo, N. Chemmangattuvalappil, D. K. Ng, R. Elyas, 2, 1573.
C.-L. Chen, R. D. Elms, H.-Y. Lee, I.-L. Chien, S. Chong, [36] J. C. Goncalves, A. E. Rodrigues, Can. J. Chem. Eng. 2015, 93,
C. H. Chong, Chemical Engineering Process Simulation, 2205.
Elsevier, Amsterdam, The Netherlands 2017. [37] P. Lutze, R. Gani, J. M. Woodley, Chem. Eng. Process. Process.
[3] G. S. Patience, Can. J. Chem. Eng. 2018, 96, 2312. Intensif. 2010, 49, 547.
[4] A. W. Westerberg, H. P. Hutchinson, R. L. Motard, P. Winter, [38] E. Barbera, F. Mantoan, A. Bertucco, F. Bezzo, Can. J. Chem.
Process Flowsheeting, Cambridge Universities Press, Cam- Eng. 2020, 98, 1893.
bridge, UK 1979. [39] K. D. Dahm, R. P. Hesketh, M. J. Savelski, Chem. Eng. Educ.
[5] R. S. Mah, R. S. Mah, Eds., Chemical Process Structures and 2002, 36, 192.
Information Flows, Butterworth-Heinemann, Stoneham, MA [40] A. M. Borreguero, J. L. Valverde, J. M. García-Vargas, L.
1990, p. 125. Sánchez-Silva, Comput. Appl. Eng. Educ. 2019, 27, 1267.
[6] G. Buzzi-Ferraris, F. Manenti, Nonlinear Systems and Optimi- [41] S. Yang, L. Yang, C. He, Process Saf. Environ. Prot. 2001,
zation for the Chemical Engineer, Wiley-VCH, Weinheim, 79, 329.
Germany 2013. [42] C. Pirola, Educ. Chem. Eng. 2019, 28, 54.
[7] Aspen Plus 2004.1 getting started using equation oriented [43] T. M. Komulainen, R. Enemark-Rasmussen, G. Sin, J. P.
modeling, https://sites.ualberta.ca/CMENG/che312/F06Ch Fletcher, D. Cameron, Educ. Chem. Eng. 2012, 7, e153.
E416/HysysDocs/AspenPlus20041GettingStartedEOModeling [44] L. Calvo, C. Prieto, Educ. Chem. Eng. 2016, 17, 65.
(accessed: August 2020). [45] A. Rafael, F. Bernardo, L. Ferreira, M. Rasteiro, J. Teixeira,
[8] H.-S. Chen, M. Stadtherr, AIChE J. 1985, 31, 1843. Educ. Chem. Eng. 2007, 2, 20.
[9] A. C. Dimian, C. S. Bildea, A. A. Kiss, Applications in Design [46] L. Fregolente, H. Venturelli, J. Rodrigues, E. Da Silva, I.
and Simulation of Sustainable Chemical Processes, Elsevier, Diniz, R. Wolf Maciel Maria, Chem. Eng. Trans. 2018,
Amsterdam, The Netherlands 2019. 69, 427.
[10] R. Raoni, A. R. Secchi, E. C. Biscaia Jr., Comput. Chem. Eng. [47] R. M. Ghoniem, H. Abas, H. Bdair, Procedia Comput. Sci.
2017, 96, 169. 2018, 135, 490.
[11] F. Zhao, X. Chen, L. Zhu, AIChE J. 2017, 63, 2764. [48] C. Foreman, M. Hilditch, N. Rockliff, H. Clarke, A Compari-
[12] H. A. J. Watson, PhD thesis, Massachusetts Institute of Tech- son of Student Perceptions of Physical and Virtual Engineering
nology, Cambridge, MA 2018. Laboratory Classes, Springer International Publishing, Cham,
[13] D. Himmelblau, Chem. Eng. Sci. 1967, 22, 883. Switzerland 2020, p. 151.
[14] H. Weinblatt, J. ACM 1972, 19, 43. [49] C. Pirola, F. Galli, C. Peretti, Comput. Chem. Eng. 2020,
[15] R. Barkley, R. Motard, Chem. Eng. J. 1972, 3, 265. 140, 106.
[16] J. H. Christensen, D. F. Rudd, AIChE J. 1969, 15, 94. [50] P. R. C. Camargo, E. de Barros Basso, FCC Training
[17] R. S. Upadhye, E. A. Grens II, AIChE J. 1975, 21, 136. Simulator–Petrobras Experience, 2nd Mercosur Congress on
[18] R. L. Motard, A. W. Westerberg, AIChE J. 1981, 27, 725. Chemical Engineering, 4th Mercosur Congress on Process
[19] J. F. Boston, S. L. Sullivan Jr., Can. J. Chem. Eng. 1974, Systems Engineering, August 2005.
52, 52. [51] J. Isimite, F. Baganz, V. C. Hass, J. Chem. Technol. Biotechnol.
[20] R. B. Saeger, P. R. Bishnoi, Can. J. Chem. Eng. 1986, 64, 759. 2018, 93, 1529.
[21] J. Jelínek, Comput. Chem. Eng. 1988, 12, 195. [52] D. S. Patle, D. Manca, S. Nazir, S. Sharma, Virtual Reality
[22] J. H. Wegstein, Commun. ACM 1958, 1, 9. 2019, 23, 293.
[23] O. Orbach, C. M. Crowe, Can. J. Chem. Eng. 1971, 49, 509. [53] J. Haydary, Reactors, John Wiley & Sons, Ltd, Hoboken, NJ
[24] C. G. Broyden, Math. Comput. 1965, 19, 577. 2019, Ch. 5, p. 101.
[25] N. J. van Eck, L. Waltman, Scientometrics 2010, 84, 523. [54] S. S. Elnashaie, Modelling, Simulation and Optimization of
[26] Web of Science Core Collection 2020, http://apps. Industrial Fixed Bed Catalytic Reactors, Vol. 7, CRC Press,
webofknowledge.com (accessed: April 2020). Boca Raton, Florida 1994.
DE TOMMASO ET AL. 2319

[55] J. Haydary, Chemical Process Design and Simulation: Aspen [86] J. J. Quiroz-Ramirez, E. Sanchez-Ramirez, S. Hernandez-
Plus and Aspen Hysys Applications, John Wiley & Sons, Castro, J. G. Segovia-Hernandez, J. M. Ponce-Ortega, ACS
Hoboken, NJ 2019. Sustainable Chem. Eng. 2017, 5, 4018.
[56] M. A. Mirzakhani, N. Tahouni, M. H. Panjeshahi, Energy [87] J. C. Eslick, D. C. Miller, Comput. Chem. Eng. 2011, 35, 1488.
2017, 130, 382. [88] S. Sharma, G. P. Rangaiah, Comput. Chem. Eng. 2013, 56, 155.
[57] M. Sardarmehni, N. Tahouni, M. H. Panjeshahi, Energy 2017, [89] B. Linnhoff, J. R. Flower, AIChE J. 1978, 24, 633.
127, 623. [90] I. C. Kemp, Pinch Analysis and Process Integration: A User
[58] M. V. Kanischev, L. M. Ulyev, R. E. Chibisov, M. A. Vasilyev, Guide on Process Integration for the Efficient Use of Energy,
Chem. Eng. Trans. 2018, 70, 1099. Elsevier, Oxford 2011.
[59] Q. Zhang, M. Yang, G. Liu, X. Feng, J. Cleaner Prod. 2016, [91] B. Linhoff, J. Flower, AIChE J. 1978, 24, 642.
112, 4799. [92] K. Brown, Process Integration Initiative: A Review of the Pro-
[60] S. S. Chauhan, S. Khanam, Energy 2019, 173, 364. cess Integration Activities Funded under the Energy Efficiency
[61] R. R. Tan, D. C. Y. Foo, K. B. Aviso, D. K. S. Ng, Appl. Energy R and D Programme, Energy Technology Support Unit
2009, 86, 605. (ETSU), Harwell 1989.
[62] G. Venkatesh, Resour. Environ. Inf. Eng. 2019, 1, 1. [93] T. Umeda, J. Itoh, K. Shiroko, Chem. Eng. Prog. 1978, 74, 70.
[63] D. I. Garnett, G. S. Patience, Chem. Eng. Prog. 1993, 89, 76. [94] P. M. Pereira, M. C. Fernandes, H. A. Matos, C. P. Nunes,
[64] R. Williams Jr., Chem. Eng. 1947, 54, 124. Appl. Therm. Eng. 2017, 114, 523.
[65] G. S. Patience, D. C. Boffito, J. Adv. Manuf. Process. 2020, 2, [95] E. Rosen, Computer Applications to Chemical Engineering
e10039. (Eds: R. G. Squires, G. V. Reklaitis), ACS Publications 1980,
[66] H. J. Lang, Chem. Eng. 1947, 54, 27. Ch 1. https://doi.org/10.1021/bk-1980-0124
[67] H. J. Lang, Chem. Eng. 1947, 54, 130. [96] N. M. Kazerooni, R. Eslamloueyan, M. Biglarkhani, Int.
[68] H. J. Lang, Chem. Eng. 1948, 55, 112. J. Ind. Chem. 2019, 10, 107.
[69] K. M. Guthrie, Chem. Eng. 1969, 76, 114. [97] Y. Xu, H. Dinh, Q. Xu, F. T. Eljack, M. M. El-Halwagi,
[70] K. Guthrie, Process Plant Estimating, Evaluation, and Control, J. Cleaner Prod. 2020, 120, 129.
Craftsman Book Company of America, Carlsbad, CA 1974. [98] P. Thomas, S. Bhattacharyya, A. Patra, G. Rao, Comput.
[71] G. D. Ulrich, P. T. Vasudevan, Chemical Engineering Process Chem. Eng. 1998, 22, 1515.
Design and Economics: A Practical Guide, Process Publishing, [99] M. Saidi, H. Kadkhodayan, J. Environ, Chem. Eng. 2020, 8,
Durham, NH 2007. 103772.
[72] G. D. Ulrich, P. T. Vasudevan, Chem. Eng. 2006, 66. [100] H. Jafari, A. Sheikhi, R. Sotudeh-Gharebagh, Can. J. Chem.
[73] R. A. Turton, J. A. Shaeiwitz, D. Bhattacharyya, W. B. Eng. 2020, 98, 1003.
Whiting, Analysis, Synthesis and Design of Chemical Pro- [101] Y. Shen, Z. Tang, D. Guo, C. Wu, D. Wang, A. Jiang, Y. Wen,
cesses, 5th ed., Pearson Education, Ann Arbor, MI 2018. Appl. Therm. Eng. 2019, 160, 113983.
[74] Y. Feng, G. P. Rangaiah, Chem. Eng. 2011, 22. [102] NIST Standard Reference Database Number 69, https://
[75] W. D. Seider, J. D. Seader, D. R. Lewin, Product & Process webbook.nist.gov/chemistry/ (accessed: August 2020).
Design Principles: Synthesis, Analysis and Evaluation, John [103] Engineering data and technical references to solve engineer-
Wiley & Sons, Hoboken, NJ 2009. ing problems, https://app.knovel.com/web/index.v?jsp=main
[76] R. Sinnott, G. Towler, Chemical Engineering Design: SI Edi- (accessed: August 2020).
tion, Butterworth-Heinemann, Oxford 2019. [104] B. E. Poling, J. M. Prausnitz, J. P. O’Connell, The Proper-
[77] N. Patel, N. Padhiyar, Chem. Eng. Res. Des. 2017, 119, 160. ties of Gases and Liquids, 5th ed., McGraw-Hill, New York
[78] Practical optimization: a gentle introduction, https://eclass. 2001.
gunet.gr/modules/document/file.php/LABGU361/PRACTICAL [105] T. Jensen, A. Fredenslund, P. Rasmussen, Ind. Eng. Chem.
%20OPTIMIZATION%20GENTLE%20INTRODUCTION% Fundam. 1981, 20, 239.
202015.pdf (accessed: August 2020). [106] W. Wakeham, M. Assael, Bulg. Chem. Commun. 2019, 51, 9.
[79] L. T. Biegler, Nonlinear Programming: Concepts, Algorithms, [107] K. C. Chao, J. D. Seader, AIChE J. 1961, 7, 598.
and Applications to Chemical Processes, Vol. 10, SIAM, [108] B. I. Lee, M. G. Kesler, AIChE J. 1975, 21, 510.
Philadelphia, PA 2010. [109] O. Redlich, J. N. S. Kwong, Chem. Rev. 1949, 44, 233.
[80] R. A. Stubbs, S. Mehrotra, Math. Program. 1999, 86, 515. [110] D.-Y. Peng, D. B. Robinson, Ind. Eng. Chem. Fundam. 1976,
[81] J. A. Nelder, R. Mead, Comput. J. 1965, 7, 308. 15, 59.
[82] C. A. C. Coello, G. B. Lamont, D. A. Van Veldhuizen, Evolu- [111] P. A. Gupte, P. Rasmussen, A. Fredenslund, Ind. Eng. Chem.
tionary Algorithms for Solving Multi-Objective Problems, Fundam. 1986, 25, 636.
Vol. 5, Springer, New York 2007. [112] G. M. Wilson, J. Am. Chem. Soc. 1964, 86, 127.
[83] J. M. Ponce-Ortega, L. G. Hernandez-Perez, Optimization of [113] H. Renon, J. M. Prausnitz, AIChE J. 1968, 14, 135.
Process Flowsheets through Metaheuristic Techniques, [114] D. S. Abrams, J. M. Prausnitz, AIChE J. 1975, 21, 116.
Springer, New York 2019. [115] A. Fredenslund, R. L. Jones, J. M. Prausnitz, AIChE J. 1975,
[84] F. Galli, C. Pirola, D. Previtali, F. Manenti, C. L. Bianchi, 21, 1086.
J. Cleaner Prod. 2018, 171, 147. [116] R. Gani, J. P. O’Connell, Comput. Chem. Eng. 1989, 13, 397.
[85] J. J. Quiroz-Ramirez, E. Sanchez-Ramirez, S. Hernandez, [117] E. C. Carlson, Chem. Eng. Prog. 1996, October, 35.
J. H. Ramirez-Prado, J. G. Segovia-Hernandez, Ind. Eng. [118] J. G. Hayden, J. P. O’Connell, Ind. Eng. Chem. Process Des.
Chem. Res. 2017, 56, 1823. Dev. 1975, 14, 209.
2320 DE TOMMASO ET AL.

[119] J. Gmehling, M. Kleiber, B. Kolbe, J. Rarey, Chemical Ther- Glassey, A. Stankiewicz, H. van den Berg, F. Galli, K. Ogden,
modynamics for Process Simulation, John Wiley & Sons, R. Santos, A. Smith-Schoettker, J. Cantin, H. Akse, G.
Hoboken, NJ 2019. Patience, J. M. Commenge, J. Harmsen, J. L. Dubois, N.
[120] S. T. Holland, R. Abbas, D. Hildebrandt, D. Glasser, Ind. Eng. Afraz, R. Bos, Educ. Chem. Eng. 2020, 32, 15.
Chem. Res. 2010, 49, 327. [125] CHEMSAFE, https://dechema.de/en/chemsafe.html (accessed:
[121] R. C. Pattison, M. Baldea, AIChE J. 2014, 60, 4104. August 2020).
[122] Chemistry effects predicting simulation software, https:// [126] A. Kummer, T. Varga, J. Loss Prev. Process Ind. 2019, 58, 22.
www.ansys.com/products/fluids/ansys-chemkin-pro [127] J. E. Coon, M. Kush, M. C. Rowland, J. R. Cunningham,
(accessed: August 2020). presented at the AIChE Spring Meeting, New Orleans, US,
[123] D. F. Rivas, D. C. Boffito, J. A. F. Albanese, J. Glassey, N. Afraz, March 1998.
H. Akse, K. Boodhoo, R. Bos, J. Cantin, Y. W. Chiang, J. M.
Commenge, J. L. Dubois, F. Galli, J. Harmsen, S. Kalra, F. Keil,
R. Morales-Mendez, F. J. Navarro-Brull, T. Noel, G. Patience, D. How to cite this article: De Tommaso J, Rossi F,
Reay, R. Santos, A. Smith-Schoettker, A. Stankiewicz, H. van
Moradi N, Pirola C, Patience GS, Galli F.
den Berg, T. V. Gerven, J. van Gestel, M. van der Stelt, M. van
de Ven, R. Weber, Educ. Chem. Eng. 2020, 32, 1.
Experimental methods in chemical engineering:
[124] D. F. Rivas, K. Boodhoo, R. Morales-Mendez, Y. W. Chiang, Process simulation. Can J Chem Eng. 2020;98:
F. Navarro, T. V. Gerven, J. van Gestel, R. Weber, S. Kalra, 2301–2320. https://doi.org/10.1002/cjce.23857
D. C. Boffito, D. Reay, J. A. F. Albanese, F. Keil, T. Noel, J.