Thermodynamic analysis of the driving force approach: Reactive systems

Computers and Chemical Engineering 0 0 0 (2019) 106509

Contents lists available at ScienceDirect

Computers and Chemical Engineering journal homepage: www.elsevier.com/locate/compchemeng

Thermodynamic analysis of the driving force approach: Reactive systems Teresa Lopez-Arenas a, Mauricio Sales-Cruz a, Rafiqul Gani b,d, Eduardo S. Pérez-Cisneros c,∗ a

Departamento de Procesos y Tecnología, Universidad Autónoma Metropolitana Cuajimalpa, Avenida Vasco de Quiroga 4871, C.P. 05348 Ciudad de México, México b PSE for SPEED Company Ltd., Skyttemosen 6, DK-3450 Allerod, Denmark c Departamento de Ingeniería de Procesos e Hidráulica, Universidad Autónoma Metropolitana-Iztapalapa, Av. San Rafael Atlixco 186, C.P. 09340 Ciudad de México, México d College of Control Science and Engineering, Zhejiang University, Hangzhou, China

a r t i c l e

i n f o

Article history: Received 28 February 2019 Revised 9 June 2019 Accepted 30 June 2019 Available online xxx Keywords: Driving force approach Thermodynamic analysis Reactive systems Element concept

a b s t r a c t A thermodynamic analysis of the driving force based approach for design of reactive systems has been performed. The thermodynamic analysis is based on the element concept for the solution of the chemical and phase equilibrium. Through the thermodynamic analysis, the fundamental relationships between the element reactive driving force and the energy involved in achieving a desired reactive separation, are highlighted. In this work, the explicit relationships of the element reactive driving force with the heat of vaporization and the reaction enthalpies of the reacting mixture are shown. Validation of the thermodynamic analysis is carried out considering three binary element reactive systems: (i) the isomerization of n-butane reaction, (ii) the dimerization of propylene reaction, both reactive systems consider an inert component, and (iii) the reactive system representing the MTBE production reaction. © 2019 Elsevier Ltd. All rights reserved.

1. Introduction The modeling and design of reactive separation processes has been thoroughly studied since the two well-known successful industrial applications: the production of methyl acetate from the esterification of acetic acid with methanol (Agreda et al., 1990) and the production of MTBE from the etherification of isobutene with methanol (Smith and Huddleston, 1982). The methods available for synthesis and design of reactive distillation columns can be divided into two main categories: optimization methods and conceptual design methods, with the latter including various graphical methods. Reactive distillation design via optimization methods is based on the main idea that the system can be represented as a special case of a multiphase reactor network. In a single reactive distillation column, reaction and distillation occur simultaneously on every stage and the flows are counter-currently contacted. Therefore, synthesis of reactive distillation columns is strongly related to the synthesis of multiphase reactor networks including reaction and separation. Ciric and Gu (1994) formulated a mixed integer non-linear programming model (MINLP), incorporating reaction kinetics, heats of reaction and liquid holdup volumes,

∗

Corresponding author. E-mail address: [email protected] (E.S. Pérez-Cisneros).

https://doi.org/10.1016/j.compchemeng.2019.06.034 0098-1354/© 2019 Elsevier Ltd. All rights reserved.

solution of which yields the optimal number of stages, the optimal feed stage location(s) and feed rates and the optimal reflux ratio at a minimal total annualized cost. The optimization methods proved to be viable methods for the design of reactive distillation systems. They are able to provide solutions which are very close to the global optimum, even for highly non-ideal systems with multiple reactions. However, for the synthesis and design of a single reactive distillation column, the optimization methods are very time consuming, and they are more appropriate to analyze a wide variety of flowsheet configurations. Till now, computers have superseded graphical techniques as the main tool for distillation design and performance evaluation. Nevertheless, graphical techniques are still widely used in modern distillation technology. They provide a means of visualizing the process and enable spotting of pinched conditions, excessive reflux, incorrect feed points, and a non-optimum thermal condition of the feed. They are powerful for analyzing computer solutions. Other applications are the screening and optimization of design options, providing initial estimates for computer calculations. When such graphical techniques and the computational procedures are combined for the design of distillation columns, the visualization and fast computation become to be a powerful hybrid method for the design and evaluation of many reactive separation processes. On the other side, it is well known that the thermodynamic analysis of a distillation column is important for synthesizing and developing energy efficient

2

T. Lopez-Arenas, M. Sales-Cruz and R. Gani et al. / Computers and Chemical Engineering 000 (2019) 106509

Nomenclature A AI AII Aji A B bβ bfeed bj Cpi E FDW G H Hj Ky Kx M NC NR n ni P Q R S T U Wj Xj xi Yj yi Z

formula matrix partition sub-matrix I (Eq. (19)) partition sub-matrix II (Eq. (19)) number of reaction invariant elements j in molecule i element A element B elemental composition vector in phase β elemental composition vector fed amount of element j heat capacity of species i matrix E (Eq. (29)) element reactive driving force Gibbs energy matrix H (Eq. (28)) element-based enthalpy of element j chemical equilibrium constant in the vapor phase chemical equilibrium constant in the vapor phase total number of elements representing the reactive system total number of component present in the reactive system total number of independent chemical reactions molar component vector moles of species i in the reactive mixture total pressure matrix Q (Eq. (30)) universal gas constant matrix S (Eq. (20)) temperature matrix U (Eq. (32)) elemental mole fraction of element j element composition j in the liquid phase mole fraction of component i in the liquid phase element composition j in the vapor phase mole fraction of component i in the vapor phase matrix Z (Eq. (20))

Greek letters α AB elemental reactive relative volatility β phase β (vapor or liquid) γi activity coefficient of species i μi chemical potential of species i λj elemental chemical potential of element j ϕi fugacity coefficient of species i AB elemental non-ideal phases contribution π parameter (T, P) for derivative of Eq. (27) Фj elemental non-ideal contribution of element j in the vapor phase Гj elemental non-ideal contribution of element j in the liquid phase

distillation processes (Bandyopadhyay, 2002). It allows the thermodynamic efficiency of the process to be quantified, regions with poor energy efficiency to be identified, and the thermodynamic targets to be defined. Several studies show the application of exergy analysis to improve the thermodynamic efficiency of a distillation column (Atkinson, 1987; Ratkje et al., 1995; Taprap and Ishida, 1996; Agrawal and Herron, 1997). Atkinson (1987) has developed a graphical representation of exergy loss in a distillation column. Ratkje et al. (1995) have analytically shown that

entropy generation for a distillation column is at a minimum when the driving force for separation is distributed uniformly along the column. Taprap and Ishida (1996) have presented different exergy losses in a distillation column on energy-utilization diagrams. These diagrams identify the amount of energy transformation and exergy loss of individual process steps. Agrawal and Herron (1997) have given equations to quantify thermodynamic efficiency of a distillation column that separates an ideal binary mixture, with constant relative volatility, into pure components. Also, Blahusiak et al. (2016) proposed a quantitative efficiency analysis based on exergy analysis. They analyzed distillation of near-ideal binary systems from a heat engine perspective (considering distillation as a process creating separation work) and investigated the influence of the feed composition and thermal properties of separated compounds on the internal efficiency of the heat engine. Specifically, for reactive systems, most of the existing work related to graphical design of RDC is based on the use of transformed composition variables proposed by Doherty and Buzad (1992). Also, the design methods based on the element concept (Perez-Cisneros et al., 1997) shown their capability for RDC design as they are similar to those design methods used for non-reactive systems, i.e., McCabe–Thiele and Ponchon–Savarit graphical methods. An extension of the element concept for the energy-efficient design of RD columns through the driving force based approach was proposed by Sanchez-Daza et al. (2003) for the production of benzene from the disproportionation of toluene reaction. Recently, Mansouri et al. (2015) have shown how the driving force concept can be used for integrated process design and control. That is, from a process design point of view, optimal/near optimal design in terms of energy consumption is obtained at the highest driving force; and from a controller design point of view, the best controller structure and set-point values for controlled and manipulated variables are obtained at this point. With an analytical analysis it is demonstrated that at the maximum driving force, the sensitivity of the controlled variables to disturbances is the lowest and at the same time, the sensitivity of controlled variables to manipulated variables (actuators) is the highest. Despite several applications of the driving force based approach for the design and control of reactive distillation columns, an explanation of the theory based on fundamental thermodynamic analysis is still not established. Therefore, the objective of this work is to extend the previous thermodynamic analysis of the driving force concept applied to non-reactive systems (Perez-Cisneros and Sales-Cruz, 2017) to also cover reactive systems where the chemical and physical equilibrium systems are represented by thermodynamically consistent elements.

2. Thermodynamic fundamentals The element concept is derived from chemical model theory, where the equations of chemical equilibrium together with any appropriate physical model yielding the chemical potentials are embedded into an element-based model, called the chemical model (Michelsen, 1995; Perez-Cisneros et al., 1997). The solution of the M chemical model equations (M being the minimum number of elements representing the reactive system) together with the condition of equilibrium (equality of the component chemical potentials in all co-existing phases) provides the element phase compositions for the reactive system. One attractive feature of this concept is its capability to handle the problem of reactive-phase equilibrium in the same manner as the case when no reactions are taking place in the system. That is, this approach reduces the chemical and physical equilibrium problem to an identical physical equilibrium problem for a mixture of elements representing the system.

T. Lopez-Arenas, M. Sales-Cruz and R. Gani et al. / Computers and Chemical Engineering 000 (2019) 106509

2.1. Chemical and physical equilibrium For the thermodynamic analysis of the driving force based approach for design of reactive systems let us start with the basic statement of the chemical and physical equilibrium problem. For a system with NP phases and NC chemical species, the fundamental thermodynamic relationship is given by the Gibbs free energy as:



β

G = G T , P, ni



(1)

where ni β (i = 1,2,…,NC; β = 1,2,…,NP) represents the number of moles of species i in phase β . The Gibbs free energy is an extensive property, proportional to the amount of material in the system. This implies that Eq. (1) must be a homogeneous function of degree one in ni β :



β

G T , p, kni





β

= kG T , p, ni



(2)

From Euler’s theorem (Tester and Modell, 1997) on homogeneous functions we have

G=

NP  NC  μβi nβi

(3)

β =1 i=1

where the chemical potential μβ i is defined as

  ∂G β β β ≡ μ = μ T , P, n i i j ∂ nβi i

(5)

β =1 i=1

subject to the M constraints: NP  NC 

β =1 i=1

β

A ji ni − b j = 0

j = 1, 2, . . . , M



NP  NC M NP  NC    β β β Lˆ = n i μi − λj A ji ni − b j

β =1 i=1



(7)

β =1 i=1

j=1

and the necessary conditions for a stationary point are M  ∂ Lˆ β = μi − A ji λ j = 0 i = 1, 2, . . . , NC; β = 1, 2, . . . , NP β ∂ ni j=1

