Newton-Raphson algorithm, equilibrium calculations

Numerical Methods and Data Structure. Both EQ3NR and EQ6 make extensive use of a combined method, using a "continued fraction" based "optimizer" algorithm, followed by the Newton-Raphson method, to make equilibrium calculations. The method uses a set of master or "basis" species to reduce the number of iteration variables. Mass action equations for the non-basis species are substituted into mass balance equations, each of which corresponds to a basis species. 

C.M. Bethke (52) has shown that significant numerical advantages in such calculations can be realized by switching into the basis set a mineral species that is in partial equilibrium with the aqueous phase. This avoids expansion of the size of the Jacobian matrix and reduces computation time. A method based on this concept is being developed for use in the 3270 version of EQ3/6. The concept appears to show promise for improvement of the "optimizer" algorithm as well as the Newton-Raphson one. 

For the bubble-T calculation in the phi-phi form, a viable alternative to Newton-Raphson is presented in Figure 11.1. This algorithm is composed of three principal parts an initialization, an outer loop that searches for the unknown T, and an inner loop that searches for the vapor-phase mole fractions y. The algorithm can be used for any number of components, but it is restricted to equilibrium between two phases. In the special case of a single component, the algorithm is equivalent to the Maxwell equal-area construction given in (8.2.22). 

The central portion of the algorithm in Figure 11.6 exactly parallels the standard Rachford-Rice procedure. First, we use (11.1.27)-(11.1.29) to compute the mole fractions for all phases, then we compute all fugacity coefficients and all activity coefficients. With those quantities we can obtain new estimates for the Cs and Ks from the phase-equilibrium relations (11.1.15) and (11.1.24). Now we use (11.1.31) and (11.1.32) to calculate values for the Rachford-Rice functions, Fj and F2, and test for convergence. If our convergence criteria are not met at iteration k, then we use the Newton-Raphson method to estimate the unknown L and V at the next iteration (fc + 1). 

Applying the Newton-Raphson method in a nonlinear finite element system will yield results only in the pre-collapse range, but it will fail to give information about the post-collapse response. To circumvent this limitation, a constraint can be added into the finite element system, which relates the load increment and the incremental displacements within each iteration (Fig. 12). This technique allows the calculation of the whole equilibrium path, even beyond the critical limit points. A number of different solution algorithms have been proposed in the literature (Riks 1979 Crisfield 1981 Ramm 1981 Bathe and Dvorkin 1983). 

---