## Solving Equations Fractional Distillation Method -sndholm Naphtali

X.N

where x, represents all of the variables of stage j and f, represents all of the functions of stagey. Because of this block form of the Jacobian, the computer time and storage is not great, possibly less than that required by a Tomich method. All of the global Newton methods have a sparse Jacobian matrix. About 95 percent of the matrix is zero's and it is quite simple to store only the nonzero elements. There are many numerical methods routines available specifically for solving sparse Jacobians. The Naphtali-Sandholm procedure is as follows

1. Set initial temperatures and total vapor and liquid rates for each stage. Calculate initial component vapor rates using the tridiagonal matrix method and find the component liquid rates by applying the absorption factor, /y- = Ay

2. Calculate all of the independent functions and the norm.

3. Fill and solve the sparse Jacobian for a new set of temperatures, Tj's, and component vapor and liquid rates, fy's and /y's, by solving the independent equations during one pass through the Newton-Raphson procedure.

4. The norm of the independent functions must be very small at the solution (see Sec. 4.2.6 for criteria). If it is not small or other solution criteria (Sec. 4.2.3) have not been met, return to step 2.

Note that there are no other side calculations between passes through the Newton-Raphson procedure. Total flow rates and duties are calculated after a solution of the column is found.

Christiansen et al. (54) applied the Naphtali-Sandholm method to natural gas mixtures. They replaced the equilibrium relationships and component vapor rates with the bubble-point equation and total liquid rate to get practically half the number of functions and variables [to N(C + 2)]. By exclusively using the Soave-Redlich-Kwong equation of state, they were able to use analytical derivatives of K-values and enthalpies with respect to composition and temperature. To improve stability in the calculation, they limited the changes in the independent variables between trials to where each change did not exceed a preset maximum. There is a Naphtali-Sandholm method in the FraChem program of OLI Systems, Florham Park, New Jersey; CHEMCAD of Coade Inc. of Houston, Texas; PRO/II of Simulation Sciences of Fullerton, California; and Distil-R of TECS Software, Houston, Texas. Variations of the Naphtali-Sandholm method are used in other methods such as the homotopy methods (Sec. 4.2.12) and the nonequilibrium methods (Sec. 4.2.13).

The Goldstein-Stanfield (45) method. This method chooses the stage temperatures, total vapor rates, total liquid rates, and liquid compo sitions of each stage as the independent variables. Different from Naphtali-Sandholm, vapor compositions are not among the independent variables. Like Naphtali-Sandholm, the independent functions include normalized energy and component balances [Eqs. (4.77) and (4.78)]. In addition, a summation equation [Eq. (4.71)] and a normalized total material balance are used, giving a total of N(C + 3) equations for the column. Also, Goldstein and Stanfield arranged the functions by equation type with the component material balances arranged by component and not by stage. Vector of independent functions is of the form:

£ = [mn.,.nilN...mcl...mCNH1...HNSl...SftM1...MN]'!: (4.83) with the corresponding set of independent variables:

X* = [*„.. .xw.. .xa ...xCNTl...TNVl...VNLl.. ,LNf (4.84)

Like the Naphtali-Sandholm, the Jacobian is mostly zeros, although it is not a block-banded matrix. The Goldstein-Stanfield method uses a series of matrix manipulations to first find the changes in the temperatures and total vapor rates, the ATj's and AV}'s. These are then used to find the changes in the liquid rates, AL/s, and the actual liquid compositions, the x/s. The steps of the Goldstein-Stanfield method follow those of Naphtali-Sandholm, except there are these unique matrix manipulations in solving the independent equations. Like Naphtali-Sandholm, the condenser and reboiler duties and side-product rates are specified. The Goldstein-Stanfield method is included in the COPE simulation program of Exxon, the MAXPLUS program of Kesler Engineering, Highland Park, New Jersey, and the DS03 program of Paul Oleson and Associates, Laguna Beach, California. Bondy (59) extended the Goldstein-Stanfield method to nonideal three-phase systems for the CHEMDIST method in the PRO/11 program of Simulation Sciences of Fullerton, California.

While Goldstein and Stanfield (45) did not include derivatives for K-values and enthalpies in their original formulation, these can be readily included. This can improve the speed of the method with reasonably ideal mixtures. Goldstein and Stanfield (45) and Boston (75) stated that the Naphtali-Sandholm approach of grouping the functions by stage is more suitable for columns with many stages and few components. Because both the component vapor and liquid rates are among the independent variables, the Naphtali-Sandholm approach should also work better for highly nonideal mixtures. The Goldstein-Stanfield approach of arranging functions by component, the tendency to ignore derivatives of if-values and enthalpies with respect to composition, the choice of independent variables, and the matrix manipu lations should be better for columns with few stages and many components, and has been shown to work well for refinery fractionators.

The Ishii and Otto method. The method of Ishii and Otto (47) follows an approach similar to the Goldstein-Stanfield method. The energy balance, summation equation for the liquid compositions, and the component balances are the independent functions. These equations are grouped by type, and as in the Goldstein-Stanfield method, adjustments to the temperatures and total vapor rates, ATj's and AV^'s, are found first after the series of matrix manipulations. The liquid compositions for each stage are then found using the total rates and temperatures. While Ishii and Otto did include derivatives of if-values and enthalpies, they were based only on the simple Chao-Seader correlation and with several approximations in the derivatives. The method may reduce computer resource needs, but since the compositions calculation is split from the temperatures and total rates calculations, and due to the approximations in the derivatives, it may have difficulty solving nonideal systems. Holmes et al. (109) successfully applied the Ishii-Otto method to amine strippers in the program TSWEET of Bryan Research and Engineering, Bryan, Texas. The Ishii-Otto method is also used by the AMSYM program of D. W. Robinson and Associates of Edmonton, Alberta,

The almost band algorithm. This method of Gallun and Holland (40) and Holland (8) is unique in its flexibility, wide combination of equations for solving different types of columns, and for the strong effort in developing sparse Jacobian matrix solution techniques. It follows the Naphtali-Sandholm approach of grouping the equations by stage in a block-banded form. It allows for mixing and matching, adding and eliminating variables and equations to fit the application, and solves for all variables simultaneously without the matrix manipulations of Goldstein-Stanfield.

Gallun and Holland bring activity coefficients into the independent functions to separate them from the calculation of the ideal mixture /¿"-values. Presenting the activity coefficients directly in the equilibrium equation emphasizes how they and their derivatives cannot be and are not ignored when considering nonideal mixtures.

The method of Gallun and Holland is the broadest application of the MESH equations in a global Newton method and may solve the widest range of columns. Formulations by Gallun and Holland (40) for distillation columns included adding the total material balance to give freedom in specifications or to substitute these for the equilibrium equations for more ideal mixtures.

Gallun and Holland (40) also presented solution techniques for the sparse matrices of the global Newton method. There are many other solution techniques available. The Jacobian matrix of these methods contain values in only 3 to 10 percent of the elements. To save storage and allow the solution of tall columns and mixtures with many components, none of the zero elements should have to be stored or manipulated.

The matrix solution techniques of the block-banded formulations of Naphtali and Sandholm (42) and of Holland (8) are generally simpler than that of the other global Newton methods. Also, the Naphtali-Sandholm and almost band methods are better suited for nonideal mixtures than other global Newton methods.

### 4.2.10 Inside-out methods

The inside-out algorithm has become one of the most popular methods because of its robustness and its ability to solve a wide variety of columns. The inside-out concept was developed by Boston and Sullivan (69) and Boston (70). Russell (72) presented an inside-out method that works well for many refinery fractionators. Jelinek (73) presented a simplified Russell method.

In the previous four groups of methods (Sees. 4.2.5 and 4.2.7 to 4.2.9), the MESH variables of temperatures, total flow rates, and component flow rates are the primary solution variables and are used to generate the A'-values and enthalpies from complex correlations. These methods update the MESH variables in an outer, though sometimes partitioned, loop with the if-values and enthalpies updated whenever the MESH variables change. The inside-out concept reverses this by using the complex it-value and enthalpy correlations to generate parameters for simple if-value and enthalpy models. These parameters are unique for each stage and become the variables for the outside loop. The inside loop consists of the MESH equations and is a variation on the BP, SR, and 2N Newton methods. Since the if-values and enthalpies are simple, the inside loop works well for a wide range of mixtures and is little affected by the nonideality of mixtures. In every step through the outside loop, the simple models are updated using MESH variables (mostly temperatures and compositions) from the inside loop. This sets up the next pass through the inside loop.

The Boston method. The original Boston and Sullivan (69) inside-out algorithm has evolved to that presented by Boston and Britt (71), Boston (70), and Trevino-Lozano et al. (74). There have been more advances to the Boston method which were not reported in the literature. The concepts of outside loop if-value and enthalpy models with the inside-loop column solution is preserved throughout.

The outer-loop A"-value model is based on the Kb method (Sec. 4.2.5):

where T* is a reference temperature for the K-value correlation. Outer-loop variables, A; and Br are generated for each stage from a reference Kbj Re{ of a composite component:

where the u>,'s are weight factors found by a method in Boston and Britt (71). The temperatures and compositions used to get the ify-;actual are the latest from the inside loop. Simple relative volatilities are among the outside-loop variables, and are used in the Kb method to calculate the temperatures and whenever if-values are needed in the inside loop

These simple relative volatilities change little over the range of temperatures that is seen on a given stage and greatly simplify tempera-:ure and composition calculations in the inside loop. For nonideal mixtures, an activity coefficient for each component accounts for composition effects in the inside loop. This activity coefficient has a simple model, similar to the Kb model:

where the new outer-loop variables, av and for each component are determined from the actual activity coefficient model at the current stage temperature and stage composition. The simple if-values used in the inside loop are easily determined from

The activity coefficient model is dropped for ideal or nearly ideal mixtures.

The enthalpy models are actually for the enthalpy departure function. Normally, the molar enthalpy of a phase is found from the ideal gas enthalpy and the enthalpy departure

where H* and h° are the ideal gas enthalpies based on the vapor- and liquid-phase compositions and AHVJ and AH^ are the enthalpy departures calculated by the enthalpy correlation. The simple enthalpy departure used in the inside loop is a straight-line function:

where T* is a base temperature. Thus, the outside-loop calculation consists of updating the terms of the simple if-value, activity, and enthalpy models. In some systems, such as ammonia-water, the enthalpy model gives too crude an approximation and is therefore inappropriate. In those systems, the enthalpy model is not used and enthalpies of the inside loop are calculated directly by the enthalpy correlations.

The inside loop consists of the actual calculation of the MESH variables and is similar to one of the BP or SR methods. For the inside loop, Boston defines two new variables, R, and ptJ:

where FXt is for component i in all feeds. The term Rj is a combination of temperature (through Kb) and the ratio of the total vapor and liquid rates. This combination frees the method from some sensitivity to wide- or narrow-boiling mixtures. From Rj and the compositions can be calculated directly:

t = 1 i - 1 The temperatures are calculated from a Kb method:

with the new temperature calculated using the K-value model:

2j (In Kty - Aj)/Bj + 1/T* Total liquid rates are also found using the Rjs:

and total vapor rates are found by total material balances. The energy balances, using the temperatures and compositions found in the pre vious steps, form the equations of a Broyden's method (Sec. 4.2.6). Their solution is used to update the Rj's.

The outer-loop variables, the a^'s A/s, B/s, ay's, 6y's, Cjs, DJs, Ejs, and Fj's, are updated after each inside-loop solution using the latest temperatures and compositions from the inside loop. Using the new KbjtRef's, the inside-loop variables, the R/s, are reinitialized and the inside-loop calculation begins again.

One major problem with all column methods is the initial values. Boston uses a set of scale factors that force a component and total flow rate distribution throughout the tower that conform to product specifications. The distribution is forced by multiplying the component flow rates by a base stripping factor for each stage:

where Sb is a scalar multiplier on the base stripping factors that is changed before the start of every column trial. This scalar forces the distribution of components and flows up or down the column in order to prevent ne gative or zero flow rates and to have the sum of the component flow rates approach a specified or estimated total flow rate. Once the flow rates are adjusted, the temperature profile can be smoothed or generated by using the flow rates to calculate dew or bubble points on selected stages.

Boston (75) added a middle loop to allow for column specifications and constraints. The arrangement of equations in the inner loop, where the solution of the MESH variables occurs, allows for only a few control or specified variables, such as fixed reflux ratio and product rates. The middle loop adjusts the control variables to meet the specifications. The middle loop can be built as an optimization method with process specification equations and economic objectives and constraints.

The calculation sequence of the Boston method is as follows:

1. Set initial temperatures and total rates for each stage.

2. Initialize the outside-loop variables; the reference base K-value, Kbj Ri(\ the relative volatilities, ay; and parameters Aj and B, of the lvalue model, parameters a,; and bi} of the activity coefficient model and parameters Cj, Dt, Ej, and F, of the enthalpy models, using the actual iT-value and enthalpy correlations and the estimated set of temperatures and compositions.

3. Adjust the component flow rates using the base stripping factor of Eq. (4.96). Update the total flow rates and temperatures with these adjusted component flows.

4. Initialize the inside-loop variables, Rj's, based on the current set of flow rates and reference A'6j Ret-'s.

The inside loop consists of steps 5, 6, 7, 8, and 9.

5. Calculate the inside loop py's, Kbj's, x,/s, and y,/s.

6. Get new stage temperatures by the Kb method.

### 7. Calculate the total liquid rates using the RjS.

8. Get a new set of Rj s by solving the energy balances as the independent functions of the quasi-Newton technique of Broyden (Sec. 4.2.6). The derivatives of the Jacobian matrix are generated numerically and must include steps 5, 6, and 7 each time the independent variables are perturbed. The Jacobian is not recalculated after the first trial in the loop.

9. Any calculation of the independent functions of step 8 must include a pass through steps 5, 6, and 7. The norm of the functions with the new RyS should be small to leave the loop (see Sec. 4.2.6 for criteria). If not, return to step 5.

10. Update the outside-loop variables; the reference base if-value, Kbj,Ref> the relative volatilities, ct^; and parameters Aj and Bj of the K- value model, parameters a(j and b^ of the activity coefficient model and parameters C]t Dp Ej, and F, of the enthalpy models, using the complex if-value and enthalpy correlations and the current set of temperatures and compositions.

11. The test for solution is that the relative volatilities and the parameters of the /f-value and enthalpy models generated in the outer loop using the complex correlations have not changed greatly between outer-loop trials. If this or other solution criteria (Sec. 4.2.3) have not been met, return to step 3.

From the equations, it may appear that the Boston method is most appropriate for narrow-boiling mixtures. However, the forcing style of the method allows it to also work well for wide-boiling mixtures. The Boston method has been found to be the better inside-out method for tall, high-purity (superfractionator) types of columns. The Boston method was extended to absorbers by Trevino-Lozano et al. (74). There, the inside loop was changed to have a base stripping factor as the independent variables in a Broyden's method solution of the energy balances. These base stripping factors were also used in a tridiagonal matrix solution of the component balances, while the Kb method for calculating the temperatures remained the same. These modifications to the inside-out method can and have been extended to distillation columns. In later variations of the Boston method, the tridiagonal matrix method of Boston and Sullivan (25) replaces Eq. (4.96; to calculate the compositions. Also, for better solution of absorb ers and strippers, a sum-rates method can replace the inside loop. Boston's method was extended to three-phase distillation by Boston and Shah (76) and to reactive distillation by Venkataraman et al. (77). A well-supported Boston method with a wide variety of features, options, and power to solve a wide range of columns is in the RADFRAC and the MULTIFRAC methods in ASPENPlus of ASPENTech of Cambridge, Massachusetts.

The Russell method. Russell (72) follows the same structure in the outer loop as the Boston method except the base temperature in the if-value model is removed from the equation:

The reference Khj is generated as in the Boston method [Eq. (4.86)] but the weight factors, wj/s, have been simplified to be based on the change in the K-values with respect to temperature:

where the y,/s, temperatures, and compositions used to get the A'y actuaJ's are the latest available in the calculation. The relative volatilities are still found from Eq. (4.87) and the a,jS, A/s, and Bj's are among the outer-loop variables. Russell also has an activity coefficient model similar to Boston, and Russell's method uses Boston's simple enthalpy models.

Russell's method differs from Boston's in the inside loop by calculating the component flow rates by the tridiagonal matrix method and the total flow rates by summing the component rates. The temperatures are calculated by Russell's version of the Kb method. As in Boston's method, one energy balance is drawn per stage, for a total of N equations per column, and these are the independent functions of a Broyden's method solution. Additional equations can join the independent functions to support specifications for product quality, stage temperatures, internal flow rates, etc. Independent variables of Russell's inside loop matrix are logarithms of base stripping and withdrawal factors:

.actual'

.actual'

In Sbl = IntKy LjiVj)

In RVJ = ln(Wy/V}) for a vapor side product In Rl, = ln(Wu/Lj) for a liquid side product