(8) and NP  NC  ∂ Lˆ β =− A ji ni + b j = 0 ∂λj β =1 i=1

j = 1, 2, . . . , M

(9)

Eqs. (8) and (9) constitute the set of equations which determine an equilibrium state and combined with any thermodynamic model that provides the chemical potentials, given T, P and b, can be solved for the NC∗ NP+M unknowns, ni β (i = 1,2,…,NC; β = 1,2,…,NP) and λj (j = 1,2,…,M). 2.2. The chemical model

is used to describe a thermodynamic system, the condition for thermodynamic equilibrium of a closed system is defined as the state for which the total Gibbs free energy attain its minimum with respect to all possible composition changes at the given T, P and coexisting phases. That is, NP  NC  β β n i μi

The solution of the constrained optimization problem given by Eqs. (5) and (6) can be obtained through the Lagrange multiplier formulation. The corresponding Lagrangian function Lˆ is defined as:

(4)

which is a homogeneous function of degree zero in nj β ; that is, μβ is an intensive property. When the Gibbs function (Eq. (3))

min G(n ) =

3

(6)

where G(n) is the total Gibbs free energy of a system containing NC species and NP phases. Eq. (6) represents the M independent element mass balances, where the coefficients Aji represent the number of reaction invariant elements j in molecule i, with the formula matrix A as a full rank matrix of M× NC elements and bj is the total number of gram-atoms (moles) of element j in the system. It should be noted that the total number of independent elements (M), which can be atoms, molecules or groups (radicals) is smaller than the number of components (NC) in the reactive system. That is, an element in the chemical model theory may be an invariant atom (i.e., carbon, hydrogen, oxygen, etc.), an invariant group or radical (i.e., OH− , H+ , etc.) or a complete molecule (i.e., isobutene, methanol, etc.). A formula matrix A can be defined with coefficients Aij , where a molecule i (columns) contains Aij invariant elements j (rows). Whether the elements should be selected as atoms, groups or compounds depend on the simplicity of the reactive system and, in order to simplify the matrix operations, the selection of elements as compounds is highly recommended. It should be pointed out that, if chemical reactions are not occurring, the number of invariant elements representing the non-reactive system is the same as that of the number of species present in the mixture and it leads to a symmetrical square formula matrix A (NC × NC).

Michelsen and Mollerup (2004) proposed a formulation of the chemical-physical equilibrium problem through the use of chemical models. This concept is derived from chemical theory, where the equations of chemical equilibrium together with any appropriate physical model yielding the chemical potentials are incorporated into an element-based model (called the chemical model). The solution of the chemical model equations together with the equilibrium condition (equality of the component chemical potentials in all co-existing phases) provides the element phase compositions for the reactive system. In this formulation the relationship between the Gibbs free energy and the Lagrange multipliers is exploited. In order to explain this approach let us consider only chemical equilibrium for a single phase, that is β =1, and for simplicity, let us drop out the superscript β in Eq. (3). By introducing Eqs. (8) and (9), that is, the stationary point conditions, into Eq. (3) the following Gibbs free energy at equilibrium is obtained:

Geq =

NC 

n i μi =

i=1

NC  i=1

 ni

M  j=1

 A ji λ j

=

M 

b jλ j

(10)

j=1

In addition, it is known that the equilibrium Gibbs energy is a homogeneous function of degree one in the element composition b (Tester and Modell, 1997), i.e.

Geq (kb ) = kGeq (b )

(11)

thus,

∂ Geq = λj ∂bj

(12)

This relation between the elemental composition vector b and the Lagrange multiplier vector λ is equivalent to the relation between the molar composition vector n and the chemical potential vector μ, and a fully consistent thermodynamic description of a chemically equilibrated phase is obtained in terms of b as the (element) composition vector and λ as the corresponding element chemical potential vector. Extending the above formulation to consider NP phases, the multiphase equilibrium conditions now become

      λ(1) b(1) = λ(2) b(2) = · · · = λ(NP ) b(NP )

(13)

4

T. Lopez-Arenas, M. Sales-Cruz and R. Gani et al. / Computers and Chemical Engineering 000 (2019) 106509

with NP 

bβ = b f eed

(14)

β =1

Substitution of the mole vector n (Eq. (21)) into Eq. (3) enable us to formulate the equilibrium problem in terms of the NR composition variables nII as

min G(nI nII ) = min G(nI (nII ), nII )

And in analogy with the conventional phase equilibrium condition,

μ (1 )



n (1 )



=

μ (2 )



n (2 )



= ··· =

μ(NP )



n(NP )



(15)

or also formulated as :

ZT μ =

NC 

Zk,i μi (nII ) = 0

for k = 1, 2, . . . , NR

(25)

i=1

and NP 

(24)

nβ = n f eed

(16)

β =1

An interesting feature of the chemical model approach should be noted. The final set of Eqs. (8), (13), and (14) are formally identical to the set of the equations solved for non-reactive phase equilibrium calculations, except for the use of elements. Therefore, the same method of solution can be employed in reactive as well as non-reactive phase equilibrium calculations. Also, it is possible to define element “mole” fractions in a similar manner as in non-reactive systems for compound mole fractions. The total element amount in any phase β is given by β

bT =

M  β bj

(17)

j=1

Thus, the element “mole” fractions are given as β

bj

Wj = β = β M bT k=1 bk

(18)



(19)

where the leading M × M block AI is assumed to be non-singular. Since A is of full rank, a non-singular leading submatrix can always be constructed by permuting the columns of A. We further select a vector n such that An = b. Following Fletcher (1981, Chapter 10), we can construct two matrices S and Z, given by



S=



A−1 I , 0



Z=

−A−1 AII I I



(20)

where the dimension of S is NC × M and that of Z is NC × (NC−M=NR). NR is the total number of independent chemical reactions. In order to compute the element chemical potential for each phase, it is necessary to solve the chemical equilibrium problem for the respective phase as indicated by Eq. (8). This can be done simply by considering the following partition for the composition vector n into two vectors



n=

nI ··· nII

(21)

where the vector nI has the first M elements and the vector nII contains the remaining NC−M =NR components. With this partition the element conservation equation (Eq. (9)) can be written as

A n = A I nI + A II nII = b

(22)

and from this expression we obtain nI as

nI = A−1 I (b − AII nII )

S j,i μi

for j = A, B, . . . , M

(26)

i=1

The computational algorithm to determine the element chemical potentials, λj , has been reported by Perez-Cisneros (1997). The derivatives with respect to the unknown variables (b in this case) are obtained by differentiating Eq. (8) with respect to a parameter, π , where π may represent temperature, T, pressure, P, or an element of the b-vector. Let π represents T or P. Then,

∂λ ∂μ = −Q ∂π ∂π

(27)

∂λ = −H ∂b 

In general the formula matrix A can be partitioned in two submatrices as follows:



M 

(28)

Explicit expressions for the submatrices E, Q, and H using the matrices Z and S from Eq. (20) are given by Fletcher (1981):

2.3. The element chemical potential (λ) and its derivatives

. A = AI .. AII

λj =

and the composition derivatives are given by

β

bj

β

The solution of the above unconstrained optimization problem yields values for the component chemical potentials μ from which the element chemical potentials can subsequently be calculated as,

(23)

E = Z Z TU Z

−1

ZT

(29)

Q T = S − EU S

(30)

H = −QU S

(31)

with

Ui j =



∂ μi ∂nj

 (32) T,P

where the dimension of matrices E, QT and H are NC × NC, NC × M and M × M, respectively. It should be clear that, in the computation of λ the evaluation of the chemical potentials for each species is required. 3. The element driving force based approach for design of reactive systems Similar to the extended equilibrium condition for non-reactive systems (Perez-Cisneros and Sales-Cruz, 2017), it is possible to write for a binary element system (A and B) an extended chemical and phase equilibrium condition in terms of the element chemical potentials as:

 ideal  ideal  nonideal  nonideal  − λlA + λvA − λlA λvA  ideal  ideal  nonideal  nonideal  − λvB − λlB + λvB − λlB =0

(33)

where A refers to the most volatile element. The connection between the element chemical potentials and the chemical potentials for the different species is given by Eq. (26) with the terms of the S matrix as:

T. Lopez-Arenas, M. Sales-Cruz and R. Gani et al. / Computers and Chemical Engineering 000 (2019) 106509

λA =

M 

Thus, the ideal element chemical potential can be written as:

SA,i μi = SA,1 μ1 + SA,2 μ2 + · · · + SA,M μM

(34)

i=1

and

λB =

M 

SB,i μi = SB,1 μ1 + SB,2 μ2 + · · · + SB,M μM

(35)

i=1

And the matrices S and Z (Eq. (20)) become for two elements and NC compounds as:

⎡

⎡

Z 1,1 ⎢ Z 2,1 ⎢ . ⎢ .. ⎢ ⎢ZM,1 Z=⎢ ⎢ 1 ⎢ ⎢ 0 ⎢ . ⎣ .. 0

⎤

SA,1 SA,2 SB,2 ⎥ ⎢ SB,1 S=⎣ ⎦, 0 0 SNC,1 = 0SNC,2 = 0

Z 1,2 Z 2,2 .. . ZM,2 0 1 .. . 0

··· ··· ··· 0 0 .. . 0

⎤

Z1,NR Z2,NR ⎥ .. ⎥ ⎥ . ⎥ ZM,NR ⎥ ⎥ 0 ⎥ ⎥ 0 ⎥ .. ⎥ ⎦ . 1 (36)

Thus, we can define the ideal contributions of the element chemical potentials as:



M 

RT ln YA = RT

SA,i ln(yi )

λoAvap =



SA,i μvi 0 = SA,1 μv10 + SA,2 μ2



M 

RT ln PT = RT



M 

RT ln XA = RT

SA,i ln (P )

M 

= RT [SA,1 ln(x1 ) + SA,2 ln(x2 )]



i=1



M 

RT ln PSA = RT

SA,i ln(PiS )

M 

RT ln YB = RT



λoB vap =

RT ln PT = RT



M 

= RT [SB,1 ln(P ) + SB,2 ln(P )]

SB,i ln(xi )



RT

= RT

M 

