16.2.3. Modified Damped Newton’s Method

The method by which Twopnt arrives at a solution to the 1-D governing equations is described here in detail. After discretization on a given mesh, we have a system of nonlinear algebraic equations that Twopnt attempts to solve by a damped Newton’s method. Newton’s method determines a sequence of iterations or approximate solutions that approach the true solution. For the sake of notational ease we call these approximate solution vectors . When any arbitrary is substituted into the finite difference analog of Equation 13–2 and Equation 13–3 , they do not equal zero as they would if the true solution were substituted; in general, they equal a residual vector that we will call . The objective is to find a vector that satisfies

(16–12)

In our case the vector is composed as follows:

(16–13)

The corresponding vector is composed of the residuals of the energy equation, the species equation, and an equation specifying that . The final equation is included for purposes of maintaining a banded Jacobian structure. The ordering of the vector corresponds to the order of the vector; it begins with the residuals of the left boundary condition , followed by the residuals at the interior mesh points, and finally the residuals at the right boundary, .

Provided the initial estimate of the solution is sufficiently good, Newton’s method produces a sequence that converges to the solution of the nonlinear equations . The purest form of the algorithm

(16–14)

is too expensive and delicate to be used in practice. On the one hand, evaluation of the Jacobian matrices is time consuming, and on the other hand, convergence usually requires a very good initial estimate . Ansys Chemkin employs the usual remedies. First, the Jacobian matrix is replaced by one, , inherited from a previous step of the algorithm. Second, the full step from to may be cut short by a damping parameter . In this way the iteration becomes

(16–15)

where, , and

(16–16)

The inverse Jacobian matrix in Equation 16–15 is not actually computed; instead a system of linear equations is solved for the undamped correction vector .

Our selection of the matrix and of the parameter is governed by a look-ahead procedure that forces certain behavior upon the algorithm. Thus, having tentatively chosen the program looks ahead to the calculation of . The criterion[121] for accepting is that the undamped steps decrease in magnitude,

(16–17)

This prevents the iteration from stepping away from the region where there is good reason to believe a solution lies. Should fail this criterion, the application rejects it and retries the step with a halved damping parameter or a new Jacobian matrix. The damping parameter is initially chosen to be as large as possible so that does not violate various bounds that are set on the solution variables. We know, for example, that the temperature and mass flow rate must be positive, and that the species mass fractions must be between zero and one. The concentrations of many species, such as fuels downwind of a flame, are close to zero and frequently threaten to place the solution out of bounds.

The Newton iteration continues until the maximum norm of the undamped correction vector is reduced to within a user-specified tolerance. Specifically, the criteria for Newton iteration convergence is when the solution correction vector satisfies

(16–18)

The relative and absolute tolerances are parameters that govern the convergence criteria for the Newton iteration. Roughly speaking the relative tolerance indicates how many significant digits the converged solution should contain, and the absolute tolerance serves to exclude solution components smaller than from the relative convergence criteria. Typically, the absolute tolerance should be smaller than the maximum mass fraction of any species of interest. The relative tolerance should be in the range of 10-3 to 10-4, which will serve to provide 3 to 4 digits of accuracy.

If damping cannot produce a suitable correction, then a new Jacobian is computed. If, after computing a new Jacobian, a damped Newton step still cannot produce a suitable correction, then the application begins to take pseudo time steps. This strategy is described further in Pseudo Time-Stepping Procedure .