## Info

In matrix language, the Newton-Raphson equation is

where f^ is a vector of the independent functions at column trial k:

The Axfe represents a vector of the changes in the independent variables that bring an improvement in the functions, f^, from trial k to trial k + 1:

where for any xk + = xk + Ax*. The objective in the Newton-Raphson is to find a set of independent variables xk = [xix2...xj (4.62)

so that all of the independent functions are very nearly zero. The MESH equations presented earlier were written so that they all equal zero in order to fulfill this requirement of the Newton-Raphson method.

Generating the derivatives is the most difficult part of the method. Writing analytical derivatives of composition-dependent K-values and enthalpies can be very difficult, especially if a method such as UNIQUAC and NRTL is being used. The alternative is numerical derivatives. For any independent variable, these can be found by:

8* 2e where e is a perturbed value for x, e.g., e = O.OOlx Numerical derivatives also have inherent problems. If the perturbed value e is too small, the derivative will be meaningless. If it is too large, the derivative may actually cause the Newton-Raphson to move away from the solution. After all, the Newton-Raphson is an approximation technique. In the derivatives, it assumes that the nonlinear MESH equations are linear over short distances. If too large a distance is used, the convergence will be very coarse.

If analytical derivatives for the if-values and enthalpies are available, they should be used. More effort will be needed when developing analytical derivatives but once tested and error-free, the rigorous method will be more reliable. Lucia and Westman (52) used a hybrid method that included a mixture of analytical and numerical derivatives and found it improved convergence.

In addition to the solution criteria in Sec. 4.2.3, the Newton-Raphson method requires an additional check on convergence. This check, termed the norm or the square root of the sum of the squares, tests that all of the functions are driven nearly to zero:

When the norm is small (i.e., <10~4), the Newton-Raphson method has found a set of variables that reasonably satisfies the independent equations. If the norm is at a small value but the other solution crite ria have not been reached, the calculation should continue with the norm driven increasingly smaller.

Produced from the manipulation of the Jacobian are the changes in the variables, i.e., the Axk vector. The variables for the next trial are calculated from x* + 1 = x* + sk Ax* (i.e., xUk + x = xltk + sk etc.). The sk scalar is generated to ensure that the norm of functions improve between trial k + 1 and trial k. Usually, sk = 1 but may have to be smaller on some trials. The Newton-Raphson method assumes that the curves of the independent functions are close to linear and the slopes will point toward the answers. The MESH equations can be far from linear and the full predicted steps, Ax*, can take the next trial well off the curves. The sk scaler helps give an improved step search or prevents overstepping the solution. Holland (8) and Broyden (119) present formulas for getting sk.

There should also be a restriction in the Newton-Raphson so that the independent variables calculated in a trial are within a tolerance. When the independent variables are temperatures, they should be in a range, Tmin < Tj < Tmax. Also, all flow rates must be greater than zero. If any variable is not within the tolerances, all steps, Ax^, should be halved.

Holland (8) recommends LU (lower-upper) factorization for calculating the Astfc vector. Older methods such as a Gauss-Jordan elimination should be avoided because they build excessive computer error. In LU factorization, the Jacobian is rearranged into two matrices; J* = L* Uk where the L^ matrix or lower triangle has elements on the diagonal and below and the U* matrix or upper triangle has elements on the diagonal and above. LU factorization has similar characteristics of the forward and back substitution of the Thomas algorithm for tri-diagonal matrices. There are many LU factorization routines available. Press et al. (121) and Kahaner et al. (122) have well-documented routines in their texts. Commercial packages such as from IMSL, Inc., in Houston, Texas, the Harwell Subroutine Library from AERE Harwell in Oxford, England, and the LINPACK system from the National Energy Software Center, Argonne National Laboratory, Argonne, Illinois, all have a more than one LU factorization and matrix solution routine.

A basic Newton-Raphson iterative calculation proceeds as follows:

1. Assume a set of independent variables x^, where k - 0 for the first trial.

2. Using x^, calculate the values for the independent functions fk.

Test-if a solution has been reached by calculating the norm.

3. Calculate the partial derivatives for the Jacobian matrix J* [Eq. (4.58)].

4. Using LU factorization, solve the Jacobian for the Ax*, vector. Add this to the existing set of independent variables, x*, to get a new set of independent variables, x* + j. Continue through the rigorous method. The next pass through the Newton-Raphson begins at step 2.

In some Newton-Raphson-based rigorous methods, a poor set of starting values can cause the calculation to never approach a solution. Occasionally, the calculation oscillates, with values swinging to either side of the solution. A good simulation program should include means for detecting and preventing this, e.g., by damping or limiting the change to the next set of variables. A Newton-Raphson method will normally take even steps toward the solution.

The quasi-Newton methods. In the Newton-Raphson method, the Jacobian is filled and then solved to get a new set of independent variables in every trial. The computer time consumed in doing this can be very high and increases dramatically with the number of stages and components. In quasi-Newton methods, recalculation of the Jacobian and its inverse or LU factors is avoided. Instead, these are updated using a formula based on the current values of the independent functions and variables. Broyden's (119) method for updating the Jacobian and its inverse is most commonly used. For LU factorization, Bennett's (120) method can be used to update the LU factors. The Bennett formula is

J*+1 = k+iU* + i = kU* + & - (1 -«*£)] C Ax, (4.65)