SB,i ln(

−



λvB − λlB

nonideal

=0

(40)

A

YB



+ RT ln

where

RT ln(Enonideal )= AB



XB PsB XA PsA

λvA − λlA



+ RT ln(Enonideal )=0 AB

nonideal

−



λvB − λlB

(41)

nonideal

(42)

i=1

i=1

i=1

i=1

 v nonideal  l nonideal ˆ A ) − RT ln(A ) − λA = RT ln(

λA M 

SA,i ln(ϕˆ i ) − RT

M 

SA,i ln(γi )

i=1

M 

SB,i ln(ϕˆ i ) − RT M 

 SA,i ln

M 

SB,i ln(γi )

i=1

ϕˆ i γi

 − RT

M  i=1

 SB,i ln

ϕˆ i γi

 (44)

It should be noted that the element mole fractions (WA v , WB v , WA l , WB l ) cannot be explicitly introduced into Eq. (41) since the reactive system has not been defined yet through the corresponding element matrix. For the non-reactive case (Gani and Bek-Pedersen, 2001), the driving force is defined as the most volatile two-phase composition difference FD =y1 −x1 , in the case of reactive systems the element driving force definition FDW for a system in chemical and physical equilibrium can be established as the two-phase element composition difference:



PiS

Y 

i=1

= RT [SB,1 ln(x1 ) + SB,2 ln(x2 )]

l0 l0 SB,i μl0 i = SB,1 μ1 + SB,2 μ2



ln PSB

RT ln

RT ln (EAB ) = RT

i=1

nonideal

i=1

i=1

λoB liq =

λvA − λlA

Notice that the element reference chemical potentials of the liquid and vapor phases are considered similar. Eq. (40) can be written in condensed form as:

= RT

i=1

M 

RT ln XB = RT

 v0

SB,i ln(P )



 v nonideal  l nonideal ˆ B ) − RT ln(B ) − λB = RT ln(

λB

= RT [SB,1 ln(y1 ) + SB,2 ln(y2 )]

SA,i μvi 0 = SB,1 μv10 + SB,2 μ2 M 

+

i=1





RT ln(YA PT ) − RT ln(XA PsA ) − RT ln(YB PT ) + RT ln(XB PsB )

= RT

i=1

and introducing Eqs. (39) into (33) we can write:

and the non-ideal element chemical potential differences are:

= RT SA,1 ln(P1S ) + SA,2 ln(P2S )

SB,i ln(yi )

(39)

(43)

i=1 M 

λ0Bliq + RT ln(XB ) + RT ln(PsB )

i=1

(37)



=

i=1

i=1

and for element B as:

λ0Bvap + RT ln(YB ) + RT ln(PT )

M M    l nonideal nonideal = RT ln(B ) = SB,i (μvi ) = RT SB,i ln(γi ) λB





=

i=1

l0 l0 SA,i μl0 i = SA,1 μ1 + SA,2 μ2

λ0Aliq + RT ln(XA ) + RT ln(PsA )

M M    v nonideal nonideal ˆ B) = = RT ln(

SB, j (μvj ) = RT SB,i ln(ϕˆ i ) λB

i=1

λoAliq =

=

M M    l nonideal nonideal = RT ln(A ) = SA,i (μli ) = RT SA,i ln(γi ) λA

= RT [SA,1 ln (P ) + SA,2 ln (P )]

SA,i ln(xi )

λ0Avap + RT ln(YA ) + RT ln(PT )

i=1

i=1

=

M M    v nonideal nonideal ˆ A) = = RT ln(

SA,i (μvi ) = RT SA,i ln(ϕˆ i ) λA

 v0

i=1

 v ideal λA  l ideal λA  v ideal λB  l ideal λB

For a ϕ −γ approach we have for the non-ideal contributions to the element chemical potentials:

= RT [SA,1 ln(y1 ) + SA,2 ln(y2 )]

i=1 M 

5

  ) = RT SB,1 ln(P1S ) + SB,2 ln(P2S )

i=1

(38)

FDW = WAv − WAl

(45)

6

T. Lopez-Arenas, M. Sales-Cruz and R. Gani et al. / Computers and Chemical Engineering 000 (2019) 106509

Thus, Eq. (41) can now be written in a general form in terms of element mole fractions as:



  XB WAl PsB YA {FDW }   RT ln + RT ln YB {FDW } XA WAl PsA 





+ RT ln(Enonideal )=0 AB

{FDW }

where YA and YB are functions of the element driving force FDW and XA {WA l } and XB {WA l } are functions of the liquid element mole fraction WA l . It should be noted that an explicit relationship of the reactive driving force based approach cannot be obtained until the reactive system is defined. Also, it should be pointed out that the driving force approach is not limited to binary component non-reactive mixtures or binary element reactive mixtures. In fact, the extended chemical and physical equilibrium condition (Eq. (33)) can be written for a multi-element reactive system as:

 ideal  ideal  nonideal  nonideal  − λlA + λvA − λlA λvA  ideal  ideal  nonideal  nonideal  − λvB − λlB + λvB − λlB  ideal  ideal  nonideal  nonideal  − λCv − λCl + λCv − λCl ··· − ···  ideal  ideal  nonideal  nonideal  − λlM + λvM − λlM =0 − λvM

=

+

V,ideal HW =

j=1

W jv H 0j V

T0

i=1



W jL H 0j V − HVap j

C pig dT i

(50)



j=1

HVap = j

NC 

Q j,i hVap i

(51)

i=1

Where Hj 0V and Hj vap are the ideal gas enthalpy and heat of vaporization of element j at temperature T, respectively. It should be noted that the Q matrix coefficients are given by Eq. (30) (see Appendix C). By subtracting Eqs. (48) and (49) and introducing Eqs. (50) and (51), the following expression is obtained: V,nonideal L,nonideal V,ideal L,ideal R E HW − HW = HW − HW + HW − HW M 

W jv H 0j V −

M 





R E W jL H 0j V − HVap + HW − HW j

(52)

j=1







R E + WBL HB0V − HBVap + HW − HW

where R HW = WAv hRA + WBv hRB ;

hRB = −RT 2

 QB,i

hEB = −RT 2

NC 

 QB,i

i=1

NC 

hRA = −RT 2

∂ ln(ϕˆ i ) ∂T

E HW = WAl hEA + WBl hEB ;

hEA = −RT 2

∂ ln(γi ) ∂T

NC 

 QA,i

i=1



 (53)

 QA,i

i=1



 ∂ ln(ϕˆ i ) ; ∂T

 ∂ ln(γi ) ; ∂T (54)

Reorganizing Eq. (53) we obtain:



V,nonideal L,nonideal HW − HW = WAv − WAL





HA0V − HB0V + WAL HAVap

R E + WBv HBVap + HW − HW

(55)

and introducing the element driving force (FDW ) definition into Eq. (55), the following relationship between the element driving force and the reactive phase enthalpy change is obtained: V,nonideal L,nonideal HW − HW



FDW =

HA0V − HB0V







−

WAL HAVap + WBv HBVap



 R  E HW − HW  −  0V 0V

HA0V − HB0V





(56)

HA − HB

(48) (49)

NC  i=1



where HW R and HW E are the element based residual and excess non-ideal contributions. The element-based ideal enthalpies are obtained using the relationship between the derivatives of the element chemical potential with respect to temperature, that is, the element-based enthalpy related to the enthalpies of the different species (see Eq. (27) and Appendix B): M 

T

V,nonideal L,nonideal HW − HW = WAv HA0V + WBv HB0V + WAL HA0V − HAVap

Similar to the non-reactive case (Thompson, 20 0 0), it is possible to define the enthalpy for a reactive non-ideal vapor–liquid mixture based on element composition as:

E HW

M 

 Q j,i

and for a binary element mixture we have from Eq. (52)

4. The element driving force and the reactive phase energy change

L,ideal HW

NC 

Q j,i h0i V =

j=1

We know that the driving force is calculated for a binary pair of compounds or elements. In a binary mixture, there is only one binary pair. In a ternary mixture, there are three binary pairs but only two adjacent pairs that need to be considered. Thus, for a NC component non-reactive mixture or M elements for a reactive mixture, there are NC-1 or M-1 adjacent binary pairs and the easiest separation is the one corresponding to the binary pair with the highest driving force. Also, it has been shown that the chemical model can be used for reactive systems considering multiple products and different reaction routes. In a previous work (Gani et al., 1998) the element-based approach for the solution of the phase and chemical equilibrium problem was applied in the modeling of an industrial gas–liquid reactor considering both, kinetically controlled and equilibrium reactions, simultaneously. The reactive system consisted of 13 components and 9 elements were used to represent the chemical reactions. Thus, the element driving force approach can be straightforward extended to reactive systems with multiple products and different reactions routes.

L,nonideal HW

L,ideal HW =

=

(47)

V,nonideal V,ideal R HW = HW + HW

NC  i=1

(46) {FDW }

H 0j V =

It is interesting to observe the terms in Eq. (56) compared with the terms of the driving force obtained through the element Gibbs energy analysis in Eq. (46). Clearly, both Eqs. (56) and (46) must render the same value of FDW , considering the ideality or not of the element binary mixture. The maximum element driving force obtained from Eq. (56) is:

dFDW = dWAL



β

0 AB

V,nonideal L,nonideal d HW − HW



−β

0 AB



dWAL

d WAL HAVap + WBv HBVap dWAL



T. Lopez-Arenas, M. Sales-Cruz and R. Gani et al. / Computers and Chemical Engineering 000 (2019) 106509



−β

0 AB

R E d HW − HW



dWAL

=0

0 = 1/ (H 0V − H 0V ) It should be noted that β 0 Where βAB and B . A AB Hj vap are functions depending only on temperature. Therefore, at the maximum element driving force, the following energy condition must be satisfied:











V,nonideal L,nonideal R E d HW − HW − WAL HAVap + WBv HBVap − HW − HW



dWAL

=0

(58) By considering the definition of the non-ideal element-based enthalpies given by Eqs. (48) and (49), the maximum element driving force condition can be written as: d



=

V,nonideal HW

d



V,ideal HW

