DISTILLATION SCIENCE

Larry Coleman
Consultant on Demand
Distillation Science

Larry Coleman
(unable to fetch text document from uri [status: 403 (Forbidden)])
TABLE OF CONTENTS
This is a ten-part series of technical articles on Distillation Science, as is currently practiced on an industrial level. It is organized as it
would be used for the process design of large-scale distillation columns, which may differ from the way topics are introduced to
students academically. As a part of distillation column design, these articles also involve the vapor-liquid equilibria (VLE) of binary
systems, with potential extension to multi-component systems.

1: OVERVIEW
While chlorosilanes and their electronic impurities are used as recurring examples, the technology developed here has great
application in other mildly polar and hydrogenated compounds, which are normally excluded from theoretical treatment in many
texts. The emphasis is on industrial level applications.

2: VAPOR PRESSURE
This article deals with the pure component VP relationships commonly found in textbooks, as well as their limitations. This article
lays the basis for subsequent articles on improving VP equations and tying in Equations of State.

3: CRITICAL PROPERTIES AND ACENTRIC FACTOR

This article deals with the tabulation of these properties for selected fluids based on globally collected data. Data analysis and
validation are discussed, as well as estimation techniques for those fluids for which data is either poor or non-existent. Critical
properties are used to convert temperature, pressure, and specific volume from conventional units used in both chemistry and
chemical engineering, to the reduced form.

4: NEW VAPOR PRESSURE EQUATION

This article shows how a new vapor pressure equation form allows practice of distillation applications at the elevated pressures more
common to industry. The culmination of this article is a thermodynamically consistent equation form that is valid between the
atmospheric boiling point and the critical point, and which allows evaluation of other required distillation properties such as saturated
phase densities and latent heat of vaporization.

5: EQUATION OF STATE
This article deals with the recommendations for Equation Of State (EOS), as well as the mathematical techniques for solving such
non-intrinsic EOS equation forms. The module explains why an EOS is needed to evaluate the pure-component physical properties
that are used in distillation science, reviews the various options and makes a recommendation. Also included are techniques for
solving the EOS cubic equations and some work-arounds near the critical point.

6: FUGACITY
This article deals with the one of the two possible departures from ideal pure-component vapor pressure in binary mixtures, as is
commonly found in practical application of distillation science.

7: LIQUID ACTIVITY COEFFICIENTS

This article deals with Liquid Activity Coefficients, the second type of departure from pure-component vapor pressure. It describes
the application and estimation of Liquid Activity Coefficients, as commonly found in practical application of distillation science.
Various activity coefficient models are reviewed along with some limited data and a recommended estimation technique given where
data is lacking.

8: VLE ANALYSIS METHODS

This article discusses the recommended methodology used when data-collecting binary systems, in order to assure that systemic errors
are minimized. This topic also deals with validation whenever data collection is done on reactive fluids, which can disproportionate or
dimerize during study.

9: PUTTING IT ALL TOGETHER

This article shows how to combine the components of the above Distillation Science topics in a practical application, as well as
illustrating why it is necessary to include departures from ideal behavior in real binary systems, as commonly encountered in
industrial practice.

10: CONVERGENCE STRATEGY

This last article shows how VP equation solutions are best obtained for additional fluids in non-intrinsic or nested-loop equation
forms, such as the recommended vapor pressure equations. While this topic is more mathematical or computer-science oriented than
chemistry, it is a necessary technique to understand when dealing with modern technology.

BACK MATTER

1 11/7/2021
INDEX
GLOSSARY

2 11/7/2021
1: Overview
This is Part I, Overview of a ten-part series of technical articles on Distillation Science, as is currently practiced on an
industrial level. It is organized as it would be used for the industrial process design of distillation columns, which may differ
from the way it is introduced to students academically. As a part of distillation column design, these articles also involve the
vapor-liquid equilibria (VLE) of binary systems, with potential extension to multi-component systems.
It is assumed that the reader has a basic understanding of distillation, as a means of separating two or more volatile fluids by
the difference between their boiling points; and also the purpose of a distillation column's various components. A variety of
texts exist that explain these concepts on a basic level. From such pre-requisites, the purpose of this series of articles is to
expand these basics to the degree necessary for commercial applications with real fluids - which typically do not follow ideal
behavior. It is also assumed that the reader is familiar with the concepts of molar units, mole fractions, the bonding of chemical
elements, vapor pressure, latent heat of vaporization, and a compound's critical point - all of which are found in college
freshman-level textbooks. In discussing mathematical relationships, it is assumed that the concepts of differentiation and
integration are known to the reader.
While chlorosilanes and their electronic impurities are used as recurring examples, the technology developed here has great
application in other mildly polar and hydrogenated compounds, which are normally excluded from theoretical treatment in
many texts. The chlorosilane homologue starts with silane (SiH4), and concludes with silicon tetrachloride (SiCl4), as the Si-H
bonds are incrementally swapped for Si-Cl bonds.
Chlorosilanes (and many of the electronic impurities) are non-naturally occurring compounds, but are key to the manufacture
of high-purity silicon (solar photovoltaics and modern electronic integrated circuits, aka “computer chips”); as well as silicon-
based chemicals such as silicones and organic/inorganic coupling agents. Applications for distillation science are found in both
bulk separations of these fluids, as well as the purification by high-reflux processing, to parts-per-billion levels.
Furthermore, this technology is applicable to a wide variety of other non-organic fluids such as refrigerants, biologic pre-
cursors, and pharmaceutical pre-cursors (i.e., not the final biological or pharmaceutical product, but rather the compounds used
as catalysts and building-block agents for assembling these specialized chemical products).
To make the discussion of Distillation Science more generic, temperature, pressure, and molar volume ( T, P, and V) are often
expressed in dimensionless reduced units, designated respectively as Tr , Pr , and Vr. Tr = T/Tc ; Pr=P/Pc ; Vr = V/Vc . By
definition, Tr =1, Pr =1, and Vr =1 simultaneously at the critical point. When temperature and pressure are given as T and P,
these are in absolute units of Kelvins and atmospheres. For those more familiar with pressures in other absolute units such as
PSIA, bar absolute, and Pascals absolute, converting those units to absolute atmospheres is easily done using internet-
accessible conversion tables and calculators.
Part II, Existing Vapor Pressure Equations deals with the pure component VP relationships commonly found in
textbooks, as well as their limitations. This article lays the basis for subsequent articles.
Part III, Critical Properties and Acentric Factor deals with the tabulation of these properties for selected fluids based
on globally collected data. Data analysis and validation are discussed, as well as estimation techniques for those fluids for
which data is either poor or non-existent. Critical properties are used to convert temperature, pressure, and specific volume
from conventional units used in both chemistry and chemical engineering, to the reduced form of Part IV. In Parts V and
VII, the value of the acentric factor is important.
Part IV, New Vapor Pressure Equation expands on the basics of Part II and the results of Part III, to show how a new
vapor pressure equation allows the practice of distillation applications at the elevated pressures more common to industry.
The culmination of this article is a thermodynamically consistent equation that is valid between the atmospheric boiling
point and the critical point, and which allows evaluation of other required distillation properties such as saturated phase
densities and latent heat of vaporization.
Part V, Equation of State deals with the recommendations for best EOS, as well as the mathematical techniques for
solving such non-intrinsic EOS equation forms.
Part VI, Fugacity deals with the departure of apparent pure-component vapor pressure in binary mixtures, as is commonly
found in practical application of distillation science. The equations for evaluating fugacity coefficents are given.
Part VII, Liquid Activity Coefficients deals with the application and estimation of Liquid Activity Coefficients, as
commonly found in practical application of distillation science. Various Liquid Activity Coefficient models are reviewed

Larry Coleman 1.1 9/26/2021 [Link]

 along with some limited data and a recommended estimation technique given where data is lacking.
Part VIII, VLE Analysis Methods discusses the recommended methodology used when data-collecting binary systems,
to assure that systemic errors are minimized. This topic also deals with validation whenever data collection is done on
reactive fluids, which can disproportionate or dimerize during study.
Part IX, Putting It All Together shows how to combine the components of the above distillation science articles in a
practical application.
Part X, Convergence Strategy illustrates how solutions are best obtained for additional fluids using non-intrinsic or
nested-loop equation methods, such as the recommended vapor pressure equations. While this topic is more mathematical
or computer-science oriented than chemistry-centered, it is a necessary technique to understand when dealing with modern
technology. No complex mathematics are required.
Equation Notation Used
Notation Usage Units Notation Usage Units

NBP Normal Boiling Point °K q Watson exponent -

Pressure, Vapor Latent heat of

P, VP atmospheres ΔHvap Cal/g-mol
Pressure Vaporization
Pc Critical Pressure atmospheres ΔHvb Latent heat @ Tb Cal/g-mol
Thek-Stiel VP
Pr Reduced Pressure - T-S ΔHvb Cal/g-mol
equation parameter
R Gas constant Cal/g-mol°K Z compressibility -
compressibility of
T Temperature °K Zv & ZL saturated vapor & -
liquid
Temperature @ (Zv-ZL) @
Tb °K ΔZvap -
Normal Boiling Point vaporization
Tc Critical Temp °K α Riedel derivative -
Tr Reduced Temp - αc α @ critical -
Vapor Pressure
Tbr Reduced Tb - ψ -
derivative
Vc Critical Volume cc/g-mol ø fugacity coefficient -
liquid activity
k H bonding parameter - γ -
coefficient
Binary Interaction
kij - ω acentric factor -
Coefficient

Units
The units used in these articles are temperatures in °K, pressures in atmospheres (absolute), and molar volumes in cc/gram-
mole. Molecular weight (MW), and compressibility (Z) are dimensionless. For the above units, the value of the gas constant
(R) is 82.057 atm-cc/mole-°K.

Larry Coleman 1.2 9/26/2021 [Link]

2: Vapor Pressure
Distillation Science (a blend of Chemistry and Chemical Engineering)
This is Part II, Vapor Pressure of a ten-part series of technical articles on Distillation Science, as is currently practiced on
an industrial level. See also Part I, Overview for introductory comments, the scope of the article series, and nomenclature.
Part II, Vapor Pressure deals with the existing pure component vapor pressure (VP) equations normally found in
textbooks. The content of this article is referred to in subsequent articles. The goal of this article is to explain the
limitations of these less sophisticated VP equations, to show how best to use them, and to set up the introduction of a
superior VP equation form in Part IV.
The original Clausius-Clapeyron equation relating VP, temperature, vaporization molar volume change (ΔVv) and latent
heat of vaporization (ΔHv) dates back to mid-19th century and is derived from thermodynamic principles. The derivation is
given in many college freshman-level texts and results from thermodynamic equilibrium between liquid and vapor phases.
First, the differential of pressure with respect to system temperature (for both vapor and liquid) is re-arranged to:
dP ΔHv
= (2.1)
dT T ΔVvap

If (assumption 1) the vaporization molar volume change (ΔVvap) is set equal to the saturated vapor volume (V) by
assuming the boiling liquid's molar volume is essentially zero; and if (assumption 2) the "Ideal Gas Law" (PV=RT) can
hold for this saturated vapor, then Equation 2.1 in differential form the relationship becomes:
d(ln P ) ΔHvap
= (2.2)
d(1/T ) RT

When integrated, this results in the simplified Clausius-Clapeyron relationship that is found in most textbooks. However,
there are limited conditions for which the above two assumptions are nearly correct: non-complex molecular structure and
very low pressure (say, below atmospheric or very near atmospheric). For most compounds, and for most pressures
encountered in normal industrial processes, neither of these assumptions holds very well and accuracy gets increasingly
worse as pressure increases beyond atmospheric.
For real compounds and pressures normally encountered in industry, a term is added to the "Ideal Gas Law" (P V = RT )
called compressibility, Z ; so then the relationship becomes:
P V = ZRT (2.3)

where compressibility Z is a function of that fluid's pressure, temperature, and other physical properties such as discussed
in the Part III. Note that Equation (2.3) holds for both vapors and liquids: with Zv being vapor compressibility, ZL being
liquid compressibility, and ΔZvap being the change in compressibility with vaporization. Now Equation (2.1) can be
transformed into a more usable relationship for all fluids and all sub-critical pressures. In differential form, it is:
d(ln P ) ΔHvap
= (2.4)
d(1/T ) ΔZvap RT

After integration between close temperatures T1 and T2 (so the ratio ΔHvap / ΔZvap can be taken as constant over that close
range), Equation (2.4) becomes:
ΔHvap
ln(P2 / P1 ) = × (1/ T1 − 1/ T2 ) (2.5)
ΔZvap R

In order for this equation to be fairly accurate, it important for T and T to be close, since both ΔHvap and ΔZvap are
1 2

actually functions of temperature. Also note that if ΔZvap is set to unity, Equation (2.5) becomes the same as the simplified
integrated Clausius-Clapeyron equation:
ΔHvap
ln(P2 / P1 ) = × (1/ T1 − 1/ T2 ) (2.6)
R

Larry Coleman 2.1 10/31/2021 [Link]

It can now be understood that ΔZvap is the measure of the net effect of removing above assumptions (1) and (2) that
allowed Equation (2.2). For most fluids near atmospheric pressure, ΔZvap is about 0.95 - 0.99, but as pressures increase
toward a fluid's critical point (where vapor and liquid merge at a singularity), ΔZvap goes to zero. ( But of course, ΔHvap
also goes to zero at the critical point, so the ratio ΔHvap / ΔZvap becomes undefined at this singularity).
The evaluation of ΔZvap normally requires use of an Equation of State (EOS), which is discussed in Part V. Using such an
EOS to assess ΔZvap would allow vapor pressure vs temperature data to exactly align with ΔHvap values. In Part IV,
Equation (2.4) will be further developed and integrated without use of an EOS, up to pressures close to critical (usually
about 95% of critical pressure in pure component systems).
Returning back to the above integrated form of Clausius-Clapeyron, Equation (2.5) suggests the commonly used empirical
means to curve-fit VP data over small temperature/pressure ranges:
Ln(P ) = A − B/T (2.7)

The term "B" does not have any exact scientific significance, but works as a curve-fitting parameter and is normally shown
with a negative sign, so as to have a positive value. Since Equation (2.7) has a limited range of application to low
pressures, it can be improved for use near ambient and slightly higher pressures with an empirical form for curve-fitting
VP vs T, called the Antoine Equation. However the constants A, B, and C have no scientific basis either and the equation
form can still only be used over modest ranges ( and low pressures).
B
ln(P ) = A − (2.8)
T +C

Additional constants D, E, F can be added to make the “extended Antoine” relationship for empirical curve-fitting;
however there is still no correlation between the constants and scientific meaning.
B
2
ln(P ) = A − +D×T +E ×T + F × Ln(T ) (2.9)
T +C

From an industrial distillation column design perspective, even the extended Antoine VP relationships found in handbooks
are not entirely adequate: they rarely reproduce the correct VP values at Tb, Tc and at Tr=0.7 (where the acentric factor is
determined). That inadequacy undermines attempts to use modern Equations of State (EOS are discussed in Part V),
forcing the use of the antiquated Van der Waals EOS. With any EOS, the extended Antoine VP relationship cannot connect
latent heats and saturated (vapor and liquid) phase densities to vapor pressure. Additionally, none of these empirical VP
equations reproduces the inflection point required by thermodynamics of a real fluid's VP vs T plot, normally occurring
between Tr= 0.7 and Tr = 0.85. For economic reasons, many industrial processes operate at these higher pressures, so using
empirical VP equations can lead to poor distillation column design. This is especially true in modern industrial
applications where complex computer programs have automated several aspects of distillation column design.
In designing such an industrial-level distillation column, the various process simulation packages (e.g., ASPEN, VMG,
etc) would then be fed with rather poor VP estimations. It is not uncommon with such simulation software to become non-
convergent or to have the problem solution get stuck on a singularity. In the vernacular of the early days of computing,
“Garbage in – Garbage out”.
The solution to this quandary is to find a better VP equation form that possesses all the criteria lacking in the VP equations
shown in this article. A solution to that is proposed in Part IV.
However, Equation (2.7) and Equation (2.8) are not without merit, as long as they are used for the purpose they were
derived: only for narrow ranges of temperature and at moderate pressures. Equation (2.8) (Antoine) is especially useful in
correlating data taken around the atmospheric boiling point, Tb, to get a more accurate value from several experimental
data points, rather than just one. It also allows data from several sources to be inter-compared, as long as they had a similar
temperature range.

While evaluating the “A, B, and C” of Equation (2.8) may seem daunting, the solution is readily managed using some
algebraic manipulation and a spreadsheet multiple regression function (e.g., MS Excel). Equation (2.8) is re-written as the
algebraically equivalent

Larry Coleman 2.2 10/31/2021 [Link]

 AC − B
ln(P ) = A + − C × Ln(P )/T (2.10)
T

Then VP vs T data (always as absolute pressure and temperature) are converted into three columns for each data set: the
dependent variable is “Ln(P)”, and the two independent variables are “1/T” and “Ln(P)/T”. When the regression is run, the
regression's intercept will equal Equation (2.8)’s “A”; and the second independent variable’s coefficient (i.e., for Ln(P)/T )
will equal the negative of Equation (2.8)’s “C”. The value of Equation (2.8)’s “B” is determined algebraically from the
first independent variable’s coefficient (i.e., for 1/T)= (AC-B). See an example of using this solution procedure in below
Table 2-1.

