To keep things simple, I’m going to start the discussion with solution of 1-D single-phase equations. The equations that follow all assume that velocity is always positive, and use an upwind (donor-cell) method to get values of thermodynamic variables at volume edges. Recall that mass and energy equations are evaluated at volume centers (integer subscripts) and the momentum equations are evaluated at volume edges (half integer subscripts). Using the notation introduced in previous classes, the semi-implicit mass equation at each cell center is:

_{}.

The energy equation at each cell center is:

_{}

The momentum equation at each cell edge is:

_{}

Remember that everything with a superscript of “n” is known and can be treated as a constant in the equations. Everything with a superscript of “n+1” is unknown. I’ll get back to the form of the density used in the momentum equation later. A donor cell average here is not always a good idea. The mess at the right hand side of the momentum equation is just a linearization of the fully implicit form of V|V|. This is required to avoid problems with stability.

The first step in solving these equations is to select independent variables. Your first reaction is probably to use density, specific internal energy, and velocity. In this case the new time pressures would be determined from an equation of state. This is a good choice for very compressible flows, but a bad one when modeling solid liquid. The equations of state used in all reactor safety codes are fairly realistic and include the compressibility of water. However, the density of liquid water is relatively insensitive to pressure. That means a small variation in density can correspond to a large variation in pressure. When solving the flow equations iteratively, you may conclude that the solution has converged, because density changes are very tiny, but be wrong, because pressure changes from one iteration to the next are still significant.

As a result of the above considerations, everybody in this business uses pressure as an independent variable.

The use of specific internal energy as an independent variable is less of a numerical issue. You will see that choice in RELAP5, and you will see the use of specific enthalpy in some other codes. However, TRACE uses temperatures rather than specific energies as independent variables. That choice is not obvious from the single phase equations. However, if you look at the two-fluid equations, there is frequent use of temperature in the wall and interfacial heat transfer terms. Thermal equilibrium involves phases, or a phase and a wall moving towards the same temperature, not the same specific internal energy. I have found that the algebra involved in setting up an iterative solution of two-phase flow equations is simpler when temperatures are chosen as independent variables.

For the Semi-Implicit equations you will also find one more trick is used in choice of variables. The momentum equation at each edge is solved at the beginning of each time step to give new time velocity as a linear function of new time pressures in the adjacent volumes. Velocity becomes a dependent variable via this relationship just as density and specific internal energy are dependent variables via the equation of state. This leaves only temperatures and pressures as independent variables for the above set of equations. For two-phase equations, new time void fraction also remains as an independent variable.

For a single phase model split into N computational volumes, you end up solving N mass equations and N energy equations (2N total equations) for N pressures and N temperatures. For a two-phase system, you end up solving 2N mass equations and 2N energy equations for N liquid temperatures, N vapor temperatures, N pressures, and N void fractions (4N equations and 4N unknowns).

In the Reactor Safety community, we’ve always used a Newton Method to solve these equations. The CFD folks have traditionally used a family of methods called SIMPLE, but seem to be coming around to our way of doing things. The only major variation within reactor safety codes is that RELAP5 only uses the first linearization of a Newton Method, while TRACE and CATHARE iterate to some measure of convergence. Last I checked, RETRAN let you pick either approach. Because a single iteration has no hope of enforcing a conservative form of the difference equations, RELAP5 has special time step controls based on monitoring mass conservation.

If you’ve completely forgotten how to do a

By now you are saying, all I want to do is run the programs, what does all this mean to me? First, it’s part of understanding your experimental apparatus. Codes vary over time, and it’s hard to predict what you will be using in the future. If you end up with some off-brand code that uses density as an independent variable, be sure to keep an eye on the convergence, and set the convergence limits about three orders of magnitude tighter than you would with a pressure based code.

There are two important characteristics of

A second important characteristic of

The primary method of setting the convergence level in TRACE is via
the maximum absolute fractional change in pressure for an iteration. This
appears as EPSO in the input. I recommend
that
you set this to 1.e-5 for most calculations. TRACE
also has some relatively loose internal checks on changes in
temperature and void fraction, and a currently undocumented method to
control convergence through tests on the residuals of the mass and
energy equations. The tests on pressure
changes
normally work quite well. However, there
are
some very rare occasions (I’ve seen 3 or 4 in 25 years) when a
reasonable test on pressure still results in a failure to properly
solve the energy equation, that in turn places the solution to the code
in an unrecoverable state. Fixes for this
involve rerunning that portion of the transient with either a tighter
convergence criterion, or smaller time step size. If
you are running TRACE as a real-time plant simulator, run time is
critical, and iteration count is an important contribution to CPU time. For the range of transients to be simulated,
you
should do sensitivity studies on EPSO. For
relatively benign transients, you will find that substantial relaxation
is possible for EPSO.

When you move to other simulation codes, be sure to read the fine
print and understand the iteration procedure. You may not be
getting
quadratic convergence, and your stategy for setting the convergence
criteria will probably need a change (tighter than you've used with
TRACE).