R − HW )−(

−

L,ideal HW



L,nonideal HW



− WAL dWAL

Component

(57)

 E



− HW − WAL

HAVap

dWAL

HAVap + WBv HBVap

+ WBv

HBVap





Element A B

n-Butane (1)

Iso-pentane (2)

Iso-butane (3)

1 0 AI

0 1

1 0 AII

And assuming ideal behavior for the vapor and liquid phases, we will determine the element driving force relationship for this reactive system. Thus, the formula matrix A for this reactive system is



A=

1 0

0 1

1 0



(62)

And the partition matrices are:





1 AI = 0



0 ; 1

1 AII = 0

(63)

the matrix product [−AI −1 AII ] is =0

(59)

The above energy condition (Eq. (59)) indicates that: at the maximum element driving force, this is, at the reactive mixture composition (Wl Amax ) corresponding to the maximum element driving force (FDWmax ), the ideal reactive energy condition prevails. Therefore, at the maximum element driving force point the minimum energy to achieve a reactive vapor–liquid separation is required (ideal reactive energy condition). It should be clear that, the minimum reactive energy statement is valid for non-ideal reactive systems, as the non-ideality of the reactive mixture is introduced by considering the reactive residual and excess terms in Eq. (59). 5. The element driving force and the definition of the binary reactive systems As in the case of non-reactive systems, where the driving force approach was defined in terms of two species, in the case of reactive-systems, it is necessary to define the driving force approach in terms of two elements (binary element system) The thermodynamic analysis of the reactive driving force approach will be explained through three binary element systems: an isomerization reaction (i.e., n-butane to iso-butane) and a dimerization reaction (i.e., propylene dimerization) in the presence of an inert species and a binary element system considering the MTBE production. For the first two cases, the phases are considered as ideal and for the third case the SRK Equation of state for the non-ideal vapor phase and UNIFAC model for activity coefficient calculation in the non-ideal liquid phase are used.



−A−1 I AII =

−1 0



(64)

and the working matrices S, Z and n are



S=

1 0 0



0 1 , 0

Z=





−1 0 , 1

nI =

5.1. Binary element system: isomerization of n-butane reactive system The isomerization reaction of n-butane to iso-butane with an inert component can be written as:

n − butane(1 ) ↔ iso − butane (3 )



M 

RT ln YA = RT

(inertelement : B )

(61)

Here, butane (1) is element A, iso-pentane (2) is element B. Therefore, a reduced formula matrix A (consisting of two partitions, AI and AII ) is written as

nII = [n3 ]

(65)

SA,i ln(yi )

= RT ln(y1 )

i=1

λoAvap =

M 

SA,i μvi 0 = μv10

i=1



M 

RT ln PT = RT



SA,i ln(P )

= RT ln(P )

i=1

M 

RT ln XA = RT

SA,i ln(xi )

= RT ln(x1 )

i=1

λoAliq =

M 

l0 SA,i μl0 i = μ1

i=1



M 

RT ln PSA = RT

SA,i ln(PiS )

= RT ln(P1S )

i=1

And for element B M 

RT ln YB = RT

SB,i ln(yi )

= RT ln(y2 )

i=1

λoB vap =

M 

SA,i μvi 0 = μv20



M 

RT ln PT = RT

In terms of elements A and B the above reaction can be written as:

n−A↔i−A

,

i=1

(60)

n1 n2

Using the values of matrix S to obtain the different terms of Eqs. (37) and (38) we obtain:



Inertcomponent : iso − pentane (2 )

7

RT ln XB = RT

SB,i ln(P )

i=1

M 

SB,i ln(xi )

i=1

λoB liq =

M  i=1

= RT ln(P )

l0 SB,i μl0 i = μ2

= RT ln(x2 )

(66)

8

T. Lopez-Arenas, M. Sales-Cruz and R. Gani et al. / Computers and Chemical Engineering 000 (2019) 106509

RT

ln PSB

M 

= RT

SB,i ln(

PiS

) = RT ln(P2S )

(67)

By adding the null term (WAl − WAl ) to the element vapor composition term of Eq. (75)



i=1

and substituting the above equalities into Eq. (46) considering ideal phases we obtain:

RT ln

Y  A

YB



+ RT ln

XB PsB XA PsA

= RT ln



y  1

y2

+ RT ln

x2 P2s x1 P1s

 + RT ln(1 ) = 0 (68)

bvA = nv1 + nv3

(69)

WAv =

bvA nv + nv = v 1 v 3 v = y1 + y3 v bT n1 + n2 + n3

WBv =

bvB nv2 = v = y2 bvT n1 + nv2 + nv3

WBl =

= =

nl1 + nl3 nl1 + nl2 + nl3 nl1 +

nl2 nl2

+ nl3

WBl

= x1 + x3





+ ln

(1 + Kx )P2s WBl P1s

WAl



=0

(77)

= x2

(70)

− μ + μ − RT ln(

(71)



) + RT ln(

x3 P3s

)=0



Ps −μ01l + μ03l x1 = 3s exp x3 P1 RT



=0

(78)

(79)

(80)

FDW + WAl WBl

− FDW

= Kyx

P1s WAl P2s WBl

(81)

Working with Eq. (81) the expression for FDW is obtained as:

KyxWAl WBl P1s − WAl WBl P2s WBl P2s + KyxWAl P1s

(82)

Dividing the numerator and denominator of the above equation over the total pressure P and defining the k1 ideal =P1 s /P and k2 ideal =P2 s /P, we have

KyxWAl (1 − WAl )kideal − WAl (1 − WAl )kideal 1 2

(1 − WAl )kideal + KyxWAl kideal 2 1

(83)

kideal 1

(84)

kideal 2

Then Eq. (83) becomes:

= Kx−1

(73) FDW =

ideal WAl (1 − WAl )αAB − WAl (1 − WAl ) ideal (1 − WAl ) + WAl αAB

(85)

Eq. (85) can be reduced to:

FDW =

ideal WAl αAB

1 + WAl



ideal −1 αAB

 − WAl

(86)

Certainly, Eq. (86) is similar to the equation given by Gani and BekPedersen (2001) for a two-component non-reactive system.

WAl XA = x1 = 1 + Kx YB = y2 = WBv = 1 − WAv

5.2. Binary element system: dimerization of propylene reactive system

XB = x2 = WBl = 1 − WAl

(74)

Introducing Eqs. (74) into (68) we obtain:



+ RT ln



(1 + Ky ) (1 + Kx )

ideal αAB = Kyx

Ky−1

WAv YA = y1 = 1 + Ky

WA 1 (1 + Ky ) WBv

(1 + Kx ) P2s WBl + ln (1 + Ky ) P1s WAl

Also, it is possible to define an ideal reactive relative volatility as:

It should be noted that in this ideal case, the chemical equilibrium constants Ky and Kx are only functions of temperature. By substituting the component mole fractions into the element fractions (Eq. (70)) we obtain:

 v



(1 + Ky ) P1s WAl (1 + Kx ) P2s WBl

=

(72)



=

− FDW

FDW =

By solving Eq. (72) for y1 and x1 we have:

−μ01v + μ03v y1 = exp y3 RT

WBl − FDW

FDW =

− μ01v + μ03v − RT ln(y1 P ) + RT ln(y3 P ) = 0 x1 P1s



it is obtained

Introducing the expressions for the component chemical potentials of species 1 and 3 into Eq. (71) we obtain:

0l 3

FDW + WAl

Kyx =

Z1,1 μl1 + Z1,2 μl2 + Z1,3 μl3 = −μl1 + 0 + μl3 = 0

RT ln

(76)

and defining

Z1,1 μv1 + Z1,2 μv2 + Z1,3 μv3 = −μv1 + 0 + μv3 = 0



FDW + WAl 1 l 1 + K ( y ) WB − FDW

FDW + WAl

Also, the chemical equilibrium condition (Eq. (25)) can be written as:

0l 1

=0

WAl

or

thus the element mole fractions can be written as:

blB blT



ln

blB = nl2

blT

P1s



And using the definition of reactive driving force FDW = WAv − WAl we have



blA = nl1 + nl3

=

(1 + Kx )P2s (1 − WAl )



And reducing the above expression the following equation is obtained

bvB = nv2

WAl



ln

But we know from the formula matrix and Eq. (22) that

blA

WAv − WAl + WAl 1 l l v 1 + K ( y ) 1 − WA + WA − WA

+ RT ln

+ RT ln(Enonideal ) AB



RT ln

(1 + Kx )P2s (1 − WAl ) P1s

WAl

The dimerization reaction of propylene to 2-methyl-2-pentene with an inert component (benzene) is given as:

 =0

(75)

2 propylene(1 ) ↔ 2 − methyl − 2 − pentene(3 ) Inertcomponent : Benzene(2 )

(87)

T. Lopez-Arenas, M. Sales-Cruz and R. Gani et al. / Computers and Chemical Engineering 000 (2019) 106509



In terms of elements A and B the above reaction can be written as:

RT ln XB = RT

(inertelement : B )

2 A ↔ A2

(88)

Here, propylene (1) is element A, benzene (2) is element B. Therefore, a reduced formula matrix A (consisting of two partitions, AI and AII ) is written as Component

Propylene (1 )

Benzene (2 )

1 0 AI

0 1

2 0 AII

1 0

0 1

2 0



(89)





1 AI = 0



0 ; 1



(90)

1 0 0

−2 = 0



(91)



0 1 , 0

Z=





−2 0 , 1

nI =



M 

RT ln YA = RT

n1 , n2

nII = [n3 ]

(92)

SA,i ln(yi )

== RT ln(y1 )

SB,i ln(PiS )

= RT ln(P2S )

(94)

λoAvap =

M  i=1



RT ln PT = RT

RT ln XA = RT

=

SA,i ln(P )

= RT ln(P )

SA,i ln(xi )

SA,i μ = μ

SA,i ln(PiS )



= RT ln(P1S )

SB,i ln(yi )

= RT ln(y2 )



RT ln PT = RT

M  i=1

 + RT ln(1 ) = 0

(95)

(96)

WAv =

bvA nv + 2nv y1 + 2y3 = v 1 v 3 v = v bT n1 + n2 + 2n3 1 + y3

WBv =

bvB nv2 y2 = v = v bT n1 + nv2 + 2nv3 1 + y3

WAl = WBl =

blA blT blB blT

= =

nl1 + 2nl3 nl1 + nl2 + 2nl3 nl2 nl1 + nl2 + 2nl3

=

x1 + 2x3 1 + x3

=

x2 1 + x3

(97)

(98)

(99)

SB,i ln(P )

( x 1 )2 (93)

x3



=P =

−1

−2μ01v + μ03v exp RT

P3s

 s 2 P1





−2μ01l + μ03l exp RT

= Ky−1

 = Kx−1

(100)

(x1 )2 = Kx−1 x3

(101)

It should be noted that in this ideal case the chemical equilibrium constants Ky and Kx are functions of temperature and pressure. By substituting the component mole fractions into the element fractions (Eq. (97)) we obtain:

SA,i μvi 0 = μv20

i=1

x2 P2s x1 P1s

(y1 )2 = Ky−1 y3

i=1

λoB vap =

y2



+ RT ln

blA = nl1 + 2nl3

y3

And for element B

M 

1

But we know from the formula matrix and Eq. (22) that

( y 1 )2

i=1

RT ln YB = RT

y 

+ RT ln(Enonideal ) AB

− 2μ01l + μ03l − 2RT ln(x1 P1s ) + RT ln(x3 P3s ) = 0

= RT ln(x1 )

l0 1



M 

= RT ln



Solving Eq. (99) for y1 and x1 we have:

i=1

RT ln PSA = RT

YB

XB PsB XA PsA

− 2μ01v + μ03v − 2RT ln(y1 P ) + RT ln(y3 P ) = 0

l0 i

M 

+ RT ln

Introducing the expressions for the component chemical potentials of species 1 and 3 considering ideal phases into Eq. (98) we obtain:

i=1 M 

A

Z1,1 μl1 + Z1,2 μl2 + Z1,3 μl3 = −2μl1 + 0 + μl3 = 0

i=1

M 



Y 

Z1,1 μv1 + Z1,2 μv2 + Z1,3 μv3 = −2μv1 + 0 + μv3 = 0

SA,i μvi 0 = μv10 M 

RT ln

Also, the chemical equilibrium condition (Eq. (25)) is given as:

i=1

λ

M 

And the element mole fractions can be written as:

Using the values of matrix S to obtain the different terms of Eqs. (37) and (38) we obtain:

o liq A



blB = nl2

and the working matrices S, Z and n are

S=

l0 SB,i μl0 i = μ2

i=1

bvB = nv2

2 AII = 0

the matrix product [−AI −1 AII ] is



M 

λoB liq =

bvA = nv1 + 2nv3

and the partition matrices are:

−A−1 I AII

= RT ln(x2 )

And substituting the above equalities into Eq. (46) considering ideal phases we obtain:

And assuming ideal behavior for the vapor and liquid phases, we will determine the element driving force relationship for this reactive system: Thus, the formula matrix A for this reactive system is

A=

SB,i ln(xi )

i=1

RT ln PSB = RT

2-methyl-2-pentene (3)

i=1

Element A B



M 

9

= RT ln(P )

WAv + y3WAv = y1 + 2y3

  (y1 )2 Ky WAv − 2 − y1 + WAv = 0

10

T. Lopez-Arenas, M. Sales-Cruz and R. Gani et al. / Computers and Chemical Engineering 000 (2019) 106509

WAl + x3WAl = x1 + 2x3





(x1 )2 Kx WAx − 2 − x1 + WAl = 0

(102)

It is clear from Eq. (102), that in order to obtain the values of y1 and x1 it is necessary to solve the quadratic expressions and, furthermore, from the chemical equilibrium expressions (104) to obtain y3 and x3 . Eqs. (93) and (94) can be written in terms of the element driving force definition (FDW =WA v −WA l ) as:



 l

YA = f1 FDW , WA = y1 =





⎡

⎛

1−











1 − 4Ky FDW + WAl − 2 FDW + WA



⎢

× ⎣1 + Ky ⎝



 l

XA = f3 WA = x1 =

XB =



f4 WAl











1 − 4Ky FDW + WAl − 2 FDW + WAl

 ⎞2

1−

 l



= x2 = 1

⎤





⎛

− WAl

⎢ ⎣1 + Kx ⎝

1−





S=

 ⎞2



1 − 4KxWAl WAl − 2





1 AI = 0



2Kx WAl − 2

⎤

⎠⎥ ⎦ (103)

Introducing Eqs. (103) into (95) we obtain:

⎛

1 AII = 1



−1 = −1

(109)

1 0 0

0 1 , 0



Z=

RT ln

Y  A

YB

= RT ln



+ RT ln

y  1

y2

blA = nl1 + nl3 ;

WAv = (104) WAl =

(105)

In terms of elements A and B the above reaction can be written as:

(106)

Here, Isobutene (1) is element A, Methanol (2) is element B and the third component MTBE (3) is AB. Therefore, a reduced formula matrix A (consisting of two partitions, AI and AII ) is written as

A B

1 0 AI

0 1

1 1 AII

(110)

+ RT ln



x2 P2s x1 P1s

+ RT ln(Enonideal ) AB

 + RT ln(Enonideal )=0 AB

(111)

bvB = nv2 + nv3 blB = nl2 + nl3 (112)

bvA y1 + y3 = ; bvT 1 + y3 blA blT

=

x1 + x3 ; 1 + x3

And considering non-ideal behavior for the vapor and liquid phases, we will find the element driving force relationship for

bv y2 + y3 WBv = vB = bT 1 + y3 WBl =

blB blT

=

x2 + x3 1 + x3

(113)

The mole fractions of component 3, this is MTBE, can be obtained using the chemical equilibrium condition given by Eq. (25) for both phases as:

The MTBE reactive system is defined by the reaction of isobutene (1) and methanol (2) to produce MTBE (3):

MTBE

XB PsB XA PsA



v v v v Zk,i μi = μ1 + μ2 − μ3 = RT ln

i=1

CH3 OH

nII = [n3 ]

and the element mole fractions at both phases are:

5.3. Binary element system: MTBE reactive system

i-C4 H8

n1 , n2

blT = nl1 + nl2 + 2nl3

NC 

Component element

 nI =

bvT = nv1 + nv2 + 2nv3

The element driving force value FDW is obtained by the solution of the non-linear Eq. (104) as a function of WA l at its corresponding equilibrium temperature and pressure.

A + B ↔ AB

−1 −1 , 1

Using the values of matrix S into Eq. (60) we obtain:

bvA = nv1 + nv3 ;



⎜ ⎟ 1− 1−4Ky (FDW +WAl −2 ) (FDW +WAl ) ⎜ ⎟ ⎜ ⎟ 2Ky (FDW +WAl −2 ) ⎡ ⎤ RT ln ⎜   2 ⎟ ⎜ ⎟ ⎜ l )(FDW +WAl ) ⎦ ⎟ ⎝ 1 − FDW − W l ⎣1 + Ky 1− 1−4Ky (FDW +WA −2 ⎠ A 2Ky (FDW +WAl −2 )

Isobutene(1 ) + Methanol(2 )+ ↔ MTBE(3 )

(108)

But we know from the formula matrix and Eq. (22) that

⎞

⎡ ⎛   2 ⎤ ⎞   1− 1−4Kx WAl (WAl −2 ) ⎦ ⎟ ⎜ 1 − WAl ⎣1 + Kx 2Kx (WAl −2 ) ⎜ ⎟ ⎜ P2s ⎟ ⎟=0  + RT ln ⎜ ⎜ P1s ⎟ 1− 1−4Kx WAl (WAl −2 ) ⎜ ⎟ ⎝ ⎠ 2Kx (WAl −2 )

(107)



0 ; 1



1 − 4KxWA WAl − 2

⎡

1 1

and the working matrices S and Z are



2Kx WAl − 2



0 1

And the partition matrices are:

−A−1 I AII

⎠⎥ ⎦

2Ky FDW + WAl − 2



1 A= 0









the matrix product [−AI −1 AII ] is

2Ky FDW + WAl − 2

YB = f2 FDW , WAl = y2 = 1 − FDW − WAl 1−

 l

this reactive system: Thus, the formula matrix A for this reactive system is

+ RT ln NC 

 PP  P

 +

 + RT ln

P1s P2s P3s



 +

ϕˆ 1 ϕˆ 2 y1 y2 ϕˆ 3 y3



μ01v + μ02v − μ03v

l l l l Zk,i μi = μ1 + μ2 − μ3 = RT ln

i=1





RT

γˆ1 γˆ2 x1 x2 γˆ3 x3





μ01v + μ02v − μ03v RT

=0

 =0

(114)

by simplifying Eq. (114) the relationship of the mole fractions for both phases is obtained as:

 0v  μ3 − μ02v − μ01v ϕˆ 3 y3 v = P −1 exp = Kϕ Ky = Keq RT ϕˆ 1 ϕˆ 2 y1 y2  s s −1  0l  μ3 − μ02l − μ01l P1 P2 γˆ3 x3 l = exp = Kγ Kx = Keq P3s RT γˆ1 γˆ2 x1 x2

(115)

T. Lopez-Arenas, M. Sales-Cruz and R. Gani et al. / Computers and Chemical Engineering 000 (2019) 106509

11

where the non-ideal contribution to the equilibrium constants are:

ϕˆ 3 ϕˆ 1 ϕˆ 2 γˆ3 Kγ = γˆ1 γˆ2 Kϕ =

(116)

By substituting the component mole fractions into the element fractions (Eq. (113)) we obtain:









WAv Kϕ Ky y21 + 1 − 2WAv Kϕ Ky + Kϕ Ky y1 − WAv = 0 WAl Kγ Kx x21 + 1 − 2WAl Kγ Kx + Kγ Kx x1 − WAl = 0 y2 =

1 − y1 ; 1 + Kϕ Ky y1

x2 =

1 − x1 1 + Kγ Kx x1

(117)

In order to obtain the values of y1 and x1 it is necessary to solve the non-linear expressions Eq. (117)). It should be pointed out that Kϕ and Kγ are functions of temperature, pressure and component mole fractions. Furthermore, from the chemical equilibrium expressions (Eq. (115)) y3 and x3 are obtained. Eqs. (37) and (38) can be written in terms of the element driving force definition (FDW =WA v −WA l ) as:









FDW + WAl Kϕ Ky y21 + 1 + Kϕ Ky − 2Kϕ Ky FDW + WAl

−



FDW + WAl





=0



y1



WAl Kγ Kx x21 + 1 + Kγ Kx − 2Kγ KxWAl x1 − WAl = 0





YA {FDW } = y1 ;

XA WAl = x1

YB {FDW } = y2 =

1 + YA {FDW } ; 1 + Kϕ KyYA {FDW }



XB WAl



= x2 =



1 + XA WAl





1 + Kγ Kx XA WAl



RT ln(Enonideal ) = RT ln AB

ϕˆ 1 γ2 ϕˆ 2 γ1



 (118)

Fig. 1. Isomerization of n-butane ideal reactive system at P = 1 atm. Inert: isopentane. (a) Reactive driving force diagram, (b) variation of the chemical equilibrium constants Ky and Kx and reactive relative volatility α AB ideal with respect to Wl A .

Introducing Eqs. (118) into (111) we obtain:

 RT ln

YA {FDW } YB {FDW }



    XB WAl P2s  l  s + RT ln(Enonideal + RT ln )=0 AB XA WA P1

(119) The element driving force value FDW for the MTBE reactive system is obtained by the solution of the non-linear Eq. (119) as a function of WA l at its corresponding equilibrium temperature and total pressure. 6. Results and discussion The reactive driving force computation for the three examples presented in the above section using their specific developed equations, is discussed in this section. Also, Appendix A gives the complete thermodynamic and thermo-chemical properties information required for the evaluation of the three reactive systems. Appendix B shows the mathematical expressions for the computation of the reactive element-based enthalpies and Appendix C presents an example of the calculation of the Q matrix needed for the computation of the element-based enthalpies of the vapor and liquid phases. 6.1. Isomerization of n-butane The catalytic isomerization of n-butane to isobutane is an industrially important reaction due to the increasing demand of isobutane in the production of synthetic rubber and in the

alkylation of isoparaffins, which are considered as an alternative for octane boosters instead of oxygenate and aromatic compounds that are subjected to strict environmental restrictions (Wittcoff et al., 2004). The isomerization reaction can be considered as a reversible reaction, and if a catalyst works at high pressure, equilibrium shifts to the reverse side. The knowledge of the equilibrium conditions is important for the development of the appropriate catalyst, the kinetic model and for the efficient operation of the catalytic reactor. Based on the driving force concept, a thermodynamic analysis of the isomerization reaction of n-butane in the presence of an inert component (isopentane) is performed. Considering the equations developed in Section 5.1 for this reactive system, Fig. 1(a) shows the reactive driving force diagram with the maximum reactive driving force FDWmax =0.336 located at Wl Amax =0.355. Eq. (86) was used for the computation of FDW . Fig. 1(b) shows the variation of the chemical equilibrium constants Ky and Kx which are only functions of temperature, with respect to Wl A . Eqs. (73) and (80) were used for the computation of such constants. It can be noted that these equilibrium relationships increase along the whole Wl A range and that Ky is greater than Kx indicating that the composition of isobutane is higher in the vapor phase. The reactive relative volatility value goes from α AB ideal =2.5 to α AB ideal =2.85 and, when Wl A is close to 1, the higher value of the reactive relative volatility indicates that the composition of the light components (n-butane and isobutane) is higher. Fig. 2(b) shows the derivative of the reactive driving force with respect to the element liquid composition indicating the point where the