where sk is the scalar found above such that norm* + : < norm* and C is a scalar found by

Just as LU factorization is more efficient and more accurate than calculating the inverse, the use of the Bennett method is less time consuming than the Broyden method.

The quasi-Newton methods reduce the computer time spent per trial by reevaluating the Jacobian but at the expense of a greater number of trials. The total computer time sometimes exceeds that of the conventional Newton-Raphson method. The quasi-Newton methods also are more sensitive to how close the initial values are to a final solution and tend to fail more readily. They may have to be avoided altogether for difficult problems or where the number of independent variables exceeds 100.

Quasi-Newton methods like Broyden's and Bennett's do work well when the curves are almost linear across the entire range. To have independent functions nearly linear, they need to be designed so that they are quite simple. The inside-out method (Sec. 4.2.10) is one of the methods that do this.

4.2.7 Sum-rates (SR) method

The SR method is suitable for modeling absorbers and strippers. For some extremely wide boiling systems, especially those with noncon-densables, it is the best method. It has been found to work very well for the side strippers of a refinery fractionator. Absorbers typically have a rich gas bottom stage feed and a lean oil top stage feed. The equations of the SR method do not allow its direct use for reboiled absorbers, absorbers with condensers, or distillation columns. For these columns, other methods like that of Tomich (Sec. 4.2.8) or Russell Sec. 4.2.10) can be used.

The SR method was first presented by Sujata (35) and McNeese (36), and was given its name by Friday and Smith (1). The application of the tridiagonal method for solving the material balances was added by Burningham and Otto (34).

In the SR method, temperatures are the dominant variables and are found by a Newton-Raphson solution of the stage energy balances. Compositions do not have as great an influence in calculating the temperatures as do heat effects or latent heats of vaporization. The component flow rates are found by the tridiagonal matrix method. These are summed to get the total rates, hence the name sum rates.

The single independent equation for each stage is its energy balance. The equation should be written in terms of the ratio of the energy entering the stage to that leaving the stage, i.e.,

Li-ihi-i

The enthalpy of a feed, such as the lean oil feed, La, on the top stage or the rich gas feed, VN + on the bottom stage will appear in the numerator. An interreboiler duty also appears in the numerator of the stage. The enthalpy for a side product or an intercondenser duty appears in the denominator of the stage.

The number of independent functions is equal to the number of stages:

The independent variables of the Newton-Raphson method are the stage temperatures:

The SR method requires that the total molal flow rate, composition, location, and thermal condition of all feeds be specified, including the lean oil, La, to the top stage and the rich gas, VN + 1( to the bottom stage. The steps of the algorithm are as follows:

1. Set initial temperatures and total vapor rates for each stage. Initial liquid rates can be found by total material balances.

2. Based on the most recent set of temperatures and total flow rates, calculate the component vapor rates using the tridiagonal matrix method. Calculate the component liquid rates by applying the absorption factor [Eq. (4.6)].

3. Calculate the total flow rates by summing the component flow rates.

4. Calculate all of the independent functions and their norm. Get a new set of temperatures T- s by solving the independent equations during one pass through the Newton-Raphson procedure.

5. The norm of the independent functions must be very small at the solution (see Sec. 4.2.6 for criteria). This should be tested in step 4 before proceeding through the calculation and solution of the Jacobian. If the norm is not small or other criteria (Sec. 4.2.3) have not been met, return to step 2,

Whenever the stage temperatures are varied, such as in the calculation of the derivatives for the Jacobian, steps 2 and 3 must be included in the calculation of the independent functions. The temperatures affect the Jf-values and therefore compositions, and these changes cascade through the column. Because of this, there are no zero elements in the Jacobian and a tridiagonal matrix technique cannot be used in its solution.

The test of the summation of the vapor component rates to see if a solution has been reached should be based on the total vapor flow rate of the previous trial and the component rates of the current trieil. If the solution is based on the independent functions alone, the SR

method may not calculate through enough trials and the solution may not material balance.

The SR method can be applied to distillation columns, but the equations of the algorithm do not allow the solution of the condenser and the reboiler with the other stages in the column. Because only the energy balances are used as independent functions, reboiler and condenser duties, reflux ratio, and the boilup ratio have to be specified. This overspecifies the column and the solution cannot be found. The condenser and the reboiler can be solved as separate unit operations in a flowsheet as demonstrated by Fonyo et al. (39). The SR method is used in the ABSBR step of FLOWTRAN of Monsanto, St. Louis, Missouri, and also in both the public release version of ASPEN and in ASPENPlus of AspenTech, Cambridge, Massachusetts.

### 4.2.8 2N Newton methods

In both the BP and SR methods (Sees. 4.2.5 and 4.2.7) the temperatures and the total flow rates are calculated separate of each other. An alternative is to calculate the temperatures and total flow rates together in a Newton-Raphson technique. The name 2N Newton comes from that there are two equations per stage for a total of 2 x N functions and variables per column for the Newton-Raphson. In all three of these approaches, the component flow rates are still calculated in an intermediate step.

The best-known presentations are by Tomich (32), Holland (8), and Orbach et al. (33). These vary in their choice of Newton-Raphson equations and independent variables and each may solve a different range of columns. These methods have been shown to work well for wide-boiling mixtures including refinery fractionators, absorber-stripper columns, and reboiled absorbers.

The Tomich method. This method has the summation equations for the vapor and liquid compositions are merged for the first independent equation:

The xy's and y^'s are found from the component balances and the equilibrium equation.

The other independent function is the stage energy balance and is the same as Eq. (4.67) of the SR method above.

In the original statement of the Tomich method, condenser and reboiler duties must be specified, but the equations can be readily adjusted to make these variable. The independent functions are then

The independent variables are the stage temperatures and total vapor flow rates xk = [TlT2...Tn_, Tn V2... Vw_, VW]T (4.73)

The calculation sequence of the Tomich method is as follows:

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

2. Find total liquid rates by stage-to-stage material balance.

3. Based on the most recent set of temperatures and total flow rates, calculate the liquid compositions using the tridiagonal matrix method and find the vapor compositions by the equilibrium equation.

4. Get a new set of temperatures and vapor rates by solving the independent equations during one pass through the Newton-Raphson procedure.

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

Though the summation equations and the energy balances are the independent functions for the Newton-Raphson, the component and total material balances must be recalculated whenever the T- s and V)'s vary. This includes when the functions are evaluated during the filling of the Jacobian and so there is a full Jacobian, i.e., no zero elements. While the derivatives well off the diagonal of the Jacobian are small, they cannot be ignored. Ignoring these is likely to cause the method to fail or to converge to the wrong answer.

Tomich recommends filling and inverting the Jacobian only once and using the Broyden method (Sec. 4.2.6) to update the inverse. This reduces computer time per trial but increases the number of column trials or passes through the procedure and for some columns may also decrease reliability of the method. It is the author's experience that for many columns, solution is easier to reach when the Jacobian is filled and inverted in each trial and Broyden's method is not used.

