Solution of the Finite Volume Equations

Assignment :  Finish HW6

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 Newton iteration, take a look at this reference.  If you are fuzzy on how to apply a Newton method to more than one equation and one unknown, take a look at this reference.  Once you’ve got that down, take a look at this description of the solution of the semi-implicit equations.

So What!

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 Newton’s Method that you need to consider when running a problem. The first is the “radius of convergence”.  If the equations that you are solving have an unpleasant non-linearity (not too unusual for us), and your initial guess at the answer is too far from the actual solution, the iteration will fail to converge. Fortunately, your initial guess is the old time (level “n”) values of the variables, and you can make that as close to the final solution as you want by lowering the time step size.  A code will do that for you automatically.  If the iteration count on one time step is too high, the size of the next step is reduced.  If an iteration fails, the time step is restarted using a smaller (usually half) value of the step size.  If a time period of a calculation has many time step backups due to iteration failures (look at error messages in trcmsg for TRACE, and plot time step size), you are getting a strong hint that you set the maximum time step size too high for that time period. Look at your plot of time step size vs. time, and select a maximum time step size somewhere below the size at which backups occur, and the minimum value after the failure.  I usually start by setting my maximum time step halfway between the largest and smallest time steps that I see in a region with frequent iteration failures.  This does not negate the need for further time step sensitivity studies.  It just gives you a better place to start.

A second important characteristic of Newton iterations is quadratic convergence.  Once the iteration is getting close to the answer, the difference between the answer from the iteration and the actual solution drops very quickly.  For example if a measure of maximum fractional error is 0.01 on one iteration, the error on the next will be about 0.0001, and on the following it will drop to about 1.0e-8. If I’m simply looking at the fractional change in pressure over an iteration (δP/P), it will follow a similar trend.  This behavior tells you two important things.  Whatever measure of convergence you are using , is based on conditions at the current iteration.  After you’ve applied the correction on the last iteration, things will be much closer to the correct answer.  For example if I stop an iteration when the maximum absolute fractional change in pressure is less than 1.e-5, then all I’ve missed from the next iteration is a further fractional correction less than 1.e-10.  The second thing that you should conclude is that if you’ve chosen a reasonable convergence criteria, and you hit 8 to 10 iterations without convergence, something is seriously wrong (probably outside the convergence radius).  It’s best to let the code backup and drop the time step.  I tell TRACE users to never set the maximum allowed number of iterations (OITMAX) greater than 10.

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).

 

Back to List of Lectures / Home

Maintained by John Mahaffy : jhm@psu.edu