12

T. Lopez-Arenas, M. Sales-Cruz and R. Gani et al. / Computers and Chemical Engineering 000 (2019) 106509

Fig. 2. Isomerization of n-butane ideal reactive system at P = 1 atm. Inert: isopentane. (a) Reactive driving force diagram, (b) reactive driving force derivative function with respect to Wl A .

derivative value is zero and the location of the maximum reactive driving force Wl A max . Eq. (59) was used to calculate the derivative of FDW. Fig. 3(b) shows the dimensionless reactive element-based enthalpies for the vapor and liquid phases HW Videal /RT, HW Lideal /RT, respectively. Eq. (B.13) given in Appendix B were used to calculate the element-based enthalpies of the vapor and liquid phases. From Fig. 3, it should be pointed out that the connection between the reactive enthalpies and the reactive driving force gives the possibility to determine the optimal thermal feed condition qW as well as the energy required for any reactive distillation process producing isobutane. In Fig. 3(b) it can be noted that, if the element feed composition is set as WAF =[Wl A ]max , the energy required for the reactive separation depends on the thermal condition of the feed. That is, if qW =1 (saturated liquid feed) or qW =0 (saturated vapor feed) two energy consumptions of a reactive distillation column are obtained since the slopes of the element enthalpy tie-lines are different. Fig. 4 shows the connections between the reactive driving force and the equilibrium temperature and the chemicalphysical equilibrium component compositions. Eq. (73) was used for the calculation of the component compositions at both phases. Fig. 4(b) shows the boiling temperature 301 K at Wl A =0 (pure isopentane) and the equilibrium temperature of 262 K for a n-butane and isobutane reactive mixture at Wl A ≈1. Fig. 4(c) shows the component distribution in the liquid and vapor phases. It can be noted in Fig. 4(c) that at the maximum driving force composition, the largest vapor–liquid component composition difference for isobutane and isopentane is obtained. Also, it should be pointed out that at Wl A ≈1 an isobutane rich mixture (xisobutane = 0.8668) is reached, considering that Wl A = xn-butane + xisobutane .

Fig. 3. Isomerization of n-butane ideal reactive system at P = 1 atm. Inert: isopentane. (a) Reactive driving force diagram, (b) dimensionless reactive element-based enthalpy graph.

6.2. Dimerization of propylene 2-Methyl-1-pentene is a dimer of propylene produced with the aid of a tri-n-propyl aluminum catalyst. This reaction is the basis for isoprene synthesis. This process was used by Goodyear and it involves the dimerization of propylene to 2-methyl-1-pentene. The reaction goes in 90% yield at about 300 °C and elevated pressures with a tri-n-propyl-aluminum catalyst. This dimer is isomerized at about 100 °C with a silica–alumina catalyst to the more stable 2-methyl-2-pentene, which in turn is demethanated to isoprene and methane at 660 °C in the presence of superheated steam and catalytic amounts of hydrogen bromide (Wittcoff et al., 2004). The reversible dimerization of propylene step is highly important since a complete dimerization-isomerization would render higher production of isoprene. Thus, the appropriate knowledge of the phase and chemical equilibrium conditions is a key issue. Considering the equations developed in Section 5.2 for this reactive system, Fig. 5(a) shows the reactive driving force diagram at P = 5 atm with benzene as inert species, with the maximum reactive driving force FDWmax =0.0794 located at Wl Amax =0.465. The non-linear Eq. (104) was solved by a Newton–Raphson method to obtain FDW . It is important to note the low values of the reactive driving force for this reactive system 0
T. Lopez-Arenas, M. Sales-Cruz and R. Gani et al. / Computers and Chemical Engineering 000 (2019) 106509

13

Fig. 5. Ideal dimerization of propylene reactive system at P = 5 atm. Inert: benzene. (a) Reactive driving force diagram, (b) variation of the chemical equilibrium constants Ky and Kx with respect to Wl A .

Fig. 4. Isomerization of n-butane ideal reactive system at P = 1 atm. Inert: isopentane. (a) Reactive driving force diagram, (b) equlibrium temperature as a function of Wl A , (c) component distribution as a function of Wl A .

range and its high values indicate that the chemical reaction is an irreversible reaction in both phases. In fact, the high value of Kx ≈ 2 × 107 means that the conversion to 2-methyl-2-pentene is far from a reversible reaction in the liquid phase. On the other hand, it should be noted that while Kx is only a function of temperature, Ky is function of both, temperature and pressure. Fig. 6(b) shows the derivative of the reactive driving force with respect to the element liquid composition indicating where the function cross the zero line and the location of the maximum reactive driving force Wl A max . Eq. (59) was used to calculate the derivative of FDW. Fig. 7(b) shows the dimensionless element reactive enthalpies for the vapor and liquid phases HW Videal /RT, HW Lideal /RT, respectively. Eq. (B.13) given in Appendix B were used to calculate the element-based enthalpies of the vapor and liquid phases. From

Fig. 7, it should be pointed out that the connection between the reactive enthalpies and the reactive driving force gives the possibility to determine the optimal thermal feed condition qW as well as the energy required for any reactive distillation producing 2methyl-2-pentene. In Fig. 7(b) it can be noted that, if the element feed composition is set as WAF =[Wl A ]max , the energy required for the reactive separation depends on the thermal condition of the feed. This is, if qW =1 (saturated liquid feed) or qW =0 (saturated vapor feed) two energy consumptions of a reactive distillation column is obtained since the slope of the element enthalpy tie-lines are different. Fig. 8(a) and (b) shows the connection between the reactive driving force and the equilibrium temperature at two different total pressures (P = 5, 8 atm). The non-linear Eq. (104) was solved by a Newton–Raphson method to obtain FDW at the two different pressures. Fig. 8(b) shows the boiling temperature of 416.6 K and 440.6 K at Wl A =0 (pure benzene) at P = 5 and P = 8 atm, respectively. It can be noted in Fig. 8(a) that by increasing the total pressure (P = 8 atm) the reactive driving force is decreased and in Fig. 8(b) it is observed an increment of the equilibrium temperature. Fig. 8(c) shows the liquid component composition as a function of Wl A at P = 5 atm. It should be pointed out that as Wl A is increased (increasing propylene composition) there is an instantaneous transformation to dimer 2-methyl-2-pentene, thus the liquid composition of propylene is close to zero, and a reactive mixture with very small amounts of propylene, benzene and the dimer 2-methyl-2-pentene is obtained. Also, it should be noted

14

T. Lopez-Arenas, M. Sales-Cruz and R. Gani et al. / Computers and Chemical Engineering 000 (2019) 106509

Fig. 6. Ideal dimerization of propylene reactive system at P = 5 atm. Inert: benzene. (a) Reactive driving force diagram, (b) reactive driving force derivative function with respect to Wl A .

that at Wl A =0.99 a 2-methyl-2-pentene rich mixture (xpropylene =0.0 0 016, x2-methyl-2-pentene = 0.9820, xbenzene =0.01783) is reached. 6.3. MTBE production Methyl tert-butyl ether (MTBE) is an important chemical used as an octane booster in gasoline to replace tetra ethyl lead. Maximum production of the MTBE can be achieved using reactive distillation (RD) process that is operated at the optimum operating conditions and column configuration. However, optimizing the column configuration such as tray or catalyst location is experimentally expensive. Therefore, a reliable model of the MTBE non-ideal chemical-physical equilibrium is important to find the optimum conditions for MTBE production. In this case, the non-ideality of the phases is considered using the SRK Equation of State for the vapor phase and UNIFAC model for the liquid phase at P = 10 atm. Considering the equations developed in Section 5.3 for this reactive system, Fig. 9(a) shows the reactive driving force diagram at P = 10 atm. The non-linear system of Eqs. (118) and (119) was solved by a Newton–Raphson method to obtain FDW . It can be observed in Fig. 9(a) that this reactive system presents two maximum reactive driving force points: FDWmax =0.2115 located at Wl Amax =0.155 and FDWmax =0.1805 located at Wl Amax =0.663. Also, it is important to note that a minimum reactive driving force is found for this reactive system FDWmin =0.117 located at Wl Amin =0.467 indicative of, by changing the total pressure, a reactive azeotrope node could exist. Fig. 9(b) shows the variation of the chemical equilibrium constants Kv eq and Kl eq which, in

Fig. 7. Ideal dimerization of propylene reactive system at P = 5 atm. Inert: benzene. (a) Reactive driving force diagram, (b) dimensionless reactive element-based enthalpy graph.

this case are functions of temperature, pressure and composition, with respect to Wl A . Eqs. (115) and (116) were used to calculate the chemical equilibrium constants. It can be noted that these equilibrium relationships increase along the whole Wl A range, but these grow exponentially after the section where the second maximum driving force is (Wl Amax =0.663). The above means that there is a higher composition of MTBE in the vapor and liquid phases compared to the reactant methanol, despite isobutene composition increases as Wl A is augmented. Fig. 9(c) shows the effect of the non-idealities on the reactive driving force. It can be noted that while Kϕ value is constant around unity indicating ideal reactive vapor phase, Kγ increases to a value of Kγ ≈ 13, signifying that the effect of the non-ideality of the liquid phase is very important in the determination of the conversion to MTBE. Also, it can be noted that the non-ideal term AB nonideal considered in the reactive driving force equation (see Eq. (119)) increased exponentially to a value of AB nonideal =17, demonstrating that while the fugacity coefficients of isobutene and methanol are unity, the activity coefficient of methanol is 17 times larger than the activity coefficient of isobutene, thus, the non-ideality of methanol in the liquid phase is a key variable of the reactive system. Fig. 10(b) shows the derivative of the reactive driving force with respect to the element liquid composition indicating the different points where the function cross the zero line and the location of the maximum reactive driving forces Wl A max and the minimum reactive driving force Wl Amin . Eq. (59) was used to calculate the derivative of FDW. Fig. 11(b) shows the dimensionless reactive

T. Lopez-Arenas, M. Sales-Cruz and R. Gani et al. / Computers and Chemical Engineering 000 (2019) 106509

Fig. 8. Ideal dimerization of propylene reactive system at P = 5 and 8 atm. Inert: benzene. (a) Reactive driving force diagrams, (b) equlibrium temperatures (P = 5, 8 atm) as a function of Wl A , (c) component distribution at P = 5 atm as a function of Wl A .

