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., T_r=T/T_c, P_r=P/P_c), 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 "Y_n" is identical to the desired target value of "Y" = "Y_T", and "X_n" 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 = X₀, and work the scientific relationship through to a first possible solution, Y₀. (If there is no rough estimation possible, start with a known “safe” initial X₀ and calculate Y₀.) Then make an extremely small ΔX step (say 0.1% change in X₀). In the next iteration, X₁=X₀+ΔX would be evaluated to get Y₁. If the change in ΔY (i.e., ΔY= Y₁-Y₀) 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=(\Delta Y/ \Delta X) = (Y_{1} -Y_{0})/(X_{1} -X_{0}) \nonumber\), which approximates the scientific 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 Y_T as the target solution value, make the next guess for X as \[X_{2} =X_{1} + K \times (Y_{T} -Y_{1})/m \label{10-1} \]
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: \[ X_{i+1}=X_{i} +K \times (Y_{T} -Y_{i}) \times (X_{i} -X_{i-1}) / (Y_{i} -Y_{i-1}) \label{10-2} \]
6. Continue iterating "n" times until the solution, Y_n, is adequately close to the target solution, Y_T. 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_{n} -Y_{n-1}) \nonumber\) or requiring an excessive number of loops to converge, that requires adjusting the relaxation 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) \alpha - \Delta H_{vap}/RT \nonumber\) using the known ambient boiling temperature, T_b @ P=1 atmosphere and a reasonable estimate for latent of vaporization,\(\Delta H_{vap} \nonumber\). Then the rough starting estimate is \(T_{0}=T_{b}- \Delta H/[R \times Ln(P_{T})] \nonumber\), where P_T is the known pressure target, in atmospheres. This will give a high estimated T₀ for P_T > atmospheric, and a low estimated T₀ for P_T < 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, P_r, since these relationships involve calculating compressibility, Z. For liquid fugacity coefficients, avoid starting out or getting too close to P_r=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(T_r, P_r) equation is best to use (i.e., Equation 6-6 or Equation 6-7): one is better for lower values of P_r 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 Z_v (the largest of the three) and which is Z_L (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 \(\gamma/ \Delta X \nonumber\) 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 (\ref{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.

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 T_r=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 T_b, T_c, P_c 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 T_c and P_c. 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 T_c and P_c , 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 T_c and P_c).

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 \(\alpha_{c} \nonumber\) parameter (a thermodynamic consistency parameter from Part IV) , which is adjusted until the VP equation gives the correct saturation temperature at atmospheric pressure, aka the normal boiling point T_b. That loop is started by using Reidel’s estimation method: calculating the rough value of \(\psi_{b}= (-35 +36/T_{br} + 42LnT_{br} -T_{br}^6 \nonumber\), and then plugging in that rough value of \(\psi_{b} \nonumber\) in the equation below:

\[ \alpha_{c} \approx (0.315 \psi_{b} + LnP_{c})/(0.838 \psi_{b} + LnT_{rb}) \label{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 P_c and T_br=T_b/T_c :

\[\omega \cong \frac{-lnP_{c}+A+B/T_{br} +C*LnT_{br} +D*T_{br}^6}{E+F/T_{br} +GLnT_{br}+H*T_{br}^6} \label{10-4} \]

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

\[ q = 0.37028 + 0.065404*ω \label{10-5} \]

The outer-most loop is the T-S heat variable, “T-S ΔH_v”, 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 T_r = 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 ΔH_vap given in chemistry handbooks and from internet sources, or based on similar compounds where no published estimate of atmospheric ΔH_vap 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 ΔH_v”. Then estimate a starting ω per Equation (\ref{10-4}); a from that estimate a starting \(\alpha_{c} \nonumber\) per Equation 10-3. Now calculate values for constants “A”, “B₀-B₃”, “c”, and “n” in Equation 4-3, and see if the inner-most loop’s convergence is satisfied. If not, iterate the value of \(\alpha_{c} \nonumber\) until the correct VP is shown at T_b of 1.0 atmosphere (i.e., the normal boiling point). This inner-most loop is best converged by approaching the VP @ T_b 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 (\ref{10-2}). Note that if the first guess of \(\alpha_{c} \nonumber\) shows a VP @ T_b of >1.0 atmospheres, just keep reducing the value of \(\alpha_{c} \nonumber\) by 10% until the convergence approach is from lower values of T_b. 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”, “B₀‑B₃”, “c”, and “n” in Equation 4-3) and determine a new value for ω based on the VP at T_r=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 ΔH_v”, 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 T_r=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.