For every side product, a withdrawal factor is added to the independent variables. With it, a specification function joins the energy balances in the independent functions. These and the energy balances are the functions of a Newton-Raphson method, and whenever they are calculated (as in filling the Jacobian), a pass must be made through the calculations of the compositions and temperatures.

The outside-loop variables; the relative volatilities, a^'s; the constants of the Kb model, Ay's and BJs; and the constants of the enthalpy models, C7 and Dj and Ej and Fj are not updated except in the outer loop. The factors, SbJ's, RVj's, and Rvs, are initialized at the start of each inside loop. The base stripping factors enter in inside-loop solution where they are needed to get the stripping factors of the component material balances:

Russell organizes the tridiagonal matrix method to calculate the component liquid rates instead of the vapor rates (as in Sec. 4.2.3) but either can be used. The component vapor rates are found by v,j = Stl l!r The total flow rates are found by summing the component flow rates:

The inside loop KbJ is found using the outside-loop relative volatilities:

2 a'MJLJ

with the new temperatures calculated using the outside-loop K model constants of the composite base component and Eq. (4.101)

The calculation sequence of the Russell method is as follows:

1. Set initial temperatures and total rates for each stage,

2. Initialize the outside-loop variables; the reference Kbj s; the relative volatilities, oty's; and parameters AJs and BJ's of the if-value model and parameters C/s, DJs, EJs, and FJs of the enthalpy models, using the complex A"-value and enthalpy correlations and the estimated set of temperatures and compositions.