element-based enthalpies for the vapor and liquid phases HW V /RT, HW L /RT, respectively. It should be noted that in this case we have non-ideal element-based reactive enthalpies. Eq. (B.13) given in Appendix B were used to calculate the element-based enthalpies of the non-ideal vapor and non-ideal liquid phases. From Fig. 11, it should be pointed out that the connection between the reactive enthalpies and the reactive driving force gives the possibility to determine the optimal thermal feed condition qW as well as the energy required for any reactive distillation producing MTBE. In Fig. 11(b) it can be noted that, if the element feed composition is set as WAF =[Wl A ]max , in this case two feeds, the energy required for the reactive separation depends on the thermal condition of the feeds. This is, if qW =1 (saturated liquid feed) or qW =0 (saturated vapor feed) two energy consumptions of a reactive section

15

Fig. 9. Non-ideal MTBE reactive system at P = 10 atm. Vapor phase SRK EoS. Liquid Phase UNIFAC model. (a) Reactive driving force diagram, (b) variation of the chemical equilibrium constants Kv eq and Kl eq with respect to Wl A , (c) effect of the non-ideality of the phases on the equilibrium constants and reactive driving force.

of a distillation column is obtained since the slopes of the reactive element-based enthalpy tie-lines are different. It should be pointed out that, at the location of the minimum reactive driving force Wl Amin , a non-reactive section should be considered as an intermediate part of the reactive distillation column in order to break the reactive azeotrope node condition (Wl Amin =0). Fig. 12(a) and (b) shows the connection between the reactive driving force and the equilibrium temperature at P = 10 atm. Fig. 12(b) shows the boiling temperature of 406 K at Wl A =0 (pure methanol) and 337 K at Wl A =1 (pure isobutene) at P = 10, respectively. Also, it can be noted in Fig. 12(b) that temperature increases at the location where the maximum liquid composition of MTBE is found. Fig. 12(c) shows the liquid component composition as a function of Wl A at P = 10 atm. The non-linear system of Eq. (117) was solved by a

16

T. Lopez-Arenas, M. Sales-Cruz and R. Gani et al. / Computers and Chemical Engineering 000 (2019) 106509

Fig. 10. Non-ideal MTBE reactive system at P = 10 atm. Vapor phase SRK EoS. Liquid phase UNIFAC model. (a) Reactive driving force diagram, (b) reactive driving force derivative function with respect to Wl A .

Newton–Raphson method to obtain the component compositions. It should be pointed out that as Wl A is increased from zero to Wl A =0.24 (increasing isobutene liquid composition) there is a transformation to MTBE preferentially in the vapor phase. However, after this point, the transformation to MTBE occurs preferentially in the liquid phase, reaching a maximum MTBE composition in the liquid phase (xMTBE =0.652) at Wl A =0.531. Finally, it should be pointed out that the connection between the reactive driving force diagram and the component composition graph, clearly leads to the design of a reactive distillation column considering two reactive sections where the maximum driving forces are located and a non-reactive section where the minimum reactive driving force appears. 7. Conclusions A fundamental thermodynamic analysis of the driving force based approach for design of reactive systems has been performed. The thermodynamic relationships between the reactive driving force definition and the energy involved in the chemical and phase equilibrium are demonstrated for ideal and non-ideal binary element mixtures. Through the thermodynamic analysis, the fundamental relationships between the element reactive driving force and the energy involved in the reactive separation problem are revealed. In this work, the explicit relationship of the element reactive driving force with the heat of vaporization and the reaction enthalpies of the reacting mixture was shown. It is clearly shown that the reactive driving force approach may be considered as a useful thermodynamic tool to identify the element or species mixture composition nodes of minimum equilibrium energy. The graphical visualization of the element driving force

Fig. 11. Non-ideal MTBE reactive system at P = 10 atm. Vapor phase SRK EoS. Liquid phase UNIFAC model. (a) Reactive driving force diagram, (b) dimensionless reactive element-based enthalpy graph.

and the element-based enthalpy diagrams are highly valuable for the design and analysis of reactive distillation columns (LópezArenas at al., 2019). Such diagrams indicate the location where the minimum energy is required for a specific reactive separation and, by using the associated temperature and component composition diagrams, it the corresponding temperature and phase component composition distribution can easily determined. Furthermore, with the element-based enthalpy diagram, it can be identified the optimal thermal feed condition (qW ) in order to reduce the energy consumption. This key information is highly valuable since it reduces the heating and cooling in any reactive distillation process. Through the validation examples it was found that the energy associated to the non-idealities and the thermo chemical properties of the different species can severely affect the energy required for the reactive separation of the mixture. Also, with the graphical visualization of the element driving force approach, it properly reveals the reactive azeotrope node location and the energy associated at this equilibrium point. It should be pointed out that the thermodynamic analysis of the element driving force approach in the present work has been performed for a binary element vapor–liquid, but this could be extended to multi-element, multi-reaction systems. Finally, the element driving force approach may be directly applied to systems considering reactive liquid– liquid equilibrium and reactive liquid–solid equilibrium. This is, the concept can be straightforward extended to consider an element driving force approach for reactive liquid–liquid extraction and reactive crystallization processes.

T. Lopez-Arenas, M. Sales-Cruz and R. Gani et al. / Computers and Chemical Engineering 000 (2019) 106509

17

Formula matrix:

1 0

0 1

1 0

Critical properties: i

Tc (K)

Pc (Bar)

w

1 2 3

425.20 460.40 408.20

3.8000D01 3.3900D01 3.6500D01

0.199 0.227 0.183

Antoine equation constants (ln Ps i =A + B/(T + C): Pa, K)): i

A

B

C

1 2 3

20.56850 20.51453 20.73320

−2154.70 −2348.67 −2032.79

−34.411 −40.051 −33.150

Ideal gas heat capacity equation coefficients (J/mol) Cp=ai + bi T+ ci T2 + di T3 (K) i

ai

bi

ci

di

1 2 3

9.487d00 −9.525d00 −1.39d00

3.313d−01 5.066d−01 3.847d−01

−1.108d−04 −2.729d−04 −1.846d−04

−2.822d−09 5.723d−08 2.895d−08

Ideal gas heat of formation and Gibbs free energy of formation (J/mol): i

( Hfi )ig

( Gfi )ig

1 2 3

−1.262d05 −1.546d05 −1.346d05

−1.610d04 −1.482d04 −2.090d04

A.2. Dimerization of propylene NC=3; M = 2 Components: PROPYLENE (1) BENZENE (2) 2-METHYL-2-PENTENE (3) Chemical reaction:

2 Propylene(2A ) <==> 2 − methyl − 2 − pentene(A2 ) Inert : Benzene(B ) Fig. 12. Non-ideal MTBE reactive system at P = 10 atm. Vapor phase SRK EoS. Liquid phase UNIFAC model. (a) Reactive driving force diagram, (b) equlibrium temperature as a function of Wl A , (c) component distribution as a function of Wl A .

Formula matrix:

1 0

0 1

2 0

Critical properties: Appendix A. Thermodynamic and thermo chemical properties A.1. Isomerization of n-butane NC=3; M = 2 Components: N-BUTANE (1) ISO-PENTANE (2) ISO-BUTANE (3) Chemical reaction:

n − Butane(A ) <==> Iso − butane(A ) Inert : Iso − pentane(B )

i

Tc (K)

Pc (Bar)

w

1 2 3

364.90 562.20 518.00

4.6000D01 4.8900D01 3.2800D01

0.144 0.212 0.229

Antoine equation constants (ln Ps i =A + B/(T + C): Pa, K)): i

A

B

C

1 2 3

20.59550 20.68970 20.83500

−1807.53 −2726.81 −2725.89

−26.150 −55.628 −47.640

Ideal gas heat capacity equation coefficients (J/mol) Cp=ai + bi T+ ci T2 + di T3 (K)

18

T. Lopez-Arenas, M. Sales-Cruz and R. Gani et al. / Computers and Chemical Engineering 000 (2019) 106509

i

ai

bi

ci

di

j,i

1

2

3_

1 2 3

3.710d00 −3.392d00 −1.475d01

2.345d−01 4.739d−01 5.669d−01

−1.160d−04 −3.017d−04 −3.341d−04

2.205d−08 7.130d−08 7.963d−08

1 2 3 4 5

2 1 0 0 0

0 0 1 0 0

3 0 0 1 1

Ideal gas heat of formation and Gibbs free energy of formation (J/mol): i

( Hfi )ig

( Gfi )ig

1 2 3

2.043d04 8.298d04 −6.653d04

6.276d04 1.297d05 7.126d04

UNIFAC binary interaction parameters Aij (cal/mol):

A.3. MTBE production

1

2

3

4

5

1 2 3 4 5

0.0 −35.36 16.51 0.0 83.36

86.02 0.0 −12.52 86.02 26.51

697.20 787.60 0.0 697.2 238.4

0.0 −35.36 16.51 0.0 83.36

251.5 214.5 −128.6 251.5 0.0

NC=3; M = 2 Components:

Appendix B. Reactive enthalpy calculation

ISOBUTENE (1) METHANOL (2) MTBE (3)

Based on the fundamental multi-component thermodynamics (Thompson, 20 0 0), similar equations can be used to calculate the element based reactive enthalpies. The element enthalpy for any phase can be obtained from the element chemical potential as:

∂ H j = −RT ∂T

Chemical reaction:



λj

(B.1)

RT

and the binary element vapor and liquid enthalpies are:

Formula matrix:

0 1



2

ISOBUTENE(1 ) + METHANOL(2 ) <==> MTBE(3 )

1 0

i,j

1 1

HAv = −RT 2

Critical properties:

HAl = −RT 2

i

Tc (K)

Pc (Bar)

w

1 2 3

417.90 512.58 497.10

3.9990D01 8.0959D01 3.4300D01

0.1893 0.5688 0.2700

Antoine equation constants (ln

Ps

i =A + B/(T + C):

i

A

B

C

1 2 3

20.64556 23.49989 20.71616

−2125.74886 −3643.31362 −2571.58460

−33.160 −33.434 −48.406

Pa, K)):

ai

bi

ci

di

1 2 3

1.605d01 2.115d01 2.534d00

2.804d−01 7.092d−02 5.136d−01

−1.091d−04 2.587d−05 −2.569d−04

9.098d−09 −2.852d−08 4.303d−08

Ideal gas heat of formation and Gibbs free energy of formation (J/mol): ( Gfi )ig

1 2 3

−1.6903d04 −2.0108d05 −2.9288d05

5.8073d04 −1.6238d05 −1.1732d05

UNIFAC model parameters: Number of groups=5 j

R

Q

1 2 3 4 5

0.9011 1.1173 1.4311 0.2195 1.1450

0.8480 0.9880 1.4320 0.0000 1.0880



λvA RT

λ

l A



HBv = −RT 2

;

 HBl = −RT 2

;

RT

v = W vHv + W vHv; HW B B A A

i

( Hfi )ig

∂ ∂T



∂ ∂T ∂ ∂T





λvB RT

λlB



 (B.2)

RT

For the binary mixture we have.

Ideal gas heat capacity equation coefficients (J/mol) Cp=ai + bi T+ ci T2 + di T3 (K)

