This section describes some of the numerical methods used by QTRAN which may be of interest to the user. These methods include QTRAN’s predictor-corrector algorithm which is used for transient calculations, the Strongly Nonlinear Point Successive Over Relaxation (SNPSOR) point iterative solution algorithm which is used to solve the systems of nonlinear equations which are posed by implicit transient integration or by steady-state calculations, convergence acceleration schemes used by the SNPSOR algorithm, table interpolation schemes, the phase change algorithm used by QTRAN, and the finite element to resistor/capacitor translation algorithm.

QTRAN uses a predictor-corrector algorithm for transient problems. This algorithm uses a first order predictor with a second order corrector equation. Details are presented in the following sections.

QTRAN currently uses simple linear extrapolation in the time domain to predict the future value of a node’s temperature.

A family of one-step algorithms has been developed by Hughes (see Ref. 1. in Appendix A) and are claimed to be unconditionally stable for nonlinear type of problems.

QTRAN’s Modified Hughes Equation:

If the heat absorbed into the capacitor is treated as merely another heat flow, the left hand side of this equation can be moved to the right side and presented as:

where the value of T[t+β*dt] is taken as:

Furthermore, since t and T[i,(t)] may be assumed known and hence constant, using the corrector equation to integrate with respect to time may be seen to be equivalent to finding a value of T[i,(t+dt)] that zeros the following function:

The corrector algorithm of the previous section uses a weighting factor β which is used to determine the amount of “implicitness” that will be used for the formula. As per Ref. 2. in Appendix A, it is possible to choose this factor β such that fourth order or higher spatial accuracy can be achieved. For a 1-D problem, the method of choosing β is detailed as follows:

QTRAN currently offers this basic scheme as an optional adaptive weighting algorithm. A significant feature of QTRAN’s adaptive explicit/implicit weighting scheme is a lower cut-off value for the Fourier Modulus. If the Fourier Modulus is such that the current time step being used by QTRAN is less than 0.35 of the stable time step for a given node, that node is then integrated with β = 0.0 (totally explicit). This can speed up calculations very significantly.

QTRAN uses a Strongly Nonlinear Point Successive Over Relaxation (SNPSOR) equation solving algorithm which is a diagonalized hybrid of Newton’s First Order Method and Newton’s Second Order Method when solving the systems of equations generated by a thermal network. Many nonlinear computer programs use a linear point iterative technique such as Gauss-Seidel coupled with successive over relaxation (SOR). Such an algorithm will frequently prove to be satisfactory. However, this algorithm can run into difficulties if the equation system is strongly nonlinear as is the case with “hot” thermal radiation problems where temperatures of 2000+ C may be encountered. When linear SOR algorithms encounter these problems, under relaxation is usually used to force convergence. However, this under relaxation will usually cause convergence to be unacceptably slow and will frequently adversely affect accuracy due to premature convergence.

An approach to the problem is to use a hybrid of Newton’s First Order Method and Newton’s Second Order Method. Newton’s First Order Method is essentially the same as using the first two terms of a Taylor series to estimate a function's zero, while Newton’s Second Order Method uses the first three terms to estimate the zero. Unfortunately, Newton’s Second Order Method requires three function evaluations per iteration if it is to be used. Furthermore, the Second Order Method is generally prohibitively costly and wasteful on conduction problems since conductive heat transfer is generally very close to being linear (First Order). It should also be noted that conduction dominated nodes are generally the most numerous in large thermal system networks. The approach used by QTRAN is a hybrid of the First Order and Second Order methods. Specifically, second derivatives are not calculated for conduction terms (i.e., conductive resistors are assumed to be linear within a given iteration and their contribution to the node’s First derivative is simply the sum of the resistor’s reciprocal values (conductances)). All other heat contributions to a node (e.g., capacitances, convective resistors, radiative resistors, mass flow resistors, and heat sources are assumed to be nonlinear and hence Second derivatives are calculated for these terms). Conductive resistor heat flows (which are generally the most numerous) are therefore calculated only once per iteration while other heat flows are calculated three times per iteration. The result is a rapid yet solidly convergent algorithm capable of both rapid computation for conduction problems and reliable convergence for the strongly nonlinear problems. For more information concerning Newton’s Second Order Method, see Ref. 3. in Appendix A.

While the hybrid Newton’s method works great for most problems, there are a few classes of problems that present pathological cases for failing any Newton’s method. These are the classic “S” shaped curves, and phase transition problems introduce exactly this type of “S” shape into the nodal heat flow rate curve. Fortunately, QTRAN can detect a phase transition and will shift to a bisection method for any node currently undergoing a phase transition. Bisection methods are relatively slow, but they are immune to the “S” shaped curves.