3. Initialize the base stripping factors, SbJs, based on the current set of flow rates and reference KhjRef's.

The inside loop consists of steps 4, 5, 6, 7, and 8.

4. Calculate the component liquid rates by the tridiagonal matrix method. Component vapor rates are found by the component (not the base) stripping factors, vtJ - StJ ltJ.

5. Get the total flow rates by summing the component flow rates.

6. Get new stage temperatures by the if^-method.

7. Get a new set of stripping factors by solving the energy balances and specification equations as the independent functions of the quasi-Newton technique of Broyden (Sec. 4.2.6). The derivatives of the Jacobian matrix are generated numerically and must include steps 4, 5, and 6 each time the independent variables are perturbed. The Jacobian is not recalculated after the first trial in the loop but is instead updated by Broyden's equation.

8. Any calculation of the independent functions must include a pass through steps 4, 5, and 6. The norm of the functions with the new stripping factors should be small (see Sec. 4.2.6 for criteria) to leave the loop. If not, return to step 4.

9. Update the outside-loop variables; the reference Kbj Ref's\ the relative volatilities, ay's; and parameters AJs and Bj s of the if-value model and parameters C/s, D/&, E/s, and F/s of the enthalpy models, using the complex if-value and enthalpy correlations and the current set of temperatures and compositions.

