- Email: [email protected]
Differential equations for the calculation of isopleths of multicomponent fluid mixtures
Accepted Manuscript Differential equations for the calculation of isopleths of multicomponent fluid mixtures Ulrich K. Deiters PII:
S0378-3812(17)30122-X
DOI:
10.1016/j.fluid.2017.03.022
Reference:
FLUID 11439
To appear in:
Fluid Phase Equilibria
Received Date: 15 November 2016 Revised Date:
25 January 2017
Accepted Date: 27 March 2017
Please cite this article as: U.K. Deiters, Differential equations for the calculation of isopleths of multicomponent fluid mixtures, Fluid Phase Equilibria (2017), doi: 10.1016/j.fluid.2017.03.022. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
x’ = const
AC C
EP
TE D
M AN U
SC
RI PT
p
T
RI PT
ACCEPTED MANUSCRIPT
Ulrich K. Deiters
SC
Differential equations for the calculation of isopleths of multicomponent fluid mixtures
AC C
EP
TE
D
M AN U
¨ Germany Institute of Physical Chemistry, University of Cologne, Luxemburger Str. 116, D-50939 Koln,
ACCEPTED MANUSCRIPT
Abstract
RI PT
An isoplethic phase envelope is a locus of two-phase equilibria where one of the phases is kept at constant composition. The thermodynamic conditions for such a phase envelope are expressed as a set of first-order ordinary differential equations in which molar concentrations appear as central variables instead of mole fractions. The resulting formalism is applicable to multicomponent mixtures, well-suited for machine calculations, and allows a very rapid calculation of phase envelopes with arbitrary equations of state.
SC
Application to the formalism to pure fluids yields simple differential equations that are analoga of Clapeyron’s equation, but permit the calculation of orthobaric densities.
1
M AN U
Keywords: fluid phase equilibria, Gibbs–Konovalov equations, algorithm, multicomponent mixtures, isochoric thermodynamics
Introduction
Thermodynamics is built on differential calculus, and therefore it is not surprising that many problems can be expressed by means of algebraic equations or differential equations. Wellknown examples are
TE
D
• the calculation of vapor-pressure curves, which can be accomplished by evaluating Maxwell’s equal-area criterion (algebraic equation) or by integrating the Clausius–Clapeyron equation (differential equation),
EP
• or the calculation of adiabatic expansion curves for ideal gases, for which probably all thermodynamics textbooks of the world present a differential equation, but which can also be accomplished by setting the entropy to a constant value and solving the resulting algebraic equation.
AC C
Likewise, phase equilibria of fluid mixtures can be calculated from the conditions that the temperatures, the pressures, and the chemical potentials of the components of coexisting phases must agree, which leads to a set of nonlinear algebraic equations. Differential equations were already developed at the end of the 19th Century by Konovalov [1,2] and Gibbs [3] and became known as Gibbs–Konovalov rules; a particularly clear and concise derivation of these rules for multicomponent mixtures can be found in the textbook of Haase [4]. The Gibbs–Konovalov rules can be used to prove some thermodynamic theorems, e.g., that, in isothermal phase diagrams of binary mixtures, azeotropes must lie at pressure extrema. Because of their inconvenient mathematical structure—an example is discussed in the next section—they never became popular for computing phase envelopes. But this is true for the classical formulation in terms of mole fractions and partial molar derivatives only. It could recently been shown that, in a density-based formulation of thermodynamics, there exist analoga to the Gibbs–Konovalov rules that can be used to develop very efficient and robust algorithms for the computation of isothermal as well as isobaric phase envelopes [5].
1
ACCEPTED MANUSCRIPT
RI PT
There is another kind of phase envelope with a paramount importance for chemical engineering, namely the isoplethic phase envelope. It represents the boundary between the singlephase and the two-phase regions for a mixture with a given composition (i.e., set of mole fractions) and comprises of states at which boiling or condensation sets in. The so-called synthetic method for experimentally determining phase equilibria, according to which samples of known overall composition a manipulated until phase transitions are detected, yields isopleths as primary results. Recently the synthetic method has even been applied to computer simulations of fluid mixtures [6]. We shall therefore extend the previous work and derive differential equations for isoplethic phase envelopes. Particular attention will be given to the behavior of these equations at critical points of mixtures.
2
Differential equations for isopleths
M AN U
SC
The isopleth of a pure fluid is its vapor pressure curve. We shall derive differential equations for the orthobaric densities—analoga of Clapeyron’s equation in the density-based formalism.
D
Fig. 1 serves to illustrate the problem. It shows the phase diagrams of the (methane + propane) system, calculated with the Peng–Robinson equation of state with standard quadratic mixing rules [7, 8]. This is a simple system with an uninterrupted critical curve connecting the critical points of the pure fluid. The diagram furthermore contains the x1 = 0.8 isopleth—a curve marking the onset of boiling of a liquid having this composition, or the onset of condensation of a vapor having this composition. The phase in equilibrium with the boiling liquid, or the condensing vapor, respectively, has a different composition, but the same pressure and temperature, and therefore cannot be shown in this diagram.
TE
The isopleth meets the critical curve of the mixture in a binary critical point. It should be noted that this critical point usually does not lie at a maximum of the isopleth. The task is to calculate the isopleth efficiently, producing not only its pT contour, but also the composition of the coexisting phase.
EP
In the following sections we shall denote coexisting phases with the superscripts 0 and 00 ; the former indicates the phase with fixed composition, and the latter the phase in equilibrium with it.
AC C
It is important to distinguish between partial derivatives and derivatives calculated along some special path. For example, • (∂ρ0i /∂T )ρj6=i is the partial derivative of the concentration of component i in phase 0 with respect to temperature, i.e., the sensitivity of ρ0i to temperature changes while all other concentrations ρj6=i are kept constant. • dρ0i /dT denotes the sensitivity of ρ0i to temperature changes while phase equilibrium is maintained.
2.1
The Gibbs–Konovalov equations
For the readers’ convenience, we briefly report the isoplethic Gibbs–Konovalov equations for N -component mixtures and their derivation.
2
ACCEPTED MANUSCRIPT
12
RI PT
10
p/MPa
8
6
SC
4
0 50
100
150
M AN U
2
200
250
300
350
400
T/K
D
Fig. 1: Phase diagram of the (methane + propane system), calculated with the Peng–Robinson equation of state. ◦: pure-fluid critical point, •: binary critical point, grey curves: —: critical curve, · · · · · · : vapor pressure curve; black curves: isopleth: —: liquid branch, - - -: vapor branch.
TE
For an equilibrium between two phases (denoted here by 0 and 00 ), the chemical potentials of all components i have to be equal, µ0i (p, T, x0 ) = µ00i (p, T, x00 ),
i = 1, . . . , N ,
(1)
AC C
EP
where xχ denotes the vector of all mole fractions of phase χ. This equilibrium condition already implies that pressure and temperature must be the same for both phases. As it must be fulfilled along a phase boundary curve, the variations of the chemical potentials must also be equal, !
dµ0i = dµ00i ,
(2)
i = 1, . . . , N .
The total differential of the chemical potential of component i in T px coordinates is dµi = −Si dT + Vi dp +
N −1 X j=1
∂µi ∂xj
dxj ,
(3)
where Si and Vi are the partial molar entropy and partial molar volume, respectively, of component i, ∂S Si ≡ ∂ni p,T,nj6=i (4) ∂V Vi ≡ . ∂ni p,T,nj6=i 3
ACCEPTED MANUSCRIPT
Here ni denotes the amount of substance of component i. Substituting Eq. (3) into Eq. (2) and rearranging the terms gives N −1 X j=1
∂µ0i ∂xj
dx0j −
N −1 X j=1
∂µ00i ∂xj
dx00j = (Si0 − Si00 ) dT − (Vi0 − Vi00 ) dp i = 1, . . . , N .
RI PT
= −∆Si dT + ∆Vi dp ,
(5)
The ∆ indicates a difference of a property between the two phases, e.g., ∆Si ≡ Si00 − Si0 .
Along an isoplethic phase envelope, the mole fractions of the 0 phase are kept constant, which means dx0i = 0 for all components i, and therefore the first term on the left-hand side vanishes. Multiplying Eq. (5) by x00i and taking the sum over all components i results in x00i
i=1
N −1 X j=1
∂µ00i ∂xj
dx00j = −
N X
N X
SC
−
N X
x00i ∆Si dT +
i=1
x00i ∆Vi dp .
(6)
i=1
j=1
M AN U
Switching the order of the summations in the left-hand side term yields # "N N N N −1 X X X X ∂µ00 i =− x00i ∆Si dT + x00i ∆Vi dp . − dx00j x00i ∂xj i=1
i=1
(7)
i=1
The term in brackets is zero because of the Gibbs–Duhem equation. Hence Eq. (6) reduces to 0=−
N X
x00i ∆Si dT +
D
i=1
N X
x00i ∆Vi dp .
(8)
i=1
EP
TE
Dividing this equation by dT yields a differential equation for the isoplethic phase envelope of the 0 phase, PN 00 xi ∆Si dp = Pi=1 . (9) N 00 dT i=1 xi ∆Vi This equation can evidently be regarded as a generalization of Clapeyron’s equation for multicomponent mixtures.
AC C
In the textbook of Haase [4, Sec. 50] this equation appears as PN 00 dp i=1 xi ∆Hi = P , 00 dT T N i=1 xi ∆Vi
(10)
where the ∆Hi denote differences pf partial molar enthalpies. Furthermore, dividing Eq. (5) by −dT (again assuming dx0i = 0 for all components) gives N −1 X j=1
∂µ00i ∂xj
dx00j dp = ∆Si − ∆Vi , dT dT
i = 1, . . . , N .
(11)
We have thus obtained a single differential equation which describes the run of the phase envelope curve in a pT diagram, and a set of linear equations that can be solved for the derivatives of the mole fractions of the 00 phase. In principle, the complete isoplethic phase envelope can therefore be obtained by integrating these differential equations. 4
ACCEPTED MANUSCRIPT
2.2
RI PT
Such an approach, however, is not practical: The partial molar entropies and volumes are derivatives calculated at constant pressure, whereas the equations of state of fluids are usually functions of density. The conversion between pressure- and density-dependent properties can, in principle, be made with Jacobi determinants, but this requires the calculation of additional derivatives. Alternatively, the partial molar derivatives could be obtained directly by numerical differentiation, but then one would have to solve the equation of state for the density within the differentiation algorithm, which is not only time-consuming, but affects the precision of the results.
Density-based formulation
M AN U
SC
The classical thermodynamic formulation of equilibrium conditions, which is based on mole fractions and either pressure or molar volumes, has the disadvantage that the mole fractions are not independent. Moreover, mole fractions and molar volumes have different dimensions, so that the xVm space has no well-defined metric, and this precludes the application of some tools from differential geometry, unless the molar volumes are made dimensionless by the introduction of an artificial normalizing factor. But there exists a more elegant formulation. In the density-based thermodynamic formalism, the composition of the system is expressed in terms of molar densities or concentrations, ρi ≡
ni , V
(12)
TE
D
where ni denotes the amount of substance of component i in a system having the total volume V . The formulation is sometimes called “isochoric thermodynamics”, as working with molar densities is equivalent to working with amounts of substance in a fixed volume. Conversion between ρi and the usual molar volume/mole fraction system is achieved by ρ=
N X
ρi =
i=1
1 , Vm
xi =
ρi . ρ
(13)
EP
xi ρi = , Vm
AC C
The central thermodynamic potential in the density-based formulation of thermodynamics [9] is the Helmholtz energy density Ψ, Ψ(ρ, T ) ≡
A , V
(14)
where A denotes the Helmholtz energy and ρ ≡ (ρ1 , ρ2 , . . . , ρN ) the vector of molar densities. The chemical potentials are given by ∂Ψ µi = i = 1, . . . , N (15) ∂ρi T,ρj6=i and the pressure by1 p = −Ψ +
N X ∂Ψ i=1
1
∂ρi
ρi .
In the previous publication [5] the factor ρi had erroneously been omitted.
5
(16)
ACCEPTED MANUSCRIPT
Using vector notation, these equations can be rewritten [9] as µ = ∇Ψ
(17)
p = −Ψ + µ · ρ ,
(18)
RI PT
and
where the dot indicates the scalar (inner) product, and ∇ differentiation with respect to all densities ρi .
··· ··· .. .
Ψ12 Ψ22 .. .
Ψ1N Ψ2N .. .
with
Ψij ≡
∂2Ψ ∂ρi ∂ρj
.
(19)
T
M AN U
Ψ≡
Ψ11 Ψ21 .. .
SC
The local curvature of the Ψ(ρ) surface is described by the Hessian matrix of the Helmholtz energy density,
ΨN 1 ΨN 2 · · · ΨN N
Ψ is symmetric and, for a stable state, positive definite.
As shown elsewhere [5], the total differential of the chemical potential can be obtained from Eq. (17) by calculating the derivatives with respect to temperature and all densities, ∂µ dµ = Ψ dρ + dT . (20) ∂T ρ
TE
D
Similarly, the total differential of the pressure can be obtained from Eq. (18), ∂p dp = (Ψρ) · dρ + dT . ∂T ρ
(21)
AC C
EP
As in the previous section, we postulate that the variations of the chemical potentials and the pressure must be equal along the phase boundary, 00 0 ! ∂µ ∂µ 00 0 00 00 0 0 dµ − dµ = Ψ dρ − Ψ dρ + − dT ∂T ρ ∂T ρ (22) ∂µ ! 00 00 0 0 = Ψ dρ − Ψ dρ + ∆ dT = 0 , ∂T ρ ∂p ∂p dp − dp = (Ψ ρ ) · dρ − (Ψ ρ ) · dρ + − ∂T ρ ∂T ρ ∂p ! = (Ψ00 ρ00 ) · dρ00 − (Ψ0 ρ0 ) · dρ0 + ∆ =0. ∂T ρ 00
0
00 00
00
0 0
0
(23)
Multiplication of Eq. (22) by ρ00 and observing the equality (Ψ dρ) · ρ = (Ψρ) · dρ, which holds for symmetric matrices Ψ [5], yields ! ∂µ 00 00 00 0 00 0 (Ψ ρ ) · dρ − (Ψ ρ ) · dρ + ∆ · ρ00 dT = 0 . (24) ∂T ρ 6
ACCEPTED MANUSCRIPT
Replacing the first term with Eq. (23) gives 0 0
0
(Ψ ρ ) · dρ − ∆
∂p ∂T
0 00
0
− (Ψ ρ ) · dρ +
∆
ρ
∂µ ∂T
!
· ρ00 dT = 0
(25)
ρ
dρ0 Ψ (ρ − ρ ) · = dT 0
00
0
∆
∂µ ∂T
!
RI PT
and, after division by dT ,
00
·ρ −∆ ρ
∂p ∂T
.
(26)
ρ
This is a system of linear equations for dρ0 /dT , the vector of derivatives of the concentrations ρ0i with respect to temperature along the isoplethic phase envelope.
Then Eq. (26) becomes 0
0 dρ
Ψ ∆ρ · x
0
dT
M AN U
dρ0 dρ0 = x0 . dT dT
=
SC
This equation can be further simplified by using the condition of an isopleth, namely that all independent mole fractions x0k , k = 1, . . . , N − 1 are constant. Hence the derivative vector can be written as
∆
∂µ ∂T
!
00
·ρ −∆
ρ
∂p ∂T
TE
D
from which the derivative is obtained immediately, ! ∂p ∂µ 00 ·ρ −∆ ∆ ∂T ρ ∂T ρ dρ0 = . dT Ψ0 ∆ρ · x0
(27)
,
(28)
ρ
(29)
AC C
EP
Once dρ0 /dT and therefore the dρ0i /dT are known, the dρ00i /dT can be obtained from Eq. (22) by writing it as 00 0 ∂µ 00 dρ 0 dρ , Ψ =Ψ −∆ (30) dT dT ∂T ρ which is a system of linear equations for the derivatives dρ00 /dT .
2.3
Infinite dilution
In principle it is not necessary to study the case of infinite dilution for isoplethic calculations: An N -component mixture in which the density of one component is zero can of course be treated as an (N −1)-component mixture in which all components have positive concentrations. But from the programmer’s viewpoint it is advantageous to have an algorithm that is executable even if one or more components have zero concentrations. This problem has been considered in a previous publication already [5]. Therefore we shall merely present the results here.
7
ACCEPTED MANUSCRIPT
For ρi → 0, the main-diagonal elements Ψχii of the Hessian matrices diverge; the leading term is RT /ρχi . In Eq. (26), however, merely products Ψ0ii ρχi occur, which remain finite. In particular, it can be shown that lim Ψ0ii ρ0i = RT 0
lim Ψ0ii ρ00i = RT φi
ρ0i →0
with φi = exp −
00 µr,∞ i
− RT
0 µr,∞ i
!
χ
,
RI PT
ρi →0
(31)
Evidently, there are no divergent terms.
M AN U
SC
denotes the residual chemical potential of component i in phase χ at infinite diluwhere µr,∞ i tion. The ith element of the Ψ0 ∆ρ product in Eq. (29) thus becomes Ψ0i1 ∆ρ1 . . .. .. Ψ0i · ∆ρ = RT = Ψ0i1 ∆ρ1 + . . . + RT (φi − 1) + . . . + Ψ0iN ∆ρN . (32) 0 · ∆ρi ρ i . .. . . . 0 ∆ρ ΨiN N
D
In order to suppress divergent behavior in Eq. (30), we multiply it from the left side by a diagonal matrix formed from a density vector, 00 dρ0 ∂µ dρ00 0 0 0 0 00 = (diagρ Ψ ) · − diagρ ∆ . (diagρ Ψ ) · (33) dT dT ∂T ρ
TE
For ρ0i , ρ00i → 0, the resulting matrix product on the left-hand side has row vectors of the form (34) (diagρ00 Ψ00 )i = 0, . . . , RT φ−1 i ,...,0 ,
EP
and the matrix product on the right-hand side of the form (diagρ0 Ψ0 )i = (0, . . . , RT, . . . , 0) .
(35)
2.4
AC C
Again, there are no divergent terms left, and therefore Eq. (33) remains solvable. Critical points of mixtures
At the critical point of a mixture the coexisting phases become identical. This implies ρ0 = ρ00 , and hence all ∆ properties in Eq. (29) become zero. Consequently, this equation breaks down and cannot be used any longer to determine the density derivatives dρ0 /dT . On the other hand, there are some regularities that are valid at critical points: • For a stable state, the Ψ(ρ) surface must be convex. The convexity is lost at a critical point. The local geometry can be described with an eigenvalue equation Ψuk = λk uk ,
k = 1, . . . , N ,
(36)
where the eigenvalues λk are measures of the local curvature. At a critical point, the lowermost eigenvalue, λ1 , is zero, i.e., Ψ(ρ) becomes “flat” in one direction. The associated 8
ACCEPTED MANUSCRIPT
normalized eigenvector u1 points into this direction, which is the direction of the phase split, i.e., u1 is aligned with ∆ρ [9,10] and, consequently, with the temperature derivative, ∆
dρ = f u1 . dT
(37)
RI PT
f is a yet to be determined scalar. • We conjecture that, in a density–density diagram, the phase split is symmetrical close to a critical point.
ρ0 = ρc +
dρ0 dρc f dT = ρc + dT − dT u1 , dT dT 2
(38)
D
which yields
M AN U
SC
For binary mixtures a simple method for obtaining the derivative vectors dρχ /dT can then be constructed from geometric considerations. The situation is depicted in Fig. 2: Fig. 2a shows three isothermal phase envelopes of the (methane + propane) system. The middle isotherm (251.6 K for the parameters used here) has its critical point at x1 = 0.8; this critical point corresponds to the one shown in Fig. 1. The isoplethic condition x01 = 0.8 is represented by a straight line. The projection of the dρ0 /dT vector lies on this line; the projection of the dρ00 /dT vector forms an angle with the other vector, and the critical curve bisects this angle. Fig. 2b shows that the isopleth (the curve of the 0 phase) lies in the x1 = 0.8 plane. The curve of the coexisting phase (00 phase) intersects the isotherms on the other side of the critical point; both curves meet at the critical point. Here the connode joining the equilibrium phases (in this diagram the ∆ρ vector) has zero length and the direction of the critical eigenvector. The density along the isopleth can evidently be described by
TE
dρ0 dρc f = − u1 . dT dT 2
(39)
EP
Here dρc /dT denotes the slope of the critical curve. Because of Eq. (27), this equation results in x0
dρ0 dρc = − f u1 . dT dT
(40)
AC C
For a binary mixture, this is a well-defined system of linear equations, which can be solved for dρ0 /dT and f . Once these are known, dρ00 dρc f = + u1 dT dT 2
(41)
yields the other derivative vector. The situation is more complicated for multicomponent mixtures. Eq. (40) remains valid in principle, but now there is no uniquely defined critical curve anymore, so that dρc /dT is no longer well-defined. Because of Gibbs’ phase rule, the critical states of an N -component mixture lie on an (N − 1)-dimensional (hyper-)surface in the (N + 1)-dimensional (ρ, T ) space. It is possible, however, to define 1-dimensional critical curves along which the mole fraction of one component, xj , varies freely, while the ratios of all other mole fractions xi6=j remain fixed: xj xi6=j
= 0 ... xj ... 1 = xi,xj =0 . . . (1 − xj )xi,xj =0 . . . 0 9
(42)
ACCEPTED MANUSCRIPT
12
246.6 K 251.6 K 256.6 K
ρ2/(mol dm–3)
8
6
2
0
2
4
6 8 ρ1/(moldm –3)
10
12
M AN U
0
SC
4
RI PT
10
ρ2
AC C
EP
TE
T
D
(a) Results obtained with the Peng–Robinson equation of state. The arrows are projections of the dρχ /dT vectors at the critical point.
ρ1
(b) Schematic three-dimensional representation. Blue curve: isopleth x01 = 0.8, red: locus of 00 phase. Fig. 2: Density–density diagram of the (methane + propane) system. Black curves: —: isothermal phase envelope: liquid branch, - - -: vapor branch, · · · · · · : connode; •: binary critical point; grey curves: —: critical curve,—: x1 = 0.8 locus.
In other words, the fixed-ratio condition defines a (hyper-)plane in (ρ, T ) space, and the critical curve is the intersection locus of that plane with the Ψ(ρc , T ) surface. Let us denote such a critical curve by T → ρc |xj l . Evidently, there are N such critical curves. We can then rewrite 10
ACCEPTED MANUSCRIPT
12
RI PT
10
lg c2
8
6
SC
4
0
–0.04
–0.02
M AN U
2
0
0.02
0.04
f/(moldm –3 K–1)
D
Fig. 3: Norm of the c vector (cf. Eqs. (43–45)) as a function of the scaling factor f for the calculation of critical points of the (methane + propane + decane) system (PPR78 equation of state). —: xc = (0.60, 0.36, 0.04) at 337.67 K, - - -: xc = (0.40, 0.50, 0.10) at 419.75 K.
TE
Eq. (40) as
N 0 c X dρ dρ = x0 cj dT dT j=1
− xj l
f u1 , 2
(43)
AC C
EP
where the cj denote a set of coefficients; they must obey the equation N X
cj = 1 .
(44)
j=1
Eqs. (43) and (44) represent N + 1 linear equations for N + 2 unknowns. A unique solution can be found by treating f as a given value and then observing that, during a variation of f , c2 must be minimal. Examples are shown in Fig. 3 for the calculation of dρ0 /dT at two critical points of the ternary mixture (methane + propane + decane) using the PPR78 method [11, 12]. The minima are evident (note the logarithmic scale!) and can be located easily. Similar to the binary case (Eq. (41)), the derivative of the coexisting phase is obtained from N dρ00 X dρc f = cj + u1 . (45) dT dT 2 j=1
xj l
Hence it is possible to determine the density derivatives at a mixture critical point, but the evaluation involves the calculation of eigenvectors and the temperature derivative of the density vector along the critical curve(s), which means a significant computational effort. 11
ACCEPTED MANUSCRIPT
For the calculation of critical states of mixtures several methods have been proposed, e.g. ˜ those of Heidemann and Khalil [13] or Michelsen [14]. Here the method of Quinones-Cisneros [10, Sec. 5.8.4] was used, which is based on the analysis of the eigenvalues of the Ψ matrix, an entity already occurring in the differential equation formulation.
Application
RI PT
2.5
If an initial two-phase equilibrium state is known, Eqs. (27) and (29) allow a straightforward calculation of the derivative vector dρ0 /dT , and then Eq. (33) represents a system of linear equations for dρ00 /dT . It is therefore possible to formulate the task of computing the densities of the coexisting phases as a system of ordinary differential equations.
SC
The solution can be obtained with adaptive integration algorithms; for this work, the method of Cash and Karp was used [15].
M AN U
Fig. 4 shows the density–temperature functions for the (methane + propane) system, which are underlying the phase diagram Fig. 1. In this representation the isopleth (0 phase) and the locus of the coexistent phase (00 phase) can be well distinguished. The density derivatives at the critical point according to Eqs. (40–41) have been indicated; they match the slopes of the phase envelopes well.
D
Even if the density derivatives at the critical point can be computed, it is not advisable to make the integration “run through” a critical point, as the numerical precision of the ∆ properties in Eq. (29) decreases in the vicinity of critical states. Moreover, the calculation of derivatives along critical curves is computationally costly. A much better alternative is extrapolating the isopleth beyond the critical point, determining the accurate phase compositions by means of an algebraic solver, and then restarting the integration from there.
AC C
EP
TE
One might wonder whether the calculation of phase envelopes from differential equations is reliable at all, for numerical errors tend to accumulate in the course of the integration procedure. Table 1 presents calculated data for a piece of the x01 = 0.8 isopleth of the (methane + propane) system, obtained either by integrating the differential equations (29) and (33) or by solving the algebraic equations of phase equilbrium as described in [9]. It turns out that the agreement between the two methods is very good over a range of over 80 K. As expected, the differential equation method becomes unreliable close to the critical point (251.6 K for this isopleth) due to the degradation of the Ψχ matrices as well as he cancellation of terms in the right-hand side of Eq. (26). This problem can be remediated by following the integration with a single postiteration step as outlined in Appendix B: then the results of the algebraic solver and the differential equation method agree fully. The derivatives of Ψ(ρ, T ) were obtained by numerical differentiation (by means of Romberg’s scheme [16]) or alternatively by analytic differentiation. Table 1 shows that the deviations are negligible, and that numerical diffrentiation can be reliable indeed. The integration of the differential equations is not affected by azeotropy, for the ρχ vectors differ even if the xχ vectors agree. In cases where there are three-phase equilibria, the algorithm would follow an isopleth through metastable and unstable portions. Checking the sign of the lowermost eigenvalues of the Hessian matrices Ψχ would reveal (diffusion) instability. But the proof that a fluid–fluid equilibrium is a stable one requires a global stability test—including tests for the formation of solid phases. Figs. 1 and 4 show that the isopleth passes through a temperature maximum, the maxcondentherm state, which is typical for systems containing a supercritical component. Hence it is not possible to compute the complete isoplethic envelope in one sweep by integration over 12
ACCEPTED MANUSCRIPT
12
RI PT
10
ρ1/(mol dm–3)
8
6
SC
4
0 50
100
150
M AN U
2
200
250
300
350
400
T/K
(a) Molar concentration of methane.
D
12
EP
ρ2/(mol dm–3)
8
TE
10
6
AC C
4
2
0 50
100
150
200
250
300
350
400
T/K
(b) Molar concentrations of propane. Fig. 4: Molar concentrations of the (methane + propane) system, calculated with the Peng–Robinson equation of state. ◦: pure-fluid critical point, •: mixture critical point; black curves: —: isopleth x1 = 0.8, liquid branch, - - -: vapor branch; grey curves: —: critical curve, —: coexisting (00 ) phase of the liquid, - - -: of the vapor. The arrows indicate the derivatives at the critical point (Eqs. (40–41)).
13
ACCEPTED MANUSCRIPT
Table 1: Equilibrium phase properties along the x01 = 0.8 isopleth of the (methane + propane) system, calculated with the Peng–Robinson equation of state. States calculated (a) by solving the algebraic equations of phase equilibrium, (d) by integrating the differential equations (using (n) numerical/(a) analytical differentiation, (p) 1 postiteration step).
200 210 220
230
240
250
*
ρ002 3 (cm /mol)
0.88131 1.33302 1.33302 1.91979 1.91978 2.65237 2.65237 3.53312 3.53312 4.54998 4.54998 5.66766 5.66765 6.81920 6.81918 6.81918 7.91044 7.91038 7.91039 8.84221 8.84191 8.84195 8.84221 9.52413 9.51547 9.51571 9.52413
0.99975 0.99940 0.99940 0.99865 0.99865 0.99707 0.99707 0.99376 0.99376 0.98697 0.98698 0.97369 0.97369 0.95025 0.95026 0.95026 0.91484 0.91485 0.91485 0.86811 0.86820 0.86819 0.86811 0.81038 0.81809 0.81795 0.81039
1.86404×10−2 1.79509×10−2 1.79508×10−2 1.72049×10−2 1.72049×10−2 1.64009×10−2 1.64009×10−2 1.55401×10−2 1.55401×10−2 1.46281×10−2 1.46281×10−2 1.36746×10−2 1.36746×10−2 1.26884×10−2 1.26884×10−2 1.26884×10−2 1.16693×10−2 1.16693×10−2 1.16693×10−2 1.06039×10−2 1.06037×10−2 1.06037×10−2 1.06039×10−2 9.46714×10−3 9.46090×10−3 9.46108×10−3 9.46714×10−3
4.66011×10−3 4.48771×10−3 4.48772×10−3 4.30122×10−3 4.30123×10−3 4.10022×10−3 4.10023×10−3 3.88503×10−3 3.88504×10−3 3.65703×10−3 3.65704×10−3 3.41864×10−3 3.41865×10−3 3.17210×10−3 3.17210×10−3 3.17210×10−3 2.91734×10−3 2.91733×10−3 2.91733×10−3 2.65097×10−3 2.65092×10−3 2.65093×10−3 2.65097×10−3 2.36678×10−3 2.36523×10−3 2.36527×10−3 2.36678×10−3
8.32530×10−4 1.25364×10−3 1.25364×10−3 1.82646×10−3 1.82646×10−3 2.59412×10−3 2.59412×10−3 3.60485×10−3 3.60484×10−3 4.88809×10−3 4.88808×10−3 6.38285×10−3 6.38282×10−3 7.83836×10−3 7.83825×10−3 7.83826×10−3 8.89470×10−3 8.89428×10−3 8.89430×10−3 9.37505×10−3 9.37226×10−3 9.37262×10−3 9.37505×10−3 9.32639×10−3 9.20286×10−3 9.20516×10−3 9.32622×10−3
2.12038×10−7 7.52267×10−7 7.52264×10−7 2.46367×10−6 2.46366×10−6 7.62456×10−6 7.62450×10−6 2.26320×10−5 2.26318×10−5 6.45077×10−5 6.45070×10−5 1.72472×10−4 1.72469×10−4 4.10332×10−4 4.10311×10−4 4.10314×10−4 8.28008×10−4 8.27877×10−4 8.27882×10−4 1.42439×10−3 1.42276×10−3 1.42297×10−3 1.42439×10−3 2.18234×10−3 2.04638×10−3 2.04873×10−3 2.18213×10−3
SC
RI PT
ρ001 3 (cm /mol)
M AN U
190
ρ02 3 (cm /mol)
D
180
ρ01 3 (cm /mol)
TE
170
a* a d a d a d a d a d a d a dn da a dn da a dn da dp a dn da dp
x001
AC C
150 160
p MPa
EP
T K
: initial state for the differential equation method
temperature. In situations where this is a disadvantage, another property than temperature can be used as “scan variable”. This problem has recently been explored by Venkatarathnam [17]. For the test system presented here, the overall density of the fixed-composition phase, ρ0 , is a good choice. Fig. 5 shows that this property is monotonic along the phase envelope and that T (ρ0 ) as well as ρ00i (ρ0 ) behave well (i.e., do not fold back or have poles). The differential
14
ACCEPTED MANUSCRIPT
40
350
35
300
30 25
T/K
250 200
20
T
15
100
SC
150
10
ρ"1 ρ"2
0
5
M AN U
50 0
ρ"(mol dm–3) i
RI PT
400
5
10
15
20
25
30
ρ’/(mol dm–3)
TE
D
Fig. 5: Temperature and molar concentrations of the coexisting phase, ρ001 and ρ002 , as function of the density of the fixed-composition phase along the x1 = 0.8 isopleth of the (methane + propane) system, calculated with the Peng– Robinson equation of state. •: mixture critical point; black curves: —: temperature along the liquid branch, - - -: vapor branch; grey curves: —: molar concentrations of coexisting phase of liquid, - - -: of the vapor.
equations to be integrated are in this case
(Ψ0 ∆ρ) · x0 from Eq. (29) ! ∂p ∂µ 00 ·ρ −∆ ∆ ∂T ρ ∂T ρ dρ00 dT ∂µ Ψ00 0 = Ψ0 x0 − 0 ∆ from Eq. (30) . dρ dρ ∂T ρ
(46)
AC C
EP
dT = dρ0
Most of the calculations presented here were made with the Peng–Robinson equation of state, using common quadratic mixing rules [7,8]. But the method is not restricted to cubic equations of state: It can be applied to any equation of state that can be represented as Helmholtz energy density Ψ(ρ, T ). It must be emphasized that, for the method proposed here, it is not necessary to calculate the molar volumes of phases for given pressure or temperature, nor is it necessary to iteratively solve nonlinear equations for the phase equilibrium conditions. This makes the new method particularly useful in connection with complicated non-cubic equations of state. The calculations reported here were made with the ThermoC software package [18], which uses numerical differentiation to obtain the required derivatives. Alternatively, analytic differentiation can be used, of course. The method employed for numerical differentiation as well as analytic expressions of the derivatives of the Peng–Robinson equation of state have been given in a previous publication [5].
15
ACCEPTED MANUSCRIPT
3
Pure fluids
3.1
The vapor pressure curve
For a pure fluid (N = 1, x1 = 1) Eq. (26) reduces to 0
dT
=
∆
∂µ1 ∂T
!
00
ρ −∆ ρ
where the only element of the Hessian matrix is
∂p ∂Vm
,
1 ρ2 κ T
=+
T
(47)
ρ
(48)
M AN U
Ψ11
1 ∂ 2 (Am ρ) =− 3 = ∂ρ2 ρ
∂p ∂T
SC
dρ Ψ011 ∆ρ
RI PT
The isopleth of a pure fluid is its vapor pressure curve. Therefore the application of the formalism of isochoric thermodynamics suggests itself.
because of (∂Am /∂Vm ) = −p. κT is the isothermal compressibility. Furthermore, the temperature derivative of the chemical potential of a pure fluid is ∂µ ∂Gm 1 ∂p = = −Sm + . (49) ∂T ρ ∂T ρ ρ ∂T ρ Inserting Eqs. (48) and (49) into Eq. (47) yields
! ∂p0 ∆Sm ρ ρ + ∆ρ ∂T ρ 0 ! ∂p ρ0 ρ00 + ∆Sm ∆ρ ∂T ρ 0 00
(50)
TE
D
ρ0 κ0 dρ0 =− T dT ∆ρ =
−ρ0 κ0T
EP
and, by switching the phase labels,
ρ0 ρ00 ∆Sm + ∆ρ
∂p00 ∂T
! .
(51)
ρ
AC C
dρ00 = −ρ00 κ00T dT
These two equations describe the variation of the orthobaric densities with temperature along the vapor pressure curve; they are the equivalents of Clapeyron’s equation in the density-based thermodynamic formalism. Incidentally, ρ0 ρ00 /∆ρ = −1/∆Vm . Nowadays the calculation of vapor pressures from an equation of state is usually accomplished by means of solving Maxwell’s criterion, ∆Am = pσ (Vm00 − Vm0 ) ,
(52)
where pσ denotes the vapor pressure and Z ∆Am = −
00 Vm
0 Vm
p(Vm , T ) dVm = RT ln
ρ00 + ρ0
Z
ρ00
ρ0
Z(ρ, T ) − 1 dρ . ρ
(53)
ρ0 and ρ00 are roots of the equation of state for the given pressure pσ , i.e., p(Vmχ , T ) = pσ . For a given temperatures, one guesses pσ , determins the roots Vmχ , evaluates Maxwell’s criterion, 16
ACCEPTED MANUSCRIPT
and then varies pσ and repeats the calculation until the criterion is fulfilled. Finding the roots is easy for so-called cubic equations of state, but may be rather complicated and time-consuming otherwise2 . Furthermore, it may happen that the equation of state does not have two suitable roots, in which case another pσ value has to be tried.
3.2
RI PT
Alternatively, one might calculate vapor pressure curves by numerically integrating the differential equations (50) and (51) and then simply using pσ = p(ρχ , T ). This way, there is no need to solve the equation of state for the phase densities, and convergence problems are eliminated.
The low-pressure region
M AN U
SC
Far away from the critical point, the vapor pressure is low and the vapor may be treated as an ideal gas. Identifying the 0 phase with the liquid and the 00 phase with the vapor, we may therefore assume ∆Sm ρ00 (∂p/∂T )ρ0 and ∆ρ ≈ −ρ0 . Then Eq. (50) reduces to 0 0 dρ0 ∂p ∂p 1 ∂Vm0 1 ∂Vm0 0 0 ≈ −ρ κT =+ 02 =− 02 = −ρ0 αp0 (54) dT ∂T ρ ∂p T ∂T Vm ∂T p Vm Vm or
dVm0 ≈ dT
∂Vm0 ∂T
= Vm0 αp0 ,
(55)
p
which means that the density variation of the liquid phase along the vapor pressure curve is approximately determined by the isobaric thermal expansivity.
D
An analogous treatment of Eq. (51), observing κ00T ≈ (RT ρ00 )−1 and (∂p00 /∂T )ρ ≈ Rρ00 , leads to
TE
dρ00 ρ00 ≈ (∆Sm − R) dT RT
(56)
EP
in the low-pressure limit. Substituting ρ00 = p/(RT ) finally gives dp ∆Sm =p . dT RT
(57)
AC C
With the usual assumption ∆Hm = T ∆Sm = const, integration yields the Clausius–Clapeyron equation.
4
Conclusion
As a complement to a previous publication, the differential equations for isoplethic phase envelopes of multicomponent mixtures were derived in the framework of isochoric thermodynamics. The calculation of phase envelopes thus becomes an initial-value problem. It is of a kind that can be integrated by means of standard adaptive methods, for instance the so-called Runge–Kutta-45 algorithm or the Cash–Karp algorithm [15, 21]; the latter was used for this work. 2
There exist algorithms which avoid the calculation of phase volumes from pressure, e.g., the method of Kretzschmar et al. [19], the so-called area method [20], or the isochoric thermodynamics algorithm [9]; the last two methods were developed for mixtures, but can be applied to pure substances, too.
17
ACCEPTED MANUSCRIPT
The core of the phase envelope algorithm is, of course, the calculation of the temperature derivatives of the orthobaric density vectors, dρχ /dT . The implementation in a computer program is straightforward, particularly if a programming language is used that supports vector and matrix operations. The necessary steps have been summarized in Alg. 1 for the readers’ convenience.
RI PT
The only thermodynamic derivatives needed are first- and second-order derivatives of the Helmholtz energy density with respect to temperature and molar densities; there is no need for differentiations at constant pressure. The differentiations can be performed numerically or analytically; derivatives for the Peng–Robinson equation have been presented in the previous publication.
SC
The computation of phase envelopes from differential equations, as proposed here, can be used with arbitrary equations of state. It is very efficient particularly in connection with noncubic equations of state, as it does not require the calculation of phase densities from given pressure and temperature.
M AN U
The differential equation method, however, requires an accurate initial state, and this has to be computed by solving the algebraic equations of phase equilibrium. Therefore the differential equation method cannot replace the algebraic method; but it can complement it by providing a rapid inter- and extrapolation of isopleths. In the vicinity of critical points, the differential equation method becomes unrealiable. On the other hand, it is not affected by azeotropy.
TE
Acknowledgments
D
Application of the isopleth differential equations to pure substances generates analoga of Clapeyron’s equation in the framework of isochoric thermodynamics. Integration yields the orthobaric densities. The differential equations may be useful for the rapid calculation of phase densities in power plant simulations.
AC C
EP
˜ We wish to thank Dr. Sergio E. Quinones Cisneros (University of Cologne and F-Thermo Services GmbH, Cologne, Germany) and Dr. Ian Bell (NIST Boulder (Colorado), USA) for helpful discussions and for critically reading the manuscript.
Symbols A BX cj
Helmholtz energy
critical amplitude of the property X weight factor for “mixing” density derivative vectors along critical curves
c
vector of coefficients, c = (c1 , c2 , . . . , cN )
ej
jth unit vector in ρ space
G
Gibbs energy
N
number of components
ni
amount of substance of component i
p
pressure
R
universal gas constant 18
ACCEPTED MANUSCRIPT
Algorithm 1: Calculation of density derivatives for isopleth calculations, dρ0 /dT and dρ00 /dT , from densities ρ0 , ρ00 and temperature T .
SC
M AN U
vector φ, b, x0 , µχT (χ =0 ,00 ) χ matrix CP , Ψχ (χ =0 ,00 ) 0 0 x = ρ /( i ρ0i ) forall ρ0i = 0 do calculate φi [Eq. (31)] for phases 0 , 00 do Ψχ = Hessian(Ψ(ρχ , T )) [note #1] µχT = ∇(∂Ψχ /∂T ) end b = Ψ0 ∆ρ forall ρ0i = 0 do bi + = RT (φi − 1) [note #2] dρ0 /dT = ((∆µT ) · ρ00 − ∆(∂p/∂T )) /(b · x0 ) [Eq. (29)] s0 = dρ0 /dT x0
RI PT
Input: fundamental equation Ψ(ρ, T ) and its parameters, density vectors ρ0 , ρ00 , temperature T Result: vectors of derivatives s0 ≡ dρ0 /dT , s00 ≡ dρ00 /dT
D
C 0 = (diag ρ0 Ψ0 ) forall ρ0i = 0 do Cii0 = RT b = C 0 s0 − diag ρ0 ∆µT C 00 = diagρ00 Ψ00 forall ρ00i = 0 do Cii00 = RT /φi solve C 00 · s00 = b [cf. Eq. (33), note #3] return s0 , s00
EP
TE
#1: If a ρχi is zero, assign a large, but finite value to the corresponding Ψχii . #2: + =: increment operator. #3: a linear problem.
partial molar entropy with respect to component i
T
temperature
uk
kth normalized eigenvector of Ψ
V Vi xi x
AC C
Si
volume
partial molar volume with respect to component i mole fraction of component i mole fraction vector
Z
compression factor, Z = p/(ρRT )
α
critical exponent of the isochoric heat capacity
αp
isobaric thermal expansivity
β
critical exponent of orthobaric densities or entropies
γ
critical exponent of the isothermal compressibility
critical exponent of the orhobaric isochoric pressure coefficients
φi
density ratio of component i, φi = ρ00i /ρ0i 19
isothermal compressibility
λk
kth eigenvalue of Ψ, λ1 = min(λ1 , . . . , λN )
µi
chemical potential of component i
µ
chemical potential vector, µ = (µ1 , µ2 , . . . , µN )
ρ
total density, ρ = Vm−1
ρi
molar density, concentration of component i, ρi = xi /Vm
ρ
density vector, ρ = (ρ1 , ρ2 , . . . , ρN )
Ψ
Helmholtz energy density, Ψ = A/V
Ψ
Hessian matrix of Ψ
Ψij
matrix element of Ψ, Ψij = (∂ 2 Ψ/(∂ρi ∂ρj ))
Ψi
ith row vector of Ψ
M AN U
SC
κT
RI PT
ACCEPTED MANUSCRIPT
Subscripts c pure-fluid critical state m
molar property
xj l
curve along which xj is varied and the ratios of all mole fractions xi6=j are fixed
D
xj = 0 value attained if xj = 0
TE
Superscripts c property on a critical curve state residual property
χ
phase indicator
0
phase indicator; along isopleths, the phase with fixed composition
00
phase indicator; along isopleths, the phase coexisting with the 0 phase
∞
infinite dilution
AC C
EP
r
References
¨ ¨ [1] D. Konowalow, Uber die Dampfspannungen der Flussigkeitsgemische, Ann. Phys. 250 (1881) 34–52, volume also cited as Wied. Ann. 14. ¨ ¨ [2] D. Konowalow, Uber die Dampfspannung von gemischten Flussigkeiten, Ann. Phys. 250 (1881) 219–226, volume also cited as Wied. Ann. 14. [3] J. W. Gibbs, On the equilibrium of heterogeneous substances, in: R. G. van Name, H. A. Bumstead (Eds.), The Scientific Papers of J. Willard Gibbs, Vol. 1: Thermodynamics, Dover Publications, New York, 1961, reprint of a book of 1906 containing articles originally published in 1876–1878. 20
ACCEPTED MANUSCRIPT
[4] R. Haase, Thermodynamik der Mischphasen, Springer, Berlin, 1956. [5] U. K. Deiters, Differential equations for the calculation of fluid phase equilibria, Fluid Phase Equilib. 428 (2016) 164–173. doi:10.1016/j.fluid.2016.04.014.
RI PT
[6] R. J. Sadus, Molecular simulation of the phase behavior of fluids and fluid mixtures using the synthetic method, J. Chem. Phys. 137 (2012) 054507:1–7. doi:10.1063/1.4739853. [7] D. Y. Peng, D. B. Robinson, A new two-constant equation of state, Ind. Eng. Chem. Fundam. 15 (1976) 59–64. [8] D. B. Robinson, D. Y. Peng, The characterization of the heptanes and heavier fractions for the GPA Peng–Robinson program, GPA Research Report RR-28 (1978) 1–36.
SC
˜ [9] S. E. Quinones-Cisneros, U. K. Deiters, An efficient algorithm for the calculation of phase envelopes of fluid mixtures, Fluid Phase Equilib. 329 (2012) 22–31. doi:10.1016/j.fluid.2012.05.023.
M AN U
[10] U. K. Deiters, T. Kraska, High-Pressure Fluid Phase Equilibria—Phenomenology and Computation, Vol. 2 of Supercritical Fluid Science and Technology (E. Kiran, ed.), Elsevier, Amsterdam, 2012. [11] J.-N. Jaubert, R. Privat, F. Mutelet, Predicting the phase equilibria of synthetic petroleum fluids with the PPR78 approach, AIChE J. 56 (2010) 3225–3235. doi:10.1002/aic.12232.
D
[12] R. Privat, J.-N. Jaubert, Classification of global fluid phase equilibrium behaviors in binary systems, Chem. Eng. Res. Des. 91 (2013) 1807–1839.
TE
[13] R. A. Heidemann, A. M. Khalil, The calculation of critical points, AIChE J. 26 (1980) 769– 779. [14] M. L. Michelsen, Calculation of critical points and phase boundaries in the critical region, Fluid Phase Equilib. 16 (1984) 57–76.
EP
[15] J. R. Cash, A. H. Karp, A variable-order Runge–Kutta method for initial-value problems with rapidly varying right-hand sides, ACM Trans. Math. Software 16 (1990) 201–222. doi:10.1145/79505.79507.
AC C
¨ Ingenieure, Bibliographisches [16] G. Jordan-Engeln, F. Reutter, Numerische Mathematik fur Institut, Mannheim, 1985. [17] G. Venkatarathnam, Density marching method for calculating phase envelopes, Ind. Eng. Chem. Res. 53 (2014) 3723–3730. doi:10.1021/ie403633d. [18] U. K. Deiters, ThermoC project homepage: http://thermoc.uni-koeln.de/index.html, accessed on 2016-02-17. [19] H.-J. Kretzschmar, T. Zschunke, J. Klinger, A. Dittmann, An alternative method for the numerical calculation of the Maxwell criterion in vapour pressure computations, in: ˇ M. P´ıchal, O. Sifner (Eds.), Proceedings of the 11th International Conference on Properties of Water and Steam, Hemisphere, New York, 1990, pp. 324–331. [20] A. E. Elhassan, R. J. B. Craven, K. M. de Reuck, The Area method for pure fluids and an analysis of the two-phase region, Fluid Phase Equilib. 130 (1997) 167–187.
21
ACCEPTED MANUSCRIPT
AC C
EP
TE
D
M AN U
SC
RI PT
[21] W. H. Press, S. A. Teukolsky, W. T. Vetterling, B. P. Flannery, Numerical Recipes in C, Cambridge University Press, Cambridge, 2002.
22
ACCEPTED MANUSCRIPT
Appendix A
Derivatives at a mixture critical point under isothermal conditions
For completeness’ sake we briefly summarize equations for the calculation of the density derivatives at a critical point in an isothermal phase diagram.
c
RI PT
• For a binary mixture, a critical point is always an extremum of the isothermal phase diagram, hence dρ0 dρ00 (58) = →∞. dp dp c
M AN U
SC
• For a ternary mixture, the derivatives can be obtained in a similar fashion as in Section 2.4. To describe the phase boundary at constant temperature and under variation of x0j (keeping the ratio of the mole fractions of the other components constant), we can set dρc f dρ0 = − u1 , (59) dp 0 dp 2 xj l
T
where dρc /dp|T is the slope of the critical curve at constant temperature and u1 the normalized critical eigenvector. The scalar f must be adjusted in such a way that the ratio of the components of dρ0i6=j /dp|x0j l has the correct value,
x0j l
1 dρ0k − 0 xk,0 dp
D
1 dρ0i x0i,0 dp
= 0,
i, k 6= j .
(60)
x0j l
EP
TE
Together, Eqs. (59) and (60) constitute a system of linear equations for f and the components of dρ0 /dp. The derivative vector of the coexisting phase (00 ) is then dρ00 dρc f = + u1 . (61) dp 0 dp 2 xj l
T
AC C
For example, if the isothermal phase envelope is calculated with the constraints T = const, x01 = const, Eq. (59) becomes dρc f dρ0 = − u1 (62) dp 0 dp 2 x1 l
T
and Eq. (60)
1 x02,x1 =0
dρ02 dp
− x01 l
1 x03,x1 =0
dρ03 dp
=0.
(63)
x01 l
• For a quaternary or higher mixture, Eq. (60) can be fulfilled by writing the derivative
23
ACCEPTED MANUSCRIPT
vector of the 0 phase as dρ0 dp
and
(64)
SC
with
0 . .. 0 ej ≡ 1 0 .. . 0
RI PT
dρ0j ej + gx0xj =0 dp x0j l x01,xj =0 .. . 0 xj−1,xj =0 , x0xj =0 ≡ 0 x0 j+1,xj =0 .. . 0 xN,xj =0 =
M AN U
where g is a yet to be determined scalar. There exist N critical curves p → ρc |T,xj l . Hence Eq. (59) becomes N X dρ0j dρc f 0 ck ej + gxxj =0 − = − u1 dp dp 2 k=1
xk l
N X
ck = 1
(65)
k=1
D
N X k=1
dT c ck dp
=0.
xk l
AC C
EP
TE
The last equation, which contains the derivatives of the temperature along the critical curves, implements the T = const condition. Eq. (65) is a set of linear equations for the unknowns dρ0j /dp, g, and c. As in Section 2.4, the physically reasonable solution is found by minimizing c2 with respect to f . Finally, the derivative vector for the 00 phase is obtained from N dρ00 X dρc f = ck + u1 . (66) dp dp 2 k=1
xk l
24
ACCEPTED MANUSCRIPT
Appendix B
Postiteration
If desired, the accuracy of isopleth calculations with the differential equation method can be improved by postiteration. The system of equations to be solved is µ00i (ρ00 , T ) = µ0i (ρ0 , T ) 00
0
i = 1, . . . , N
(67)
0
p (ρ , T ) = p (ρ , T ) . The vector of residuals is therefore
N p00
−
N p0
(68)
SC
µ001 − µ01 .. . y= . 00 0 µ −µ
RI PT
00
M AN U
As ρ0 = ρ0 x0 and x0 is fixed during an isopleth calculation, the vector of unknowns can be defined as 00 ρ1 .. z= . . (69) 00 ρ N
ρ0
D
A single step of a Newton–Raphson iteration can then be described by z := z − δz ,
TE
(70)
where the correction vector δz is obtained by solving the linear system of equations (71)
EP
B δz = y
AC C
and the matrix B = {Bij } is defined by 00 ∂µi = Ψ00ij ∂ρ00 00 j ρ ,T 0 k6=j 0 0 i − ∂µ ∂ρ0 x0 ,T = −Ψi · x Bij = 00 ∂p 00 00 ∂ρ00j ρ ,T = Ψj · ρ k6 = j − ∂p0 = −ρ0 x0 · Ψ0 x0 ∂ρ0
x0 ,T
i = 1, . . . N ∧ j = 1, . . . N i = 1, . . . , N ∧ j = N + 1 . i = N + 1 ∧ j = 1, . . . , N
(72)
i=N +1 ∧ j =N +1
Usually a single iteration step is sufficient. As the derivatives required for the postiteration also appear in the differential equations, the additional programming effort is negligible.
25
ACCEPTED MANUSCRIPT
Highlights: • differential equations of isoplethic phase envelopes • applicable to multicomponent systems • can be rapidly integrated; avoids convergence problems or internal volume calculations
AC C
EP
TE D
M AN U
SC
RI PT
• for pure fluids: density-based Clapeyron equation (gives orthobaric densities)
1
S0378-3812(17)30122-X
DOI:
10.1016/j.fluid.2017.03.022
Reference:
FLUID 11439
To appear in:
Fluid Phase Equilibria
Received Date: 15 November 2016 Revised Date:
25 January 2017
Accepted Date: 27 March 2017
Please cite this article as: U.K. Deiters, Differential equations for the calculation of isopleths of multicomponent fluid mixtures, Fluid Phase Equilibria (2017), doi: 10.1016/j.fluid.2017.03.022. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
x’ = const
AC C
EP
TE D
M AN U
SC
RI PT
p
T
RI PT
ACCEPTED MANUSCRIPT
Ulrich K. Deiters
SC
Differential equations for the calculation of isopleths of multicomponent fluid mixtures
AC C
EP
TE
D
M AN U
¨ Germany Institute of Physical Chemistry, University of Cologne, Luxemburger Str. 116, D-50939 Koln,
ACCEPTED MANUSCRIPT
Abstract
RI PT
An isoplethic phase envelope is a locus of two-phase equilibria where one of the phases is kept at constant composition. The thermodynamic conditions for such a phase envelope are expressed as a set of first-order ordinary differential equations in which molar concentrations appear as central variables instead of mole fractions. The resulting formalism is applicable to multicomponent mixtures, well-suited for machine calculations, and allows a very rapid calculation of phase envelopes with arbitrary equations of state.
SC
Application to the formalism to pure fluids yields simple differential equations that are analoga of Clapeyron’s equation, but permit the calculation of orthobaric densities.
1
M AN U
Keywords: fluid phase equilibria, Gibbs–Konovalov equations, algorithm, multicomponent mixtures, isochoric thermodynamics
Introduction
Thermodynamics is built on differential calculus, and therefore it is not surprising that many problems can be expressed by means of algebraic equations or differential equations. Wellknown examples are
TE
D
• the calculation of vapor-pressure curves, which can be accomplished by evaluating Maxwell’s equal-area criterion (algebraic equation) or by integrating the Clausius–Clapeyron equation (differential equation),
EP
• or the calculation of adiabatic expansion curves for ideal gases, for which probably all thermodynamics textbooks of the world present a differential equation, but which can also be accomplished by setting the entropy to a constant value and solving the resulting algebraic equation.
AC C
Likewise, phase equilibria of fluid mixtures can be calculated from the conditions that the temperatures, the pressures, and the chemical potentials of the components of coexisting phases must agree, which leads to a set of nonlinear algebraic equations. Differential equations were already developed at the end of the 19th Century by Konovalov [1,2] and Gibbs [3] and became known as Gibbs–Konovalov rules; a particularly clear and concise derivation of these rules for multicomponent mixtures can be found in the textbook of Haase [4]. The Gibbs–Konovalov rules can be used to prove some thermodynamic theorems, e.g., that, in isothermal phase diagrams of binary mixtures, azeotropes must lie at pressure extrema. Because of their inconvenient mathematical structure—an example is discussed in the next section—they never became popular for computing phase envelopes. But this is true for the classical formulation in terms of mole fractions and partial molar derivatives only. It could recently been shown that, in a density-based formulation of thermodynamics, there exist analoga to the Gibbs–Konovalov rules that can be used to develop very efficient and robust algorithms for the computation of isothermal as well as isobaric phase envelopes [5].
1
ACCEPTED MANUSCRIPT
RI PT
There is another kind of phase envelope with a paramount importance for chemical engineering, namely the isoplethic phase envelope. It represents the boundary between the singlephase and the two-phase regions for a mixture with a given composition (i.e., set of mole fractions) and comprises of states at which boiling or condensation sets in. The so-called synthetic method for experimentally determining phase equilibria, according to which samples of known overall composition a manipulated until phase transitions are detected, yields isopleths as primary results. Recently the synthetic method has even been applied to computer simulations of fluid mixtures [6]. We shall therefore extend the previous work and derive differential equations for isoplethic phase envelopes. Particular attention will be given to the behavior of these equations at critical points of mixtures.
2
Differential equations for isopleths
M AN U
SC
The isopleth of a pure fluid is its vapor pressure curve. We shall derive differential equations for the orthobaric densities—analoga of Clapeyron’s equation in the density-based formalism.
D
Fig. 1 serves to illustrate the problem. It shows the phase diagrams of the (methane + propane) system, calculated with the Peng–Robinson equation of state with standard quadratic mixing rules [7, 8]. This is a simple system with an uninterrupted critical curve connecting the critical points of the pure fluid. The diagram furthermore contains the x1 = 0.8 isopleth—a curve marking the onset of boiling of a liquid having this composition, or the onset of condensation of a vapor having this composition. The phase in equilibrium with the boiling liquid, or the condensing vapor, respectively, has a different composition, but the same pressure and temperature, and therefore cannot be shown in this diagram.
TE
The isopleth meets the critical curve of the mixture in a binary critical point. It should be noted that this critical point usually does not lie at a maximum of the isopleth. The task is to calculate the isopleth efficiently, producing not only its pT contour, but also the composition of the coexisting phase.
EP
In the following sections we shall denote coexisting phases with the superscripts 0 and 00 ; the former indicates the phase with fixed composition, and the latter the phase in equilibrium with it.
AC C
It is important to distinguish between partial derivatives and derivatives calculated along some special path. For example, • (∂ρ0i /∂T )ρj6=i is the partial derivative of the concentration of component i in phase 0 with respect to temperature, i.e., the sensitivity of ρ0i to temperature changes while all other concentrations ρj6=i are kept constant. • dρ0i /dT denotes the sensitivity of ρ0i to temperature changes while phase equilibrium is maintained.
2.1
The Gibbs–Konovalov equations
For the readers’ convenience, we briefly report the isoplethic Gibbs–Konovalov equations for N -component mixtures and their derivation.
2
ACCEPTED MANUSCRIPT
12
RI PT
10
p/MPa
8
6
SC
4
0 50
100
150
M AN U
2
200
250
300
350
400
T/K
D
Fig. 1: Phase diagram of the (methane + propane system), calculated with the Peng–Robinson equation of state. ◦: pure-fluid critical point, •: binary critical point, grey curves: —: critical curve, · · · · · · : vapor pressure curve; black curves: isopleth: —: liquid branch, - - -: vapor branch.
TE
For an equilibrium between two phases (denoted here by 0 and 00 ), the chemical potentials of all components i have to be equal, µ0i (p, T, x0 ) = µ00i (p, T, x00 ),
i = 1, . . . , N ,
(1)
AC C
EP
where xχ denotes the vector of all mole fractions of phase χ. This equilibrium condition already implies that pressure and temperature must be the same for both phases. As it must be fulfilled along a phase boundary curve, the variations of the chemical potentials must also be equal, !
dµ0i = dµ00i ,
(2)
i = 1, . . . , N .
The total differential of the chemical potential of component i in T px coordinates is dµi = −Si dT + Vi dp +
N −1 X j=1
∂µi ∂xj
dxj ,
(3)
where Si and Vi are the partial molar entropy and partial molar volume, respectively, of component i, ∂S Si ≡ ∂ni p,T,nj6=i (4) ∂V Vi ≡ . ∂ni p,T,nj6=i 3
ACCEPTED MANUSCRIPT
Here ni denotes the amount of substance of component i. Substituting Eq. (3) into Eq. (2) and rearranging the terms gives N −1 X j=1
∂µ0i ∂xj
dx0j −
N −1 X j=1
∂µ00i ∂xj
dx00j = (Si0 − Si00 ) dT − (Vi0 − Vi00 ) dp i = 1, . . . , N .
RI PT
= −∆Si dT + ∆Vi dp ,
(5)
The ∆ indicates a difference of a property between the two phases, e.g., ∆Si ≡ Si00 − Si0 .
Along an isoplethic phase envelope, the mole fractions of the 0 phase are kept constant, which means dx0i = 0 for all components i, and therefore the first term on the left-hand side vanishes. Multiplying Eq. (5) by x00i and taking the sum over all components i results in x00i
i=1
N −1 X j=1
∂µ00i ∂xj
dx00j = −
N X
N X
SC
−
N X
x00i ∆Si dT +
i=1
x00i ∆Vi dp .
(6)
i=1
j=1
M AN U
Switching the order of the summations in the left-hand side term yields # "N N N N −1 X X X X ∂µ00 i =− x00i ∆Si dT + x00i ∆Vi dp . − dx00j x00i ∂xj i=1
i=1
(7)
i=1
The term in brackets is zero because of the Gibbs–Duhem equation. Hence Eq. (6) reduces to 0=−
N X
x00i ∆Si dT +
D
i=1
N X
x00i ∆Vi dp .
(8)
i=1
EP
TE
Dividing this equation by dT yields a differential equation for the isoplethic phase envelope of the 0 phase, PN 00 xi ∆Si dp = Pi=1 . (9) N 00 dT i=1 xi ∆Vi This equation can evidently be regarded as a generalization of Clapeyron’s equation for multicomponent mixtures.
AC C
In the textbook of Haase [4, Sec. 50] this equation appears as PN 00 dp i=1 xi ∆Hi = P , 00 dT T N i=1 xi ∆Vi
(10)
where the ∆Hi denote differences pf partial molar enthalpies. Furthermore, dividing Eq. (5) by −dT (again assuming dx0i = 0 for all components) gives N −1 X j=1
∂µ00i ∂xj
dx00j dp = ∆Si − ∆Vi , dT dT
i = 1, . . . , N .
(11)
We have thus obtained a single differential equation which describes the run of the phase envelope curve in a pT diagram, and a set of linear equations that can be solved for the derivatives of the mole fractions of the 00 phase. In principle, the complete isoplethic phase envelope can therefore be obtained by integrating these differential equations. 4
ACCEPTED MANUSCRIPT
2.2
RI PT
Such an approach, however, is not practical: The partial molar entropies and volumes are derivatives calculated at constant pressure, whereas the equations of state of fluids are usually functions of density. The conversion between pressure- and density-dependent properties can, in principle, be made with Jacobi determinants, but this requires the calculation of additional derivatives. Alternatively, the partial molar derivatives could be obtained directly by numerical differentiation, but then one would have to solve the equation of state for the density within the differentiation algorithm, which is not only time-consuming, but affects the precision of the results.
Density-based formulation
M AN U
SC
The classical thermodynamic formulation of equilibrium conditions, which is based on mole fractions and either pressure or molar volumes, has the disadvantage that the mole fractions are not independent. Moreover, mole fractions and molar volumes have different dimensions, so that the xVm space has no well-defined metric, and this precludes the application of some tools from differential geometry, unless the molar volumes are made dimensionless by the introduction of an artificial normalizing factor. But there exists a more elegant formulation. In the density-based thermodynamic formalism, the composition of the system is expressed in terms of molar densities or concentrations, ρi ≡
ni , V
(12)
TE
D
where ni denotes the amount of substance of component i in a system having the total volume V . The formulation is sometimes called “isochoric thermodynamics”, as working with molar densities is equivalent to working with amounts of substance in a fixed volume. Conversion between ρi and the usual molar volume/mole fraction system is achieved by ρ=
N X
ρi =
i=1
1 , Vm
xi =
ρi . ρ
(13)
EP
xi ρi = , Vm
AC C
The central thermodynamic potential in the density-based formulation of thermodynamics [9] is the Helmholtz energy density Ψ, Ψ(ρ, T ) ≡
A , V
(14)
where A denotes the Helmholtz energy and ρ ≡ (ρ1 , ρ2 , . . . , ρN ) the vector of molar densities. The chemical potentials are given by ∂Ψ µi = i = 1, . . . , N (15) ∂ρi T,ρj6=i and the pressure by1 p = −Ψ +
N X ∂Ψ i=1
1
∂ρi
ρi .
In the previous publication [5] the factor ρi had erroneously been omitted.
5
(16)
ACCEPTED MANUSCRIPT
Using vector notation, these equations can be rewritten [9] as µ = ∇Ψ
(17)
p = −Ψ + µ · ρ ,
(18)
RI PT
and
where the dot indicates the scalar (inner) product, and ∇ differentiation with respect to all densities ρi .
··· ··· .. .
Ψ12 Ψ22 .. .
Ψ1N Ψ2N .. .
with
Ψij ≡
∂2Ψ ∂ρi ∂ρj
.
(19)
T
M AN U
Ψ≡
Ψ11 Ψ21 .. .
SC
The local curvature of the Ψ(ρ) surface is described by the Hessian matrix of the Helmholtz energy density,
ΨN 1 ΨN 2 · · · ΨN N
Ψ is symmetric and, for a stable state, positive definite.
As shown elsewhere [5], the total differential of the chemical potential can be obtained from Eq. (17) by calculating the derivatives with respect to temperature and all densities, ∂µ dµ = Ψ dρ + dT . (20) ∂T ρ
TE
D
Similarly, the total differential of the pressure can be obtained from Eq. (18), ∂p dp = (Ψρ) · dρ + dT . ∂T ρ
(21)
AC C
EP
As in the previous section, we postulate that the variations of the chemical potentials and the pressure must be equal along the phase boundary, 00 0 ! ∂µ ∂µ 00 0 00 00 0 0 dµ − dµ = Ψ dρ − Ψ dρ + − dT ∂T ρ ∂T ρ (22) ∂µ ! 00 00 0 0 = Ψ dρ − Ψ dρ + ∆ dT = 0 , ∂T ρ ∂p ∂p dp − dp = (Ψ ρ ) · dρ − (Ψ ρ ) · dρ + − ∂T ρ ∂T ρ ∂p ! = (Ψ00 ρ00 ) · dρ00 − (Ψ0 ρ0 ) · dρ0 + ∆ =0. ∂T ρ 00
0
00 00
00
0 0
0
(23)
Multiplication of Eq. (22) by ρ00 and observing the equality (Ψ dρ) · ρ = (Ψρ) · dρ, which holds for symmetric matrices Ψ [5], yields ! ∂µ 00 00 00 0 00 0 (Ψ ρ ) · dρ − (Ψ ρ ) · dρ + ∆ · ρ00 dT = 0 . (24) ∂T ρ 6
ACCEPTED MANUSCRIPT
Replacing the first term with Eq. (23) gives 0 0
0
(Ψ ρ ) · dρ − ∆
∂p ∂T
0 00
0
− (Ψ ρ ) · dρ +
∆
ρ
∂µ ∂T
!
· ρ00 dT = 0
(25)
ρ
dρ0 Ψ (ρ − ρ ) · = dT 0
00
0
∆
∂µ ∂T
!
RI PT
and, after division by dT ,
00
·ρ −∆ ρ
∂p ∂T
.
(26)
ρ
This is a system of linear equations for dρ0 /dT , the vector of derivatives of the concentrations ρ0i with respect to temperature along the isoplethic phase envelope.
Then Eq. (26) becomes 0
0 dρ
Ψ ∆ρ · x
0
dT
M AN U
dρ0 dρ0 = x0 . dT dT
=
SC
This equation can be further simplified by using the condition of an isopleth, namely that all independent mole fractions x0k , k = 1, . . . , N − 1 are constant. Hence the derivative vector can be written as
∆
∂µ ∂T
!
00
·ρ −∆
ρ
∂p ∂T
TE
D
from which the derivative is obtained immediately, ! ∂p ∂µ 00 ·ρ −∆ ∆ ∂T ρ ∂T ρ dρ0 = . dT Ψ0 ∆ρ · x0
(27)
,
(28)
ρ
(29)
AC C
EP
Once dρ0 /dT and therefore the dρ0i /dT are known, the dρ00i /dT can be obtained from Eq. (22) by writing it as 00 0 ∂µ 00 dρ 0 dρ , Ψ =Ψ −∆ (30) dT dT ∂T ρ which is a system of linear equations for the derivatives dρ00 /dT .
2.3
Infinite dilution
In principle it is not necessary to study the case of infinite dilution for isoplethic calculations: An N -component mixture in which the density of one component is zero can of course be treated as an (N −1)-component mixture in which all components have positive concentrations. But from the programmer’s viewpoint it is advantageous to have an algorithm that is executable even if one or more components have zero concentrations. This problem has been considered in a previous publication already [5]. Therefore we shall merely present the results here.
7
ACCEPTED MANUSCRIPT
For ρi → 0, the main-diagonal elements Ψχii of the Hessian matrices diverge; the leading term is RT /ρχi . In Eq. (26), however, merely products Ψ0ii ρχi occur, which remain finite. In particular, it can be shown that lim Ψ0ii ρ0i = RT 0
lim Ψ0ii ρ00i = RT φi
ρ0i →0
with φi = exp −
00 µr,∞ i
− RT
0 µr,∞ i
!
χ
,
RI PT
ρi →0
(31)
Evidently, there are no divergent terms.
M AN U
SC
denotes the residual chemical potential of component i in phase χ at infinite diluwhere µr,∞ i tion. The ith element of the Ψ0 ∆ρ product in Eq. (29) thus becomes Ψ0i1 ∆ρ1 . . .. .. Ψ0i · ∆ρ = RT = Ψ0i1 ∆ρ1 + . . . + RT (φi − 1) + . . . + Ψ0iN ∆ρN . (32) 0 · ∆ρi ρ i . .. . . . 0 ∆ρ ΨiN N
D
In order to suppress divergent behavior in Eq. (30), we multiply it from the left side by a diagonal matrix formed from a density vector, 00 dρ0 ∂µ dρ00 0 0 0 0 00 = (diagρ Ψ ) · − diagρ ∆ . (diagρ Ψ ) · (33) dT dT ∂T ρ
TE
For ρ0i , ρ00i → 0, the resulting matrix product on the left-hand side has row vectors of the form (34) (diagρ00 Ψ00 )i = 0, . . . , RT φ−1 i ,...,0 ,
EP
and the matrix product on the right-hand side of the form (diagρ0 Ψ0 )i = (0, . . . , RT, . . . , 0) .
(35)
2.4
AC C
Again, there are no divergent terms left, and therefore Eq. (33) remains solvable. Critical points of mixtures
At the critical point of a mixture the coexisting phases become identical. This implies ρ0 = ρ00 , and hence all ∆ properties in Eq. (29) become zero. Consequently, this equation breaks down and cannot be used any longer to determine the density derivatives dρ0 /dT . On the other hand, there are some regularities that are valid at critical points: • For a stable state, the Ψ(ρ) surface must be convex. The convexity is lost at a critical point. The local geometry can be described with an eigenvalue equation Ψuk = λk uk ,
k = 1, . . . , N ,
(36)
where the eigenvalues λk are measures of the local curvature. At a critical point, the lowermost eigenvalue, λ1 , is zero, i.e., Ψ(ρ) becomes “flat” in one direction. The associated 8
ACCEPTED MANUSCRIPT
normalized eigenvector u1 points into this direction, which is the direction of the phase split, i.e., u1 is aligned with ∆ρ [9,10] and, consequently, with the temperature derivative, ∆
dρ = f u1 . dT
(37)
RI PT
f is a yet to be determined scalar. • We conjecture that, in a density–density diagram, the phase split is symmetrical close to a critical point.
ρ0 = ρc +
dρ0 dρc f dT = ρc + dT − dT u1 , dT dT 2
(38)
D
which yields
M AN U
SC
For binary mixtures a simple method for obtaining the derivative vectors dρχ /dT can then be constructed from geometric considerations. The situation is depicted in Fig. 2: Fig. 2a shows three isothermal phase envelopes of the (methane + propane) system. The middle isotherm (251.6 K for the parameters used here) has its critical point at x1 = 0.8; this critical point corresponds to the one shown in Fig. 1. The isoplethic condition x01 = 0.8 is represented by a straight line. The projection of the dρ0 /dT vector lies on this line; the projection of the dρ00 /dT vector forms an angle with the other vector, and the critical curve bisects this angle. Fig. 2b shows that the isopleth (the curve of the 0 phase) lies in the x1 = 0.8 plane. The curve of the coexisting phase (00 phase) intersects the isotherms on the other side of the critical point; both curves meet at the critical point. Here the connode joining the equilibrium phases (in this diagram the ∆ρ vector) has zero length and the direction of the critical eigenvector. The density along the isopleth can evidently be described by
TE
dρ0 dρc f = − u1 . dT dT 2
(39)
EP
Here dρc /dT denotes the slope of the critical curve. Because of Eq. (27), this equation results in x0
dρ0 dρc = − f u1 . dT dT
(40)
AC C
For a binary mixture, this is a well-defined system of linear equations, which can be solved for dρ0 /dT and f . Once these are known, dρ00 dρc f = + u1 dT dT 2
(41)
yields the other derivative vector. The situation is more complicated for multicomponent mixtures. Eq. (40) remains valid in principle, but now there is no uniquely defined critical curve anymore, so that dρc /dT is no longer well-defined. Because of Gibbs’ phase rule, the critical states of an N -component mixture lie on an (N − 1)-dimensional (hyper-)surface in the (N + 1)-dimensional (ρ, T ) space. It is possible, however, to define 1-dimensional critical curves along which the mole fraction of one component, xj , varies freely, while the ratios of all other mole fractions xi6=j remain fixed: xj xi6=j
= 0 ... xj ... 1 = xi,xj =0 . . . (1 − xj )xi,xj =0 . . . 0 9
(42)
ACCEPTED MANUSCRIPT
12
246.6 K 251.6 K 256.6 K
ρ2/(mol dm–3)
8
6
2
0
2
4
6 8 ρ1/(moldm –3)
10
12
M AN U
0
SC
4
RI PT
10
ρ2
AC C
EP
TE
T
D
(a) Results obtained with the Peng–Robinson equation of state. The arrows are projections of the dρχ /dT vectors at the critical point.
ρ1
(b) Schematic three-dimensional representation. Blue curve: isopleth x01 = 0.8, red: locus of 00 phase. Fig. 2: Density–density diagram of the (methane + propane) system. Black curves: —: isothermal phase envelope: liquid branch, - - -: vapor branch, · · · · · · : connode; •: binary critical point; grey curves: —: critical curve,—: x1 = 0.8 locus.
In other words, the fixed-ratio condition defines a (hyper-)plane in (ρ, T ) space, and the critical curve is the intersection locus of that plane with the Ψ(ρc , T ) surface. Let us denote such a critical curve by T → ρc |xj l . Evidently, there are N such critical curves. We can then rewrite 10
ACCEPTED MANUSCRIPT
12
RI PT
10
lg c2
8
6
SC
4
0
–0.04
–0.02
M AN U
2
0
0.02
0.04
f/(moldm –3 K–1)
D
Fig. 3: Norm of the c vector (cf. Eqs. (43–45)) as a function of the scaling factor f for the calculation of critical points of the (methane + propane + decane) system (PPR78 equation of state). —: xc = (0.60, 0.36, 0.04) at 337.67 K, - - -: xc = (0.40, 0.50, 0.10) at 419.75 K.
TE
Eq. (40) as
N 0 c X dρ dρ = x0 cj dT dT j=1
− xj l
f u1 , 2
(43)
AC C
EP
where the cj denote a set of coefficients; they must obey the equation N X
cj = 1 .
(44)
j=1
Eqs. (43) and (44) represent N + 1 linear equations for N + 2 unknowns. A unique solution can be found by treating f as a given value and then observing that, during a variation of f , c2 must be minimal. Examples are shown in Fig. 3 for the calculation of dρ0 /dT at two critical points of the ternary mixture (methane + propane + decane) using the PPR78 method [11, 12]. The minima are evident (note the logarithmic scale!) and can be located easily. Similar to the binary case (Eq. (41)), the derivative of the coexisting phase is obtained from N dρ00 X dρc f = cj + u1 . (45) dT dT 2 j=1
xj l
Hence it is possible to determine the density derivatives at a mixture critical point, but the evaluation involves the calculation of eigenvectors and the temperature derivative of the density vector along the critical curve(s), which means a significant computational effort. 11
ACCEPTED MANUSCRIPT
For the calculation of critical states of mixtures several methods have been proposed, e.g. ˜ those of Heidemann and Khalil [13] or Michelsen [14]. Here the method of Quinones-Cisneros [10, Sec. 5.8.4] was used, which is based on the analysis of the eigenvalues of the Ψ matrix, an entity already occurring in the differential equation formulation.
Application
RI PT
2.5
If an initial two-phase equilibrium state is known, Eqs. (27) and (29) allow a straightforward calculation of the derivative vector dρ0 /dT , and then Eq. (33) represents a system of linear equations for dρ00 /dT . It is therefore possible to formulate the task of computing the densities of the coexisting phases as a system of ordinary differential equations.
SC
The solution can be obtained with adaptive integration algorithms; for this work, the method of Cash and Karp was used [15].
M AN U
Fig. 4 shows the density–temperature functions for the (methane + propane) system, which are underlying the phase diagram Fig. 1. In this representation the isopleth (0 phase) and the locus of the coexistent phase (00 phase) can be well distinguished. The density derivatives at the critical point according to Eqs. (40–41) have been indicated; they match the slopes of the phase envelopes well.
D
Even if the density derivatives at the critical point can be computed, it is not advisable to make the integration “run through” a critical point, as the numerical precision of the ∆ properties in Eq. (29) decreases in the vicinity of critical states. Moreover, the calculation of derivatives along critical curves is computationally costly. A much better alternative is extrapolating the isopleth beyond the critical point, determining the accurate phase compositions by means of an algebraic solver, and then restarting the integration from there.
AC C
EP
TE
One might wonder whether the calculation of phase envelopes from differential equations is reliable at all, for numerical errors tend to accumulate in the course of the integration procedure. Table 1 presents calculated data for a piece of the x01 = 0.8 isopleth of the (methane + propane) system, obtained either by integrating the differential equations (29) and (33) or by solving the algebraic equations of phase equilbrium as described in [9]. It turns out that the agreement between the two methods is very good over a range of over 80 K. As expected, the differential equation method becomes unreliable close to the critical point (251.6 K for this isopleth) due to the degradation of the Ψχ matrices as well as he cancellation of terms in the right-hand side of Eq. (26). This problem can be remediated by following the integration with a single postiteration step as outlined in Appendix B: then the results of the algebraic solver and the differential equation method agree fully. The derivatives of Ψ(ρ, T ) were obtained by numerical differentiation (by means of Romberg’s scheme [16]) or alternatively by analytic differentiation. Table 1 shows that the deviations are negligible, and that numerical diffrentiation can be reliable indeed. The integration of the differential equations is not affected by azeotropy, for the ρχ vectors differ even if the xχ vectors agree. In cases where there are three-phase equilibria, the algorithm would follow an isopleth through metastable and unstable portions. Checking the sign of the lowermost eigenvalues of the Hessian matrices Ψχ would reveal (diffusion) instability. But the proof that a fluid–fluid equilibrium is a stable one requires a global stability test—including tests for the formation of solid phases. Figs. 1 and 4 show that the isopleth passes through a temperature maximum, the maxcondentherm state, which is typical for systems containing a supercritical component. Hence it is not possible to compute the complete isoplethic envelope in one sweep by integration over 12
ACCEPTED MANUSCRIPT
12
RI PT
10
ρ1/(mol dm–3)
8
6
SC
4
0 50
100
150
M AN U
2
200
250
300
350
400
T/K
(a) Molar concentration of methane.
D
12
EP
ρ2/(mol dm–3)
8
TE
10
6
AC C
4
2
0 50
100
150
200
250
300
350
400
T/K
(b) Molar concentrations of propane. Fig. 4: Molar concentrations of the (methane + propane) system, calculated with the Peng–Robinson equation of state. ◦: pure-fluid critical point, •: mixture critical point; black curves: —: isopleth x1 = 0.8, liquid branch, - - -: vapor branch; grey curves: —: critical curve, —: coexisting (00 ) phase of the liquid, - - -: of the vapor. The arrows indicate the derivatives at the critical point (Eqs. (40–41)).
13
ACCEPTED MANUSCRIPT
Table 1: Equilibrium phase properties along the x01 = 0.8 isopleth of the (methane + propane) system, calculated with the Peng–Robinson equation of state. States calculated (a) by solving the algebraic equations of phase equilibrium, (d) by integrating the differential equations (using (n) numerical/(a) analytical differentiation, (p) 1 postiteration step).
200 210 220
230
240
250
*
ρ002 3 (cm /mol)
0.88131 1.33302 1.33302 1.91979 1.91978 2.65237 2.65237 3.53312 3.53312 4.54998 4.54998 5.66766 5.66765 6.81920 6.81918 6.81918 7.91044 7.91038 7.91039 8.84221 8.84191 8.84195 8.84221 9.52413 9.51547 9.51571 9.52413
0.99975 0.99940 0.99940 0.99865 0.99865 0.99707 0.99707 0.99376 0.99376 0.98697 0.98698 0.97369 0.97369 0.95025 0.95026 0.95026 0.91484 0.91485 0.91485 0.86811 0.86820 0.86819 0.86811 0.81038 0.81809 0.81795 0.81039
1.86404×10−2 1.79509×10−2 1.79508×10−2 1.72049×10−2 1.72049×10−2 1.64009×10−2 1.64009×10−2 1.55401×10−2 1.55401×10−2 1.46281×10−2 1.46281×10−2 1.36746×10−2 1.36746×10−2 1.26884×10−2 1.26884×10−2 1.26884×10−2 1.16693×10−2 1.16693×10−2 1.16693×10−2 1.06039×10−2 1.06037×10−2 1.06037×10−2 1.06039×10−2 9.46714×10−3 9.46090×10−3 9.46108×10−3 9.46714×10−3
4.66011×10−3 4.48771×10−3 4.48772×10−3 4.30122×10−3 4.30123×10−3 4.10022×10−3 4.10023×10−3 3.88503×10−3 3.88504×10−3 3.65703×10−3 3.65704×10−3 3.41864×10−3 3.41865×10−3 3.17210×10−3 3.17210×10−3 3.17210×10−3 2.91734×10−3 2.91733×10−3 2.91733×10−3 2.65097×10−3 2.65092×10−3 2.65093×10−3 2.65097×10−3 2.36678×10−3 2.36523×10−3 2.36527×10−3 2.36678×10−3
8.32530×10−4 1.25364×10−3 1.25364×10−3 1.82646×10−3 1.82646×10−3 2.59412×10−3 2.59412×10−3 3.60485×10−3 3.60484×10−3 4.88809×10−3 4.88808×10−3 6.38285×10−3 6.38282×10−3 7.83836×10−3 7.83825×10−3 7.83826×10−3 8.89470×10−3 8.89428×10−3 8.89430×10−3 9.37505×10−3 9.37226×10−3 9.37262×10−3 9.37505×10−3 9.32639×10−3 9.20286×10−3 9.20516×10−3 9.32622×10−3
2.12038×10−7 7.52267×10−7 7.52264×10−7 2.46367×10−6 2.46366×10−6 7.62456×10−6 7.62450×10−6 2.26320×10−5 2.26318×10−5 6.45077×10−5 6.45070×10−5 1.72472×10−4 1.72469×10−4 4.10332×10−4 4.10311×10−4 4.10314×10−4 8.28008×10−4 8.27877×10−4 8.27882×10−4 1.42439×10−3 1.42276×10−3 1.42297×10−3 1.42439×10−3 2.18234×10−3 2.04638×10−3 2.04873×10−3 2.18213×10−3
SC
RI PT
ρ001 3 (cm /mol)
M AN U
190
ρ02 3 (cm /mol)
D
180
ρ01 3 (cm /mol)
TE
170
a* a d a d a d a d a d a d a dn da a dn da a dn da dp a dn da dp
x001
AC C
150 160
p MPa
EP
T K
: initial state for the differential equation method
temperature. In situations where this is a disadvantage, another property than temperature can be used as “scan variable”. This problem has recently been explored by Venkatarathnam [17]. For the test system presented here, the overall density of the fixed-composition phase, ρ0 , is a good choice. Fig. 5 shows that this property is monotonic along the phase envelope and that T (ρ0 ) as well as ρ00i (ρ0 ) behave well (i.e., do not fold back or have poles). The differential
14
ACCEPTED MANUSCRIPT
40
350
35
300
30 25
T/K
250 200
20
T
15
100
SC
150
10
ρ"1 ρ"2
0
5
M AN U
50 0
ρ"(mol dm–3) i
RI PT
400
5
10
15
20
25
30
ρ’/(mol dm–3)
TE
D
Fig. 5: Temperature and molar concentrations of the coexisting phase, ρ001 and ρ002 , as function of the density of the fixed-composition phase along the x1 = 0.8 isopleth of the (methane + propane) system, calculated with the Peng– Robinson equation of state. •: mixture critical point; black curves: —: temperature along the liquid branch, - - -: vapor branch; grey curves: —: molar concentrations of coexisting phase of liquid, - - -: of the vapor.
equations to be integrated are in this case
(Ψ0 ∆ρ) · x0 from Eq. (29) ! ∂p ∂µ 00 ·ρ −∆ ∆ ∂T ρ ∂T ρ dρ00 dT ∂µ Ψ00 0 = Ψ0 x0 − 0 ∆ from Eq. (30) . dρ dρ ∂T ρ
(46)
AC C
EP
dT = dρ0
Most of the calculations presented here were made with the Peng–Robinson equation of state, using common quadratic mixing rules [7,8]. But the method is not restricted to cubic equations of state: It can be applied to any equation of state that can be represented as Helmholtz energy density Ψ(ρ, T ). It must be emphasized that, for the method proposed here, it is not necessary to calculate the molar volumes of phases for given pressure or temperature, nor is it necessary to iteratively solve nonlinear equations for the phase equilibrium conditions. This makes the new method particularly useful in connection with complicated non-cubic equations of state. The calculations reported here were made with the ThermoC software package [18], which uses numerical differentiation to obtain the required derivatives. Alternatively, analytic differentiation can be used, of course. The method employed for numerical differentiation as well as analytic expressions of the derivatives of the Peng–Robinson equation of state have been given in a previous publication [5].
15
ACCEPTED MANUSCRIPT
3
Pure fluids
3.1
The vapor pressure curve
For a pure fluid (N = 1, x1 = 1) Eq. (26) reduces to 0
dT
=
∆
∂µ1 ∂T
!
00
ρ −∆ ρ
where the only element of the Hessian matrix is
∂p ∂Vm
,
1 ρ2 κ T
=+
T
(47)
ρ
(48)
M AN U
Ψ11
1 ∂ 2 (Am ρ) =− 3 = ∂ρ2 ρ
∂p ∂T
SC
dρ Ψ011 ∆ρ
RI PT
The isopleth of a pure fluid is its vapor pressure curve. Therefore the application of the formalism of isochoric thermodynamics suggests itself.
because of (∂Am /∂Vm ) = −p. κT is the isothermal compressibility. Furthermore, the temperature derivative of the chemical potential of a pure fluid is ∂µ ∂Gm 1 ∂p = = −Sm + . (49) ∂T ρ ∂T ρ ρ ∂T ρ Inserting Eqs. (48) and (49) into Eq. (47) yields
! ∂p0 ∆Sm ρ ρ + ∆ρ ∂T ρ 0 ! ∂p ρ0 ρ00 + ∆Sm ∆ρ ∂T ρ 0 00
(50)
TE
D
ρ0 κ0 dρ0 =− T dT ∆ρ =
−ρ0 κ0T
EP
and, by switching the phase labels,
ρ0 ρ00 ∆Sm + ∆ρ
∂p00 ∂T
! .
(51)
ρ
AC C
dρ00 = −ρ00 κ00T dT
These two equations describe the variation of the orthobaric densities with temperature along the vapor pressure curve; they are the equivalents of Clapeyron’s equation in the density-based thermodynamic formalism. Incidentally, ρ0 ρ00 /∆ρ = −1/∆Vm . Nowadays the calculation of vapor pressures from an equation of state is usually accomplished by means of solving Maxwell’s criterion, ∆Am = pσ (Vm00 − Vm0 ) ,
(52)
where pσ denotes the vapor pressure and Z ∆Am = −
00 Vm
0 Vm
p(Vm , T ) dVm = RT ln
ρ00 + ρ0
Z
ρ00
ρ0
Z(ρ, T ) − 1 dρ . ρ
(53)
ρ0 and ρ00 are roots of the equation of state for the given pressure pσ , i.e., p(Vmχ , T ) = pσ . For a given temperatures, one guesses pσ , determins the roots Vmχ , evaluates Maxwell’s criterion, 16
ACCEPTED MANUSCRIPT
and then varies pσ and repeats the calculation until the criterion is fulfilled. Finding the roots is easy for so-called cubic equations of state, but may be rather complicated and time-consuming otherwise2 . Furthermore, it may happen that the equation of state does not have two suitable roots, in which case another pσ value has to be tried.
3.2
RI PT
Alternatively, one might calculate vapor pressure curves by numerically integrating the differential equations (50) and (51) and then simply using pσ = p(ρχ , T ). This way, there is no need to solve the equation of state for the phase densities, and convergence problems are eliminated.
The low-pressure region
M AN U
SC
Far away from the critical point, the vapor pressure is low and the vapor may be treated as an ideal gas. Identifying the 0 phase with the liquid and the 00 phase with the vapor, we may therefore assume ∆Sm ρ00 (∂p/∂T )ρ0 and ∆ρ ≈ −ρ0 . Then Eq. (50) reduces to 0 0 dρ0 ∂p ∂p 1 ∂Vm0 1 ∂Vm0 0 0 ≈ −ρ κT =+ 02 =− 02 = −ρ0 αp0 (54) dT ∂T ρ ∂p T ∂T Vm ∂T p Vm Vm or
dVm0 ≈ dT
∂Vm0 ∂T
= Vm0 αp0 ,
(55)
p
which means that the density variation of the liquid phase along the vapor pressure curve is approximately determined by the isobaric thermal expansivity.
D
An analogous treatment of Eq. (51), observing κ00T ≈ (RT ρ00 )−1 and (∂p00 /∂T )ρ ≈ Rρ00 , leads to
TE
dρ00 ρ00 ≈ (∆Sm − R) dT RT
(56)
EP
in the low-pressure limit. Substituting ρ00 = p/(RT ) finally gives dp ∆Sm =p . dT RT
(57)
AC C
With the usual assumption ∆Hm = T ∆Sm = const, integration yields the Clausius–Clapeyron equation.
4
Conclusion
As a complement to a previous publication, the differential equations for isoplethic phase envelopes of multicomponent mixtures were derived in the framework of isochoric thermodynamics. The calculation of phase envelopes thus becomes an initial-value problem. It is of a kind that can be integrated by means of standard adaptive methods, for instance the so-called Runge–Kutta-45 algorithm or the Cash–Karp algorithm [15, 21]; the latter was used for this work. 2
There exist algorithms which avoid the calculation of phase volumes from pressure, e.g., the method of Kretzschmar et al. [19], the so-called area method [20], or the isochoric thermodynamics algorithm [9]; the last two methods were developed for mixtures, but can be applied to pure substances, too.
17
ACCEPTED MANUSCRIPT
The core of the phase envelope algorithm is, of course, the calculation of the temperature derivatives of the orthobaric density vectors, dρχ /dT . The implementation in a computer program is straightforward, particularly if a programming language is used that supports vector and matrix operations. The necessary steps have been summarized in Alg. 1 for the readers’ convenience.
RI PT
The only thermodynamic derivatives needed are first- and second-order derivatives of the Helmholtz energy density with respect to temperature and molar densities; there is no need for differentiations at constant pressure. The differentiations can be performed numerically or analytically; derivatives for the Peng–Robinson equation have been presented in the previous publication.
SC
The computation of phase envelopes from differential equations, as proposed here, can be used with arbitrary equations of state. It is very efficient particularly in connection with noncubic equations of state, as it does not require the calculation of phase densities from given pressure and temperature.
M AN U
The differential equation method, however, requires an accurate initial state, and this has to be computed by solving the algebraic equations of phase equilibrium. Therefore the differential equation method cannot replace the algebraic method; but it can complement it by providing a rapid inter- and extrapolation of isopleths. In the vicinity of critical points, the differential equation method becomes unrealiable. On the other hand, it is not affected by azeotropy.
TE
Acknowledgments
D
Application of the isopleth differential equations to pure substances generates analoga of Clapeyron’s equation in the framework of isochoric thermodynamics. Integration yields the orthobaric densities. The differential equations may be useful for the rapid calculation of phase densities in power plant simulations.
AC C
EP
˜ We wish to thank Dr. Sergio E. Quinones Cisneros (University of Cologne and F-Thermo Services GmbH, Cologne, Germany) and Dr. Ian Bell (NIST Boulder (Colorado), USA) for helpful discussions and for critically reading the manuscript.
Symbols A BX cj
Helmholtz energy
critical amplitude of the property X weight factor for “mixing” density derivative vectors along critical curves
c
vector of coefficients, c = (c1 , c2 , . . . , cN )
ej
jth unit vector in ρ space
G
Gibbs energy
N
number of components
ni
amount of substance of component i
p
pressure
R
universal gas constant 18
ACCEPTED MANUSCRIPT
Algorithm 1: Calculation of density derivatives for isopleth calculations, dρ0 /dT and dρ00 /dT , from densities ρ0 , ρ00 and temperature T .
SC
M AN U
vector φ, b, x0 , µχT (χ =0 ,00 ) χ matrix CP , Ψχ (χ =0 ,00 ) 0 0 x = ρ /( i ρ0i ) forall ρ0i = 0 do calculate φi [Eq. (31)] for phases 0 , 00 do Ψχ = Hessian(Ψ(ρχ , T )) [note #1] µχT = ∇(∂Ψχ /∂T ) end b = Ψ0 ∆ρ forall ρ0i = 0 do bi + = RT (φi − 1) [note #2] dρ0 /dT = ((∆µT ) · ρ00 − ∆(∂p/∂T )) /(b · x0 ) [Eq. (29)] s0 = dρ0 /dT x0
RI PT
Input: fundamental equation Ψ(ρ, T ) and its parameters, density vectors ρ0 , ρ00 , temperature T Result: vectors of derivatives s0 ≡ dρ0 /dT , s00 ≡ dρ00 /dT
D
C 0 = (diag ρ0 Ψ0 ) forall ρ0i = 0 do Cii0 = RT b = C 0 s0 − diag ρ0 ∆µT C 00 = diagρ00 Ψ00 forall ρ00i = 0 do Cii00 = RT /φi solve C 00 · s00 = b [cf. Eq. (33), note #3] return s0 , s00
EP
TE
#1: If a ρχi is zero, assign a large, but finite value to the corresponding Ψχii . #2: + =: increment operator. #3: a linear problem.
partial molar entropy with respect to component i
T
temperature
uk
kth normalized eigenvector of Ψ
V Vi xi x
AC C
Si
volume
partial molar volume with respect to component i mole fraction of component i mole fraction vector
Z
compression factor, Z = p/(ρRT )
α
critical exponent of the isochoric heat capacity
αp
isobaric thermal expansivity
β
critical exponent of orthobaric densities or entropies
γ
critical exponent of the isothermal compressibility
critical exponent of the orhobaric isochoric pressure coefficients
φi
density ratio of component i, φi = ρ00i /ρ0i 19
isothermal compressibility
λk
kth eigenvalue of Ψ, λ1 = min(λ1 , . . . , λN )
µi
chemical potential of component i
µ
chemical potential vector, µ = (µ1 , µ2 , . . . , µN )
ρ
total density, ρ = Vm−1
ρi
molar density, concentration of component i, ρi = xi /Vm
ρ
density vector, ρ = (ρ1 , ρ2 , . . . , ρN )
Ψ
Helmholtz energy density, Ψ = A/V
Ψ
Hessian matrix of Ψ
Ψij
matrix element of Ψ, Ψij = (∂ 2 Ψ/(∂ρi ∂ρj ))
Ψi
ith row vector of Ψ
M AN U
SC
κT
RI PT
ACCEPTED MANUSCRIPT
Subscripts c pure-fluid critical state m
molar property
xj l
curve along which xj is varied and the ratios of all mole fractions xi6=j are fixed
D
xj = 0 value attained if xj = 0
TE
Superscripts c property on a critical curve state residual property
χ
phase indicator
0
phase indicator; along isopleths, the phase with fixed composition
00
phase indicator; along isopleths, the phase coexisting with the 0 phase
∞
infinite dilution
AC C
EP
r
References
¨ ¨ [1] D. Konowalow, Uber die Dampfspannungen der Flussigkeitsgemische, Ann. Phys. 250 (1881) 34–52, volume also cited as Wied. Ann. 14. ¨ ¨ [2] D. Konowalow, Uber die Dampfspannung von gemischten Flussigkeiten, Ann. Phys. 250 (1881) 219–226, volume also cited as Wied. Ann. 14. [3] J. W. Gibbs, On the equilibrium of heterogeneous substances, in: R. G. van Name, H. A. Bumstead (Eds.), The Scientific Papers of J. Willard Gibbs, Vol. 1: Thermodynamics, Dover Publications, New York, 1961, reprint of a book of 1906 containing articles originally published in 1876–1878. 20
ACCEPTED MANUSCRIPT
[4] R. Haase, Thermodynamik der Mischphasen, Springer, Berlin, 1956. [5] U. K. Deiters, Differential equations for the calculation of fluid phase equilibria, Fluid Phase Equilib. 428 (2016) 164–173. doi:10.1016/j.fluid.2016.04.014.
RI PT
[6] R. J. Sadus, Molecular simulation of the phase behavior of fluids and fluid mixtures using the synthetic method, J. Chem. Phys. 137 (2012) 054507:1–7. doi:10.1063/1.4739853. [7] D. Y. Peng, D. B. Robinson, A new two-constant equation of state, Ind. Eng. Chem. Fundam. 15 (1976) 59–64. [8] D. B. Robinson, D. Y. Peng, The characterization of the heptanes and heavier fractions for the GPA Peng–Robinson program, GPA Research Report RR-28 (1978) 1–36.
SC
˜ [9] S. E. Quinones-Cisneros, U. K. Deiters, An efficient algorithm for the calculation of phase envelopes of fluid mixtures, Fluid Phase Equilib. 329 (2012) 22–31. doi:10.1016/j.fluid.2012.05.023.
M AN U
[10] U. K. Deiters, T. Kraska, High-Pressure Fluid Phase Equilibria—Phenomenology and Computation, Vol. 2 of Supercritical Fluid Science and Technology (E. Kiran, ed.), Elsevier, Amsterdam, 2012. [11] J.-N. Jaubert, R. Privat, F. Mutelet, Predicting the phase equilibria of synthetic petroleum fluids with the PPR78 approach, AIChE J. 56 (2010) 3225–3235. doi:10.1002/aic.12232.
D
[12] R. Privat, J.-N. Jaubert, Classification of global fluid phase equilibrium behaviors in binary systems, Chem. Eng. Res. Des. 91 (2013) 1807–1839.
TE
[13] R. A. Heidemann, A. M. Khalil, The calculation of critical points, AIChE J. 26 (1980) 769– 779. [14] M. L. Michelsen, Calculation of critical points and phase boundaries in the critical region, Fluid Phase Equilib. 16 (1984) 57–76.
EP
[15] J. R. Cash, A. H. Karp, A variable-order Runge–Kutta method for initial-value problems with rapidly varying right-hand sides, ACM Trans. Math. Software 16 (1990) 201–222. doi:10.1145/79505.79507.
AC C
¨ Ingenieure, Bibliographisches [16] G. Jordan-Engeln, F. Reutter, Numerische Mathematik fur Institut, Mannheim, 1985. [17] G. Venkatarathnam, Density marching method for calculating phase envelopes, Ind. Eng. Chem. Res. 53 (2014) 3723–3730. doi:10.1021/ie403633d. [18] U. K. Deiters, ThermoC project homepage: http://thermoc.uni-koeln.de/index.html, accessed on 2016-02-17. [19] H.-J. Kretzschmar, T. Zschunke, J. Klinger, A. Dittmann, An alternative method for the numerical calculation of the Maxwell criterion in vapour pressure computations, in: ˇ M. P´ıchal, O. Sifner (Eds.), Proceedings of the 11th International Conference on Properties of Water and Steam, Hemisphere, New York, 1990, pp. 324–331. [20] A. E. Elhassan, R. J. B. Craven, K. M. de Reuck, The Area method for pure fluids and an analysis of the two-phase region, Fluid Phase Equilib. 130 (1997) 167–187.
21
ACCEPTED MANUSCRIPT
AC C
EP
TE
D
M AN U
SC
RI PT
[21] W. H. Press, S. A. Teukolsky, W. T. Vetterling, B. P. Flannery, Numerical Recipes in C, Cambridge University Press, Cambridge, 2002.
22
ACCEPTED MANUSCRIPT
Appendix A
Derivatives at a mixture critical point under isothermal conditions
For completeness’ sake we briefly summarize equations for the calculation of the density derivatives at a critical point in an isothermal phase diagram.
c
RI PT
• For a binary mixture, a critical point is always an extremum of the isothermal phase diagram, hence dρ0 dρ00 (58) = →∞. dp dp c
M AN U
SC
• For a ternary mixture, the derivatives can be obtained in a similar fashion as in Section 2.4. To describe the phase boundary at constant temperature and under variation of x0j (keeping the ratio of the mole fractions of the other components constant), we can set dρc f dρ0 = − u1 , (59) dp 0 dp 2 xj l
T
where dρc /dp|T is the slope of the critical curve at constant temperature and u1 the normalized critical eigenvector. The scalar f must be adjusted in such a way that the ratio of the components of dρ0i6=j /dp|x0j l has the correct value,
x0j l
1 dρ0k − 0 xk,0 dp
D
1 dρ0i x0i,0 dp
= 0,
i, k 6= j .
(60)
x0j l
EP
TE
Together, Eqs. (59) and (60) constitute a system of linear equations for f and the components of dρ0 /dp. The derivative vector of the coexisting phase (00 ) is then dρ00 dρc f = + u1 . (61) dp 0 dp 2 xj l
T
AC C
For example, if the isothermal phase envelope is calculated with the constraints T = const, x01 = const, Eq. (59) becomes dρc f dρ0 = − u1 (62) dp 0 dp 2 x1 l
T
and Eq. (60)
1 x02,x1 =0
dρ02 dp
− x01 l
1 x03,x1 =0
dρ03 dp
=0.
(63)
x01 l
• For a quaternary or higher mixture, Eq. (60) can be fulfilled by writing the derivative
23
ACCEPTED MANUSCRIPT
vector of the 0 phase as dρ0 dp
and
(64)
SC
with
0 . .. 0 ej ≡ 1 0 .. . 0
RI PT
dρ0j ej + gx0xj =0 dp x0j l x01,xj =0 .. . 0 xj−1,xj =0 , x0xj =0 ≡ 0 x0 j+1,xj =0 .. . 0 xN,xj =0 =
M AN U
where g is a yet to be determined scalar. There exist N critical curves p → ρc |T,xj l . Hence Eq. (59) becomes N X dρ0j dρc f 0 ck ej + gxxj =0 − = − u1 dp dp 2 k=1
xk l
N X
ck = 1
(65)
k=1
D
N X k=1
dT c ck dp
=0.
xk l
AC C
EP
TE
The last equation, which contains the derivatives of the temperature along the critical curves, implements the T = const condition. Eq. (65) is a set of linear equations for the unknowns dρ0j /dp, g, and c. As in Section 2.4, the physically reasonable solution is found by minimizing c2 with respect to f . Finally, the derivative vector for the 00 phase is obtained from N dρ00 X dρc f = ck + u1 . (66) dp dp 2 k=1
xk l
24
ACCEPTED MANUSCRIPT
Appendix B
Postiteration
If desired, the accuracy of isopleth calculations with the differential equation method can be improved by postiteration. The system of equations to be solved is µ00i (ρ00 , T ) = µ0i (ρ0 , T ) 00
0
i = 1, . . . , N
(67)
0
p (ρ , T ) = p (ρ , T ) . The vector of residuals is therefore
N p00
−
N p0
(68)
SC
µ001 − µ01 .. . y= . 00 0 µ −µ
RI PT
00
M AN U
As ρ0 = ρ0 x0 and x0 is fixed during an isopleth calculation, the vector of unknowns can be defined as 00 ρ1 .. z= . . (69) 00 ρ N
ρ0
D
A single step of a Newton–Raphson iteration can then be described by z := z − δz ,
TE
(70)
where the correction vector δz is obtained by solving the linear system of equations (71)
EP
B δz = y
AC C
and the matrix B = {Bij } is defined by 00 ∂µi = Ψ00ij ∂ρ00 00 j ρ ,T 0 k6=j 0 0 i − ∂µ ∂ρ0 x0 ,T = −Ψi · x Bij = 00 ∂p 00 00 ∂ρ00j ρ ,T = Ψj · ρ k6 = j − ∂p0 = −ρ0 x0 · Ψ0 x0 ∂ρ0
x0 ,T
i = 1, . . . N ∧ j = 1, . . . N i = 1, . . . , N ∧ j = N + 1 . i = N + 1 ∧ j = 1, . . . , N
(72)
i=N +1 ∧ j =N +1
Usually a single iteration step is sufficient. As the derivatives required for the postiteration also appear in the differential equations, the additional programming effort is negligible.
25
ACCEPTED MANUSCRIPT
Highlights: • differential equations of isoplethic phase envelopes • applicable to multicomponent systems • can be rapidly integrated; avoids convergence problems or internal volume calculations
AC C
EP
TE D
M AN U
SC
RI PT
• for pure fluids: density-based Clapeyron equation (gives orthobaric densities)
1