QTRAN currently offers the user the option of selecting either a fixed or an adaptive over-relaxation scheme. In addition, a single over-relaxation parameter may be applied to the entire system of equations, or groups of relaxation parameters may be applied to groups of nodes dependent on the type of node and boundary conditions present, or individual relaxation parameters can be applied to each individual node. The following discussion applies to which-ever method of application is selected. Over-relaxation schemes can speed the solution of elliptic problems by a factor of over 30 for many problems which are commonly encountered. Numerical experiments by the author with QTRAN’s adaptive relaxation algorithm (i.e., QTRAN estimates and continuously updates the relaxation parameter for the adaptive algorithm) have shown that the adaptive algorithm is over twice as fast as the optimal fixed over- relaxation algorithm for many transient problems. QTRAN’s adaptive relaxation algorithm has also proven to be quite successful in estimating the optimal relaxation parameter for steady- state problems. Special note should be made that over-relaxation (not under relaxation) can be used even with the “hot” thermal radiation problems. QTRAN’s SNPSOR algorithm does not (in general) force the user to resort to under-relaxation to gain a sufficient radius of convergence (with the possible exception of a few sublinear problems), whereas codes relying on linear SOR solvers frequently do.

QTRAN’s relaxation parameter estimation algorithm is based on the convergence rate R of the iterative solution, where R is calculated from:

E[i-1] is the largest iterative delta (positive or negative) at the i-1'th iteration, and E[i] is the largest iterative delta (positive or negative) at the i'th iteration. In the limit as the iteration count i goes to infinity, the following formula (see eq. 3-83, 3-113, and 3-116 of Ref. 10. in Appendix A) may be used to estimate the relaxation parameter (assuming that the eigenvalues of the iteration matrix based on the current relaxation parameter have not yet become complex).

if

else

endif

RELAXmaximum is the maximum allowed relaxation parameter which can be set by the user. The theoretical maximum value is 2.0 and the code limits it to a value of 1.999. RELAXdamp is a damper on the rate of increase of the relaxation parameter. As with other applications of numerical analysis, this parameter does not have a theoretical basis, but numerical license allows for its interjection. This parameter tends to prevent the over-relaxation parameter from going beyond its optimum value.

As stated previously, this algorithm works fairly well as long as i approaches infinity and the eigenvalues of the iteration matrix have not become complex. However, i never approaches infinity in practice and eigenvalues frequently become complex.

QTRAN approaches the i going to infinity problem by recognizing that what one is really after is to simply get a “good” estimate of the convergence rate R in some averaged sense. Specifically, R will frequently change values drastically in early iterations for a given time step or for a steady-state problem, and also after the relaxation parameter has been updated. What one wants to do is simply wait for the R value to stabilize before updating the estimate of R. This can be easily done by simply requiring that, for a finite number of iterations (say (3) iterations): (1) R be less than 1.0, (2) that the sign of the iterative delta’s are constant, and (3) that the value of RELAXcurrent has not been updated more recently than (4) iterations ago. If these three conditions are met, you can assume that the value of R[i] is approaching the asymptotic limit of R[infinity].

One minor complication to all of this is that the above algorithm for calculating RELAXnew is valid only as long as the eigenvalues of the iteration matrix do not become complex (this was alluded to earlier). For more information, see Ref. 9. in Appendix A. This brings up two major problems: (1) how does one know when the eigenvalues have gone complex, and (2) what to do then, since there is no theory for predicting the optimal relaxation parameter for a system with complex eigenvalues. The answer lies in purely empirical numerical engineering.

First, one observable fact is that the largest iterative delta’s E[i] have a tendency to oscillate in sign (+ and -) once the eigenvalues have gone complex (they can change sign at other times too, but the E[i] will frequently oscillate continuously with complex eigenvalues). Second, the cause of complex eigenvalues is usually that the value of RELAXcurrent is greater than RELAXoptimum for the equation system being solved (although you have no way of knowing how much greater). However, since an oscillating sign is an indication that RELAXcurrent is too large, QTRAN simply subtracts 0.01 from the value of RELAXcurrent each time that the sign of E[i-1]/E[i] changes. Eventually, the value of RELAXcurrent is reduced back to near RELAXoptimum. It is also worth knowing that from a practical consideration, RELAXcurrent needs only to be within about 0.002 of RELAXoptimum in order to be “exact”.

Finally, the value of RELAXcurrent is limited so that 0.01 <= RELAXcurrent <= 1.999.

QTRAN does not base its convergence tests solely upon the iterative delta E[i] as do most other commercial codes. Instead, the test performed is based upon the E[i] value plus the historical convergence rate R. The formula for the convergence test is as follows:

if

else

end if

This basic formula can be derived from the equations presented in Convergence Acceleration Schemes, 218 along with a formula used to compute the sum of an infinite series.

Energy balance is often used as a convergence criteria but it is not used as a criteria in QTRAN even though it is printed out.