The 2Af Newton-Raphson method. This method of Holland (8) is similar to the Tomich method but adds two innovations. The summation equa tion for each stage is replaced by either a dew-point or bubble-point equation as an independent function and the total vapor rate is replaced by a multiplier 9, on the ratio of the total vapor to liquid rates LjfVj. Since the dew- and bubble-point equations are combinations of both the summation equation and equilibrium equation, both now directly come into play in the independent functions. The other independent functions are still the stage energy balances. The independent functions are represented by

The independent variables are the temperatures and the 9, multipliers:

The multiplier 0, is defined as follows:

where (Lj/Vj)ca uses the most recent value of the total flow rates. The corrected Lj/Vj ratio is used in the absorption and stripping factors of the component balances and in the total material balances.

The fy multiplier, by being a correction to the ratio of total flow rates, adds a measure of stability to the solution by preventing the existence of negative flow rates. The steps of the 2N Newton-Raphson method follow those of Tomich.

Application. Both the Tomich and the 2N Newton-Raphson methods are proven methods and have been applied often. The Tomich method was part of the GMB system of The Badger Company, Cambridge, Massachusetts, and is in many in-house simulators. Both methods are best for wide- or middle-boiling separations. Because one of the equations in the 2N Newton-Raphson method is a dew- or bubble-point equation, it may work better for middle or more narrow-boiling mixtures than the Tomich method. Both methods have also been applied to absorber-strippers, though an SR method is still the best method for the most wide-boiling absorber-strippers. Because of the full Jacobian more numbers to manipulate), for columns over 50 stages these methods will use excessive computer time and memory. Also, the solution of the Jacobian is prone to failure when the number of stages is high, and so these methods are not recommended for tall columns.

### 4.2.9 Global Newton methods

One group of methods that is very popular are the global Newton methods, also called the simultaneous correction fSC) methods. The most applied of these are that of Naphtali and Sandholm (42) and Goldstein and Stanfield (45). Other global Newton methods have been presented by Ishii and Otto (47) and Gallun and Holland (40). The general principles of a global Newton method are common to all these methods. The Naphtali and Sandholm method will be discussed first to gain an understanding of the global Newton methods. Advantages and differences of the others listed above will follow. Other methods of merit but not presented here are that of Ferraris (41), Ricker and Grens (57), and Kubicek, Hlavacek, and Prochaska (46). Global Newton methods have been extended to include additional equations and variables for solving three-phase distillation columns as done by Block and Hegner (50), Wu and Bishnoi (51) and Bondy (59), and for reactive distillation columns by Holland (8) and Rafal et al. (58).

In the global Newton methods, all of the equations are solved together in a Newton-Raphson technique. The methods vary in their choice of variables and MESH equations for the Newton-Raphson calculation, but none of the MESH equations is solved in any separate step outside of the Newton-Raphson calculation (e.g., the component balances in the 2N Newton methods). In the BP, SR, and 2N Newton methods, the compositions lag the other MESH calculations (since K-values and enthalpies are generated using the compositions from the previous trial) and compositions of each component are calculated independently of the others (due to the use of the tridiagonal matrix). These are major disadvantages with highly nonideal systems, where AT-values (especially activity coefficients) and enthalpies are highly composition dependent and where the composition of one component cannot be readily decoupled from those of others. The global Newton method includes the component balances among the Newton-Raphson independent functions and compositions join other MESH variables as the independent variables.

The global Newton methods are the most sensitive of the rigorous methods to the quality of the initial values and often require initial values near the answer. While they are the most powerful in solving nonideal mixtures, it may be necessary to use another rigorous method, such as a BP or SR method, to develop initial values approaching the solution before the global Newton method takes over the calculation.

The Naphtali-Sandholm (42) method. This method chooses the stage temperatures and component vapor and liquid rates from the MESH variables as the independent variables of the Newton-Raphson calcu lation. The MESH equations serving as independent functions and grouped by stage are the energy balance, a component balance for each component and an equilibrium equation for each component. There are 2C + 1 equations per stage with N(2C + 1) equations for the whole column. The energy balance is the normalized version presented for the SR or Tomich methods [Eq. (4.67)] except that the total flow rates are not independent variables. Instead, they are calculated by summing the component flow rates to give c c fy+i 2 vo+i + hj-1 2 h-i —^c-^--1 = 0 W.77)