This curve-fit parameter solution technique does not work with Equation (2.9), so often the “Extended Antoine” equation
is stated differently as
B 2
 Ln(P ) = A − + C × Ln(T ) + DT + E T (2.11)
T

requiring a multiple linear regression (such as with MS Excel) be done with the four “independent” variables: 1/T, Ln(T),
T, and T2. Assuming the temperature range is narrow, the curve-fit quality is almost always statistically better with
Equation (2.8) than with Equation (2.11), since there are two fewer constants. If the range is broader, such as spanning
from atmospheric to the half-way point of critical pressure, then Equation (2.11) will give the better curve-fit.
When using experimental data to determine the normal boiling point (Tb), Equation (2.8) seems to work best. When using
experimental data to determine the VP at Tr= 0.7 (i.e., evaluating the acentric factor), or inter-comparing VP measurements
over a similar broader range, Equation (2.11) is preferred. Equation (2.9) is rarely used for design or data comparison, and
is just mentioned for historical purposes.

Example 2.1

PCl3 is an important impurity to remove in producing high quality polysilicon for use in solar arrays and electronic
integrated circuits ( like computer chips). Part of the purification process involves distillation, and so it is desirable to
know the NBP of this compound with high confidence. Researching the NIST database shows an Antoine expression (
Equation (2.8) ), but the constants given are recalculated from a 1947 paper by Dan Stull published in I&EC. In
making that paper's VP tables, he took data from several early-to-mid 20th century sources and plotted them on a Cox
Chart (a graphical approximation method from the 1920's), then read "best fit" values from the Cox Chart at selected
pressures. So this data source is questionable, with a possibility of data being "overly massaged". However from
DeChema and Infotherm online global databases, experimental data is available from four more recent sources, with
some mild disagreement in the pressure range of 1.0 atmospheres = 760 Torr = 101.325 kPa. It is decided to download
the original data and develop an Antoine curve-fit, to best determine a most-likely NBP for PCl3.
Table 2-1 shows the experimental data, sorted by temperature. DeChema and Infotherm report vapor pressures in kPa,
so that pressure unit is used in the below calculation. To each data point row, columns for "1/T" and "Ln(P)/T" are
added. The last two columns show the VP predicted by the Antoine regression and the difference between data point
and predicted VP. Under the tabulated data, the results of the regression are shown, as well as the calculation of
Antoine constants A, B, and C. And finally the best value for PCl3's NBP based on data along with an assessment of
accuracy.
Table 2-1 Curve-fitting PCl3 VP data to the Antoine Equation, plus regression results
Ln(VP), Ln(VP)/T,
1/T, independent Antoine Prediction
T, K VP, kPa dependent independent
variable 1 predicted VP, kPa Difference
variable variable 2

361.85 151.988 5.024 2.7636E-03 1.3884E-02 150.362 1.626

361.75 151.988 5.024 2.7643E-03 1.3887E-02 149.964 2.024

358.55 136.389 4.916 2.7890E-03 1.3709E-02 137.544 -1.155
355.75 126.790 4.843 2.8110E-03 1.3612E-02 127.191 -0.402
352.95 116.524 4.758 2.8333E-03 1.3481E-02 117.318 -0.794

Larry Coleman 2.3 10/31/2021 [Link]

 Ln(VP), Ln(VP)/T,
1/T, independent Antoine Prediction
T, K VP, kPa dependent independent
variable 1 predicted VP, kPa Difference
variable variable 2

349.15 104.525 4.649 2.8641E-03 1.3316E-02 104.680 -0.155

348.55 101.325 4.618 2.8690E-03 1.3250E-02 102.764 -1.439
348.45 101.325 4.618 2.8699E-03 1.3250E-02 102.447 -1.122
347.95 101.325 4.618 2.8740E-03 1.3273E-02 100.870 0.455
346.15 95.192 4.556 2.8889E-03 1.3162E-02 95.317 -0.125
343.75 88.259 4.480 2.9091E-03 1.3034E-02 88.213 0.046
343.15 86.526 4.460 2.9142E-03 1.2999E-02 86.490 0.036
338.15 73.461 4.297 2.9573E-03 1.2707E-02 72.952 0.508
337.55 72.127 4.278 2.9625E-03 1.2675E-02 71.425 0.703

The regression results using MS Excel's data analysis function show an R2

= 0.9971 with 14 data points
Intercept = 9.7954; Variable 1 coefficient = -2640.04; Variable 2 coefficient
=181.437
Calculate Antoine "A", "B", and "C" per the above procedure:
Intercept = 9.7954 = Antoine "A"
Variable 1 coeff. = -2640.04 so (- Intercept x Variable 2 coeff. - Variable 1 coeff. ) = Antoine "B" = 862.795
Variable 2 coeff. = 181.437 = -1 x Antoine "C", so Antoine "C" = -181.437

The Antoine Equation for PCl3 is determined as

862.795
ln(P ) = 9.7954 − (2.12)
T − 181.437

Solving Equation (2.12) for the atmospheric pressure of 101.325 kPa gives a Tb of 348.09°K = 74.94°C ⇒ 74.9°C .
The average absolute difference of ( data -predicted) VP = the table's last column is 0.756 kPa
Note that if the NIST webbook values were used for "A", "B" and "C", a value of 348.34°K ⇒ 75.2°C would result. In
this case, the NIST listed results were pretty close the experimental values' regression. Also note that a quick search on
Wikipedia's website would have offered a Tb value of 76.1°C, or about a degree higher than actual results. But now
that the best value of Tb is known, and there can also be a reasonable evaluation of value's accuracy based on the
regression's average error. A 0.756 kPa error at 101.325 kPa pressure equates to a change in calculated Tb of 0.24°K,
so the scientific answer to the question is 74.9+ 0.2°C.

Larry Coleman 2.4 10/31/2021 [Link]

3: Critical Properties and Acentric Factor
Distillation Science (a blend of Chemistry and Chemical Engineering)
This is Part III, Critical Properties and Acentric Factor of a ten-part series of technical articles on Distillation Science,
as is currently practiced on an industrial level. See also Part I, Overview for introductory comments, the scope of the
article series, and nomenclature. The goal of this article is to explain how the critical properties and acentric factor are
used, and how they may be estimated when reliable values are not available.
As discussed in Part II, the more basic vapor pressure equations have limitations that prevent them from being used in
many modern industrial distillation systems, specifically at higher pressures. The concepts in this article will be used in
Part IV, to introduce an improved vapor pressure equation, which removes the limitations of basic vapor pressure
equations.
To implement the more sophisticated VP relationships and inter-compare the constants in Part IV, experimentally
measured P vs T data needs to be put into reduced form. This is done only with the knowledge of the critical point, which
is that singularity on the saturation line where the liquid and vapor phases become one. The critical temperature and
pressure are identified as Tc and Pc. At the critical point, the specific molar volume is Vc.

Tr = T / Tc (3.1)

Pr = P / Pc

Vr = V / Vc

(note that reduced properties have the advantage of being dimensionless)

Using reduced variables Tr, Pr and Vr from Equation (3.1) also allows introduction to the concept of the "Law of
Corresponding States" ( which is really not a scientific law per se, but more of a generally followed relationship). This
"law" expresses the generalization that those properties dependent on intermolecular forces are related to the critical
properties in the same way for all fluids. This concept underlies the development of several later articles of this series on
Distillation Science, including Parts IV through VII.
The critical compressibility is Zc, as defined as:
Pc Vc
Zc = (3.2)
RTc

which is also dimensionless.

Two other properties are important to note:
Tb is defined as the atmospheric pressure boiling point of a fluid, with T as the reduced value of the atmospheric
br

boiling point
ω is the acentric factor, which is a measure of molecular complexity as defined as:

ω = − log10 (V P ) − 1 (3.3)

where the VP (in atmospheres) is evaluated at Tr= 0.7

The term “acentric factor” comes from Kenneth Pfizer's observation in 1955 that compact and near-spherical molecules
have close to zero values when their ω is calculated. For example, neon, argon and krypton have ω values of -0.04, 0.00
and 0.00, respectively (helium has a more slightly negative ω, the value depending on isotope). Methane (CH4) has a ω of
0.011; and the molecularly larger silane (SiH4) has an ω of 0.099. An even larger and more complex molecule, like carbon
tetrachloride (CCl4) has an ω of 0.193; whereas silicon tetrachloride (SiCl4) has an ω of 0.248. So valuation of the acentric
factor from VP data and knowledge of the critical point is a good validation check against the known molecular size and
shape.
For Equations of State (see Part V) that are more advanced than Van der Waals, the acentric factor is a required property,
along with critical temperature and pressure. It is also used along with critical properties in the estimation of binary
interaction parameters in Part VII.

Larry Coleman 3.1 10/31/2021 [Link]

As mentioned in Part II, the basic VP relationships have acceptable accuracy to obtain good values of Tb from available
data. Experimental data is often taken near atmospheric pressure, but not exactly at one atmosphere absolute (i.e, 760
mmHg = 760 Torr = 101.325 Pa). Instead of blindly accepting a value of Tb from a handbook or single website, the
practicing scientist or engineer should consider the data’s source and reliability. There is ready internet availability of
atmospheric boiling points, critical property data, and acentric factors, especially on the NIST WebBook and global
websites like Dechema, Infotherm (recently acquired by John Wiley and Sons from FIZ Chemie Berlin).
A recommended way to sort through the dizzying amount of critical property data (i.e., Tc , Pc , and Vc) is to use the
correlation techniques of Lyderson ("Estimation of Critical Properties of Organic Compounds", University of Wisconsin
College of Engineering, 1955) and others, which are specifically geared toward organic compounds. With some algebraic
manipulation and adjustment of Lyderson’s parameters (since the usage intent here is for polar compounds that are not
organic, but are not ionic either), and based on data across homologues, the following relationships seem to hold and are
recommended for both validating questionable data or filling in “holes” of no data:
For Tb, the general correlating relationship within a homologue ( e.g., from silane to silicon tetrachloride, or from
phosphine to trichlorophosphine) :
n
(Tb × M W ) = A + B × MW (3.4)

where MW is the molecular weight. For chlorosilanes, chloromethanes, chlorophosphines, etc, this relationship shows a
best fit with an exponent “n” of 0.8320. The constants A and B for the homologue are typically set by the hydride and
chloride, since that data is usually more reliable. Where the core atom is not just carbon or silicon, but a blend (i.e, a
methyl silane), an exponent of 0.8455 fits the data a bit better. This contrasts with organic compounds, whose data tends to
fit best with somewhat lower exponents that are closer to 0.80.
For Tc (critical temperature), the best correlation within a homologue is found to be:
Tc − Tb
= A + B × MW (3.5)
Tc

where MW is the molecular weight and the constants A and B for the homologue are typically set by the hydride and
chloride, as long as there is no organic content and the fluids are monomeric.
For Pc (critical temperature), use
MW
n
( ) = A + B × MW (3.6)
Pc

where MW is the molecular weight, exponent "n" = 0.5672, and constants A and B for the homologue are typically set by
the hydride and chloride (as long as there is no organic content and the fluids are monomeric). For a pure organic
compound, the exponent "n" should be Lyderson’s suggested 0.5000; for methyl silanes (e.g, dimethyl silane) the best
exponent is 0.4550; and for Group III (dimeric bridge-bonded) compounds the best exponent fit is 1.03-1.08 ( 1.03 for di-
gallanes and 1.08 for diboranes).
For Vc (critical molar volume), Lyderson’s rule of
Vc = A + B × M W (3.7)

is probably still the best, where MW is the molecular weight and constants A and B for the homologue are typically set by
the hydride and chloride. For some homologues, there is a possibility that the Vc vs MW relationship is not exactly linear,
but rather has a slight concave quadratic nature; but the data is rarely good enough to determine such. Of all the critical
properties, Vc is the hardest to experimentally measure, and is frequently estimated.
There is a way to get around some of the uncertainty with Vc, and that is to validate (or make minor Vc adjustments) based
on the pattern of Zc values in the homologue where
Pc × Vc
Zc = (3.8)
R × Tc

Within each homologue there is a characteristic “check-mark” pattern to Zc values, when plotted against the number of
hydrogen to chlorine swaps of the homologue (i.e., the homologue’s hydride has zero swaps, the monochloride one swap,

Larry Coleman 3.2 10/31/2021 [Link]

the dichloride two swaps, etc). See below example data plot Figure 3-1 for the chlorosilane and chlorophosphine
homologues.
The homologue’s hydride (which tends to be a fairly spherical shaped molecule with minimal dipole moment) typically has
a Zc in the range of 0.27-0.29. The first swap replaces a small “H” atom with a larger “Cl” atom, and re-orients the
molecular shape to have a significant dipole, dropping the Zc value by about 10%. Then the next swaps reduce the dipole,
until the molecular shape returns closer to spherical, although significantly larger in diameter. Typically the Zc of the
completely chlorinated homologue fluid has a slightly larger Zc than the homologue’s hydride.

Figure 3-1, Comparison of Zc pattern within two homologues

To the extent that this “checkmark” shape is not observed in plotting calculated Zc values per Equation 3-8, the “culprit” is
almost always the value of Vc, which allows for some correction. The value of Zc is known to have a significant
relationship to molecular shape and complexity, although the stronger relationship between Zc and the acentric factor, ω, is
not so easily generalized as in organic compounds.
The importance of accurate values of Tb , Tc , Pc (and ω to an extent) is in calculating vapor pressure using the advanced
relationship of Part IV. Having an accurate value of ω is more important, along with Tc, Pc in calculating the Equation of
State for a fluid, in Part V. Having good values for ω and Vc are important to the estimation of binary interaction
parameters in Part VII.
When two different types of swaps are made in a slightly polar molecule, an alternate technique is used to best correlate
fluid properties (Tb , Tc , Pc , Vc , and ω), with a good example being methyl chlorosilanes. With methyl chlorosilanes, one
or more methyl groups are substituted for hydrogen atoms (i.e., Si-H swapped to Si-CH3), with the remaining Si-H
bonding possibly swapped to Si-Cl. Another example would be swapping C-H for C-Cl and C-F bonding in chloro-fluoro
methanes. There is no concise formula that governs such a dual substitution. Instead the concept of neural networks is
used.
In this technique, the two base curves are laid out from the same starting compound (silane in the example below), and on
each base curve the pertinent physical property values shown for each substituted homologue (chlorosilanes and methyl
silanes being the base curves in the example below) vs molecular weight. For each amount of combined substitution
(100%, 75% and 50% substitution = 0%, 25%, and 50% Si-H remaining) the data points are connected to their respective
base curves. For example, on the 75% substitution amount, one would curve-connect trichlorosilane, methyl
dichlorosilane, dimethyl monochlorosilane, and trimethyl silane.
Additional curve-connections could be made for values of constant amount of chlorine or methyl substitution (e.g., curve-
connecting monochlorosilane, methyl monochlorosilane, dimethyl monochlorosilane, and trimethyl monochlorosilane).
Using this method, the data on the dual-substituted fluids can be validated and adjustments made so that the curve-
connections are smooth.
In the example in below Figure 3-2, the values of Pc for methyl monochlorosilane (sic, me-mono) and dimethyl
monochlorosilane (sic, di-me mono) are hardened up from rough measurements.

Larry Coleman 3.3 10/31/2021 [Link]

 Figure 3-2, Neural Network for methyl chlorosilanes’ Pcvalues
The value of the neural network technique is to make adjustments where data is poor or the estimation is of questionable
accuracy. The presumption of neural networks is that for a class of compound, the continuum of structural changes should
be uniform and internally consistent (although not necessarily linear).
Neural networks were constructed for all the pertinent properties of methyl chlorosilanes, and the technique found to be
rather simple to deal with graphically using MS Excel. The resulting values were included in Table 3-3 below, for
chlorosilanes and their common impurities, whose properties were validated using the correlating techniques of Equation (
3.4) through Equation (3.8).

Table 3-3, Key properties, arranged by Group and Homologue

Fluid MW Tb Tc Pc Vc Zc ω
Group III(A)
B2H6 27.670 180.54 289.70 39.58 173.10 0.2882 0.1254
B2H5Cl 62.119ǂ 216.37 346.33 37.13 189.92 0.2481 0.1283
BH2Cl
42.284 237.72 379.56 37.10 206.73 0.2463 0.1315
as monomer
BHCl2
82.722 266.07 422.71 37.63 240.37 0.2608 0.1389
as monomer
B2HCl5
199.893ǂ 276.67 438.47 37.92 257.18 0.2711 0.1433
hypothetic
BCl3 117.169 285.88 451.95 38.20 274.00 0.2822 0.1468
AlCl3
133.341 466.86* 625.70 26.00 261.80 0.1326 0.3474
as monomer
Ga2H4Cl2
214.383ǂ 350.86 545.93 40.94 231.9 0.2120 0.2726
as dimer
Ga2H2Cl4
283.273ǂ 421.79 639.69 41.02 247.5 0.1943 0.3525
as dimer
Ga2Cl6
352.162ǂ 473.49 694.00 37.70 263.0 0.1741 0.4504
as dimer
Group IV(A)
SiH4 32.117 161.75 269.65 47.99 130.07 0.2821 0.09860
SiH3Cl 66.562 242.75 396.65 47.82 169.32 0.2488 0.1252

Larry Coleman 3.4 10/31/2021 [Link]

 Fluid MW Tb Tc Pc Vc Zc ω
SiH2Cl2 101.007 281.45 449.45 44.83 215.05 0.2614 0.1589
SiHCl3 135.452 306.15 479.15 41.15 267.28 0.2797 0.2090
SiCl4 169.896 330.72 506.95 36.50 326.00 0.2860 0.2482
SiH3(CH3) 46.144 216.48 348.35 41.53 185.20 0.2691 0.1264
SiH2Cl(CH3) 80.589 278.82 439.11 41.16 233.90 0.2672 0.1793
SiHCl2(CH3) 115.034 314.17 489.18 39.25 287.81 0.2814 0.2269
SiCl3(CH3) 149.479 339.72 517.61 35.86 345.70 0.2919 0.2655
SiH2(CH3)2 60.169 252.86 399.17 36.05 245.84 0.2706 0.1604
SiHCl(CH3)2 94.615 305.30 471.35 35.94 300.32 0.2791 0.2264
SiCl2(CH3)2 129.061 342.89 519.21 34.26 358.20 0.2880 0.2740
GeH4 76.642 184.93 307.98 54.77 128.28 0.2780 0.1270
GeH3Cl 111.087 302.91 495.10 49.95 180.17 0.2215 0.1476
GeH2Cl2 145.532 334.60 536.93 45.36 236.71 0.2437 0.1721
GeHCl3 179.976 343.54 541.41 41.42 283.96 0.2647 0.2012
GeCl4 214.421 357.28 553.16 38.10 335.86 0.2819 0.2334
SnH4 122.742 221.07 360.20 51.70 152.90 0.2674 0.1619
SnCl4 260.521 387.21 591.85 36.95 351.20 0.2672 0.2625
Group V(A)
PH3 33.998 185.41 324.75 64.51 113.33 0.2743 0.03052
PH2Cl 68.443 273.09 463.73 62.29 153.63 0.2515 0.0750
PHCl2 102.888 318.44 524.74 55.82 202.52 0.2625 0.1362
PCl3 137.333 349.25 558.95 50.00 260.00 0.2834 0.2117
POCl3 153.331 379.00 605.21 47.59 276.00 0.2645 0.1993

Table 3-3, Continued

AsH3 77.945 210.73 373.00 65.12 132.50 0.2819 0.01341

AsH2Cl 112.390 300.41 515.99 64.33 174.85 0.2657 0.0573

AsHCl2 146.835 359.66 600.00 61.57 218.29 0.2730 0.1153

AsCl3 181.281 403.30 654.00 58.35 259.56 0.2822 0.1875

SbH3 124.781 256.09 446.20 66.61 157.20 0.2860 0.01659

SbCl3 228.115 794.05 794.05 68.85 268.00 0.2832 0.3309

n-pentane 72.149 309.16 470.05 32.86 310.0 0.2641 0.2393

iso-pentane 72.149 300.82 460.56 33.17 307.1 .2695 0.2147

In above Table 3-3, those entries that are marked with a ǂ have their molecular weight given as the dimer. For AlCl3, the
asterisk on the Tb entry is for the best correlating value for the Part IV VP equation, even though it is below the triple
point. Group IIIA compounds can feature “bridge-bonding” and can be either monomeric or dimeric. The two pentane
entries at the end of the table are included because of their significance as electronic impurities, and to show that the
techniques can be extended to organics.
With Group IIIA, some compounds are not included in a homologue simply because they are impossible structurally or are
unstable. An example would be B2H3Cl3, which would put too much strain on the B-B bridge-bond. Another is AlH3
which cannot exist as a stand-alone liquid with vapor pressure at normal processing conditions, but tends to form stable

Larry Coleman 3.5 10/31/2021 [Link]

hydrides with Group I metals such as LiAlH4. However, B2HCl5 is included (but noted as hypothetic), since some
researchers claim to have seen its presence in low levels of TCS.

Larry Coleman 3.6 10/31/2021 [Link]

4: New Vapor Pressure Equation
Distillation Science (a blend of Chemistry and Chemical Engineering)
This is Part IV, New Vapor Pressure Equation of a ten-part series of technical articles on Distillation Science, as is currently
practiced on an industrial level. See Part I, Overview for introductory comments, the scope of the article series, and nomenclature.
Parts II and III are pre-requisite to Part IV.
Previous Part II, Vapor Pressure deals with the pure component equations normally found in textbooks, but which have limitations in
industrial application. Part III, Critical Properties and Acentric Factor develops those parameters used to inter-relate conventional
units of temperature, pressure and volume to their reduced property equivalents, as well as the Pfizer acentric factor that is used in
modern Equations of State.
Several thermodynamically consistent vapor pressure (VP) equation forms have been developed in the past purposely for organic fluids,
and typically have the stipulation that they do not work with either polar or heavily hydrogen-bonded fluids. Being both slightly polar
in nature as well as hydrogen-bonded, chlorosilanes and their impurities fall exactly in this exclusion, and therefore make good
continuing examples.
In this article, temperature and pressure are used in their reduced form, since this allows better inter-comparison between fluids, is
dimensionless, and allows a more global future development. The reduced form is obtained by simply dividing conventional
temperature, T, by the fluid’s critical temperature=> Tr = T/Tc. Similarly, reduced pressure is conventional pressure divided by the
fluid’s critical pressure => Pr= P/Pc. Temperature and pressure must be in absolute units ( e.g., ºK not ºC)
Referring back to Part II, Equation 2-4, the Clausius-Clapeyron equation is given in reduced format, using ψ as notation for the
thermodynamic-based derivative of the natural log of (reduced) vapor pressure with respect to the inverse of (reduced) absolute
temperature:
−d(ln Pr ) ΔHv
ψ = = (4.1)
d(1/ Tr ) ΔZRTc

In order to have an improved VP equation form that can be used for the wide range between Tb and Tc (i.e., from atmospheric pressure
to the critical point), that relationship should meet the following criteria:
1. The equation form should exactly have a reduced VP value of 1/Pc at Tb/Tc (i.e., yields the known atmospheric boiling point at
atmospheric pressure), which is a well-measured lab value and is readily estimated based on structural considerations (see Part III
article).
2. The equation form should exactly have a reduced VP value of unity at the critical temperature (i.e., match the measured value of the
critical point), where the liquid meniscus disappears. This property can also be readily estimated based on structural considerations
(see Part III article).
3. Give a reasonable value for reduced VP at Tr=0.7, corresponding to the definition of the Pfizer acentric factor, which has a basis in
molecular complexity (see Part III article).
4. Meet the two thermodynamic consistency Reidel tests, as the critical point is approached ( that "α" monotonically drop to its lowest
value at critical, and that the temperature derivative of α = zero at critical). See the left-most part of Equation (4.4) for the definition
of "α".
Such a VP equation form would allow evaluation of modern Equations of State (see Part V article), Fugacity considerations (see Part
VI article), and Binary Interaction Parameters (see Part VII article). Preferably the VP equation form would reasonably fit most fluids
(both naturally occurring as well as synthetic), including polar compounds and those with substantial hydrogen bonding. Note that the
VP equations given in Part II work over only narrow ranges; however above criteria #1 and #2 require the VP equation form to span a
very broad range. In the text "The Properties of Gases and Liquids", by Prausnitz, et al, several alternate VP equation forms are
discussed, all of which meet the above four criteria.
In his doctoral thesis at Syracuse University in 1965, under the guidance of Leonard Stiel, Richard Thek proposed exactly such an
equation form to fill the technology gap. It was subsequently published by the AIChE in 1966. Unfortunately the computing power at
that time was very limited, so equation forms with iterative solutions could not be used; and the internet did not exist to allow ready
access to global data taken over the last many decades. Using modern computer processing speed and programming, these limitations
have been resolved, and the constants determined for all chlorosilanes and their electronic impurities.
As given by Equation (4.2) below, the original Thek-Stiel VP equation in its differential form has two terms: the first that governs at
lower pressures close to atmospheric, and the second that takes over as the critical point is approached. While Thek set variables “q”
and “k” equal to constants in order to avoid iterative solutions, with modern programming these can now be evaluated as true variables.
The derivation of the final Modified Thek-Stiel VP Equation is available on request, but is not reproduced in these articles for the sake
of brevity.

Larry Coleman 4.1 10/24/2021 [Link]

 d(ln Pr ) A c
2 3 4 n
= [1 − B1 Tr + B2 Tr − B3 Tr + B4 Tr . . . ] + (Tr − k) (4.2)
2 2
d(1/ Tr ) Tr Tr

The several “B” constants result from the binomial expansion of the general Watson relationship for ΔHv. After dropping the B4 and
higher expansion terms as non-significant mathematically and integrating, the resulting reduced vapor pressure equation is:
n−1
Tr −1
−1 2 −1
ln(Pr ) = A[ B0 − Tr − B1 ln Tr + B2 Tr − (1/2)B3 Tr ] + c [ + k(Tr − 1)] (4.3)
n−1

ΔHvb
A =
q
RTc (1 − Tbr )
1
B0 = 1 − B2 + ( )B3
2

B1 = q

q(q − 1) q(q − 1)
B2 = =
2! 2
q(q − 1)(q − 2) q(q − 1)(q − 2)
B3 = =
3! 6
αc − A(1 − B1 + B2 − B3 )
c =
1 −k
A
n = (1 − k) + (1 − B2 + 2 B3 )
c

The valuation of equation parameters ΔHvb, q, k, and αc are done iteratively using the techniques in Part X, Convergence Strategy,
based on vapor pressure data, the molecule’s structure, and the value of Tbr (the reduced atmospheric boiling point, Tb/Tc= Tbr). Pr must
be unity at the critical point (Tr= 1), which allows valuing the integration constant of Equation 4-3. The two important thermodynamic
derivative functions are:
d(ln Pr ) k
−1 2 n−1
= α = A(Tr − B1 + B2 Tr − B3 Tr ) + c(Tr − ) (4.4)
d(ln Tr ) Tr

−d(ln Pr ) ΔHv
2 3 n
= = ψ = A(1 − B1 Tr + B2 Tr − B3 Tr ) + c(Tr − k) (4.5)
d(1/ Tr ) ΔZRTc

The "A" and "B" constants of Equations (4.2), (4.3), and (4.4) are the same as Equation (4.5). The evaluation of the Reidel derivative
function "α" in Equation (4.4) is used to establish thermodynamic consistency, and allows valuing an equation parameter. Note that the
Clapeyron derivative of Equation (4.5), ψ , is mathematically identical to that in Equation (4.1). ψ will be used Part V to establish the
latent heat, ΔHv, as well as the saturated vapor and liquid density of the fluids at the various conditions in the distillation column
modeling.
Using the results of Part III, the Modified Thek-Stiel VP Equation constants are given below in Table 4-1, for chlorosilanes and
electronic impurities found in the manufacture of high purity silicon. For brevity, the values of acentric factor, ω, and the critical Riedel
derivative, α , are not given in the table. They can be calculated from the other equation constants. In some instances, an electronic
c

impurity can have either monomeric or dimeric form, but the form is noted in the table. Where a compound is not stable, the reader will
note a “hole” in the table (e.g., AlH3 and Ga2H6 are not stable as liquids or vapors).
Because the solution to the Modified Thek-Stiel VP Equation is necessarily iterative, the application of this VP relationship requires the
use of some modern computing tools to perform such iterative calculations (aka “nested loops”). Depending on preference, this could
be done as a macro in a spreadsheet program like MS Excel, or a stand-alone program developed in BASIC (or MatLab or Fortran). As
mentioned above, some tips on convergence strategy are given in Part X.
Table 4-1 is rather expansive, to illustrate how general the new VP relationship is. Table organization is primarily done by Periodic
Table Group (e.g., Group III, Group IV, and Group V) of the molecule’s core atom, and then by the homologue from hydride to
chloride. Methyl chlorosilanes are included as examples of good application to fluids that form the border between classically organic
and inorganic. Two of the pentanes are included, which happen to be of industrial concern as electronic impurities, to illustrate that the
new VP equation is useful for organic fluids as well.
While volatile Group II, Group VI and Group VII fluids seem to follow the new VP relationship as well, there are no entries given in
Table 4-1 since these fluids are generally not of concern in the production of electronic materials. The reader is encouraged to follow
the evaluation techniques and stratagems given in these articles to further develop use of the new VP equations.
Table 4-1, VP solution constants, by Group and Homologue

Fluid MW A B0 B1 B2 B3 c n k

Group III(A)

Larry Coleman 4.2 10/24/2021 [Link]

B2H6 27.670 8.9910 1.1494 0.37848 -0.11762 0.063573 2.8701 4.7796 0.1197

B2H5Cl 62.119ǂ 8.1915 1.1494 0.37867 -0.11764 0.063578 3.4560 3.8431 0.1073

BH2Cl as monomer 42.284 8.3872 1.1495 0.37888 -0.11766 0.063583 3.2767 4.0925 0.0938

BHCl2 as monomer 82.722 8.9910 1.1495 0.37937 -0.11772 0.063596 2.7577 4.9953 0.0636

BCl3 117.169 9.5652 1.1496 0.37988 -0.11779 0.063609 2.2357 6.2975 0.0291

AlCl3 as monomer 133.341 11.9183 1.1512 0.39300 -0.11928 0.063892 2.2700 7.5185 0.0291

Ga2H4Cl2
214.383ǂ 7.0714 1.1506 0.38810 -0.11874 0.063799 5.3346 2.5583 0.0938
as dimer

Ga2H2Cl4
283.273ǂ 10.8499 1.1513 0.39333 -0.11931 0.063898 3.3656 4.9567 0.0636
as dimer

Ga2Cl6
352.162ǂ 11.1760 1.1520 0.39974 0.11997 0.063996 3.8566 4.5874 0.0291
as dimer

Group IV(A)

SiH4 32.117 7.4652 1.1492 0.37673 -0.11740 0.063525 3.6247 3.4716 0.0914

SiH3Cl 66.562 8.8830 1.1494 0.37847 -0.11761 0.063572 2.7614 4.9286 0.0756

SiH2Cl2 101.007 9.9221 1.1497 0.38068 -0.11788 0.063629 2.1589 6.6627 0.0598

SiHCl3 135.452 9.9796 1.1501 0.38395 -0.11827 0.063708 2.5687 5.7950 0.0446

SiCl4 169.896 10.0240 1.1504 0.38651 -0.11856 0.063765 2.8465 5.3590 0.0291

SiH3(CH3) 46.144 9.0535 1.1494 0.37855 -0.11762 0.063574 2.6380 5.1963 0.0758

SiH2Cl(CH3) 80.589 10.1402 1.1499 0.38201 -0.11804 0.063662 2.1796 6.7338 0.0600

SiHCl2(CH3) 115.034 9.1805 1.1503 0.38512 -0.11840 0.063735 3.3346 4.3854 0.0446

SiCl3(CH3) 149.479 10.1187 1.1506 0.38512 -0.11840 0.063789 2.9357 5.2665 0.0291

SiH2(CH3)2 60.169 8.2747 1.1497 0.38077 -0.11789 0.063631 3.4554 3.9214 0.0603

SiHCl(CH3)2 94.615 9.1914 1.1503 0.38509 -0.11840 0.063734 3.3232 4.4012 0.0447

SiCl2(CH3)2 129.061 10.1007 1.1507 0.38820 -0.11875 0.063801 3.0279 5.1285 0.0291

GeH4 76.642 8.1986 1.1494 0.37858 -0.11763 0.063575 3.3615 3.9445 0.0914

GeH3Cl 111.087 8.6130 1.1496 0.37993 -0.11779 0.063610 3.1743 4.3025 0.0756

GeH2Cl2 145.532 8.8966 1.1498 0.38154 -0.11798 0.063650 3.1184 4.4929 0.0598

GeHCl3 179.976 9.1230 1.1501 0.38344 -0.11821 0.063696 3.1490 4.5641 0.0446

GeCl4 214.421 9.4041 1.1503 0.38554 -0.11845 0.063744 3.1654 4.6724 0.0291

SnH4 122.742 9.0298 1.1497 0.38087 -0.11790 0.063634 3.0671 4.5745 0.0914

SnCl4 260.521 10.1994 1.1506 0.38745 -0.11867 0.063744 2.8471 5.4354 0.0291

Group V(A)

PH3 33.998 7.8041 1.1485 0.37228 -0.11684 0.063396 2.8050 4.3644 0.1009

PH2Cl 68.443 8.6412 1.1490 0.37518 -0.11721 0.063482 2.4920 5.2267 0.0766

PHCl2 102.888 8.5412 1.1495 0.37919 -0.11770 0.063591 3.0719 4.3813 0.0527

PCl3 137.333 7.2366 1.1501 0.38413 -0.11829 0.063712 4.4289 2.9789 0.0291

POCl3 153.331 9.7746 1.1500 0.38331 -0.11819 0.063693 2.4586 5.9576 0

Table 4-2, VP solutions, by Group and Homologue, continued

Fluid MW A B0 B1 B2 B3 c n k

Group V(A),cont’d

AsH3 77.945 7.4208 1.1484 0.37116 -0.11670 0.063362 2.9474 4.0298 0.1009

AsH2Cl 112.390 8.5242 1.1488 0.37403 -0.11707 0.063448 2.3972 5.3468 0.0766

AsHCl2 146.835 9.3548 1.1493 0.37782 -0.11754 0.063555 2.1821 6.2831 0.0527

Larry Coleman 4.3 10/24/2021 [Link]

AsCl3 181.281 10.0281 1.1499 0.38254 -0.11810 0.063675 2.2519 6.5171 0.0291

SbH3 124.781 8.5409 1.1484 0.37136 -0.11673 0.063368 2.0491 6.0821 0.1009

SbCl3 228.115 7.1200 1.1511 0.39192 -0.11916 0.063873 5.3455 2.6317 0.0291

n-pentane 72.149 10.0508 1.1504 0.38593 -0.11849 0.063753 2.6269 5.7673 0

iso-pentane 72.149 9.9797 1.1502 0.38432 -0.11831 0.063717 2.4514 6.0715 0

Two plots are offered in Figure 4-3 to validate the principles of this Distillation Science article. These plots are for the ψ (Psi) function
of Equation (4.5) and the α (alpha) function of Equation (4.4), for the chlorosilane homologue. The derivative function ϕ goes through
the required minimum, which produces the expected inflection in the vapor pressure curve, typically around Tr= 0.80-0.86. Silane’s
minimum ϕ occurs at Tr=0.67 and DCS’s minimum ϕ occurs at Tr= 0.89. The fact that two fluids of this homologue are outliers
illustrates why chlorosilane vapor pressures are so poorly fit by VP expressions that are intended for use with hydrocarbons. However,
it can be seen how the Reidel α function does fit the thermodynamic consistency requirements: for all homologue fluids, the α function
asymptotes to a constant value (α ) as the critical point is approached, with the slope of the α function going to zero as critical is
c

approached.

Figure 4-3 Derivative functions ψ and α , for chlorosilanes

Another validation check is shown in Figure 4-4, comparing the vapor pressures of chlorosilanes and chlorophosphines, which are
distillation separation applications in all commercial chlorosilane purification.

Larry Coleman 4.4 10/24/2021 [Link]

 Figure 4-4, Cross-Plot of Chlorophosphines vs Chlorosilanes

As expected, the chlorophosphine homologue’s VP curves “nest” around those of the chlorosilanes. There are no cross-overs of the
curves, and the location of chlorophosphine vapor pressure curves lie exactly as they are experienced in commercial practice of
chlorosilane purification distillation columns.
Since the math in Equation (4.3) is somewhat involved, an example is provided with calculations.
Example: At 3.50 atmospheres pressure what is the boiling point of trichlorosilane (TCS)?
From Part III, Table 3-3, the pertinent property values for TCS are: Tc=479.15°K and Pc=41.15 atm. MW, Tb , Vc and Zc are not
needed for this example. From Part IV, the VP solution constants for TCS are: A=9.9796; B0=1.1501; B1=0.38395; B2=-0.11827;
B3=0.063708; c=2.5687; n=5.7950; and k=0.0446.
As all advanced VP equations are, determining the TCS boiling point at 3.50 atmospheres is iterative; so it is nice to have a VP plot to
get a good first guess. Using above Figure 4‑4, for
\ln (VP) =\ln (3.50) = 1.2528, it looks like 1/T= 2.89E-3 might be a good first guess => T=346°K; so first guess Tr=346/479.15
=0.7221. Plugging that first guess value of Tr into the T-S VP equation with the constants above:
n−1
Tr −1
ln(Pr ) = A[ B0 − Tr
−1 2
− B1 ln Tr + B2 Tr − (1/2)B3 Tr ] + c[
n−1
−1
+ k(Tr
− 1)] yields a reduced VP value of Pr = 0.08273; so
VP=0.08273*41.15 = 3.404 atm, which is just a bit less than the desired 3.50 atm. So T needs to increase a little from the 346°K initial
first guess.
Re-doing the TCS VP equation with T= 347°K yields a VP of 3.50 atm (note: 347.05°K is the answer, if the solution is worked out to
two decimal places of temperature). So the boiling temperature at 3.50 atmospheres is 347°K or 73.9°C.

Larry Coleman 4.5 10/24/2021 [Link]

5: Equation of State
Distillation Science (a blend of Chemistry and Chemical Engineering)
This is Part V, Equation of State of a ten-part series of technical articles on Distillation Science, as is currently practiced on
an industrial level. See also Part I, Overview for introductory comments, the scope of the article series, and nomenclature.
Part V, Equation of State explains why an Equation of State (EOS) is needed to evaluate the pure-component physical
properties that are used in distillation science, reviews the various options and makes a recommendation. Also included are
techniques for solving the EOS cubic equations via spreadsheet (i.e., MS Excel) and some work-arounds near the critical point.
This article uses the critical constants and acentric factor of Part III. The EOS solutions of this article will feed values to Part
VI, Fugacity.
The general PVT relationship for a fluid is
P V = ZRT (5.1)

where the fluid could be a sub-cooled liquid, a saturated liquid (i.e, at its boiling point), a saturated vapor (i.e., at its dewpoint),
or a superheated vapor (i.e., a gas). In the ideal situation ( ambient pressure and temperature), ZL is very close to zero for either
a sub-cooled or saturated liquid, and Zv is close to unity for either a saturated or superheated vapor. The Ideal Gas Law has
ZV = 1 , so ΔZ = (Z − Z ) ≈ 1
v L

ΔZ is an important parameter to Equation 2-2 and 4-1 in previous articles Part II and Part IV. For calculations in this and
other Parts, temperature and pressure are in absolute units (°K and atmospheres), and volume in molar units (e.g., cc/mole). In
that case, the gas constant R has units of 82.057 atm-cc/mole-K.
Van der Waals EOS and more recent improvements
For a real fluid, Z is never 0 and almost never =1, but rather a value between 0 and 1, except for gases at high temperatures
and pressures (where Zv >1). While the Ideal Gas Law is a nice introductory concept, it is rarely used in distillation science.
Instead, a more accurate expression is needed, and so Equation (5.1) is re-written as
PV
Z = (5.2)
RT

and different models are used to give results: Z as a function of T & P . The oldest is the Van der Waals EOS from the 1873
based on a model that molecules are hard bodies that take up space and have inter-particle forces. Note that an Equation of
State covers all non-solid phases, so for saturated systems (i.e., fluids at their boiling/condensation point) it will have two real
solutions: one for the vapor and one for the liquid: Z and Z . For a non-saturated system (i.e., a super-heated vapor or sub-
V L

cooled liquid) it will have only one real solution.

Written in terms of Z, the Van der Waals EOS is
a
Z = (P + ) × (V − b)/RT (5.3)
2
V

where a and b are constants derived from just the fluid’s critical properties. While a good general approximation at lower
pressures (better than the Ideal Gas Law), this expression becomes less accurate as the critical point is approached, and it
always yields 0.375 as Zc. There are very few fluids that have such a high Zc - most are closer to 0.30. Solutions for saturated
phase values, Z and Z , require manipulation of a cubic equation and solving its roots. However inaccurate it may be, it does
v L

allow estimation of the saturated vapor and liquid molar densities: ρ = 1/V and ρ = 1/V , as well as ΔZ = Z − Z .
ν ν L L v L

In Part IV, Equation 4‑1 showed how the latent heat of vaporization, ΔH , is calculated from the differential Clausius-
v

Clapyeron equation, once a good vapor pressure relationship is known and ΔZ is closely determined. And in Part VI, it will
be seen how fugacity coefficients are calculated from an EOS, to partially account for some departures to ideal behavior in
mixtures of fluids.
In 1949, the Redlich-Kwong EOS improved on Van der Waals, and in 1972 Soave made a modification to yield the S-R-K
EOS, by adding in the effect of the acentric factor, ω. Further improvement was done in a general-application EOS by Peng-
Robinson in 1976. Since then, there have been many additional EOS proposed for specialized applications which give better
results, but require additional data the additional constants. A full history of various EOS’s is given in Wikipedia.

Larry Coleman 5.1 10/2/2021 [Link]

For purposes of practical chlorosilane distillation design of fugacity (see Part VI), the Peng-Robinson (P-R) EOS appears to
give satisfactory results. P-R does not appear to be accurate enough to extend to closely estimating saturated phase densities of
chlorosilanes, or latent heats. If enough data were available, it might be possible to determine which of the newer, more
specific EOS forms is accurate enough for those derivative properties (i.e., superior to Peng-Robinson). Unfortunately, that
data does not yet exist (although a few key data points show the failings of P-R).
Regardless of the fluid parameters chosen, the P-R EOS yields a Zc of 0.307 at the critical point. While this Zcvalue is still
higher than that of chlorosilanes, the accuracy of Zv and ZL values are adequate for fugacity estimation, through 90% of
critical pressure (Pc). Since commercial applications only rarely approach 90% of Pc for reasons of process and equipment
design, the P-R EOS is an acceptable compromise. It must be pointed out that the P-R EOS should not be used in any
calculations that involve differential or integral calculus manipulation of Z, since erroneous (or impossible) values result.
Peng-Robinson Equation of State
RT aα
p = − (5.4)
2 2
Vm − b Vm + 2b Vm − b

with
2 2
0.45724R Tc
a =
pc

0.07780RTc
b =
pc
0.5 2
α = (1 + κ(1 − Tr )

κ = 0.37464 + 1.54226ω − 0.26992ω

2
when ω < or = 0.49
+ 0.016666 ω when ω>0.49
2 3
κ = 0.379642 + 1.48503ω − 0.164423 ω

T
Tr = noting that ω, κ, α are dimensionless
Tc

In more compact and dimensionless form, the Peng-Robinson EOS can be re-stated as a cubic in Z using
αap bp
A =
2 2
and B =
R T RT

3 2 2 2 3
Z − (1 − B)Z + (A − 2B − 3 B )Z − (AB − B −B ) = 0 (5.5)

At a given temperature, T, and for fluid physical properties Tc, Pc, and ω (from Part III), and where P is the vapor pressure
calculated for that fluid by the VP equation of Part IV, dimensionless parameters “A” and “B” are determined for Equation (
5.5). This sets up a cubic equation in Z, which has three real roots (as opposed to being imaginary). The largest root is the

value of Zv; the smallest root is the value of ZL; and the third root is discarded as having no physical meaning.
Then knowing Zv and ZL, the value of ΔZ can be calculated, from which the latent heat of vaporization (at that T & P) could
be calculated using Equation 4-1. Also from Zv and ZL, the saturated phase densities can be calculated as estimates. But again,
the P-R EOS will not hold up under this much mathematical manipulation to give any more than rough estimates of these
distillation design properties; and so it is recommended to only use the P-R EOS to estimate fugacities (unless rough estimates
of latent heat and phase densities are acceptable) . The P-R EOS is offered to lay the groundwork for future improvement.
A cautionary statement is in order regarding the use of an EOS (whether Van der Waals, SRK, P-R or others) with liquid
mixtures. In some hydrocarbon applications, VLE convergence calculations attempt to combine the critical properties of a
mixture's components, as a "psuedo-critical" mixture property, and use these "psuedo-criticals" in EOS calculations along with
various mixing models. While there has been some success in this approach for certain hydrocarbons ( mostly in the oil and
gas industry), such approach rarely gives good results with inorganic polar or mildly polar mixtures, such as chlorosilanes.
Various papers have been issued over the recent years promoting such approach, suggesting the use of "binary interaction
parameters" to estimate an EOS for a mixture. There are simply too many differences in polarity, acentric factor, and
individual fluid critical properties for this to work. Therefore, the reader is discouraged from such practice. Sometimes
calculations must be worked out, checked along the way, and results validated - rather than taking an "easy way out". Certainly
with computer automation, there is no reason not to do the full calculation procedure.

Larry Coleman 5.2 10/2/2021 [Link]

One stumbling block to many is that solving the EOS equation for Z means solving a cubic equation. Currently there are no
MS Excel functions for solving the roots of a cubic equation. However a ready solution is found in an older CRC “Standard
Math Tables” book, which can be solved with a spreadsheet or programmed into a macro.

Cubic Equations

Method for solving for the roots of a cubic equation in the variable "y", with roots y1, y2 and y3
3 2
y + py + qy + r = 0 (5.6)

p
may be transformed to the format , x 3
+ ax + b = 0 by substituting for y the value, x − . Where a =
3
1

3
2
(3q − p ) and
1 3
b = (2 p − 9pq + 27r)
27

If p, q, and r are real, then

2 3

If b

4
+
a

27
>0 , there will be one real root and two conjugate imaginary roots.
2 3

If b

4
+
a

27
=0 , there will be three real roots of which at least two are equal.
2 3

If b

4
+
a

27
<0 , there will be three real roots and unequal roots.
For saturated liquids and vapors, the last case is always true so a trigonometric solution is useful for solving EOS
equations. Compute the value of the angle ϕ in the expression
−−−−−−
3
−b −a
cos ϕ = ÷ √( ) (5.7)
2 27

Then the three transformed roots x1, x2, and x3 will have the following values:
−
−−
ϕ p
x1 = 2 √
−a

3
× cos(
3
) and so cubic root y 1 =−
3
+ x1

−
−−
ϕ p
x2 = 2 √
−a

3
× cos(
3
+ 120 ) º and cubic root y 2 =−
3
+ x2

−
−−
ϕ p
x3 = 2 √
−a

3
× cos(
3
+ 240 ) º and cubic root y 3 =−
3
+ x3

There is an alternate algebraic solution to solving a cubic equation (which is frequently needed) to evaluate Z for superheated
vapors. But for liquids and vapors near saturation, this alternate solution frequently breaks down close to 80-85% of critical
pressure using MS Excel, due to software problems. Hence, the recommendation of the more stable trigonometric method,
even though its math is seemingly odd. The trigonometric root-solving technique can break down somewhere near 95% of
critical, using MS Excel because of numerical precision; but since the P-R EOS is only valid to 90% of critical, that is not a
problem. Note that there are still other cubic root solving algorithms, but these also require even higher levels of numerical
precision and are often found to break down (e.g., division by zero or requiring math with imaginary numbers).
When doing the above calculations, note that cos-1(ø) is only valid over the range between -1 and 1, returning values of ø
between 1 to 0 radians, respectively. When doing automated calculations (via a spreadsheet like MS Excel), it is a good idea to
−−−−−−
check to see if the quantity −b/(2√−a /27) is between -1 and 1, in order to avoid an error message.
3

In Equation (5.5) above, the cubic equation’s Z3 term is unity, the Z2 is [- (1-B)] = B-1, the Z term is [A‑2B‑3B3] and the
constant is [- (AB-B2-B3)]. Then the solution of that cubic equation, using Equation (5.6):
p = B−1
2
q = A − 2B − 3B
2 3
r = −(AB − B −B )
2
a = (3q − p )/3
3
b = (2 p − 9pq + 27r)/27

2 3

Intermediate parameters “p” and “r” will usually be negative, and parameter “q” will be positive; and the quantity +
b

4
a

27

will always be negative for liquids and vapors at or near saturation: so there will be three unequal roots. Note that for
2 3

superheated vapors or vapors above critical pressure (aka gases), the quantity + can be >0, in which case there is only
b

4
a

27

one real root for Z, and an alternate algebraic solution method will be needed ( i.e., the trigonometric method breaks down) . In
such case, Z will be >1.

Larry Coleman 5.3 10/2/2021 [Link]

To use the trigonometric method to get Z and ZL , first solve for the angle ø (in radians) of Equation (5.6):
v

−−−−−−
ϕ = cos
−1 3
(−b/(2 √−a /27) and so the three roots, y1, y2, and y3 will be:
−−−−− −−−−−
y1 = −p/3 + 2 √−a/3 × cos(ϕ/3) , y2 = −p/3 + 2 √−a/3 × cos(ϕ/3 + 2π/3) and
−−−−−
y3 = −p/3 + 2 √−a/3 × cos(ϕ/3 + 4π/3)

The largest valued root is Zv , and the smallest valued root is ZL, with the intermediate valued root being discarded as
meaningless.
Since the math in Part V of both the EOS and the cubic solution might seem daunting, an example is provided with
calculations.

Example 5.1

Continuing the example of Part IV, of TCS vaporizing at 3.50 atm, and 347.05°K = 73.9°C, determine the values of Zv ,
ZL , and of ΔZ = Zv - ZL. From the values of Zv and ZL determine the saturated vapor and liquid densities and then from
the value of ΔZ determine the value of the latent heat of vaporization at 3.50 atm.
Solution
From Part III, Table 3-3, the pertinent property values for TCS are: MW = 135.452; Tc=479.15°K; Pc=41.15 atm ; and
ω=0.2090. Tb , Vc and Zc are not needed for this example.
For use in the P-R EOS, at the examples temperature of 347.05°K, T r= 347.05/479.15 = 0.7243.
Using Equation (5.5) above for the P-R EOS and plugging in the values for T, P, Tr, Tc, Pc, and ω:
2 2
a = 0.45724 × 82.057 × 479.15 /41.15 = 1.7177E7

b = 0.07780 × 82.057 × 479.15/41.15 = 74.336

2
κ = 0.37464 + 1.54226 × 0.02090 − 0.26992 × 0.2090 = .68518

0.5 2
α = (1 + 0.68518 × (1 − 0.7243 )) = 1.2145@ : 347.05°K

A = 1.2145 × 1.7177E7 × 3.50/(82.057 × 347.05 )

2 2
=0.09003 @: 347.05°K, 3.50 atm
B = (74.336 × 3.50)/(82.057 × 437.05) = 9.136E − 3 @: 347.05°K, 3.50 atm
Then the cubic equation to solve (using Equation (5.6)) is Z 3
+ pZ
2
+ qZ + r = 0 , with terms:
p = B−1
2
q = A − 2B − 3B
2 3
r = −(AB − B −B )
2
a = (3q − p )/3
3
b = (2 p − 9pq + 27r)/27

p = 0.0091361 − 1 = −0.99086
2
q = 0.09003 − 2 × 0.00913161 − 3 × 0.009136 = 0.071516
2 3
r = −(0.09003 × 0.009136 − 0.009136 − 0.009136 ) = −0.00073819

so for the trigonometric solution a and b are:

2
a = (3 × 0.071516 − (−0.99086 ) )/3 = −.25575

3
b = (2(−0.99086 ) − 9 × (−0.99086) × 0.071516 + 27 × (−0.00073819))/27 = −0.049179

2 3
b /4 +a /27 = -1.49326E-5, so three real roots for Z are confirmed
−−−−−−
ϕ = cos
−1 3
(−b/(2 √−a /27) = 0.15588 radians and so the three roots are:
−−−−−
Z1 = −p/3 + 2 √−a/3 × cos(ϕ/3) = 0.91345
−−−−−
Z2 = −p/3 + 2 √−a/3 × cos(ϕ/3 + 2π/3) = 0.012439
−−−−−
Z3 = −p/3 + 2 √−a/3 × cos(ϕ/3 + 4π/3) = 0.064968

Larry Coleman 5.4 10/2/2021 [Link]

 Zv is the largest root = 0.91345 and ZLis the smallest root = 0.01244 , with the Z3 root being discarded as meaningless. Zv
and ZL will be used in Part VI to give estimates of fugacity coefficients, which partially account for departures from ideal
mixture behavior.
ΔZ = Zv - ZL= 0.91345-0.01244 = 0.90101 (note that this is less than the unity value of Part II, Equation 2-3 and 2-4 ,
which led to simplified VP equations).
If Equation (5.1) is re-written as 1

V
=
P

ZRT
that gives the density in molar units,
and then multiply by molecular weight to get density in more familiar mass units.
Predicted vapor density @ 3.50 atm, 73.9°C = 3.50/(0.91345 × 82.057 × 347.05) = 1.3455E-4 g-mole/cc . In mass
units = (135.45 x 1.3455E-4) = 1.822E-2 g/cc
Predicted liquid density @3.50 atm,73.9°C = 3.50/(0.012439 × 82.057 × 347.05) = 9.8804E-3 g-mole/cc . In mass
units = (135.45 x 9.8804E-3) = 1.338 g/cc
Looking at some TCS saturated liquid phase density data, it can be seen that the predicted saturated liquid density by the
P-R EOS is high by about 8%, with a somewhat greater inaccuracy for sub-cooled liquid TCS. So unless there is better
data available, the P-R EOS should not be thought of as better than ±10% on saturated phase density. Comparable TCS
saturated vapor phase density data does not seem to exist, so no accuracy estimates can be made.
In order to estimate the latent heat of vaporization at these conditions, it is necessary to go back to Part IV, Equation 4.5
and the example values: Tr= 0.7243; A=9.9796; B1=0.38395; B2= -0.11827; B3=0.063708; c=2.5687; n=5.7950; and
k=0.0446.
ΔHν
2 3 n
ψ = = A(1 − B1 Tr + B2 Tr − B3 Tr + c(Tr − k))
ΔZRT

and substituting in the values from above, but using a value for R=1.9859 to get ΔHv in energy units of calories/gram-
mole;
2
ΔHν /(0.90101 × 1.9859 × 347.05) = 9.9796 × (1 − 0.38395 × 0.7243 − 0.11827 × 0.7243 − 0.063708

3
× 0.7243 ) + 2.5687(0.72435.7950 − 0.0446)

So ΔHv is predicted at = 4,114 cal/g-mole. Multiplying by the MW, the latent heat of vaporization is predicted as 30.4
cal/gram at 3.50 atm and 72.9°C. Better estimates of the TCS latent heat at these conditions by previous work done under
NASA contract indicated a value closer to 50 cal/gram (but these were still estimates, not direct measurements, and are
possibly in error as well).

This illustrates the potential power of using an EOS with a thermodynamically consistent VP equation, but also shows the
degree of inaccuracy of pairing the P-R EOS with the VP derivative of Part IV, Equation 4.5. In this case, it is unknown how
much of the inaccuracy is caused by the EOS, and how much caused by VP equation’s derivative.
As will be seen in Part VI, the best use of the P-R EOS is just to estimate fugacity coefficients, which are relative quantities
rather than absolute.

Larry Coleman 5.5 10/2/2021 [Link]

6: Fugacity
Distillation Science (a blend of Chemistry and Chemical Engineering)
This is Part VI, Fugacity of a ten-part series of technical articles on Distillation Science, as is currently practiced on an
industrial level. See also Part I, Overview for introductory comments, the scope of the article series, and nomenclature.
Part VI, Fugacity deals with the expected departure of apparent pure-component vapor pressure in the vapor-liquid
equilibria (VLE) of binary mixtures, as is commonly found in practical application of distillation science. Part VII,
Liquid Activity Coefficients builds on this for mixtures that have greater-than-expected departures from ideal behavior.
This article uses the Vapor Pressure of Part IV and Equations of State of Part V.
The ideal vapor-liquid behavior of volatile fluid mixtures is governed by the combination of Raoult’s Law and Dalton’s
Law. Neither of these is a scientific law but rather a representation of ideal behavior, closely approximated when the fluids
in question have very similar properties. Raoult’s Law of liquid partial pressures ( PPi= Xi* VPi) states that the partial
pressure of a fluid is the liquid mole fraction times the pure-component vapor pressure. Dalton’s Law (Yi= PPi / Σ{PPi} )
states that the vapor mole fraction of a volatile component is that component’s partial pressure divided by the sum of all
partial pressures. If either of these two relationships is extended, with either Xi or Yi being zero or unity, a tautology
results. The nomenclature used here is: Xi and Yi are the liquid and vapor mole fractions of the ith component in the
mixture; VPi and PPi are the vapor pressure and partial pressure of the ith component.
Combining Raoult & Dalton, for a binary mixture:
X1 × V P1
Y1 = (6.1)
X1 × V P1 + X2 × V P2

X2 × V P2
Y2 =
X1 × V P1 + X2 × V P2

However, in actual practice this never exactly works out. The more volatile fluid exerts an effect on the less volatile fluid,
seemingly boosting that fluid’s vapor pressure slightly. Conversely, the less volatile fluid slightly reduces the more volatile
fluid’s vapor pressure. There are also interactions in the vapor phase due to compressibility (Zv, from Part V), which are
sometimes significant. These departures from ideal behavior are due to two factors: (1) the fluid’s phase change is at a
temperature that is above or below its pure components’ boiling points; (2) interactions occur between molecules with
different properties. From a distillation column design perspective, these departures are important, since they make the
component separations more difficult.
The technical term for the first factor of departure is fugacity, which has the same units as the fluid’s vapor pressure,
including the slight positive or negative secondary effect. The ratio between a liquid’s fugacity and its vapor pressure is the
liquid-phase fugacity coefficient; between the vapor’s fugacity and its partial pressure is the vapor-phase fugacity
coefficient. For a fluid by itself, fugacity has no meaning.
But in a mixture of two volatile fluids, it is the vapor and liquid fugacities that are in equilibrium with real fluids ( as
opposed to a mixture of ideal fluids, such as with Raoult/Dalton).
L L
F =ϕ × V Pi (6.2)
i i

where øLi denotes the liquid fugacity coefficient

V V
F =ϕ × P Pi (6.3)
i i

where øVi denotes the vapor fugacity coefficient

Consider a mixture of "fluid 1" and "fluid 2". At a given temperature, "fluid 1" has a pure-component vapor pressure of
VP1. But the presence of "fluid 2" alters that somewhat, so "fluid 1's" fugacity (F1) is different from its pure-component
vapor pressure. Conversely the presence of "fluid 1" somewhat alters "fluid 2's" fugacity (F2), away from it's pure-
component vapor pressure of VP2. The higher the mole fraction of "fluid 2", the more "fluid 1's " fugacity is altered away
from its pure-component vapor pressure. And vice-versa.

Larry Coleman 6.1 11/7/2021 [Link]

This is how fluid mixtures behave in the real world: Raoult/Dalton is just a introductory concept that would apply in an
ideal world.
(The technical term for the second factor of departure is the Liquid Activity Coefficient, which is covered extensively in
Part VII, but initially mentioned here briefly for a more complete understanding.)
Fugacity effects are always present whenever two fluids of different volatility are mixed and vaporized or condensed. The
Liquid Activity Coefficient effects are always present, but more noticed when the two fluids have significantly different
properties, and even more so when the fluids are polar.
This is best illustrated by an example, partially continued from Part IV and Part V. At 73.90°C = 347.05°K, TCS has a
pure component vapor pressure of 3.50 atmospheres. At that same temperature, STC has a pure component vapor pressure
of 1.65 atmospheres. In this example, the liquid mixture is 40% mole fraction TCS, so X1= 0.40 . Therefore STC’s liquid
mole fraction, X2= 0.60. Note: it is conventional notation for the most volatile fluid to have the lowest subscript.
Per Raoult, the partial pressures PP1 and PP2 should be 0.40 × 3.50 = 1.40 and 0.60 × 1.65 = 0.991 atmospheres,
respectively. This liquid mixture at 73.90°C would be expected to boil at 0.40 + 1.99 = 2.39 atmospheres total system
pressure; with the TCS mole fraction of the first bubble’s vapor (Y1) being 1.40/2.39 = 0.586= 0.586 (so Y2 is 0.414).
Instead, it is found that first bubble’s vapor has a TCS mole fraction (Y1) of only 0.573 (with Y2 being 0.427) indicating
that the TCS was just a little less volatile than expected and the STC was just a little more. While there is just over 2%
difference in this case between Raoult/Dalton’s expected and the actual resulting vapor mole fractions, the slight difference
is truly present. If the difference in vapor pressures between the two fluids increased, and the overall pressure increased
towards critical, the departure from the Raoult & Dalton relationship would then increase.
Two departure-from-ideality factors were purposefully shown in the above example: Fugacity and Liquid Activity
Coefficients. The vapor-phase and liquid-phase fugacity coefficients (ϕ & ϕ ) are never exactly unity (as implicit with
V
i
L
i

Raoult &Dalton), but they do often approach unity (especially at low pressures and for molecularly similar fluids).
Likewise, the Liquid Activity Coefficients (γ ) are never exactly unity for real fluid mixtures (see Part VII), but are often
i

set to unity because the fluids are non-interactive. The point to make here is that non-unity fugacity coefficients are just a
result of fluids having different volatilities (i.e., they are not fluid-dependent); whereas Liquid Activity Coefficients are
specific fluid-system dependent (e.g., TCS-STC has a different set of parameters than DCS-TCS). (There is no analog to
Liquid Activity Coefficients in the vapor phase - vapors just interact due to fugacity effects.)
In the above example, the overall deviation of 2.3% (0.586 vs 0.573) on TCS mole fraction is almost all due to fugacity
(both liquid and vapor), with a lesser amount due to Liquid Activity Coefficients. When designing a distillation for high-
reflux purification, the relative importance between these two departure factors often reverses, as well as the overall
magnitude.
The general equation (including the γ covered in Part VII), is that:
i

L V
Yi × π = P Pi = (ϕ /ϕ ) × γi × V Pi × Xi (6.4)
i i

where π is the total pressure

Setting aside further discussion of Liquid Activity Coefficients until Part VII, the remainder of this article deals with how
to determine the fugacity coefficients,ϕ & ϕ . That determination is dependent on the constants of Part III, and the
V
i
L
i

calculation of the compressibility factors Zv and ZL from Part V.

For either phase, the equation for fugacity coefficient, as derived from thermodynamics is:
Fortunately, the Peng-Robinson Equation of State (P-R EOS) is good enough for most distillation design to evaluate
fugacity coefficients based on some simplifications (see Part V). Admittedly the P-R EOS does not exactly fit
chlorosilanes (and similar fluids), and a better EOS would be preferable, if it exists. However P-R is adequate to the task,
especially with the “work-arounds” given below.
For ϕ there is an almost linear relationship between Zv and VP until close to Pr=0.4. So that allows Equation (??? ) to be
V
i

reduced to
T ,P sat
Ln(ϕ) = (Z −Z ) (6.5)
V V

Larry Coleman 6.2 11/7/2021 [Link]

 T ,P
where Z V
and Z are respectively the component’s vapor-phase compressibility at the temperature and pressure of the
sat
V

boiling mixture, and at the pure-component saturation condition (T & VP). For Pr>0.4, Equation 6-6 tends to give low
values of ϕ .ν

Since many industrial processes operate at higher pressures (closer to critical than Pr=0.4) another work-around is needed.
Between 0.3< Pr<0.90, a cubic polynomial fits well for Zv as a function of Pr, which then integrates (Zv-1)/Pr of Equation (
??? ) to the function:

ν 3 3 2 2
Ln(ϕ ) = a3 (P −P ) + a2 (P −P ) + a1 (Pr2 − Pr1 ) + a0 [(LnPr2 ) − Ln(Pr1 )] (6.6)
r2 r1 r2 r1

where Pr2 and Pr1 are the reduced pressures of the mixture boiling pressure and the component’s saturation pressure. The
constants for Equation (6.6), a3 through a0, are -0.175248; 0.393866; -0.887714; and -0.0141409, respectively.
For ϕ , the evaluation of Equation (??? ) is simpler, since ZL is fairly insensitive to pressure, simplifying the relationship
L
i

for the ith component to be:

L
Ln(ϕ ) = [ ZL × (π − V P )/V P ] (6.7)

where π is the system total pressure, VP is the pure component’s vapor pressure and Z is the compressibility at the
L

component’s saturation temperature/pressure.

To the extent that the system pressure is greater than the ith component’s vapor pressure (i.e., that liquid component is
forcibly sub-cooled), the value of ϕ will be greater than unity. To the extent that the system pressure is less than the ith
L
i

component’s vapor pressure (i.e., that liquid component is forcibly super-heated), the value of ϕ will be less than unity.
L
i

The pattern for ϕ is usually the reverse as a result of compressibility changes, so the value of ϕ /ϕ in Equation (6.4)
V
i
L
i
V
i

above shows how the distinct difference in Raoult/Dalton departure results. In the example above (40% molar TCS in STC
@ 73.9°C saturation), the values are given in Table 6-1:
Table 6-1: Comparing TCS and STC Fugacity Coefficients
TCS STC

ϕ
L
0.9976 1.0039

ϕ
V
1.0488 0.9562
ϕ
L
/ϕ
V
0.9512 1.0498

The careful reader will note an apparent discrepancy between the overall 2.3% ideality departures reported in the example
above and the fugacity ratio of Table 6-1. Just the fugacity effects should make the TCS mole fraction departure about
4.8% lower than Raoult/Dalton expectation, but the Liquid Activity Coefficient effects restores about half of the departure.
Thus it can be seen that fugacity effects and Liquid Activity Coefficient effects can be counteractive.
The “take-away” from this article is:
Raoult & Dalton describe ideal behavior of vaporizing liquid mixtures, but do not completely describe the behavior of
real fluids
All departure factors must be considered: liquid fugacity, vapor fugacity and Liquid Activity Coefficients

Larry Coleman 6.3 11/7/2021 [Link]

7: Liquid Activity Coefficients
Distillation Science (a blend of Chemistry and Chemical Engineering)
This is Part VII, Liquid Activity Coefficients of a ten-part series of technical articles on Distillation Science, as is currently
practiced on an industrial level. See also Part I, Overview for introductory comments, the scope of the article series, and
nomenclature.
Part VII, Liquid Activity Coefficients builds on Part VI, Fugacity regarding the departure of vapor-liquid equilibria (VLE) of
binary mixtures from ideal behavior, as is commonly found in practical application of distillation science. In conjunction with
previous articles, the goal of this article to to complete the explanation of equilibrium behavior of binary systems; such that Part
IX can illustrate an example of a distillation process design.
In Part VI the combination of Raoult’s Law and Dalton’s Law was introduced governing the ideal vapor-liquid behavior of volatile
fluid mixtures, with Equation 6-1 repeated for convenience.
Y1 = (X1 × V P1 )/(X1 × V P1 + X2 × V P2 ) Y2 = (X2 × V P2 )/(X1 × V P1 + X2 × V P2 ) repeat Equation 6-1
Also in Part VI the real-world departure from ideal behavior was discussed as far as the fugacity coefficients, but stopped short of
Liquid Activity Coefficients. Equation 6-4 repeated below gave the complete relationship between the partial pressure of the ith
component, PPi, and the pure-component vapor pressure VPi , for liquid and vapor mole fractions Xi and Yi, respectively. The total
system pressure (sum of the partial pressures) is denoted as π; (ϕ & ϕ ) are the liquid and vapor fugacity coefficients. γ is the ith
V
i
L
i i

component’s Liquid Activity Coefficient for that specific set of components.

Yi × π = P Pi = (ϕ
L
i
V
/ϕ
i
) × γi × V Pi × Xi repeat Equation 6-4
An example was given in Part VI of a vaporizing TCS-STC mixture, contrasting the vapor-liquid equilibria (VLE) expected from
Raoult/Dalton; from what actually occurs in the real world as a result of the ideality departures. Summarizing the previous article,
the more volatile fluid exerts an effect on the less volatile fluid, to increase its apparent vapor pressure; and the less volatile fluid
reduces the more volatile fluid’s vapor pressure. The net effect is to make the component separations more difficult via distillation
(e.g., more theoretical trays or increased amount of reflux than expected). For detailed discussion of fugacity coefficients, see Part
VI.
As a “house-keeping” note, it must be emphasized that Parts VI and VII treat fugacity coefficients and Liquid Activity
Coefficients separately. In some texts, the two topics are compressed together using specific models of mixing and Equations of
State (EOS). However, the two ideality departure factors are differently based, and their combination often gives erroneous or
nonsense results when used with polar fluids, or those that have a high degree of hydrogen bonding (such as the chlorosilanes that
are used as continuing examples). Fugacity coefficients naturally result from differences in fluid volatility. Liquid Activity
Coefficients result from differences in fluid properties, including critical temperature (Tc), critical pressure (Pc), critical volume
(Vc) and acentric factor (ω), as well as the entropy & enthalpy effects of mixing (referred to in thermodynamic texts as excess
Gibbs Free Energy).
For a given temperature and combination of fluids, the value of γ is a function of the liquid mole fraction, rising asymptotically
from unity at the pure component condition. Typically the values of γ for a system are plotted as a function of the more volatile
component’s mole fraction. An example of such a plot is given in Figure 7-1, for the TCS-STC binary at 73.9°C (the temperature of
the continuing example of previous articles).

Larry Coleman 7.1 10/31/2021 [Link]

 Figure 7-1, TCS-STC Liquid Activity Coefficients @ 73.9°C
There are several models for Liquid Activity Coefficients that are thermodynamically consistent (i.e., follow the Gibbs-Duhem
Rule of binary system thermodynamics). The two-constant Van Laar model is the easiest to manipulate, but is limited to binary
systems. The Margules model is not that different from the Van Laar model, and also limited to binary systems. The Wilson model
is more complex mathematically, but can be applied to ternary or greater-numbered component systems. Fortunately the Wilson
model parameters can be calculated from the Van Laar model. Even more complex Liquid Activity Coefficient models are known
(e.g., NRTL and UNIQUAC), but these require more data points to effectively evaluate. In the case of the above plot, the Van Laar
model is shown. In Part VIII, the techniques are discussed to evaluate experimental data on Liquid Activity Coefficients and fit
the data to a model. In most cases with chlorosilanes (and their impurities) the data quality is not that great, so the simpler Van Laar
is typically used to correlate data and for evaluating binary systems, and the Wilson model is used for ternary systems and beyond.
For a more in-depth discussion of Liquid Activity Coefficients, the reader is referred to the text " The Properties of Gases and
Liquids", by Reid, Prausnitz and Sherwood, (McGraw-Hill). The third edition is more readable, but the fifth edition is more
updated.
The Van Laar model is:
−2
Ln(γ1 ) = A1 × [1 + (A1 × X1 )/(A2 × X2 )] (7.1)

−2
Ln(γ2 ) = A2 × [1 + (A2 × X2 )/(A1 × X1 )]

with A1 and A2 being the Van Laar constants for the more volatile and less volatile component, respectively.
Conveniently,

EXP (A1 ) = γ1 @X2 → 0 (7.2)

EXP (A2 ) = γ2 @X1 → 0

which simplifies evaluation. Note from Figure 7-1 that the Liquid Activity Coefficients plot is not always symmetrical (i.e., A1≠
A2). In the plot of Figure 7-1, A1 (TCS) = 0.1752, so the TCS curve asymptotes at γ= 1.191; and A2 (STC) = 0.2086, so the TCS
curve asymptotes at γ= 1.232 .
Completing the on-going example from Part VI, for a TCS-STC liquid mixture, with 40% molar TCS, at 73.9°C, to determine the
vapor mole fraction in equilibrium (but now including the Liquid Activity Coefficient effects):
In Part VI, the pure-component vapor pressures of TCS and STC respectively were calculated from the VP equation of Part IV to
be 3.500 and 1.651 atmospheres, respectively at 73.9°C. Also from Part V, the liquid and vapor fugacities were calculated as per
Table 6-1, repeated below, using the Peng-Robinson EOS (Part V).
Repeat Table 6-1: Comparing TCS and STC Fugacity Coefficients
TCS STC

ϕ
L
0.9976 1.0039

ϕ
V
1.0488 0.9562
ϕ
L
/ϕ
V
0.9512 1.0498

Larry Coleman 7.2 10/31/2021 [Link]

The values for γ are calculated for TCS and STC @ 0.40 mf TCS, per Equation (7.1) using the values for A1 and A2 above (or just
read from Figure 7-1), to be 1.075 and 1.027 respectively. Plugging all these values into Equation 6-4 (repeated above), the partial
pressures of TCS and STC @ 73.9°C are:
P PT C S = 0.9512 × 1.075 × 3.500 × 0.40 = 1.432 atmospheres, and
P PST C = 1.0498 × 1.027 × 1.651 × 0.60 = 1.068 atmospheres.
So the value of YT CS = 1.432/(1.432 + 1.068) = .573 , which was given in Part VI, with the STC mole fraction being 0.427.
The careful reader will see that not only are the partial pressures of TCS and STC different from Raoult/Dalton, but so is the total
pressure. In the ideal behavior of Raoult/Dalton, the partial pressures were given in Part VI as 1.400 and 0.991 atmospheres, for a
total pressure of 2.39 atmospheres. But in reality, at 73.9°C and 0.40 mf TCS, the total pressure would instead be 2.50 atmospheres.
So while the ideality departures of TCS made the vapor mole fraction smaller (0.573 vs 0.586), both the TCS and STC exerted
slightly more partial pressure than Raoult would indicate. But the TCS mole fraction still went down, because the ideality departure
of STC was greater than that of TCS.
The take-away from this example is that while Raoult/Dalton will give you a “quick & dirty” value for Y/X, it is easy to get the
overall wrong answer for determining the requirements of a distillation column design. In almost every instance of industrial
practice, the use of Raoult/Dalton will undersize the number of trays needed to make a separation (or for the same number of trays,
undersize the reflux required). For this reason, it is important to know how to work out the better answer.
For many binary systems, VLE data exists – but typically not at the temperature/pressure needed for industrial design. Most
industrial designs use higher pressures for economy of size, as well as to accommodate available energy sources for driving the
column reboiler and condenser. Yet the practicalities of data collection in the laboratory usually require pressures slightly above
ambient, or at slight vacuums. This is especially true of chlorosilanes, since the normal materials of construction for higher
pressures can catalyze slow side reactions, such as disproportionation and dimerization.
Lab data is normally collected based on small amounts of fluid that are used repeatedly at somewhat different combinations of
temperature and composition. In the case of electronic impurities, some of these compounds are not very stable outside of a
chlorosilane matrix. So a correlating method must be used that allows a certain degree of extrapolation.
In 1981 Chung-Ton Lin and Thomas Daubert from Penn State University developed such a Liquid Activity Coefficient model and
published same in Industrial & Engineering Chemistry Process Design and Development, basing their work on non-polar
hydrocarbon mixtures. I have used their general method, but modified two of the constants so as to better fit available VLE data on
chlorosilanes. Using these re-evaluated constants, the revised “Lin & Daubert” method (detailed below as Equation (7.3) was used
to check against industrially obtained data on the distillation purification of TCS and silane, with good results. Obviously more data
would be preferable, especially on ternary mixtures and low-concentration electronic impurities in chlorosilanes. However until
such becomes available, Equation (7.3) below is recommended.
To develop a generic expression for the Van Laar parameters (A & A ), Lin & Daubert go back to Van Laar’s original two-term
i j

expression regarding the excess Gibbs energy of mixing, using the SRK Equation of State for partial molar compressibility and
simple mixing rules; and expand it into the following:
Ai = c × Fi + d × kij × Gi (7.3)

Aj = c × Fj + d × kij × Gj

0.5 0.5 0.5 0.5 2

Fi = 1/[ Tri Pci ] × (Ri × (1 + Mj ) ×P − Rj × (1 + Mj ) ×P )
ci cj

0.5 0.5 0.5 0.5 2

Fj = 1/[ Trj Pcj ] × (Rj × (1 + Mi ) ×P − Ri × (1 + Mi ) ×P )
cj ci
−1 0.5 0.5 0.5
Gi = T (Pcj / Pci ) × Ri × (1 + Mi ) × Rj (1 + Mi )
ri
−1 0.5 0.5 0.5
Gj = T (Pci / Pcj ) × Rj × (1 + Mj ) × Ri (1 + Mj )
rj

2 2
Mi = 0.480 + 1.574 ωi − .176 ω Mj = 0.480 + 1.574 ωj − .176 ω
i j

0/5 0.5 0/5 0.5

Ri = [1 + Mi × (1 − T )] Rj = [1 + Mj × (1 − T )]
ri rj

1/3 1/3 0.5 1/3 1/3

kij = 1 − 2 × [ V ×V ] /[ V +V ]
ci cj ci cj

Using an improved set of curve-fits for “c”,

3 2
0 < Δω < 0.03 ⇒ Ln(c) = −8.0637E + 04(Δω) + 8.2649E + 03(Δω) − 3.2891E + 02(Δω) + 5.5051

0.03 < Δω < 0.30 ⇒

4 3 2
Ln(c) = 6.4645E + 02(Δω) − 7.9857E + 02(Δω) + 3.4768E + 02(Δω) − 6.5099E + 01(Δω) + 2.5667

Larry Coleman 7.3 10/31/2021 [Link]

Reducing available VLE data to Van Laar parameters, and adjusting for vapor and liquid fugacity coefficients, the best value of “d”
is found to be 112.1 for chlorosilanes and similar fluids.
Lin & Daubert suggested that the “c” function be capped at 0.005 < Δω = |ω − ω | < 0.15 , based on their data; but using the
i j

above improved set of polynomials for “c” avoids singularities during iterative calculations, with a minimum “c” of 0.135 (as “Δω”
→∞)
Note that the "c" variable in the first term of Equation 7-3 is based on the F values and the differences between the two fluid's
i

acentric factor, ω ( i.e., their molecular complexity); whereas the second term of Equation 7-3 is based on the G values and the
i

differences between the two fluid's molecular volumes at critical (Vc).

Returning to the continuing previous TCS-STC example, and using the two components’ critical temperature (Tc), critical pressure
(Pc), critical volume (Vc) and acentric factor (ω), as well as basing the reduced temperature (Tr) on 73.9°C = 347.05°K, Equation (
7.3) gives the values of A1 and A2 used earlier in this article, of 0.1752 and 0.2086 respectively.

Frequently in evaluating different distillation designs for a process, the designer is faced with the need to make process temperature
changes and determine the effect on Y/X using Equation 6-4, repeated above. That causes no problem in the evaluation of the
fugacity ratio,ϕ /ϕ , since those factors are derived from the EOS ( see Part VI). However, if using experimentally obtained
L
i
V
i

values for γ that were obtained at a different temperature, there needs to be further adjustment. The general “rule-of-thumb” seen in
some texts is that Ln(γ) is proportional to 1/T (absolute temperature) for small temperature changes. In reality such a “rule-of-
thumb” is tantamount to making an assumption that the mixing effects of the multi-component fluid are closer to an isothermal
model than an isenthalpic one. In most cases, the temperature dependency on activity coefficient is a blend of two theoretical
models, which leads to the relationship:
b
Ln(γi ) = a + (7.4)
T

where “a” and “b” are the constants of a linear relationship (“a” = 0 leading to the isenthalpic model, or “b” = 0 leading to the
isothermal model). The use of Equation (7.3) solves the problem, since it predicts such temperature changes along the lines of
Equation (7.4) , and so can be used in a relative manner. For example, if Equation (7.3) is found to over-predict the value of Ln(γ)
by 10% of an experimentally obtained activity coefficient value, then calculate how much Equation (7.3) would predict for the
different temperature, and apply that 10% on the Ln(γ).
Note that to all intents, the Liquid Activity Coefficient is not a direct function of pressure, other than increasing the boiling
temperature increases pressure. However, fugacities are a function of both temperature and pressure, and that is why they should be
included in Y/X calculations (many texts advocate setting fugacities to unity, which is only close to being valid at very low
pressures).
To demonstrate how Equation (7.3) performs with changing temperature (from 330°K to 380°K), Figure 7-2 shows the dependency
for the asymptotic BIP’s of the TCS-STC mixture in the example above. Note that the slope of the STC Liquid Activity Coefficient
line is about 20% steeper than that of the TCS Liquid Activity Coefficient line (i.e., in the TCS-STC binary system, increasing the
temperature has a greater effect on the STC's γ than it does on the TCS's γ).

Figure 7-2, Temperature dependency of Ln(ϒ) at the asymptotes

Larry Coleman 7.4 10/31/2021 [Link]

8: VLE Analysis Methods
Distillation Science (a blend of Chemistry and Chemical Engineering)
This is Part VIII, VLE Analysis Methods of a ten-part series of technical articles on Distillation Science, as is currently practiced
on an industrial level. See also Part I, Overview for introductory comments, the scope of the article series, and nomenclature.
Part VIII, VLE Analysis Methods recommends methodology used when data-collecting the Liquid Activity Coefficients of binary
systems (see Part VII). This article also deals with validation, especially when data collection is done on reactive fluids like
chlorosilanes which can disproportionate or dimerize during study.
This article uses the information of Parts III through VII.
In Part VII a correlation technique was given for estimating Liquid Activity Coefficients. As mentioned in the previous article it is
common to collect such vapor-liquid equilibria (VLE) data in the lab, at conditions close to ambient pressure. Yet many industrial
applications frequently need to use VLE results at higher process pressure/temperatures or different compositions. To avoid notation
confusion between the liquid activity coefficient (γ ) and the vapor mole fraction (Yi), the mole fractions X and Y are shown as
i

bolded.
Recall from Parts VI and VII that fugacity coefficients (ϕ & ϕ ) are functions of temperature and pressure. However, Liquid
V
i
L
i

Activity Coefficient s ( γ ) are functions of only temperature and mole fraction. In an experimental situation, the calculated values of
i

γ (from the values of Yi and Xi = vapor-phase mole fraction and liquid-phase mole fraction) are easiest to correlate when system
i

temperature is held constant, and system pressure is allowed to vary. Stated otherwise, PXY data @ constant T is far easier to collect
and correlate than TXY @ constant P. However, in commercial applications the distillation column is typically controlled at constant
pressure. So a common industrial communication error is to request the lab to collect VLE data as TXY@P, but instead get PXY@T
instead. There is a way to convert one data set to the other relationship, but it is somewhat awkward and requires choosing activity
coefficient model.
In Figure 7-1 of Part VII, the expected profile of γ vs X @ constant T is shown, and in Figure 7-2 of Part VII the expected profile
of γ vs T @ constant X is shown, at the X=0 and X=1 axes. The reason for generating Liquid Activity Coefficients experimentally is
normally to use the lab results for distillation column design, with the confidence that γ can be accurately known at any value of X or
T, so as to establish the Y/X relationship up and down the column design. Using the principles of Parts VI and VII, these problems
can be resolved.
Having established the experimental protocol, it is necessary to consider the impact of data collection equipment on the fluids being
analyzed and to validate the quality of the sample fluids. For reactive fluids such as chlorosilanes, that means using pressure-capable
glassware is the best choice (or nickel-lined steel as a second choice), since some of the alloying elements of stainless steel will
slowly catalyze disproportionation reactions, which alter the sample compositions.
Also, it is important to validate the purity of the fluid samples, as opposed to blindly accepting the analysis of any accompanying
supplier information. Again, with chlorosilanes (as typical of reactive fluids), the shipping container can catalyze side reactions (i.e.,
the supplier’s COA was perhaps accurate when the sample was loaded into the container, but purity degraded during shipment). It is
common for 99.99% pure TCS samples (under argon inerting) from reputable suppliers to end up having several percent DCS and
several percent STC, along with a few tenths percent hexachlorodisilane and a few tenths percent hydrogen gas, when used just a few
weeks later. A good practice before loading sample fluids into the lab equipment is to double-distill samples (discarding the “lights”
and “heavies” fractions, and assuring that only the “heart’s cut” fraction is used: whose boiling pressure is identical to the expected
pure component vapor pressure).
With the above precautions, data is collected at 15-20 (or more) X,Y points, with some duplicates later in the run to establish
experimental analysis accuracy and confirm the absence of systemic error (such as contamination). At least two pairs of data points
should be collected close to X= zero and X= unity. To help confirm temperature dependency of activity coefficients, at least three sets
of data (each at a different constant temperature or pressure) should be collected.
After data collection, the lab work is suspended, but then "number-crunching" is needed to confirm the data’s validity before
reporting or using the results. If data validity is not confirmed, the error must be found and lab work repeated. First, plot out the data
set and check that both the P-X curve and P-Y curve , or T-X curve and T-Y curve show an identical intersection, at both zero and
unity mole fraction axes; and that the axis intersections are exactly representative of each pure-component’s known vapor pressure. If
the X and Y curves do not intersect at the correct value on each axis, there is some systematic error in the data set that must be
resolved before any more evaluation work is done. A corrupted dataset is likely to be of only minimal value in establishing activity
coefficients; and is more often a reason to toss it all out and re-do the work, after the systemic error is resolved.

Larry Coleman 8.1 10/31/2021 [Link]

Using an algebraic modification of previous Equation 6-4, activity coefficients, the γi i are calculated for each P-X and P-Y data
point.
V L
γi = (ϕ /ϕ ) × (Yi × π)/(V Pi × Xi ) (8.1)
i i

where π is the total pressure, V P is the vapor pressure of the ith component, and (ϕ
i
V
i
/ϕ
L
i
) is the vapor/liquid fugacity ratio of that
component.

Figure 8-1, TXY/P VLE data on TCS-STC system

In Figure 8-1 the twelve-point TXY/P data set of Zanta et al, from Chemicky Prumysl (Czech Chemical Industry) is shown, using a
spreadsheet plot. The researchers did collect a few data points near the STC axis (i.e., TCS= zero). However it would have been
preferable to have a few more collected close to the other (i.e., TCS= unity) axis; and to have that all data reported with more
precision. The researchers chose to collect data at a constant pressure of 740 mm Hg (as opposed to the more preferable PXY/T
method that avoids the need for adjustment of calculated Liquid Activity Coefficient γ for changing temperature), for reasons of their
lab equipment and simplifying their procedure; and to only collect one TXY/P data set.
If only the raw data is considered, it would appear that the T-X and T-Y curves might not mutually intersect at the STC axis, or that
the intersection is higher than the expected 56.2°C temperature representative of STC vaporization @ 740 mm Hg. The reason is that
individual data point precision is low in the region of the STC axis.
However, with data-smoothing and extrapolation of the T-X and T-Y curves via some spreadsheet curve-fitting, the two curves
appear to quite likely both intersect the STC axis near the expected 56.2°C temperature. At the TCS axis however, curve-fitting the T-
X and T-Y data appears to both result in their intersection on the TCS axis at 30.0°C, as opposed to the expected 32.3°C for pure
TCS. The most likely explanation for this data set discrepancy is that the TCS sample supplied was either a disproportionated mix of
DCS/TCS/STC as-received, and/or the fluids in the equipment disproportionated as a result of the equipment materials of
construction. A TCS axis temperature of 30.0°C is representative of a 5% DCS/95% TCS liquid mixture. Surprisingly, the
researchers’ analytical procedure did not pick up the DCS peak, which would have revealed the sample corruption ( but that is a
common error using that the researcher’s type of analyzer).
Now that the data set has been discredited, it remains to be seen what - if any- useful information can be extracted. Other than
demonstrating technique, there is little use in reducing the data that is high in “bad TCS” to get the γ asymptote intersecting the TCS
axis. That is because the XTCS* mole fractions are changing along the STC activity coefficient curve (about a high 1:19 ratio of
DCS/TCS, dropping in ratio to very little DCS in the ternary mixture, close to 56.2°C at the STC axis). It would be impossible to
guess what property values are best to use for fugacity coefficients, given a “wild-card” DCS content). The data closest to
intersecting the STC axis could be interpolated to make up a few “pseudo-points” close to XTCS=0, in order to get a rough value of
the TCS γ asymptote (based on the assumption that there is little DCS in the ternary mix at that point). However, even small amounts
of DCS near the XTCS=0 area could be expected to exert a significant Liquid Activity Coefficient effect, and the fugacity values
might be somewhat in error by assuming only the property values of TCS.
However, to illustrate the data reduction technique for TXY/P data (as opposed to the preferred PXY/T form), a few “pseudo-points”
are evaluated near the axis.
Using a fourth order curve-fit, with the STC axis temperature set to 56.20°C, the T-X and T-Y curve-fit equations are, respectively,
via spreadsheet curve-fitting functions:
T-X T = 18.949 x
4
− 64.825 x
3
+ 75.709 x
2
− 56.044x + 56.20

T-Y T = −22.243 y
4
+ 33.158 y
3
− 18.028 y
2
− 19.075y + 56.20

From those two curve-fits, the following “pseudo-points” are calculated, by solving the quartic equations for T-X and T-Y:
Table 8-1, three “pseudo-points” made from curve-fits near the STC axis
P, mm Hg T, °C XTCS* γ TCS*

Larry Coleman 8.2 10/31/2021 [Link]

 P, mm Hg T, °C XTCS* γ TCS*

740 56.2 0.0000 0.0000

740 55.0 0.00359 0.0104

740 54.0 0.0415 0.1065
740 53.0 0.0620 0.1515

In Table 8-1, XTCS* and YTCS* are subscripted as "TCS*" to denote that while TCS properties are used for vapor pressure and
fugacity estimation, the more volatile fluid is really a DCS/TCS mix. The data-point at X, Y = 0.0000, 0.0000 is given just to show
the STC axis intersection, but cannot be used to calculate the γ asymptote since that would result in “division by zero” in Equation 7-
1.
Using the properties of TCS, the values for VP, øL and øv are calculated as per Parts IV, V, and VI, and the results shown in Table 8-
2 for calculated γ TCS* at 740mm Hg (for each “pseudo-point” temperature).
Table 8-2 for calculated γ TCS* at “pseudo-point” temperatures, from PXY/T data
P, mm Hg T, °C VPTCS* øL øv γ TCS*

740 56.2 2.130 0.9983 1.0321 to be determined

740 55.0 2.055 0.9984 1.0316 1.418

740 54.0 1.994 0.9984 1.0311 1.294
740 53.0 1.934 0.9984 1.0307 1.270

The calculated values for γ TCS* in Table 8-2 are based on the X,Y “pseudo-points” of Table 8-1 and using Equation (8.1) above, to
demonstrate the methodology of experimental data reduction. The trending on fugacity coefficients and activity coefficient is correct:
øL values increase toward unity as XTCS* increases toward unity; øv values decrease toward unity as XTCS* increases toward unity; γ
TCS* values decrease toward unity as XTCS* increases toward unity.

The asymptote at the STC axis is determined by extrapolating value of γ TCS* from the three “pseudo-points”. (That curve would
have a different slope with respect to composition if it were TXY/P data, since temperature is changing as well as XTCS* in Table 7-
2.) This technique identifies the γ TCS* asymptote as 1.437 (the average value of γ TCS* obtained by extrapolating curve-fits of γ TCS*
vs XTCS* to a zero value of XTCS*, and γ TCS* vs YTCS* to a zero value of XTCS*). Compared to other researcher’s data, the calculated
asymptote has a high value. It also seems high per expectation from Equation 7-3 from Part VII, which would indicate a value closer
to 1.21 @ 56.2°C, or 1.23 @ 32.3°C (the expected TCS boiling point at 740 mm Hg). Possibly there is error from the DCS content,
which would tend to make the calculated γ TCS* values high by increasing the non-STC mole fraction in the vapor. The curve-fitting
technique could be “off” since there was low precision in the two data pairs close to the XTCS* =0 axis.
This demonstrates that there is simply no good way to “fix” data that has systematic error in it: virtually all texts suggest that when
corrupted data is encountered, it is pointless to continue with data analysis. In addition to determining reasonable asymptotic values
at either axis, Reid, Prausnitz, and Sherwood suggest in “The Properties of Gases and Liquids”, that all γ calculated values be first
adjusted to a common temperature basis (using the approximation of Ln(γ)× = constant ), then use the Gibbs-Duhem Law to
establish thermodynamic consistency. While technically the best way to validate data for consistency, it requires a significant amount
of high-quality data, and is suggested only if advanced Liquid Activity Coefficient models are to be considered (e.g., NRTL or
UNIQUAC). For Van Laar, Margules, Wilson or similar two-constant models, 15-20 data points should be sufficient.
It is always better to have fewer data points, but have the dataset internally consistent, rather than a large number of possibly
corrupted data points.
This article suggests using the Van Laar model for activity coefficients based on reasonable results with data analysis on
chlorosilanes, and the success in predicting the Van Laar constants via Equation 7-3 in Part VII. However, the reader may want to
explore other models for a better fit to experimental data, after such data has been validated. In “The Properties of Gases and
Liquids”, the authors give an exhaustive list of other activity coefficient models to consider, with some notes on “pro & con”.

Larry Coleman 8.3 10/31/2021 [Link]

9: Putting It All Together
Distillation Science (a blend of Chemistry and Chemical Engineering)
This is Part IX, Putting It All Together of a ten-part series of technical articles on Distillation Science, as is currently
practiced on an industrial level. See also Part I, Overview for introductory comments, the scope of the article series, and
nomenclature.
Part IX, Putting It All Together uses the information from Parts III through VII, showing how to combine them in a
practical example for bulk separation; and how using the techniques detailed in previous articles give answers that differ
from the ideals of basic Raoult/Dalton Law application.
See previous articles per the following detail:
Part III for critical properties and acentric factor, especially Table 3-3
Part IV for the recommended vapor pressure equation, especially Equation 4-3 and Table 4-1
Part V for Equation of State, especially Equations 5-3 and 5-4
Part VI, especially Equation 6-4 for partial pressure and Equations 6-7 and 6-8 for fugacity coefficients
Part VII for Liquid Activity Coefficients, especially Equation 7-3
In this article a practical distillation example is developed, showing the Y/X relationship with temperature for the DCS-
TCS binary system with pressure is held constant at 11.0 atmospheres, and showing the results as a TXY/P plot as well as
in tabular form. From such data, a distillation column designer would consider the best arrangement for column feed,
“tops” and “bottoms” recovery, tray count and external reflux ratio ( and from that, the energy required for the reboiler and
condenser, as well as hydraulic loading per unit mass flow of feedstock). Two variations of this example are worked-up,
showing the application of distillation science: one which only uses ideal Raoult/Dalton, and one that uses the full set of
ideality departure coefficients. Table output and plots below are done using MS Excel. Then at the end of the article, the
table results are graphically demonstrated using a McCabe-Theile plot to step off the trays for an appropriate reflux rate
and show the best feed tray.
The table of Figure 9-1 is set up to cover the range of DCS liquid mole fraction from unity to zero, in increments of 0.05
mole fraction, with 0.01/0.99 mole fraction added to illustrate asymptotic effects. For each non-zero table row, a
temperature is first assumed, and then the vapor mole fractions calculated per Raoult/Dalton, with p=11.0 atmospheres;
then that row’s temperature is iteratively adjusted until the total pressure = 11.0 atmospheres.
In the plot of Figure 9-1 shows the calculation results in TXY/P form, with the value of XDCS=1.0 representing the tightest
possible “tops” product at the distillation column condenser; and a value of XDCS=0.0 representing the tightest possible
“bottoms” product at column reboiler, all at a total reflux condition. For a finite column, the mole fraction of the “tops”
product would be a cropping of the upper XDCS values of the TXY/P curves; and the mole fraction of the “bottoms”
product a would be a cropping of the lesser XDCS values of the TXY/P curves. In such instance, the column feed tray
would be identified as that point in the column where the tray temperature is identical to the feed mixture’s boiling point.
The accompanying plot of solution shows what the TXY/P curves would look like, for such a Raoult/Dalton solution to the
DCS-TCS binary at 11.0 atmospheres.
Y1 = (X1 × V P1 )/(X1 × V P1 + X2 × V P2 ) Y2 = (X2 × V P2 )/(X1 × V P1 + X2 × V P2 ) repeat Equation 6-1

Larry Coleman 9.1 10/31/2021 [Link]

 Figure 9-1, TXY calculation table and plot, Raoult/Dalton
If the data plot of Figure 9-1 plot is carefully examined, one can see that the T-X and T-Y curves are exactly symmetrical
when “folded” along a line running from abscissa/ordinate points of (0,124.76) to (1,92.40); or “folded” orthogonally to
that same line. In other words, there are no departures from ideality in either vapor or liquid phase, either with DCS or
TCS.
Now X and Y are re-calculated including the departures from ideality, both in fugacity (ø) as well as Liquid Activity
Coefficients ( γ ).

Figure 9-2, TXY/P calculation table, all departure factors included

The table in above Figure 9-2 is likewise set up and calculations performed row-wise for each value of XDCS, but the liquid
and vapor-phase fugacity coefficients (øL and øV) and activity coefficients ( γ ) are included in the calculation of vapor
i

mole fractions (γDC S ). For each row, the temperature is iterated until the sum of partial pressures = 11.0 atmospheres.
Note that since the activity coefficients are mild functions of temperature, each row has its own Van Laar constants
determined and the value of γ calculated for that row’s liquid mole fraction. Thus all departure-from-ideality factors are
included.
L V
Yi = (ϕ × γi × V Pi × Xi )/(ϕ × π) (9.1)
i i

which is an algebraic re-arrangement of Equation 6-4 from previous Part VI.

To better show how the several ideality departure factors impact the TXY/P plot (as compared to Figure 9‑1 above), the
table in Figure 9-2 is plotted. The symmetry of the Figure 9-1 Raoult/Dalton plot is no longer there, neither along the
“fold” line running from abscissa/ordinate points of (0,124.76) to (1, 92.40); or “folded” orthogonally to that same line.
Also note that the T-X and T-Y curves are closer together in Figure 9‑3, as opposed to Figure 9-1.

Larry Coleman 9.2 10/31/2021 [Link]

 Figure 9-3, TXY/P plot, all departure factors included
For a very simple approximation of tray count, for given distillation column “tops” product , “bottoms” product and feed
mixture mole fractions, and for a given reflux rate, a McCabe-Thiele plot is constructed in Figure 9-4 to graphically step
off the number of trays above and below the feed tray. For almost all industrial applications, this tray-by-tray analysis is
done via computer modeling, such as ASPEN, VMG, HYSYS, ChemCad, etc. To apply the technique to a computer
model, the user would specify the critical constants ( as per Part III), the vapor pressure equation to be used (as per Part
IV), the EOS to be used to calculate fugacity coefficients (as per Part V), and the Liquid Activity Coefficient model to be
used (as per Part VI).
To construct a McCabe-Thiele plot, the X-Y equilibrium (in red below) is plotted on a set of axes where X= liquid mole
fraction of the most volatile component = XDCS ; and Y= vapor mole fraction of the most volatile component = YDCS . The
total reflux line (in black below) connects (X,Y) = (0,0) and (1,1); representing the maximum reflux possible. The lines are
added (in dashed blue below) for the feed (XF), and distillate product (XD) and the bottoms product (XB), from the X-axis
to the total reflux line. The feed line is extended toward the equilibrium curve, and from that intersection a line is drawn to
the end of the distillate product line. The slope of that drawn line represents the minimum internal reflux needed to make
the separation (i.e., with an infinite number of trays). The upper column operating line (in magenta below) is then drawn
connecting the end of the distillate product line with the feed line, with a slope somewhat greater than the minimum reflux
slope. In normal practice the upper column operating line slope is usually about 1.3 times the minimum slope, depending
on the relative value of operating energy vs capital cost. Then the lower column operating line is drawn, intersecting the
bottoms product line with the upper column operating line (and the feed line). The number of theoretical trays (aka stages)
can then be stepped off ( in green below) from XD to XB, and best tray determined for the feed.

Larry Coleman 9.3 10/31/2021 [Link]

 Figure 9-4, McCabe Thiele solution to example calculations
To illustrate an example of using a McCabe-Thiele plot, assume the mixture to be separated is a 35% DCS molar liquid,
pre-heated to its boiling point of 108.61°C (as shown on the 35% of Figure 9-2). Further, assume that the column “tops”
distillate product is to be 98% molar DCS and the column “bottoms” product is to be 98% molar TCS. Per the text above,
the minimum reflux (i.e., the slope of the line between XD = XDCS of 0.98, XB = XDCS of 0.02) is determined to be 3.9:1,
reflux-to-distillate. To allow a reasonable number of trays, a reflux rate of 5:1 is chosen (about 1.3 times the minimum), so
the slope of the upper column operating line is R/(R+1) = 5/(5+1) =0.8333. Stepping off trays, it is determined that the
total number of theoretical trays required is 34, with the feed tray at #22. Conventional tray counting starts at the column
condenser = tray #1.
As can be seen by examining Figure 9-4, as the column’s reflux rate is increased, upper column (magenta) operating line’s
slope increases, and the gap opens up between the (red) X-Y Equilibrium Line and the (black) Total Reflux Line. So
increasing the reflux decreases the number trays required to make the specified separation.
In a normal industrial process application, all of the calculations are done via computer simulators, allowing for variations
in operating line slope; feed quality condition (e.g., saturated vapor, saturated liquid, sub-cooled liquid, etc); changes in the
internal liquid/vapor flow ratios due to physical properties; ambient heat effects; tray efficiencies, etc; as well as
integrating the distillation column into the rest of the process design. So the McCabe-Thiele method is just used for
graphical explanation.
However, the primary consideration of the distillation column design is the X-Y equilibrium relationship. If this is
considerably in error – such as using only Raoult/Dalton – the column design will be a failure. Given the cost of industrial-
scale distillation systems being in the millions of dollars (tens of $MM for larger ones), there is every reason to get the X-
Y equilibrium relationship correct.
In the example of this article, the basic Raoult/Dalton model would result in the X-Y equilibrium relationship in the above
example having more curvature (than shown in Figure 9-4) and therefore more open gap between the equilibrium curve
and operating lines. The resulting reflux determined would be about 2/3 of required, and number of trays required would
also be about 2/3 of needs - so a major error in design. Given what is at stake financially, there is no question that the use
of proper distillation science is well worth it, even if the concepts are complex and the calculations are tedious.

Larry Coleman 9.4 10/31/2021 [Link]

10: Convergence Strategy
Distillation Science (a blend of Chemistry and Chemical Engineering)
This is Part X, Convergence Strategy of a ten-part series of technical articles on Distillation Science, as is currently
practiced on an industrial level. See also Part I, Overview for introductory comments, the scope of the article series, and
nomenclature.
This last article discusses how to approach solving non-intrinsic equations and equation systems, since these are part of the
practice of distillation science. While it can be argued that this topic is technically more related to math or computer
science, the techniques are indispensable tools to obtaining workably accurate solutions.
Most relationships in Distillation Science feature equations that are written with temperature as the starting variable, and
then one solves for the other parameters. For example, almost all vapor pressure relationships are written intrinsically with
VP as a function of temperature. But frequently the problem statement begins with a fixed pressure, and so there needs to
be an iterative procedure to solve for temperature.
More complex relationships involving two starting variables are written in temperature and pressure, with the equation’s
dependent variable written as a function of (absolute) temperature and (absolute) pressure. An example is an Equation of
State, where Z is written as f(T,P), with molar volume (V) being back-solved from the defining relationship PV=ZRT.
Some scientific relationships are written in dimensionless “reduced” variables, where T, P, and V are divided by their
critical point value (e.g., Tr=T/Tc, Pr=P/Pc), which allows for more global equation forms. The Peng-Robinson Equation of
State is an example of using reduced properties, although Equation 5-4 partially masks that by breaking the equation into
pieces.
The most basic technique for solution convergence in a non-intrinsic relationship is based on “gaming strategy”: (1) pick a
possible independent variable’s value and work it through to a solution, comparing it to the desired final dependent
variable; then (2) move the value choice a little bit and see how that helps. This is just a methodical form of “trial and
error”, and will eventually lead to a solution. However, such strategy is almost impossible to automate for easy
programming, as in a spreadsheet macro or a stand-alone BASIC program. Moreover, it can be hard to get started, may
lead to an unstable situation that cycles around a solution rather than stably converging to the required accuracy.

Newton-Raphson Method

A better convergence option is a modification of Newton-Raphson, aka Newton’s Method. In the instructions below,
"X" is the relationship's independent variable and "Y" is the dependent variable. The problem is that you have a
desired target value for "Y" and have to iteratively seek the correct value for "X", and there is no way to re-write the
relationship so that "Y" is the independent variable. The increasing subscripts denote the progression at improving the
solution, such that "Yn" is identical to the desired target value of "Y" = "YT", and "Xn" is the value for "X" that yields
that solution.
1. If there is a way to get a very rough estimate of solution, start with that as the first trial = X0, and work the
scientific relationship through to a first possible solution, Y0. (If there is no rough estimation possible, start with a
known “safe” initial X0 and calculate Y0.) Then make an extremely small ΔX step (say 0.1% change in X0). In the
next iteration, X1=X0+ΔX would be evaluated to get Y1. If the change in ΔY (i.e., ΔY= Y1-Y0) is too
inconsequential to be seen, increase the value of the ΔX step, until there is some noticeable change in ΔY. (This is
an important consideration since some relationships are complex, so it is hard to pick that first ΔX step to get a
reasonable change in Y). It does not matter whether the change to X (i.e., ΔX) is positive or negative: that will get
worked out later.
2. Calculate the quantity m = (ΔY /ΔX) = (Y − Y )/(X − X ) , which approximates the scientific
1 0 1 0

relationship’s first derivative. “m” can be positive or negative (which is why the sign in ΔX didn’t matter). A
positive value for “m” indicates that Y increases with increased X; or a negative value for “m” indicates the
opposite effect. Take notice of the sign of “m”, to see if its sign makes common sense.
3. With YT as the target solution value, make the next guess for X as
X2 = X1 + K × (YT − Y1 )/m (10.1)

10.1 10/31/2021 [Link]

 4. The variable “K” is the relaxation factor, with values of “K” closer to unity giving faster convergence, but at the
peril of overshoot or convergence instability; and “K” values closer to zero giving slow but stable convergence.
Until you know something about the scientific relationship, it is recommended using a “K” value of about 0.2-0.3
at first, especially when logarithmic relationships are involved. You can always change the “K” value once you see
how sensitive the scientific relationship is.
5. For the “ith” iteration, the iteration equation becomes:
Xi+1 = Xi + K × (YT − Yi ) × (Xi − Xi−1 )/(Yi − Yi−1 ) (10.2)

6. Continue iterating "n" times until the solution, Yn, is adequately close to the target solution, YT. It is suggested that
loop iteration accuracy should be one decimal place greater than that of the final solution’s desired precision, to
allow for round-off error and computer calculation precision.

When doing automated computer calculations, either via spreadsheet macro or stand-alone BASIC program, it is a good
idea to error-check each iteration for division by zero, and to limit the number loop iterations (say to a high number like
100). Otherwise, the iteration program can get “hung up” in an infinite loop. If you are getting essentially zero value for
the quantity (Y − Y ) or requiring an excessive number of loops to converge, that requires adjusting the relaxation
n n−1

factor, “K”, up or down to adjust convergence sensitivity.

Some ideas for making that first rough estimate, to get the iteration started, are:
For finding the saturation temperature for a known pressure (i.e., inverting a vapor pressure relationship), start with
using a basic VP relationship (inaccurate, but simple), such as Ln(P )α − ΔH /RT using the known ambient
vap

boiling temperature, Tb @ P=1 atmosphere and a reasonable estimate for latent of vaporization,ΔH . Then the rough
vap

starting estimate is T = T − ΔH /[R × Ln(P )] , where PT is the known pressure target, in atmospheres. This will
0 b T

give a high estimated T0 for PT > atmospheric, and a low estimated T0 for PT < atmospheric. It is more important to
start with a "safe" estimate than close first iteration, to avoid the pitfalls of automated solutions mentioned above.
Vapor pressure charts and tables are good ways to get started, although sometimes hard to program into an automated
sequence.
When iterating for the saturation temperature of a mixture, start with the assumption of a linear relationship between
the two fluids’ saturation temperatures at that pressure; then interpolate between those two pure-component
temperatures per the mixture’s mole fraction. There is always curvature between two fluids’ saturation temperatures, so
assuming linearity isn’t correct; but it gets the iteration loop started.
Starting an iteration loop for fugacity coefficients should only be done after a cursory look at the reduced pressure, Pr,
since these relationships involve calculating compressibility, Z. For liquid fugacity coefficients, avoid starting out or
getting too close to Pr=0.9, because the automated solution to cubic equations gets too sensitive, and one of the roots
may become an “imaginary number” (i.e., a number times the square root of one). For vapor fugacity coefficients, do
an initial check to see which Z = f(Tr, Pr) equation is best to use (i.e., Equation 6-6 or Equation 6-7): one is better for
lower values of Pr and one better for higher values. If you attempt to cross over between equations in the middle of a
solution, that discontinuity might cause an automated iteration to “blow up”.
Resist the temptation to use Newton-Raphson to get the roots for the cubic relationship Z: there are three roots, and
only two are scientifically meaningful. Besides, you need to know all three roots in order to judge which is Zv (the
largest of the three) and which is ZL (the smallest of the three). Instead, go through the procedure in Part V and solve
the cubic equation for Z, using the exact trigonometric method of Equation 5-4.
Starting the loop for activity coefficients should never be done at either extreme of mole fractions (zero or unity). The
Liquid Activity Coefficient's relationships are complex-shaped and can have γ/ΔX values that are almost zero (near
X=1) or very high (near X=0). For bulk separations, start somewhere between 0.3<X<0.7, and use a moderate value for
relaxation factor, “K” (say 1/3) in Equation (10.2) . For trace impurity distillation systems, either start closer to X=0.1
or at X=0.9, plus use a conservative relaxation parameter, say K=0.1.
Convergence Strategy for modified Thek-Stiel VP Equations
The “sting” of working out nested loop iterations has been removed by giving the solutions for the modified Thek-Stiel VP
equation (Equation 4-3) in tabular form (Table 4-1). Obtaining vapor pressure solutions for fluids other than chlorosilanes
and their impurities (see Part IV) require more complex convergence strategies than Newton-Raphson.

10.2 10/31/2021 [Link]

These VP solutions involve a system of nested-loops: the value of the Watson coefficient “q” depends on evaluating the
candidate vapor pressure equation at Tr=0.7 to get the loop’s estimate of acentric factor; and obtaining the loop’s estimate
of “c” and “n” depends on finding a value of αc that gives the correct saturation temperature at ambient pressure (aka, the
normal boiling point). The following solution strategy is advised below, for fluids other than those of Table 4-1 (or in case
one wants to check the table’s values).
First, be sure all of the values for Tb, Tc, Pc are well-validated, and the degree of hydrogen bonding is known, per the
definitions in Part IV. If the new fluid to be investigated is of a periodical table group other than III, IV, or V, the solution
will involve a more intricate solution than below, since “k” (i.e., the degree of hydrogen bonding) is not well determined
for other than Groups III-V, and a double iteration will be needed. Also be sure the vapor pressure data is also well-
validated. A good check on that is to use a simplistic extrapolation of the known VP data, to see if the rough value of the
acentric factor matches with the values of Tc and Pc. See Part III for more details on acentric factor, and especially
Equation 3-3 For example, if an extrapolation of VP data using the Antoine VP relationship (Equation 2-3) in conjunction
with values for Tc and Pc , ends up giving an acentric factor value over 0.3 for a compact molecule, or an acentric factor
under 0.1 for a complex molecule, something is wrong with the data (either the VP data or the values of Tc and Pc).
For Equation 4-3, the value of the hydrogen bonding parameter, k, is established by experimental data for Group III, IV,
and V hydrides, chlorides, oxychlorides, intermediate hydrochlorides, and several alphatic hydrocarbons. The value of k
ranges from 0.12 to zero, and is tabulated in Part IV for many fluids. The rules for determining "k" are a bit complex, and
are available upon request from the author. In brief, for Group III, IV, and V hydrides, k=0.1197; 0.0914; and 0.1009
respectively. For all fully chlorinated or methylated compounds of Groups III, IV, and V , k = 0.0291. For aliphatic
hydrocarbons, k = 0. For other hydrochloride intermediates, a quadratic function has been determined to interpolate
between all-hydride and all-chloride. These rules for "k" values are used in the entries of Table 4-3.
Having a rough estimate of the ambient pressure latent heat of vaporization is handy (but not necessary) for the first guess
at the T-S heat variable. If no estimate is available, consult one of the many handbooks of chemical compounds, such as
Yaw’s Chemical Property Handbook, or Lange’s Handbook of Chemistry and look for a latent heat value for a compound
of similar complexity. It is always better to start a new fluid’s solution to the Thek-Stiel VP equation with a low estimate
for the T-S heat variable, and then increase.
The overall solution convergence strategy to the modified Thek-Stiel VP equation involves a triple-nested iteration loop,
with the logic flow sheet of the strategy shown on Figure 10-1. The solution is obtained by improving candidate solutions
for the T-S heat variable, the acentric factor (ω ), and the αc parameter, until all loops are converged.
The inner-most iteration loop is the α parameter (a thermodynamic consistency parameter from Part IV) , which is
c

adjusted until the VP equation gives the correct saturation temperature at atmospheric pressure, aka the normal boiling
point Tb. That loop is started by using Reidel’s estimation method: calculating the rough value of
ψ = (−35 + 36/ T
b + 42LnT
br −T
br , and then plugging in that rough value of ψ in the equation below:
6
br b

αc ≈ (0.315 ψb + LnPc )/(0.838 ψb + LnTrb ) (10.3)

The middle of the nested loops is the acentric factor, which can be conveniently started using the Lee-Kister acentric factor
(ω) approximation (which is not very accurate, but a good loop-starting estimate):
The Lee-Kister approximation for acentric factor (ω) conveniently uses only Pc and Tbr=Tb/Tc :
6
−lnPc + A + B/ Tbr + C ∗ LnTbr + D ∗ T
br
ω ≅ (10.4)
6
E + F / Tbr + GLnTbr + H ∗ T
br

A = -5.92714 E = 15.2518
B = 6.09648 F = -15.6875
C = 1.28862 G = -13.4721
D = -0.169347 H=0.43577
Based on this rough estimated value of ω, a value of Watson’s coefficient “q” (from Equation 4-3) can be determined per
the relationship

10.3 10/31/2021 [Link]

 q = 0.37028 + 0.065404 ∗ ω (10.5)

The outer-most loop is the T-S heat variable, “T-S ΔHv”, which doesn't really have much scientific meaning, but is a part
of the Thek-Stiel relationship. This parameter is iterated based on the error function comparing the VP equation’s output
against the various known ( or estimated) VP data points. To get the best overall fit, I suggest including data near Tr = 0.7
(where the acentric factor is evaluated). In getting the various T-S VP Equation solutions shown on Table 4-1, I chose to
not automate this outer loop, but to rely on the more basic “gaming strategy” ⇒ always starting low, so as to monitor the
convergence better. For computer automation, I used a BASIC program (actual software used was Liberty Basic version
4.03, selected due to prior familiarity). The starting guess for the outer loop iteration was based on the atmospheric ΔHvap
given in chemistry handbooks and from internet sources, or based on similar compounds where no published estimate of
atmospheric ΔHvap was available. This choice of starting guess was based on the "safe guess" criterion, and opposed to a
scientific estimation.
Figure 10-1 below schematically gives the convergence strategy previously used with success. After insuring that all input
parameter and data are reasonable, start with the initial guess for the T-S heat variable, “T-S ΔHv”. Then estimate a starting
ω per Equation (10.4); a from that estimate a starting α per Equation 10-3. Now calculate values for constants “A”, “B0-
c

B3”, “c”, and “n” in Equation 4-3, and see if the inner-most loop’s convergence is satisfied. If not, iterate the value of α c

until the correct VP is shown at Tb of 1.0 atmosphere (i.e., the normal boiling point). This inner-most loop is best
converged by approaching the VP @ Tb from P<1 atmosphere, with convergence to ± 1.0E-4 atmospheres being seen
within about ten iterations via Newton-Raphson, using a relaxation factor of 0.4 in Equation (10.2). Note that if the first
guess of α shows a VP @ Tb of >1.0 atmospheres, just keep reducing the value of α by 10% until the convergence
c c

approach is from lower values of Tb. Trying to converge from the “high side” can be unstable.
Once the inner-most loop is converged, start changing the value of the acentric factor (ω) to converge the middle loop.
This loop converges quite easily with just simple continuous substitution. Newton-Raphson is not required. Just put in the
last value of ω, calculate the candidate VP constants (“A”, “B0‑B3”, “c”, and “n” in Equation 4-3) and determine a new
value for ω based on the VP at Tr=0.7 (see Equation 3-3 for the definition of ω).
The outer loop is converged by picking that value of the T-S heat variable, “T-S ΔHv”, that best satisfies the error function
of the VP data (i.e., predicted VP's as compared to the known data point VP's). I found the best overall fit used a least
squares approach which emphasized the data near Tr=0.7, and which compared the VP equation error on a relative level
(as opposed to an absolute level). If absolute VP equation errors are used (as opposed to relative ones), that over-
emphasizes data at higher temperatures.

10.4 10/31/2021 [Link]

10.5 10/31/2021 [Link]
Index
A E L
acentric factor Equation of State liquid activity coefficients
5: Equation of State 9: Putting It All Together
3: Critical Properties and Acentric Factor
Antoine relationship
2: Vapor Pressure
F V
fugacity vapor pressure
C 6: Fugacity 2: Vapor Pressure

compressibility Fugacity coefficients

5: Equation of State
6: Fugacity W
Watson’s coefficient
10: Convergence Strategy
Glossary

Sample Word 1 | Sample Definition 1