10. The test for solution is that the relative volatilities and the parameters of the lf-value and enthalpy models generated in the outer loop using the complex correlations have not changed greatly between outer-loop trials. If this and the other solution criteria (Sec. 4.2.3) have not been met, return to step 3.

Russell's method is insensitive to the quality of initial values. The user does not have to provide additional stage temperatures and flow rates beyond what might be estimated by the method. The scaling factor of Eq. (4.100) in Boston's method is used only in the first trial to force a better material distribution in the column.

Russell's method is one of the methods used in the MULTIFRAC option for multiple columns in ASPENPlus and is the column method in HYSIM of Hyprotech of Calgary, Alberta. Russell's method as written by Richard Russell is in the PD + Plus system available from Deer-haven Technical Software in Deerhaven, Massachusetts. The HYSIM and PD + Plus versions of Russell's method have been found to work well for refinery fractionators with sidestrippers and other similar columns including columns for which a version of Boston's method failed.

Outlook. The inside-out methods are now the methods of choice for mainstream column simulation. While other methods still have their place, the inside-out methods have displaced the BP, SR, and 2N Newton methods and their application should continue to grow. Other products that have insideyout methods include PRO/II of Simulation Sciences of Fullerton, California, and ChemStation of COADE of Houston, Texas.

### 4.2.11 Relaxation methods

One of the earliest successful solution methods was the relaxation method. The first statement of this method was by Rose, Sweeny, and Schrodt (60). Others have been by Ball (63); Jelinek, Hlavacek, and Kubicek (61); Ketchum (64); and Mori et al. (62).