i

∂ ∂T

l HW = WAl HAl + WBl HBl

(B.3)

It is know from the chemical model approach (see Eq. (31)) that the relationship between the element chemical potentials λj and the species chemical potentials μi for the vapor and liquid phases respectively, is given as:

∂ ∂T ∂ ∂T



λvj



=−

RT



λlj

i=1

 =−

RT

NC 

NC 

∂ j,i ∂T



Qv

Q lj,i

i=1

∂ ∂T

μvi



;

RT



μli

 (B.4)

RT

and for a ϕ −γ approach we have:

μvi RT

=

μ0i v RT





+ ln P yi ϕˆ i ;

μli RT

=

μ0i l RT



+ ln Pis xi γi



(B.5)

Performing the derivatives of Eq. (B.4) with respect to the temperature we obtain:

  ∂ μvi RT ∂ T RT i=1

   NC   ∂  ∂ μ0i v =− Q vj,i + ln ϕˆ i ∂ T RT ∂T i=1  l  l NC  ∂ μi ∂ λj =− Q lj,i ∂ T RT ∂ T RT i=1

 0l   NC  ∂  s ∂ ∂ μi =− Q lj,i + ln Pi + (ln γi ) ∂ T RT ∂T ∂T i=1 ∂ ∂T



λvj



=−

NC 

Q vj,i

(B.6)

T. Lopez-Arenas, M. Sales-Cruz and R. Gani et al. / Computers and Chemical Engineering 000 (2019) 106509

the matrix product [−AI −1 AII ] is

Also, we know that:

μ0i v RT

∂ ∂T

=





Gigf i

RT

;



μ0i v

=

RT

∂ ∂T





−A−1 I AII =

 ig

G f i RT

=−

Hi0

(B.7)

RT

C pig = ai + bi T + ci T 2 + di T 3 i

(B.8)

then

(B.9)

(B.10)

Introducing the above equations into Eq. (B.6), the ideal contributions to the element-based enthalpies are:

HBvideal = − HAlideal = −





=−

NC 

l QB,i





Hi0

vap

− Hi

hRB = −RT 2

 QB,i

i=1 E HAB = WAl hEA + WBl hEB ;

hEB = −RT 2

NC 

 QB,i

i=1

hRA = −RT 2

∂ ln(ϕˆ i ) ∂T



hEA = −RT 2

∂ ln(γi ) ∂T



NC 



QA,i

i=1

NC  i=1

 QA,i



∂ ln(ϕˆ i ) ; ∂T

l E HW = WAl HAlideal + WBl HBlideal + HAB

1 0

1 AI = 0



0 ; 1

−1

(C.4)

ZT

(C.6) (C.7)

∂ μi ∂ nj



(C.8) T,P

(B.13)

−1 1 −1 y2 −1

−1 −1 ⎦; 1 −1 y3

−1 −1 −1

−1 1 −1 x2 −1

−1 −1 ⎦ 1 −1 x3

1 nlT

⎣

x1

⎤

(C.9)

where nv T and nl T are the total moles at the vapor and liquid phases. Working with the above matrices the following relationships can be obtained:



Z T Uv Z T

l

Z UZ

−1 −1



=

=

1 nvT 1

 y + y −1 1 3 y1 y3

 x + x −1 1 3

nlT



−1



−1

x1 x3

−1 = Y13 ;

−1 = X13

(C.10)





⎡

y3 y1 +y3

Ev Uv S = ⎣ 0

1 0 −1

0 0 0

−1 0 1

⎤

⎡

0 0⎦; 0

The QT matrices are:

Q

x3 x1 +x3

El Ul S = ⎣ 0

−x3 x1 +x3

⎡

[Q v ] = S − E v U v S = ⎣

 l T (C.2)

−1 0 ; 1



1−

T

(C.1)

0 0 0

−1 Z T = Y13

−1 Z T = X13

1 0 −1

(C.11)

With the above equations it is possible to obtain the matrices Ev Uv S and El Ul S as:

−y3 y1 +y3



1 AII = 0

⎤

1 −1 1 y1 v ⎣ −1 U = v nT −1

E l = Z ZT U l Z

and the partition matrices are:



Ui j =

(B.12)

Considering the isomerization of n-butane, we have the following formula matrix:

0 1

nII = [n3 ]

(C.5)

E v = Z ZT U vZ

Appendix C. Calculation of the Q matrix: ideal example

1 0

n1 , n2

and the matrices Ev and El are:

 ∂ ln(γi ) ; ∂T

v = W v H videal + W v H videal + H R HW B B A A AB

A=





The final expressions for the calculation of the element-based reactive enthalpies for the reactive vapor-liquid mixture are:



 nI =

where the elements of matrix U are:

(B.11)

i=1

NC 



E = Z Z TU Z

Ul =



And the non-ideal contributions to the element-based enthalpies in terms of the fugacity and activity coefficients of the different species are given as: R HAB = WAv hRA + WBv hRB ;

−1 Z 0 , 1

The QT matrix is given as:

⎡1

l QA,i Hi0 − Hivap ;

i=1

HBlideal



0 1 , 0

⎡

v H 0 QB,i i

i=1 NC 

1 0 0

By performing the respective derivations of Eq. (C.8) for the NC=3 components and M = 2 elements, the following matrices for the ideal vapor and liquid phases are obtained:

v H 0 ; QA,i i

i=1 NC 

(C.3)

H = −QU S

∂  s  Hivap ln Pi = ∂T RT NC 



with

Also

HAvideal = −



Q T = S − EU S

 bi  2 = + ai (T − T0 ) + T − T02 RT 2  di  4  ci  3 3 + T − T0 + T − T04 3 4 H igf i

RT

−1 0

and the working matrices S, Z and n are

S=

if the heat capacity of species i is given as:

Hi0

19

⎡

1

= S − E lU l S = ⎣

y3 y1 +y3

0 y3 y1 +y3 − x1x+3x3

0 x3 x1 +x3

⎤

0 0⎦ 0

(C.12)

⎤

0 1⎦; 0

⎤

0 1⎦ 0

(C.13)

20

T. Lopez-Arenas, M. Sales-Cruz and R. Gani et al. / Computers and Chemical Engineering 000 (2019) 106509

Finally the Q matrices obtained are:

Qv =

1−

y3 y1 +y3

0

0 1

y3 y1 +y3

0



;

Ql =

1−

x3 x1 +x3

0

0 1

x3 x1 +x3



0 (C.14)

References Agrawal, R., Herron, D.M., 1997. =ptimal thermodynamic feed conditions for distillation of ideal binary mixtures. AIChE J. 43, 2984. Agreda, V.H., Partin, L.R., Heise, W.H., 1990. High-purity methyl acetate via reactive distillation. Chem. Eng. Prog. 86, 40. Atkinson, T.D., 1987. Second law analysis of cryogenic processes using a three term model. Gas Sep. Purif. 1 (2), 84. Bandyopadhyay, S., 2002. Effect of feed on optimal thermodynamic performance of a distillation column. Chem. Eng. J. 88 (1–3), 175. Blahusiak, M., Kiss, A.A., Kersten, S.R.A., Schuur, B., 2016. Quick assessment of binary distillation efficiency using a heat engine perspective. Energy 116, 20–31. Ciric, A.R., Gu, D., 1994. Synthesis of nonequilibrium reactive distillation by MINLP optimization. AIChE J. 40 (9), 1479–1487. Doherty, M.F., Buzad, G., 1992. Reactive distillation by design. Trans. Inst. Chem. Eng. 70A, 448. Fletcher, R., 1981. Practical Methods of Optimization. Volume 2: Constrained Optimization. John Wiley, New York, USA. Gani, R., Bek-Pedersen, E, 2001. Simple new algorithm for distillation column design. AIChE J. 46 (6), 1271–1274. doi:10.1002/aic.690460619. Gani, R., Jepsen, T.S., Perez-Cisneros, E.S., 1998. A generalized reactive separation unit model. Model. Simul. Asp. Comput. Chem. Eng. 22, S363–S370. Lopez-Arenas, T., Mansouri, S.S., Sales-Cruz, M., Gani, R., Pérez-Cisneros, E.S., 2019. A Gibbs energy-driving force method for the optimal design of non-reactive and reactive distillation columns. Comput. Chem. Eng. 128, 53–68 2019.

Mansouri, S.S., Sales-Cruz, M., Huusom, J.K., Woodley, J.M., Gani, R., 2015. Integrated process design and control of reactive distillation processes. IFAC-PapersOnLine 48 (8), 1120–1125. doi:10.1016/j.ifacol.2015.09.118. Michelsen, M.L., 1995. Basic concepts in the calculation of chemical equilibrium. In: Proceedings of the Nordic Conference on Reactive Separation Systems. Lyngby, Denmark. Technical University of Denmark April 28. Michelsen, M.L., Mollerup, J., 2004. Thermodynamic Models: Fundamentals and Computational Aspects, second ed. Tie-Line Publications ISBN: 87-989961-3-4. Pérez-Cisneros, E.S., 1997. Modelling, Design and Analysis of Reactive Separation Processes Ph.D. Thesis. Technical University of Denmark, Kongens Lyngby. Pérez-Cisneros, E.S., Gani, R., Michelsen, M.L., 1997. Reactive separation systems—I. Computation of physical and chemical equilibrium. Chem. Eng. Sci. 52 (4), 527– 543. doi:10.1016/S0 0 09-2509(96)0 0424-1. Pérez-Cisneros, E.S, Sales-Cruz, M., 2017. Thermodynamic analysis of the driving force approach: non-reactive systems. Comput. Chem. Eng. doi:10.1016/j. compchemeng.2017.06.021. Ratkje, S.K., Sauar, E., Hansen, E.M., Lien, K.M., Hafskjold, B., 1995. Analysis of entropy production rates for design of distillation columns. Ind. Eng. Chem. Res. 34 (9), 3001. Sánchez-Daza, O., Pérez-Cisneros, E.S., Bek-Pedersen, E., Gani, R., 2003. Graphical and stage-to-stage methods for reactive distillation column design. AIChE J. 49 (11), 2822–2841. doi:10.1002/aic.690491115. Smith, L.A., Huddleston, M.N., 1982. New MTBE design now commercial. Hydrocarb. Process. 3, 121. Taprap, R., Ishida, M., 1996. Graphic exergy analysis of processes in distillation column by energy-utilization diagram. AIChE J. 42 (6), 1633. Tester, J.W., Modell, M., 1997. Thermodynamics and its Applications. Prentice Hall. Thompson, E.V., 20 0 0. A Unified Introduction to Chemical Engineering Thermodynamics: The Laws of Thermodynamics, Material and Energy Balances, Chemical Thermodynamics. Stillwater Press, Orono, Maine 20 0 0. Wittcoff, H.A., Reuben, B.G., Plotkin, J.S., 2004. Industrial Organic Chemicals. John Wiley & Sons, Inc., Hoboken, New Jersey.