The energy balance in the QTRAN output is equal to the energy stored during a time step for transient problems and is zero for steady-state problems. If advection is present, then the energy balance will include the net energy carried into and out of the problem by the flowing fluid. When examining the energy balance for the system, it should be compared to the total energy flow in the system. If the system energy balance was 1100 when a zero was expected, it would be a large number; however, if it were compared to the total system energy which may be 1.1 x 1012, then the 1100 for the system energy balance becomes a small number.

The System Energy Balance (QSB) reported by QOUT.DAT is defined as:

where, QNF is the combined net heat flows on all nodes including the boundary nodes and QS is the heat fluxes by applied Heating or macro functions. In cases where the heat flow is driven by boundary nodes, the system energy balance is simply the sum of all net heat flow rates on every node (including the boundary nodes) as reported in QOUT.DAT. The system energy balance may be higher than the control volume heat balance (heat entering the system minus heat leaving the system) since it also reflects the numerical tolerances on the calculated nodes. As such, it provides a more realistic picture of the uncertainty.

QTRAN currently uses three distinct table interpolation schemes. These schemes are detailed in the following sections.

This is the simplest table algorithm used by QTRAN and involves simple linear interpolation between table data points. Linear extrapolation is used when data is requested that lies outside of the X-Y data provided by the user.

This is a slightly more complex algorithm that resembles cubic spline interpolation in nature. The exact computational process involves the determination of the tabular function’s first derivatives at table interval boundaries. Using the first two derivative values and the function’s values at an interval’s boundaries, a cubic equation can be derived and used for interpolation within the interval. Since neighboring intervals also use first derivatives when calculating their cubic equations, it is seen that these tabular functions maintain continuity of slope (first derivatives) from interval to interval as well as maintaining continuity of value (linear interpolation maintains continuity of value from interval to interval, but not continuity of slope). Hermite Polynomial Interpolation can thus be used to provide a “smooth” tabular function if the user desires. Quadratic extrapolation is performed if the value lies outside of the tabular data.

QTRAN also offers Linear Computed Interval (LCI) interpolation for library material property evaluations. This linear table algorithm is an extraordinarily rapid linear interpolation scheme. A table property value is returned using only two multiplication and two addition operations. No interval search is required when LCI is used. Extrapolation is performed if the value lies outside of the tabular data.

The phase change algorithm used by QTRAN accurately models latent heats by adding an artificial temperature-dependence function into the capacitor’s specific heat curve. This has the effect of taking the latent heat and spreading it out over a very narrow (user specified) temperature region. This has the advantage of maintaining a capacitor’s energy content as a function of temperature and is done for numerical and coding convenience. It is similar in function to putting a “spike” into the specific heat curve to account for latent heat effects. However, unlike some other codes which use such “spikes”, QTRAN knows exactly where these “spikes” are in the temperature/enthalpy curve and there is no possibility that QTRAN can jump over them or miss a part of the spike (and hence erroneously fail to calculate the latent heat effects properly) as some other codes can. This technique has no restrictions whatsoever in terms of whether the user runs explicit, implicit, or some weighted average. Many codes add restrictions of this type for phase change models, but not QTRAN.

There is currently no coded limit to the number of phase regions that may be assigned to any capacitor in the system, nor is there any limit to the number of times that a given capacitor may transition from one phase to another. This also includes phase direction reversal prior to the complete melting, freezing, etc., of a node. There is also no limit to the number of nodes which may be simultaneously undergoing a phase transition.

Material properties may be varied in a near step-function manner at phase transition temperatures. While QTRAN does not perform convective calculations when a material has changed phase, it does assume that conduction continues to be a heat transfer mechanism. It is thus seen that QTRAN is able to continue with heat transfer calculations after a phase transition has occurred, subject to the assumption that an appropriate value for the thermal conductivity of the newphase region is used

The QTRAN thermal analysis package includes a routine in PATQ to translate finite element data to resistor/capacitor data for the following finite elements:

PATQ program reads in finite element data for these element types and then constructs a mathematically exact resistor/capacitor representation. The network is NOT an approximation, but instead includes all finite element cross-derivative terms and any other increased accuracy that is inherent to a finite element discretization scheme using a lumped capacitance matrix. Such a translation program is highly desirable from the standpoint of being able to take advantage of the excellent finite element solids modeling packages that exist, and also from the standpoint of accuracy for skewed meshes (many finite difference schemes introduce an intolerable 0th order error for Laplacian difference operators when the mesh being used is skewed). It is also highly desirable from the standpoint that finite element formulations inherently allow for arbitrarily oriented anisotropic materials (e.g., composites), while many finite difference codes cannot perform analyses on anisotropic materials unless the material axes are by chance lined up with the mesh directions.

While the finite element based approach generates a mathematically exact resistor/capacitor representation, the resulting areas and lengths have no physical significance. The lumped capacitance network solver used in QTRAN is able to solve both finite difference and finite element based resistor/capacitor networks either separately or simultaneously. Finite difference based resistor/capacitors may be appended to the PATQ translated model by including them in the appropriate input files (i.e., CONDUCDAT, CAPDAT, etc.).