A relaxation method finds a steady-state solution of a column as if it were an operating column changing with time. The column is initialized using some realistic condition such as startup, with the liquid on every stage having the feed composition at its bubble point. The column is carried to the steady-state conditions by successive approximations of the unsteady-state distillation equations. These unsteady-state equations are modifications to the MESH equations to include changes in the MESH variables with respect to time. This process could be said to be mimicking the physical startup of the column, but the objective is not to follow the dynamic operation but to seek the steady-state solution. The earliest relaxation methods used an unsteady-state, time-dependent form of the component balances to get the compositions instead of the tridiagonal matrix method presented earlier. Later methods added calculating changes in temperature with respect to time.

The unsteady-state form of the component balance for a conventional tray written with respect to time t is dxn

(y; + iy,j + i + Lj-i x,,-! ~ Vj y,j ~ LjXij)/Uj = (4.109)

The term Uj is the material holdup on the tray. The vapor holdup is very small and the liquid holdup here is assumed not to change with time. The composition for the next trial, k + 1, is calculated using Euler's equation:

Euler's equation will suffice as long as the time step, At, is kept small. Jelinek et al. (61) and Mori et al. (62) developed techniques for calculating the time step and the composition derivative and restated the component balance equations accordingly.

In relaxation equations, the total flow rates and if-values do not change from one time step to the next. Once the compositions are found for a step, the if-values can be updated and the remaining MESH equations can be solved by a conventional column solution method. Since the relaxation equations yield a calculated set of product compositions, Ball (63) uses the theta method (Sec. 4.2.5) to correct these when initiating a new trial. As in the theta method, temperatures are updated by the Kb method and total flow rates by energy balance. Jelinek et al. (61) use the theta method for distillation columns and a sum-rates method (Sec. 4.2.7) for absorbers.

The success of relaxation is dependent on a choice of the step size, It. The method can fail if too large a step is chosen or be very slow if too small a step is chosen. Since computers are getting faster, a small step size, unless excessive, should not be a great problem. Mori et al. i,62) used the inverse of an internal vapor rate, such as the reboiler boilup, Vjv, as a step size in order to take into account both feed rate and reflux rate.

A relaxation method will make a stable, stepwise movement toward the steady-state values. Problems arising from poor initial values are absent. A BP or SR method (Sees. 4.2.5 and 4.2.7) can be converted to a relaxation method by substituting the relaxation composition calculation for the one used in the method. Brierley and Smith (106) found the relaxation method to be reliable for a wide range of columns, including stripper-absorbers, and their method is part of DISTPACK of ICI. Jelinek et al. (61) also applied their method to nonideal systems and extended BP and SR methods outside of their normal range of systems. Since global Newton methods (Sec. 4.2.9) need good initial values, a relaxation method can be used to bring the column near the solution before switching to the global Newton method. Such a technique, using a Naphtali-Sandholm-like method, should work well for solving nonideal systems. This is the method of Ketchum (64) and Brierley and Smith (106).

The Ketchum method. Ketchum (64), in his review of relaxation methods, pointed out three inherent weaknesses:

1. The compositions, .r,-/s, are corrected separate of the temperatures and total flow rates, which are calculated by a more conventional BP or SR method. The temperature calculation can be unstable or sluggish.

2. There is no estimate of the change of temperature with time. This could be especially important for wide-boiling mixtures. All of the rigorous methods are a balance between the composition and temperature calculations but ¿f-values and enthalpies are assumed to be constant with respect to time.

3. All relaxation methods tend to move slowly, especially near the solution.

Any of the global Newton methods can be converted to a relaxation form in Ketchum's method by making both the temperatures and the liquid compositions time dependent and by having the time step increase as the solution is approached. The relaxation technique should be applied to difficult-to-solve systems and the method of Naphtali and Sandholm (42) is best-suited for nonideal mixtures since both the liquid and vapor compositions are included in the independent variables. Drew and Franks (65) presented a Naphtali-Sandholm method for the dynamic simulation of a reactive distillation column but also stated that this method could be used for finding a steady-state solution.

The change in the temperature with time is expressed as dTj

The independent functions are transient versions of the Naphtali-Sandholm functions, including the component balances [Eq. (4.109)] and the energy balance of each stage:

where change in the enthalpy of the holdup, Ur with respect to temperature is the heat capacity of the liquid phase, Cpur The equilibrium equation [Eq. (4.79)] does not have a time-dependent form and remains the same.

The time step, Ai, is used to switch the method from being a relaxation method to a global Newton method. When the time step is small, e.g., At = 0.1, then the changes in the independent variables are small. The method performs like a damped Newton-Raphson method, where the steps are small but in the direction of the solution and without any oscillation. When the value of At is large, i.e., Ai = 1000, the method performs like a Newton-Raphson method. The value of At at each column trial determines the speed and stability of the method. The units of the time step are the same as the flows to and from the column. The calculation sequence of the Ketchum method is as follows:

1. Develop a set of initial total flow rates and temperatures. Also set liquid compositions for trial k = 1 and the time step size, Ai.

2. Calculate all of the independent functions and their norm. Get a new set of all of the independent variables, including the temperatures, Tjs, and liquid composition, xt -& by solving the independent equations during one pass through the Newton-Raphson procedure.

3. Besides the usual solution criteria, there is one additional criterion for the relaxation method. The maximum difference between a composition and temperature on any stage calculated in the last trial and the previous should be very small, i.e., and (4.113)

The equations will still form the same block-banded sparse matrix as in the Naphtali-Sandholm method. No matter what size the time step, the same matrix solution technique can be used to calculate the next set of independent variables.

Ketchum's relaxation method is stable and fast because the time £tep, At, acts as switch to a global Newton method. The time step can be increased as changes in the temperatures and liquid compositions become slight, i.e., the above solution criterion is within the same magnitude as the convergence limit. The global Newton methods may be the only choice for solving nonideal systems, but they suffer from needing very good starting values. A Ketchum relaxation technique applied to a global Newton method can bring the composition and temperature profiles to a point where a solution can be found. For very nonlinear, nonideal if-value and enthalpy curves, the relaxation form will step slowly along the curves of the MESH equations and not make large changes in independent variables that could place the calculation well off. The method has proved highly reliable for difficult problems such as reactive distillation [Komatsu (67)] and extractive distillation [Ishikawa and Hirata (68)]. Many dynamic simulation methods can be adapted as relaxation methods. Agreda et al. (112) used the dynamic methods in the DYFLOW system from DuPont to find the steady-state solution of a reactive distillation column. A complete study of dynamic simulation methods was presented by Holland and Liapis (10) ancj the DYFLOW system is described by Franks (66).

### 4.2.12 Homotopy-continuation methods

Homotopy or continuation methods have recently been applied to difficult-to-solve columns. The beauty of these methods is that they are a simple means of forcing a solution. Tom Wayburn's analogy is

"They are like a blind man using a rope to cross a room." The MESH equations can be difficult to solve, either due to the nature of the column (many feeds or side draws, side strippers, near-minimum reflux, etc.) or due to the nonidealities of the if-values or enthalpies, if-values and enthalpy methods for nonideal mixtures, dependent on both temperature and composition, are far from linear, can have discontinuities, and may even change slope. With global Newton methods like the Naphtali-Sandholm (Sec. 4.2.9), there are problems in estimating initial values for the column from which the desired solution can be reached. For three-phase systems, azeotropic systems, or systems of columns with two or more feed-recycle stream combinations, there may be more than one calculated solution. The method must be forced to reach the desired solution. Homotopy methods begin with a known solution of the column and from there follow a path (or integrate) to the desired solution. The known solution can be at different conditions or with much simpler if-value and enthalpy methods and stepped changes are made from there, solving the column at each step, until the final solution is reached.

The homotopy methods can be divided into two general classes, mathematical homotopies and physical or parametric homotopies. The mathematical homotopies are conventions without a physical relationship to the MESH equations and this occasionally causes problems. The physical homotopies have a basis in the MESH equations and these will be emphasized. Taylor, Wayburn, and Vickery (80) state that the physical homotopies should outperform the mathematical homotopies and are easier to implement.

The homotopy principle. A homotopy function is a continuous blending of two other functions so that where F(x) is the difficult solution of a column and G(x) is a simpler problem or the column at known conditions. The homotopy parameter t travels between values of 0 and 1 and is the term shown in Fig. 4.2

SOs)

Figur« 4.2 The homotopy path.