Modelling of two-phase natural circulation instabilities

From the modelling point of view, two types of instabilities are important. These are the static and dynamic instability. Static instability can be fully described by the steady state governing equations. Examples of this type of instability are Ledinegg (flow excursion), flow pattern transition instability, flashing instability, etc. For the dynamic instability, the feedback effects are important and hence it requires the solution of the time dependent conservation equations of mass, momentum and energy. Typical example is the density wave oscillations observed in two-phase systems. Often, thermalhydraulic systems exhibit compound instability that is usually a combination of two mechanisms. Here, the instability may begin with one of the known mechanism for static instability that may trigger some secondary mechanisms resulting in an oscillatory behaviour. Typical examples are parallel channel instability and coupled neutronic-thermalhydraulic instability of BWRs.

In generally, a two-phase natural circulation system can be designed to avoid all the different types of instability discussed above. This is done by carrying out a stability analysis, which is usually phenomena-specific and therefore cannot be described in detail in a brief note as the present one. However, the objectives of most stability analyses can be listed as:

(a) To generate stability maps which delineate the stable and unstable zones of operation;

(b) To specify adequate stability margin for design purposes; and

(c) To predict the nature of the unstable oscillatory behaviour of flow, power and temperature (in time domain) if the system passes through an unstable operating zone during power/pressure raising.

Most static instability analyses are carried out with the prime objective of generating a stability map. For the dynamic analysis, however, all the three objectives become important. Both the linear and nonlinear methods are used for the analysis of this type of instability. Mathematically, the dynamic behaviour of any system can be considered to be linear for small perturbations around the steady state operating condition. An analytical solution is sought for the perturbed equations to obtain the characteristic equation for the stability of the system. The characteristic equation is solved numerically to generate the stability map (locus of all neutrally stable points) of the system. The computer codes based on the linear stability method (e. g. RAMONA, NUFREQ, etc.) are very useful to predict the stability maps without consuming much computer time. However, it is beyond its capability to predict the non-linear effects that come into play after the neutral threshold (to achieve the third objective). Hence, direct numerical solution of the non-linear governing equations is carried out using finite difference methods. An additional objective of the non-linear computer codes is to provide an understanding of the basic physical mechanisms involved in BWR behaviour beyond neutral stability. In principle, system codes like RELAP5, ATHLET, CATHARE, etc. can be used for this purpose, although this is not the current practice (generally). If the system codes are used, this will require the use of qualified fine nodalisation as those used for the normal transient analysis (normally a coarse nodalisation) may not even predict the existence of instability (as the numerical solution technique adopted in these codes are very robust which introduce numerical damping through artificial/numerical diffusion).

For simulating coupled neutronic-thermalhydraulic instability, usually a point kinetics model is used in addition to a fuel heat transfer model (which basically solves the conduction equation in the fuel, fuel-clad gap and clad to obtain the feedback effects of fuel). Point kinetics equations are first order ordinary differential equations that can be easily integrated by using a conventional numerical technique. However, the problem is the introduction of feedback constants in them, which are spatially averaged over the whole core. So uncertainty is therefore introduced at this level in addition to the inaccuracy in neglecting the changes in core power distribution. Similarly, the way in which the neutronics and the thermal hydraulics interact during the calculation has an important effect on the capability to simulate multi­dimensional core dynamics. Grouping of core channels in thermal hydraulic regions having similar geometric and operating parameters can have strong influence on the out-of-phase mode oscillations.