Open access

The Coolant Channel Module CCM — A Basic Element for the Construction of Thermal-Hydraulic Models and Codes

Written By

Alois Hoeld

Submitted: March 6th, 2012 Published: February 13th, 2013

DOI: 10.5772/53372

Chapter metrics overview

1,792 Chapter Downloads

View Full Metrics

1. Introduction

The development of LWR Nuclear Power Plants (NPP) and the question after their safety behaviour have enhanced the need for adequate efficient theoretical descriptions of these plants. Thus thermal-hydraulic models and, based on them, effective computer codes played already very early an important role within the field of NPP safety research. Their objective is to describe both the steady state and transient behaviour of characteristic key parameters of a single- or two-phase fluid flowing along corresponding loops of such a plant and thus also along any type of heated or non-heated coolant channels being a part of these loops in an adequate way.

Due to the presence of discontinuities in the first principle of mass conservation of a two-phase flow model, caused at the transition from single- to two-phase flow and vice versa, it turned out that the direct solution of the basic conservation equations for mixture fluid along such a coolant channel gets very complicated. Obviously many discussions have and will continue to take place among experts as to which type of theoretical approach should be chosen for the correct description of thermal-hydraulic two-phase problems when looking at the wide range of applications. What is thus the most appropriate way to deal with such a special thermal-hydraulic problem?

With the introduction of a ‘Separate-Phase Model Concept’ already very early an efficient way has been found how to circumvent these upcoming difficulties. Thereby a solution method has been proposed with the intention to separate the two-phases of such a mixture-flow in parts of the basic equations or even completely from each other. This yields a system of 4-, 5- or sometimes even 6-equations by splitting each of the conservation equations into two so-called ‘field equations’. Hence, compared to the four independent parameters characterising the mixture fluid, the separate-phase systems demand a much higher number of additional variables and special assumptions. This has the additional consequence that a number of speculative relations had to be incorporated into the theoretical description of such a module and an enormous amount of CPU-time has to be expended for the solution of the resulting sets of differential and analytical equations in a computer code. It is also clear that, based on such assumptions, the interfacial relations both between the (heated or cooled) wall but also between each of the two phases are completely rearranged. This raises the difficult question of how to describe in a realistic way the direct heat input into and between the phases and the movement resp. the friction of the phases between them. In such an approach this problem is solved by introducing corresponding exchange (=closure) terms between the equations based on special transfer (= closure) laws. Since they can, however, not be based on fundamental laws or at least on experimental measurements this approach requires a significant effort to find a correct formulation of the exchange terms between the phases. It must therefore be recognised that the quality of these basic equations (and especially their boundary conditions) will be intimately related to the (rather artificial and possibly speculative) assumptions adopted if comparing them with the original conservation laws of the 3-equation system and their constitutive equations as well. The problem of a correct description of the interfacial reaction between the phases and the wall remains. Hence, very often no consistency between different separate-phase models due to their underlying assumptions can be stated. Another problem arises from the fact that special methods have to be foreseen to describe the moving boiling boundary or mixture level (or at least to estimate their ‘condensed’ levels) in such a mixture fluid (see, for example, the ‘Level Tracking’ method in TRAC). Additionally, these methods show often deficiencies in describing extreme situations such as the treatment of single- and two-phase flow at the ceasing of natural circulation, the power situations if decreasing to zero etc. The codes are sometimes very inflexible, especially if they have to provide to a very complex physical system also elements which belong not to the usual class of ‘thermal-hydraulic coolant channels’. These can, for example, be nuclear kinetic considerations, heat transfer out of a fuel rod or through a tube wall, pressure build-up within a compartment, time delay during the movement of an enthalpy front along a downcomer, natural circulation along a closed loop, parallel channels, inner loops etc.

However, despite of these difficulties the ‘Separate-Phase Models’ have become increasingly fashionable and dominant in the last decades of thermal-hydraulics as demonstrated by the widely-used codes TRAC (Lilles et al.,1988, US-NRC, 2001a), CATHENA (Hanna, 1998), RELAP (US-NRC,2001b, Shultz,2003), CATHARE (Bestion,1990), ATHLET (Austregesilo et al., 2003, Lerchl et al., 2009).

Within the scope of reactor safety research very early activities at the Gesellschaft für Anlagen- und Reaktorsicherheit (GRS) at Garching/Munich have been started too, developing thermal-hydraulic models and digital codes which could have the potential to describe in a detailed way the overall transient and accidental behaviour of fluids flowing along a reactor core but also the main components of different Nuclear Power Plant (NPP) types. For one of these components, namely the natural circulation U-tube steam generator together with its feedwater and main steam system, an own theoretical model has been derived. The resulting digital code UTSG could be used both in a stand-alone way but also as part of more comprehensive transient codes, such as the thermal-hydraulic GRS system code ATHLET. Together with a high level simulation language GCSM (General Control Simulation Module) it could be taken care of a manifold of balance-of-plant (BOP) actions too. Based on the experience of many years of application both at the GRS and a number of other institutes in different countries but also due to the rising demands coming from the safety-related research studies this UTSG theory and code has been continuously extended, yielding finally a very satisfactory and mature code version UTSG-2.

During the research work for the development of an enhanced version of the code UTSG-2 it arose finally the idea to establish an own basic element which is able to simulate the thermal-hydraulic mixture-fluid situation within any type of cooled or heated channel in an as general as possible way. It should have the aim to be applicable for any modular construction of complex thermal-hydraulic assemblies of pipes and junctions. Thereby, in contrast to the above mentioned class of ‘separate-phase’ modular codes, instead of separating the phases of a mixture fluid within the entire coolant channel an alternative theoretical approach has been proposed, differing both in its form of application but also in its theoretical background. To circumvent the above mentioned difficulties due to discontinuities resulting from the spatial discretization of a coolant channel, resulting eventually in nodes where a transition from single- to two-phase flow and vice versa can take place, a special and unique concept has been proposed. Thereby it has been assumed that each coolant channel can be seen as a (basic) channel (BC) which can, according to their different flow regimes, be subdivided into a number of sub-channels (SC-s). It is clear that each of these SC-s can consist of only two types of flow regimes. A SC with just a single-phase fluid, containing exclusively either sub-cooled water, superheated steam or supercritical fluid, or a SC with a two-phase mixture. The theoretical considerations of this ‘Separate-Region Approach’ can then (within the class of mixture-fluid models) be restricted to only these two regimes. Hence, for each SC type, the ‘classical’ 3 conservation equations for mass, energy and momentum can be treated in a direct way. In case of a sub-channel with mixture flow these basic equations had to be supported by a drift flux correlation (which can take care also of stagnant or counter-current flow situations), yielding an additional relation for the appearing fourth variable, namely the steam mass flow.

The main problem of the application of such an approach lies in the fact that now also varying SC entrance and outlet boundaries (marking the time-varying phase boundary positions) have to be considered with the additional difficulty that along a channel such a SC can even disappear or be created anew. This means that after an appropriate nodalization of such a BC (and thus also it’s SC-s) a 'modified finite volume method' (among others based on the Leibniz Integration Rule) had to be derived for the spatial discretization of the fundamental partial differential equations (PDE-s) which represent the basic conservation equations of thermal-hydraulics for each SC. Furthermore, to link within this procedure the resulting mean nodal with their nodal boundary function values an adequate quadratic polygon approximation method (PAX) had to be established. The procedure should yield finally for each SC type (and thus also the complete BC) a set of non-linear ordinary differential equations of 1st order (ODE-s).

It has to be noted that besides the suggestion to separate a (basic) channel into regions of different flow types this special PAX method represents, together with the very thoroughly tested packages for drift flux and single- and two-phase friction factors, the central part of the here presented ‘Separate - Region Approach’. An adequate way to solve this essential problem could be found and a corresponding procedure established. As a result of these theoretical considerations an universally applicable 1D thermal-hydraulic drift-flux based separate-region coolant channel model and module CCM could be established. This module allows to calculate automatically the steady state and transient behaviour of the main characteristic parameters of a single- and two-phase fluid flowing within the entire coolant channel. It represents thus a valuable tool for the establishment of complex thermal-hydraulic computer codes. Even in the case of complicated single- and mixture fluid systems consisting of a number of different types of (basic) coolant channels an overall set of equations by determining automatically the nodal non-linear differential and corresponding constitutive equations needed for each of these sub- and thus basic channels can be presented. This direct method can thus be seen as a real counterpart to the currently preferred and dominant ‘separate-phase models’.

To check the performance and validity of the code package CCM and to verify it the digital code UTSG-2 has been extended to a new and advanced version, called UTSG-3. It has been based, similarly as in the previous code UTSG-2, on the same U-tube, main steam and downcomer (with feedwater injection) system layout, but now, among other essential improvements, the three characteristic channel elements of the code UTSG-2 (i.e. the primary and secondary side of the heat exchange region and the riser region) have been replaced by adequate CCM modules.

It is obvious that such a theoretical ‘separate-region’ approach can disclose a new way in describing thermal-hydraulic problems. The resulting ‘mixture-fluid’ technique can be regarded as a very appropriate way to circumvent the uncertainties apparent from the separation of the phases in a mixture flow. The starting equations are the direct consequence of the original fundamental physical laws for the conservation of mass, energy and momentum, supported by well-tested heat transfer and single- and two-phase friction correlation packages (and thus avoiding also the sometimes very speculative derivation of the ‘closure’ terms). In a very comprehensive study by (Hoeld, 2004b) a variety of arguments for the here presented type of approach is given, some of which will be discussed in the conclusions of chapter 6.

The very successful application of the code combination UTSG-3/CCM demonstrates the ability to find an exact and direct solution for the basic equations of a 'non-homogeneous drift-flux based thermal-hydraulic mixture-fluid coolant channel model’. The theoretical background of CCM will be described in very detail in the following chapters.

For the establishment of the corresponding (digital) module CCM, based on this theoretical model very specific methods had to be achieved. Thereby the following points had to be taken into account:

  • The code has to be easily applicable, demanding only a limited amount of directly available input data. It should make it possible to simulate the thermal-hydraulic mixture-fluid situation along any cooled or heated channel in an as general as possible way and thus describe any modular construction of complex thermal-hydraulic assemblies of pipes and junctions. Such an universally applicable tool can then be taken for calculating the steady state and transient behaviour of all the characteristic parameters of each of the appearing coolant channels and thus be a valuable element for the construction of complex computer codes. It should yield as output all the necessary time-derivatives and constitutive parameters of the coolant channels required for the establishment of an overall thermal-hydraulic code.

  • It was the intention of CCM that it should act as a complete system in its own right, requiring only BC (and not SC) related, and thus easily available input parameters (geometry data, initial and boundary conditions, parameters resulting from the integration etc.). The partitioning of BC-s into SC-s is done at the beginning of each recursion or time-step automatically within CCM, so no special actions are required of the user.

  • The quality of such a model is very much dependent on the method by which the problem of the varying SC entrance and outlet boundaries can be solved. Especially if they cross BC node boundaries during their movement along a channel. For this purpose a special ‘modified finite element-method’ has been developed which takes advantage of the ‘Leibniz’ rule for integration (see eq.(15)).

  • For the support of the nodalized differential equations along different SC-s a ‘quadratic polygon approximation’ procedure (PAX) was constructed in order to interrelate the mean nodal with the nodal boundary functions. Additionally, due to the possibility of varying SC entrance and outlet boundaries, nodal entrance gradients are required too (See section 3.3).

  • Several correlation packages such as, for example, packages for the thermodynamic properties of water and steam, heat transfer coefficients, drift flux correlations and single- and two-phase friction coefficients had to be developed and implemented (See sections 2.2.1 to 2.2.4).

  • Knowing the characteristic parameters at all SC nodes (within a BC) then the single- and two-phase parameters at all node boundaries of the entire BC can be determined. And also the corresponding time-derivatives of the characteristic averaged parameters of coolant temperatures resp. void fraction over these nodes. This yields a final set of ODE-s and constitutive equations.

  • In order to be able to describe also thermodynamic non-equilibrium situations it can be assumed that each phase is represented by an own with each other interacting BC. For these purpose in the model the possibility of a variable cross flow area along the entire channel had to be considered as well.

  • Within the CCM procedure two further aspects play an important role. These are, however, not essential for the development of mixture-fluid models but can help enormously to enhance the computational speed and applicability of the resulting code when simulating a complex net of coolant pipes:

  • The solution of the energy and mass balance equations at each intermediate time step will be performed independently from momentum balance considerations. Hence the heavy CPU-time consuming solution of stiff equations can be avoided (Section 3.6).

  • This decoupling allows then also the introduction of an ‘open’ and ‘closed channel’ concept (see section 3.11). Such a special method can be very helpful in describing complex physical systems with eventually inner loops. As an example the simulation of a 3D compartment by parallel channels can be named (Jewer et al., 2005).

The application of a direct mixture-fluid technique follows a long tradition of research efforts. Ishii (1990), a pioneer of two-fluid modelling, states with respect to the application of effective drift-flux correlation packages in thermal-hydraulic models: ‘In view of the limited data base presently available and difficulties associated with detailed measurements in two-phase flow, an advanced mixture-fluid model is probably the most reliable and accurate tool for standard two-phase flow problems’. There is no new knowledge available to indicate that this view is invalid.

Generally, the mixture-fluid approach is in line with (Fabic, 1996) who names three strong points arguing in favour of this type of drift-flux based mixture-fluid models:

  • They are supported by a wealth of test data,

  • they do not require unknown or untested closure relations concerning mass, energy and momentum exchange between phases (thus influencing the reliability of the codes),

  • they are much simpler to apply,

and, it can be added,

  • discontinuities during phase changes can be avoided by deriving special solution procedures for the simulation of the movement of these phase boundaries,

  • the possibility to circumvent a set of ‘stiff’ ODE-s saves an enormous amount of CPU time which means that the other parts of the code can be treated in much more detail.

A documentation of the theoretical background of CCM will be given in very condensed form in the different chapters of this article. For the establishment of the corresponding (digital) module CCM, based on this theoretical model, very specific methods had to be achieved.

The here presented article is an advanced and very condensed version of a paper being already published in a first Open Access Book of this INTECH series (Hoeld, 2011a). It is updated to the newest status in this field of research. An example for an application of this module within the UTSG-3 steam generator code is given in (Hoeld 2011b).


2. Thermal-hydraulic drift-flux based mixture fluid approach

2.1. Thermal-hydraulic conservation equations

Thermal-hydraulic single-phase or mixture-fluid models for coolant channels or, as presented here, for each of the sub-channels are generally based on a number of fundamental physical laws, i.e., they obey genuine conservation equations for mass, energy and momentum. And they are supported by adequate constitutive equations (packages for thermo-dynamic and transport properties of water and steam, for heat transfer coefficients, for drift flux, for single- and two-phase friction coefficients etc.).

In view of possible applications as an element in complex thermal-hydraulic ensembles outside of CCM eventually a fourth and fifth conservation law has to be considered too. The fourth law, namely the volume balance, allows then to calculate the transient behaviour of the overall absolute system pressure. Together with the local pressure differences then the absolute pressure profile along the BC can be determined. The fifth physical law is based on the (trivial) fact that the sum of all pressure decrease terms along a closed loop must be zero. This is the basis for the treatment of the thermal-hydraulics of a channel according to a ‘closed channel concept’ (See section 3.11). It refers to one of the channels within the closed loop where the BC entrance and outlet pressure terms have to be assumed to be fixed. Due to this concept then the necessary entrance mass flow term has be determined in order to fulfil the demand from momentum balance.

2.1.1. Mass balance (Single- and two-phase flow)

t{A[(1-α)rW+αrS]}+zG=0 E1

Containing the density terms ρW and ρS for sub-cooled or saturated water and saturated or superheated steam, the void fraction α and the cross flow area A which can eventually be changing along the coolant channel. It determines, after a nodalization, the total mass flow G=GW+GS at each node outlet in dependence of its node entrance value.

2.1.2. Energy balance (Single- and two-phase flow)

t{A[(1-α)rWhW+αrShS-P]}+z[GWhW+GShS] = qTWL= UTWqTWF= A qD              E2

Containing the enthalpy terms hW and hS for sub-cooled or saturated water and saturated or superheated steam. As boundary values either the ‘linear power qTWL’, the ‘heat flux qTWF’ along the heated (or cooled) tube wall (with its heated perimeter UTW) or the local ‘power density term qD’ (being transferred into the coolant channel with its cross section A) are demanded to be known (See also sections 2.2.4 and 3.5). They are assumed to be directed into the coolant (then having a positive sign).

2.1.3. Momentum balance (Single- and two-phase flow)


describing either the pressure differences (at steady state) or (in the transient case) the change in the total mass flux (GF =G/A) along a channel (See chapter 3.10).

The general pressure gradient (Pz) can be determined in dependence of

  • the mass acceleration


with vS and vW denoting steam and water velocities given by the eqs.(19) and (20),

  • the static head

(Pz)s=- cos(ΦZG) gCρS+(1-α)ρw]E5

with ΦZG representing the angle between z-axis and flow direction. Hence

cos(ΦZG)= ±ΔzEL/ΔzL, with ΔzL denoting the nodal length and ΔzEL the nodal elevation height (having a positive sign at upwards flow).

  • the single- and/or two-phase friction term

(Pz)F= - fRGF|GF|2dHWρE6

with a friction factor derived from corresponding constitutive equations (section 2.2.2)

and finally

  • the direct perturbations (P/z)X from outside, arising either by starting an external pump or considering a pressure adjustment due to mass exchange between parallel channels.

2.2. Constitutive equations

For the exact description of the steady state and the transient behaviour of single- or two-phase fluids a number of mostly empirical constitutive correlations are, besides the above mentioned conservation equations, demanded. To bring a structure into the manifold of existing correlations established by various authors, to find the best fitting correlations for the different fields of application and to get a smooth transfer from one to another of them special and effective correlation packages had to be developed. Their validities can be and has been tested out-of-pile by means of adequate driver codes. Obviously, my means of this method improved correlations can easily be incorporated into the existing theory.

A short characterization of the main packages being applied within CCM is given below. For more details see (Hoeld, 2011a).

2.2.1. Thermodynamic and transport properties of water and steam

The different thermodynamic properties for water and steam (together with their derivatives with respect to P and T, but also P and h) demanded by the conservation and constitutive equations have to be determined by applying adequate water/steam tables. This is, for light-water systems, realized in the code package MPP (Hoeld, 1996 and 2011a).

Then the time-derivatives of these thermodynamic properties which respect to their independent local parameters (for example of an enthalpy term h) can be represented as

ddth(z,t) = ddth[T(z,t),P(z,t)] =  hTddtTMn(z,t) + hPddtPMn(z,t)  E7

Additionally, corresponding thermodynamic transport properties such as ‘dynamic viscosity’ and ‘thermal heat conductivity’ (and thus the ‘Prantl number’) are asked from some constitutive equations too as this can be stated, for example, for the code packages MPPWS and MPPETA (Hoeld, 1996). All of them have been derived on the basis of tables given by (Schmidt and Grigull,1982) and (Haar et al., 1988).

Obviously, the CCM method is also applicable for other coolant systems (heavy water, gas) if adequate thermodynamic tables for this type of fluids are available.

2.2.2. Single and two-phase friction factors

The friction factor fR needed in eq.(6) can in case of single-phase flow be set, as proposed by (Moody, 1994), equal to the Darcy-Weisbach single-phase friction factor.

The corresponding coefficient for two-phase flow has to be extended by means of a two-phase multiplier Φ 2PF2as recommended by (Martinelli-Nelson, 1948).

For more details see again (Hoeld, 2011a).

2.2.3. Drift flux correlation

Usually, the three conservation equations (1), (2) and (3) demand for single-phase flow the three parameters G, P and T as independent variables. In case of two-phase flow, they are, however, dependent on four of them, namely G, P, α and GS. This means, the set has to be completed by an additional relation. This can be achieved by any two-phase correlation, acting thereby as a ‘bridge’ between GS and α. For example, by a slip correlation. However, to take care of stagnant or counter-current flow situations too an effective drift-flux correlation seemed here to be more appropriate. For this purpose an own package has been established, named MDS (Hoeld, 2001 and 2002a). Due to the different requirements in the application of CCM it turned out that it has a number of advantages if choosing the ‘flooding-based full-range’ Sonnenburg correlation (Sonnenburg, 1989) as basis for MDS. This correlation combines the common drift-flux procedure being formulated by (Zuber-Findlay, 1965) and expanded by (Ishii-Mishima, 1980) and (Ishii, 1990) etc. with the modern envelope theory. The correlation in the final package MDS had, however, to be rearranged in such a way that also the special cases of α 0 or α 1 are included and that, besides their absolute values and corresponding slopes, also the gradients of the approximation function can be made available for CCM. Additionally, an inverse form had to be installed (needed, for example, for the steady state conditions) and, eventually, also considerations with respect to possible entrainment effects be taken care.

For the case of a vertical channel this correlation can be represented as

vD= 1.5 vWLIMC0CVD[(1+CVD2)3/2-(1.5+CVD2)CVD] with  vDvD0=916C0vWLIM    if α0E8

where the coefficient CVD is given by


with (in case of a heated or non-heated channel) C0 0 or 1 if α 0 resp. C0 1 if α 1 and corresponding drift velocity terms vD according to eq.(8.)

The resulting package MDS yields in combination with an adequate correlation for the phase distribution parameter C0 relations for the limit velocities vSLIM and vWLIM and thus (independently of the total mass flow G which is important for the theory below) relations for the drift velocity vD with respect to the void fraction α. All of them are dependent on the given 'system pressure P', the 'hydraulic diameter dHY' (with respect to the wetted surface ATW and its inclination angle ΦZG), on specifications about the geometry type (LGTYPE) and, for low void fractions, the information whether the channel is heated or not.

The drift flux theory can thus be expressed in dependence of a (now already on G dependent) steam mass flow (or flux) term

GS=ρ//ρ/αCGC(C0G+Aρ/vD) = AGFS αCGCE10

with the coefficient

CGC= 1-(1-ρ//ρ/)αC01ifαρ//ρ/ifα1E11

Knowing now the fourth variable then, by starting from their definition equations, relations for all the other characteristic two-phase parameters can be established. Such two-phase parameters could be the ‘phase distribution parameter C0’, the ‘water and steam mass flows GW and GS’, ‘drift, water, steam and relative velocities vD, vW,vS and vR’ (with special values for vS if α -> 0 and vW if α -> 1) and eventually the ‘steam quality X’. Their interrelations are shown, for example, in the tables of (Hoeld, 2001 and 2002a). Especially the determination of the steam mass flow gradient

GSα)GS0α=ρ//ρ/(C00G+Ar/vD0) =Ar//vS0    or =0    ifα0 and LHEATD= 0 or 1GS1α=Aρ//ρ/(1+C01(a)) (G -r/ /vSLIM) = Ar/vW1  ifα1E12

will play (as shown, for example, in eq.(52)) an important part, if looking to the special situation that the entrance or outlet position of a SC is crossing a BC node boundary (α 0 or 1). This possibility makes the drift-flux package MDS to an indispensable part in the nodalization procedure of the mixture-fluid mass and energy balance.

At a steady state situation as a result of the solution of the basic (algebraic) set of equations the steam mass flow term GS acts as an independent variable (and not the void fraction α). The same is the case after an injection of a two-phase mixture coming from a ‘porous’ channel or an abrupt change in steam mass flux GFS, as this can take place after a change in total mass flow or in the cross flow area at the entrance of a following BC. Then the total and the steam mass flow terms G and GS have to be taken as the basis for further two-phase considerations. The void fraction α and other two-phase parameters (vD, C0) can now be determined from an inverse (INV) form of this drift-flux correlation (with GS now as input).

2.2.4. Heat transfer coefficients

The nodal BC heat power terms QBMk into the coolant are needed (as explained in section 3.5) as boundary condition for the energy balance equation (2). If they are not directly available (as this is the case for electrically heated loops) they have to be determined by solving an adequate Fourier heat conduction equation, demanding as boundary condition


Such a procedure is, for example, presented in (Hoeld, 2002b, 2011) for the case of heat conduction through a U-tube wall (See also section 3.5).

For this purpose adequate heat transfer coefficients are demanded. This means a method had to be found for getting these coefficients αTW along a coolant channel at different flow regimes. In connection with the development of the UTSG code (and thus also of CCM) an own very comprehensive heat transfer coefficient package, called HETRAC (Hoeld 1988a), has been established.

This classic method is different to the ‘separate-phase’ models where it has to be taken into account that the heat is transferred both directly from the wall to each of the two possible phases but also exchanged between them. There arises then the question how the corresponding heat transfer coefficients for each phase should look like.


3. Coolant channel module CCM

3.1. Channel geometry and finite-difference nodalization

The theoretical considerations take advantage of the fact that, as sketched in fig.1, a ‘basic’ coolant channel (BC) can, as already pointed-out, according to their flow regimes (characterized by the logical LFTYPE = 0, 1, 2 or 3) be subdivided into a number (NSCT) of sub-channels (SC-s), each distinguished by their characteristic key numbers (NSC). Obviously, it has to be taken into account that their entrance and outlet SC-s can now have variable entrance and/or outlet positions.

The entire BC, with its total length zBT = zBA-zBE, can then, for discretization purposes, be also subdivided into a number of (not necessarily equidistant) NBT nodes. Their nodal positions are zBE, zBk (with k=1,NBT), the elevation heights zELBE, zELk, the nodal length ΔzBk=zBk-zBk-1, the nodal elevations ΔzELBk=zELBk-zELBk-1, with eventually also locally varying cross flow and average areas ABk and ABMk=0.5(ABk+ABk-1) and their slopes A Bkz= (ABk-ABk-1)/ΔzBk, a hydraulic diameter dHYBk and corresponding nodal volumes VBMk = ΔzBkABMk. All of them can be assumed to be known from input.

As a consequence, each of the sub-channels (SC-s) is then subdivided too, now into a number of NCT SC nodes with geometry data being identical to the corresponding BC values, except, of course, at their entrance and outlet positions. The SC entrance position zCE and their function fCE are either identical with the BC entrance values zBE and fBE or equal to the outlet values of the SC before. The SC outlet position (zCA) is either limited by the BC outlet (zBA) or characterized by the fact that the corresponding outlet function has reached an upper or lower limit (fLIMCA). This the term represents either a function at the boiling boundary, a mixture level or the start position of a supercritical flow. Such a function follows from the given BC limit values and will, in the case of single-phase flow, be equal to the saturation temperature TSATCA or saturation enthalpies (h/ or h// if LFTYPE=1 or 2). In the case of two-phase flow (LFTYPE=0) it has to be equal to a void fraction of α = 1 or = 0. The moving SC inlet and outlet positions zCE and zCA can (together with their corresponding BC nodes NBCE and NBCA = NBCE+NCT) be determined according to the conditions (zBNk-1 ≤ zCE < zBNk at k = NBCE) and (zBNk-1 ≤ zCA < zBNk at k = NBCA). Then also the total number of SC nodes (NCT=NBCA-NBCE) is given, the connection between n and k (n=k-NBCE with n=1, NCT), the corresponding positions (zNn, zELCE, zELNn), their lengths (ΔzNn=zNn-zNn-1), elevations (ΔzELNn=zELNn-zELNn-1) and volumes (VMn=zNnAMn) and nodal boundary and mean nodal flow areas (ANn, AMn).

Figure 1.

Subdivision of a ‘basic channel (BC)’ into ‘sub-channels (SC-s)’ according to their flow regimes and their discretization

3.2. Spatial discretization of PDE-s of 1-st order (Modified finite element method)

Based on this nodalization the spatial discretization of the fundamental eqs.(1) to (3) can be performed by means of a ’modified finite element method’. This means that if a partial differential equation (PDE) of 1-st order having the general form with respect to a general solution function f(z,t)

tf(z,t) +zH[f(z,t)] = R[f(z,t)]E14

is integrated over the length of a SC node three types of discretization elements can be expected:

  • Integrating a function f(z,t) over a SC node n yields the nodal mean function values fMn,

  • integrating over the gradient of a function f(z,t) yields a difference of functions values (fNn – fNn-1) at their node boundaries

and, finally,

  • integrating over a time-derivative of a function (by applying the 'Leibniz' rule) yields

zNn1(t)zNn(t)tf(z,t)dz =ΔzNn(t)ddtfMn(t) -[fNn(t) -fMn(t)]ddtzNn(t)- [fMn(t)-fNn-1(t)]ddtzNn-1(t)             (n=1, NCT)E15

This last rule plays for the here presented ‘separate-region mixture-fluid approach’ an outstanding part. It allows (together with PAX) to determine in a direct way the time-derivatives of parameters which represent either a boiling boundary, mixture or a supercritical level. This procedure differs considerably from some of the 'separate-phase methods' where, as already pointed out, very often only the collapsed levels of a mixture fluid can be calculated.

3.3. Quadratic polygon approximation procedure PAX

According to the above described three different types of possible discretization elements the solution of the set of algebraic equations will in the steady state case (as shown later-on) yield function values (fNn) at the node boundaries (zNn), the also needed mean nodal functions (fMn) will then have to be determined on the basis of fNn. On the other hand, the solution of the set of ordinary differential equations will in the transient case now yield the mean nodal functions fMn as a result, the also needed nodal boundary values fNn will have to be estimated on the basis of fMn.

It is thus obvious that appropriate methods had to be developed which can help to establish relations between such mean nodal (fMn) and node boundary (fNn) function values. Different to the ‘separate-phase’ models where mostly a method is applied (called ‘upwind or donor cell differencing scheme’) with the mean parameter values to be shifted (in flow direction) to the node boundaries in CCM a more detailed mixture-fluid approach is asked. This is also demanded because, as to be seen later-on from the relations of the sections 3.7 to 3.9, not only absolute nodal SC boundary or mean nodal function values are required but as well also their nodal slopes f Nnsand f Mnstogether with their gradients f Nnzsince according to this approach the length of SC nodes can tend also to zero.

fNns=(fNn-fNn-1)ΔzNnfCEIz(at n=1) =input  or fNn-1s(at n=NCT>1) if DzNn0E16
fMns=2(fMn-fNn-1)ΔzNnfCEIz(at n=1) =input orfNn-1z (at n=NCT>1) if ΔzNn0E17

Hence, for this purpose a special ‘quadratic polygon approximation’ procedure, named 'PAX', had to be developed. It plays (together with the Leibniz rule presented above) an outstanding part in the development of the here presented ‘mixture-fluid model’ and helps, in particular, to solve the difficult task of how to take care of varying SC boundaries (which can eventually cross BC node boundaries) in an appropriate and exact way.

3.3.1. Establishment of an effective and adequate approximation function

The PAX procedure is based on the assumption that the solution function f(z) of a PDE (for example temperature or void fraction) is split into a number of NCT nodal SC functions fn(z,t). Each of them has then to be approximated by a specially constructed quadratic polygon which have to fulfil the following requirements:

  • The node entrance functions (fNn-1) must be either equal to the SC entrance function (fNn-1 = fCE) (if n = 1) or to the outlet function of the node before (if n > 1). This is obviously not demanded for gradients of the nodal entrance functions (except for the last node at n = NCT).

  • The mean function values fMn over all SC nodes have to be preserved (otherwise the balance equations could be hurt).

  • With the objective to guarantee stable behaviour of the approximated functions (for example by excluding 'saw tooth-like behaviour’) it will, in an additional assumption, be demanded that the outlet gradients of the first NCT -1 nodes should be set equal to the slopes between their neighbour mean function values. The entrance gradient of the last node (n= NCT) should be either equal to the outlet gradient of the node before (if n = NCT > 1) or equal to a given SC input gradient (for the special case n = NCT =1). Thus

fNnz=(fz)Nn = 2fMn+1-fMnΔzNn+1+ΔzNn(n=1, NCT-1, if NCT>1)=2ΔzNn(2fNn- 3fMn+ fNn-1)fNn-1(z)ifΔzCA0(n = NCT, if NCT> 1)E18
fNn-1z=fCEz=fCEIz=2ΔzNn(3fMn–fCA- 2fCE)        (n = NCT=1)=Nnz(of the node before)                         (n = NCT> 1)E19

This means, the corresponding approximation function reaches not only over the node n. Its next higher one (n+1) has to be included into the considerations too (except, of course, for the last node). This assumption makes the PAX procedure very effective (and stable). It is a conclusive onset in this method since it helps to smooth the curve, guarantees that the gradients at the upper or lower SC boundary do not show abrupt changes if these boundaries cross a BC node boundary and has the effect that perturbations at channel entrance do not directly affect corresponding parameters of the upper BC nodes.

For the special case of a SC having shrunk to a single node (n=NCT=1) the quadratic approximation demands as an additional input to PAX (instead of the now not available term fMn) the gradient f CEIzat SC entrance. It represents thereby the gradient of either the coolant temperature T CEIzor void fraction α CEIz(in case of single- or two-phase flow entrance conditions). If this parameter is not directly available it can, for example, be estimated by combining the mass and energy balance equations at SC entrance in an adequate way (See Hoeld, 2005). This procedure allows to take care not only of SC-s consisting of only one single node but also of situations where during a transient either the first or last SC of a BC starts to disappear or to be created anew (i.e. zCA zBE or zCE zBA), since now the nodal mean value fMn at n = NCT (for both NCT = 1 or > 1) is no longer or not yet known.

3.3.2. Resulting nodal parameters due to PAX

It can be expected that in the steady state case after having solved the basic set of non-linear algebraic equations (as presented later-on in the sections 3.7, 3.8 and 3.9) as input to PAX the following parameters will be available:

  • Geometry data such as the SC entrance (zCE) and node positions (zNn) (and thus also the SC outlet boundary position zCA as explained in section 3.9) determining then in PAX the number of SC nodes (NCT),

  • the SC entrance function fN0 = fCE and (at least for the special case n=NCT=1) its gradient f CEI(z)

and finally

  • the nodal boundary functions fNn (n=1,NCT) with fCA = fNn at n = NCT and fCA= fLIMCA if zCA < zBA.

Based on these inputs PAX yields then the nodal mean function values fMn (at n=1,NCT)

fMn = (ΔzNn+1+ΔzNn)(2fNn+fNn-1)-ΔzNnfMn+13ΔzNn+1+2ΔzNn(n =1, NCT-1, NCT> 1 ifzCA=zBAor n =1, NCT-2, NCT> 2 ifzCA<zBA)=13(fCA+2 fCE)+16DzCAfCEI(z)(n = NCT= 1)=13(fCA+2 fCE)+16ΔzCAf(fCA- fNn-2)    (fCA- fNn-2)               E20

functions which are needed as initial values for the transient case.

In the transient case it can be expected that after having integrated the set of non-linear ordinary differential equations (ODE-s) (as to be shown again in the sections 3.7, 3.8 and 3.9) the SC outlet position zCA (=zNn) and thus also the total number NCT of SC nodes will now be directly available from the integration. Thereby, as input data needed for the use in the PAX procedure it has now to be distinguished between two cases:

  • if the now known SC outlet position is identical with the BC outlet (zCA=zBA) the mean nodal function values fMn for all NCT nodes (n=1,NCT) should be provided as input


  • if the SC outlet position moves still within the BC (zCA < zBA) the mean nodal function values fMn of only NCT-1 nodes (n=1, NCT-1) are demanded. Additionally, the nodal function limit values fLIMNn (usually saturation temperature values at single-phase flow resp. void fraction values limited by 1 or 0 at mixture flow conditions) are asked. Then the missing mean nodal function fMn of the last SC is not any longer needed, since this function can be determined from eq.(20) in dependence of the known outlet function fNn = fCA = fLIMCA at zNn = zCA.

These nodal input function values can then be taken (together with its input parameter fCE and the nodal positions zBE and zBn at n=1, NCT) as basic points for the PAX procedure, yielding, after rearranging the definition equations of the approximation function in an adequate way, all the other not directly known nodal function parameters of the SC.

Hence, it follows for the special situation of a SC being the last one within the BC (zCA = zBA)

  fNn=12(3fMn- fNn-1) +12ΔzNnΔzNn+1+ΔzNn(fMn+1- fMn)  (n=1, NCT-1 with NCT>1 if zCA= zBA)=3fMn- 2fCE-12ΔzCAfCEIz(n = NCT = 1   if zCA= zBA)=fCA = 2(fMn– fMn-1) + fNn-2                                                (n = NCT > 1  if zCA= zBA)E21

resp. for the case zCA < zBA

  fNn=12(3fMn- fNn-1) +12ΔzNnΔzNn+1+ΔzNn(fMn+1- fMn) (n =1, NCT–2 with NCT> 2 if zCA< zBA)=12(3fMn- fNn-1) +14ΔzNnΔzNn+1+ΔzNn(fLIMCA- fNn-1)   (n =NCT-1 with NCT> 1 if zCA< zBA)=fCA = fLIMCA                                                                             (n = NCT       if zCA< zBA)E22

The last mean nodal function fMn (at n=NCT) which is, for this case (zCA < zBA), not available from the integration, follows now directly from eq.(21) if replacing there fCA by fLIMCA.

Finally, from the eqs.(16), (17) and (19) the slopes and gradients can be determined.

The corresponding time-derivative of the last mean node function which is needed for the determination of the SC boundary time-derivative (see section 3.9) follows (as long as zCA < zBA) by differentiating the relation above

ddtfMn  = fPAXCAt+fPAXCAzddtzCA  (n =NCT if zCA< zBA)  E23

with the coefficients

fPAXCAt=13(ddtfLIMCA+2ddtfCE) -16(fCEIzddtzCE-DzCAddtfCEIz)(n =NCT=1 if zCA<zBA)=ddtfMn-1+12ddtfLIMCA -12ddtfNn-2                               (n =NCT> 1 if zCA< zBA)fPAXCAz=16fCEI(z)or = 0                                                    (n = NCT= 1 or > 1 if zCA< zBA)E24

The differentials ddtfMn-1, ddtzCA, ddtfLIMCA are directly available from CCM and, if NCT=2, the term ddtfNn-2 = ddtfCE from input too. For the case that a SC contains more than two nodes only their corresponding mean values are known, the needed term ddtfNn-2 has thus to be estimated by establishing the time-derivatives of all the boundary functions at the nodes below NCT < 2.

3.3.3. Code package PAX

Based on the above established set of equations a routine PAX had to be developed. Its objective was to calculate automatically either the nodal mean or nodal boundary values (in case of an either steady state or transient situation). The procedure should allow also determining the gradients and slopes at SC entrance and outlet (and thus also outlet values characterizing the entrance parameters of an eventually subsequent SC). Additionally, contributions needed for the calculation of the time-derivatives of the boiling boundary or mixture level can be gained (See later-on the eqs.(66) and (67)).

Figure 2.

Approximation function f(z) along a SC for both steady state and transient conditions after applying PAX (Example)

Before incorporating the subroutine into the overall coolant channel module the validity of the presented PAX procedure has been thoroughly tested. By means of a special driver code (PAXDRI) different characteristic and extreme cases have been calculated. The resulting curves of such a characteristic example are plotted in fig.2. It represents the two approximation curves of an artificially constructed void fraction distribution f(z) = α(z) along a SC with two-phase flow both for the steady state but also transient situation. Both curves (on the basis of fMn and fNn) should be (and are) identical.

3.4. Needed input parameters

3.4.1. Initial conditions

For the start of the transient calculations adequate steady state parameters have to be available as initial conditions.

3.4.2. Boundary conditions

As boundary conditions for both the steady state and especially for transient calculations the following input parameters are expected to be known. Thereby demanding only easily available BC values (They will, within CCM, then be automatically translated into the corresponding SC values):

  • Power profile along the entire BC. This means that either the nodal heat flux terms qFBE and qFBk (at BC entrance and each node k=1,NBT) or qFBE and the nodal power terms QBMk are expected to be known, either directly from input or (as explained in section 2.2.4) by solving the appropriate ‘Fourier heat conduction equation’. From the relation

QBMk12ΔzBk(qLBk+qLBk-1) =12ΔzBk(UTWBkqFBk+UTWBk-1qFBk-1) = VBMkqBMkwith qLBE= UTWBEqFBEand qDBE=UTWBEABEqFBE                         (k=1,NBT)E25

then the other BC nodal terms (qFBk or QBMk, qLBk and qDBk) can be determined too. (Hoeld, 2002b, 2004a and 2011).

  • For normalization purposes at the starting calculation (i.e., at the steady state situation) as an additional parameter the total nominal (steady state) heat power QNOM,0 is asked.

  • Channel entrance temperature TBEIN (or enthalpy hBEIN)

  • System pressure PSYS and its time-derivative (dPSYS/dt), situated at a fixed position either along the BC (entrance, outlet) or even outside of the ensemble. Due to the fast pressure wave propagation each local pressure time-derivative can then be set equal to the change in system pressure (as described in section 3.6).

  • Total mass flow GBEIN at BC entrance together with pressure terms at BC entrance PBEIN and outlet PBAIN. These three parameters are needed for steady state considerations (and partially used for normalization purposes). In the transient case only two of them are demanded as input. The third one will be determined automatically by the model. These allows then to distinguish between the situation of an ‘open’ or ‘closed channel’ concept as this will be explained in more detail in section 3.11.

  • Steam mass flow GSBEIN at BC entrance (=0 or = GBEIN at single- or 0 < GSBEIN < GBEIN at two-phase flow conditions). The corresponding entrance void fraction αBE will then be determined automatically within the code by applying the inverse drift-flux correlation.

Eventually needed time-derivatives of such entrance functions can either be expected to be known directly from input or be estimated from their absolute values.

By choosing adequate boundary conditions then also thermal-hydraulic conditions of other situations can be simulated. For example that of several channel assembles (nuclear power plants, test loops etc.) which can consist of a complex web of pipes and branches (represented by different BC-s, all of them distinguished by their key numbers KEYBC). Even the case of an ensemble consisting of inner loops (for example describing parallel channels) can be treated in an adequate way according to the concept of a ‘closed’ channel (see section 3.11).

3.4.3. Solution vector

The characteristic steady state parameters are determined in a direct way, i.e. calculated by solving the non-linear set of algebraic equations for SC-s (as being presented in the chapters 3.7, 3.8 and 3.9). Thereby, due to the nonlinearities in the set of the (steady state) constitutive equations a recursive procedure in combination with and controlled by the main program has to be applied until a certain convergence in the solution vector can be stated. The results are then combined to BC parameters and transferred again back to the main (= calling) program.

For the transient case, as a result of the integration (performed within the calling program and thus outside of CCM) the solution parameters of the set of ODE-s are transferred after each intermediate time step to CCM. These are (as described in detail also in chapter 4) mainly the mean nodal SC and thus BC coolant temperatures, mean nodal void fractions and the resulting boiling or superheating boundaries. These last two parameters allow then to subdivide the BC into SC-s yielding the corresponding constitutive parameters and the total and nodal length (zNn and ΔzNn) of these SC-s and thus also their total number (NCT) of SC nodes. Finally, the needed SC (and thus BC) time-derivatives can be determined within CCM (as described in the sections 3.7, 3.8 and 3.9) and then transmitted again to the calling program where the integration for the next time step can take place.

3.5. SC power profile

Knowing (as explained in section 3.4.2) the nodal BC power QBMk together with the linear power qLBE at BC entrance the corresponding SC power profile can now be determined too.

Hence, since linear behaviour of the linear nodal power terms within the corresponding BC nodes can be assumed it follows for the ‘linear SC power’ term

qLNn= qLCE= qLBE(=BC entr.) or = (qLCA)of the last node of the SC before(n=0 if zCE= or > zBE)= qLBk                                                                  (n=1,NCTand k=n+NBCEif LFTYPE=2)                          = qLBk                 (n=1,NCT-1and, if zCA= zBA, n= NCTwith k=n+NBCE)=qLCA= qLBk-1(qLBk– qLBk-1)ΔzCAΔzBk(n=NCTand k=NBCAif zCA< zBA)E26

for the ‘nodal SC power term’

QMn=12ΔzNn(qLNn+qLNn-1)(n=1,NCT)= QBMk– (QMCA)of the last node of the SC before   (n=1  with k=1+NBCE)= QBMk                                                                 (n=2, NCT-1and k=n+NBCE)= QMCA =DzCA[qLBk-1+12 (qLBk– qLBk-1) (qLBk– qLBk-1)ΔzCAΔzBk](n=NCTand k=NBCAif zCA< zBA)E27

and the ‘mean nodal’ and ‘nodal boundary SC power density’ terms (qMn and qNn)

qMn=QMnVMnand (independently of the node length) =12AMn(qLNn-1+qLNn) (n =1,NCT)qNn = qCE= qBE or  = (qCA)of the last node of the SC before                 (n=0 if zCE= or > zBE)              = 2qMn- qNn-1                                                                                         (n =1,NCT)E28

These last parameters are, since independent of ΔzNn, very useful for equations where it has to be expected that ΔzNn -> 0 (as this can be seen later-on by the eqs.(32) and (49)).

The SC length zCA (and thus that of its last node ΔzCA) are (for the case zCA,0 < zBA) only known in the transient case (as a result of the integration procedure). For steady state conditions the term QMCA,0 follows from energy balance considerations (eqs.(40) and (59)). Then, in a reverse manner, ΔzCA,,0 can be calculated (See eq.(65) in section 3.9).

3.6. Decoupling of mass and energy balance from momentum balance equations

Treating the conservation equations in a direct way produces due to elements with fast pressure wave propagation (which are responsible for very small time constants) a set of ‘stiff’ ODE-s. This has the consequence that their solution turns out to be enormously CPU-time consuming. Hence, to avoid this costly procedure CCM has been developed with the aim to decouple the mass and energy from their momentum balance equations. This can be achieved by determining the thermodynamic properties of water and steam in the energy and mass balance equations on the basis of an estimated pressure profile P(z,t). Thereby the pressure difference terms from a recursive (or a prior computational time step) will be added to an eventually time-varying system pressure PSYS(t), known from boundary conditions. After having solved the two conservation equations for mass and energy (now separately from and not simultaneously with the momentum balance) the different nodal pressure gradient terms can (by the then following momentum balance considerations) be determined according to the eqs.(4), (5) and (6).

It can additionally be assumed that according to the very fast (acoustical) pressure wave propagation along a coolant channel all the local pressure time-derivatives can be replaced by a given external system pressure time-derivative, i.e.,


By applying the above explained ‘intelligent’ (since physically justified) simplification in CCM the small, practically negligible, error in establishing the thermodynamic properties on the basis of such an estimated pressure profile can be outweighed by the enormous benefit substantiated by two facts:

  • Avoidance of the very time-consuming solution of stiff equations,

  • the calculation of the mass flow distribution into different channels resulting from pressure balance considerations can, in a recursive way, be adapted already within each integration time step, i.e. there is no need to solve the entire set of differential equations for this purpose (See ‘closed channel’ concept in section 3.11).

3.7. Thermal-hydraulics of a SC with single-phase flow (LFTYPE > 0)

The spatial integration of the two PDE-s of the conservation eqs.(1) and (2) over a (single-phase) SC node n yields (by taking into account the rules from section 3.2, the relations from the eqs.(7) and (29), the possibility of a locally changing nodal cross flow area along the BC and the fact that eventually VMn -> 0 ) for the transient case the relation

  • a relation for the total nodal mass flow

GNn= GNn-1- VMnMnTddtTMnMnPddtPSYS) + (rNn-rMn)ANnPSYS) + (rNn-rMn)ANnddtzNn+ (rMn-rNn-1)ANn-1ddtzNn-1                   (n=1,NCT), LFTYPE> 0E30


  • the time-derivative for the mean nodal coolant temperature (if eliminating the term GNn in the resulting equation by inserting from the equation above):

ddtTMn = TTnt+TTCAzddtzNn with ddtzNn=  ddtzCAor = 0 if n =NCT or <NCT         (n=1, NCTand zCA=zBA) or (n=1, NCT-1 and zCA< zBA), LFTYPE> 0E31

with the coefficients

TTntqMn-qGn+qPnρMnhMnTCTMn+ TddtzNn-1  withddtzNn-1=ddtzCE or = 0 if n = 1 or > 1=1VMnρMnhMnTCTMn[QMn-GNn-1(hNn- hNn-1)+VMnqPn]+ T+ TTCEzddtzNn-1E32
CTMn = 1 -ρMnTρMnhNn-hMnhMnT=1-ρMnTρMn(TNn-TMn)E33
qGn =GNn-1AMnhNnsE34
qPn = CRHPnddtPSYS  and CRHPn= 1-rMnhMnP+ρMnP(hNn- hMn)E35
 qZn =12AMn{ANnrMnhNn(s)ddtzNn+ANn-1[ρMnhMn(s)-2(ρMn-ρNn-1) hNn(s))]ddtzNn-1}E36
TTCEz=12[TMns-2(1-ρCEρMn)TNns]ACEAMnCTMnor  = 0  at n = 1 and LFTYPE > 0(if either zCE > zBE or  zCE= zBE)E37
TTCAz= 0  or  =12ANnAMnCTMnTNns(n < NCT and zCA< zBAor n = NCTand zCA= zBA), LFTYPE> 0E38

In the transient case the mean nodal coolant temperature value TMn is at the begin of each (intermediate) time step known. This either, at the first time step, from steady state considerations (in combination with PAX) or as a result of the integration procedure. Hence, additional parameters needed in the relations above can be determined too. From the PAX procedure it follow also the SC nodal terms TNn, T Nn(s)and the slopes T Mn(s)resp., for the case that ΔzNn 0, their gradients. Finally, using the water/steam tables (Hoeld, 1996), also their nodal enthalpies are fixed.

From the integration procedure then also the SC outlet position zCA is provided, allowing to determine the total number of SC nodes (NCT) too. The situation that zCA = zBT means the SC nodal boundary temperature values have, within the entire BC, not yet reached their limit values (TLIMNn=TSATNn). Thus NCT= NBCA with NBCA=NBT–NBCE. Otherwise, if zCA < zBT this limit is reached (at node n), then NCT = n. Obviously, the procedure above yields also the time-derivative of the SC outlet position moving within this channel (As described in section 3.9).

The steady state part of the total nodal mass flow (charaterized by the index 0) follows from the basic non-linear algebraic equation (30) if setting there the time-derivative equal to 0:

GNn,0= GCA,0 = GCE,0 = GBA,0 = GBE,0           (n=1,NCT), LFTYPE>0E39

Treating eq.(31) in a similar way and multiplying the resulting relation by VMn,0 yields the steady state nodal temperature resp. enthalpy terms

hNn,0= hNn-1,0+QMn,0GBE,0hNn,0/orhNn,0//(with hNn-1,0= hCE,0at n=1)(if LFTYPE = 1 or = 2 at n = 1, NBT-NBCE)E40

restricted by their saturation values. Then, with regard to eq.(27) which results from energy balance considerations, the nodal power term for the last SC node has to obey the relation

QMCA,,0= QMn,,0= (hCA,0/- hNn-1,0) GBE,0        (if n =NCT< NBCAat LFTYPE= 1)= (hNn-1,0–hCA,0//) GBE,0        (if n =NCT< NBCAat LFTYPE= 2)E41

Thus, for the steady state case, NCT is fixed too:

 NCT= n < NBCA= NBT- NBCE+1 and zCA,0(=zNn,0) < zBA (if hNn,0= hNn,0/at LFTYPE= 1)NCT= n = NBCA and zCA,0= zBA                                          (if hNn,0< hNn,0/at LFTYPE= 1)E42

with similar relations for the case LFTYPE = 2.

From the resulting steady state enthalpy value hNn,0 follows then (by using the thermodynamic water/steam tables) the corresponding coolant temperature value TNn,0 (with TNn,0 = TSATNn,0 if n = NCT and zCA < zBA) and, by applying the PAX procedure (see section 3.3), their mean nodal temperature and enthalpy values TMn,0 and hMn,0, parameters which are needed as start values for the transient calculations. Obviously, due to the non-linearity of the basic steady state equations, this procedure has to be done in a recursive way.

It can additionally be stated that both the steady state and transient two-phase mass flow parameters get the trivial form

GSNn = GSNn,0  = 0 resp. GWNn= GCE and GWNn,0= GBE,0    (n=1, NCT, if LFTYPE= 1) GSNn  = GCE  and GSNn,0= GBE,0  resp.  GWNn= GWNn,0= 0   (n=1, NCT, if LFTYPE= 2)E43


ddtαMn= 0 , αNn = αNn,0 = αMn,0 = 0   or  = 1 and XNn,0=hNn,0-hNn,0/hSWNn,0(n=1, NCT, if LFTYPE=1 or = 2)E44

3.8. Thermal-hydraulics of a SC with two-phase flow (LFTYPE = 0)

The spatial integration of the two PDE-s of the conservation) eqs.(1) and (2) (now over the mixture-phase SC nodes n) can be performed by again taking into account the rules from section 3.2, the relations from the eqs.(7) and (29), by considering the possibility of locally changing nodal cross flow areas along the BC) and the fact that eventually VMn -> 0. This yields then relations for

  • the total mass flow term

GNn = GNn-1+VMn(ρ’-ρ’’)Mn(ddtαMnGPntGZnt)(n=1,NCT, LFTYPE= 0)E45

with, if neglecting thereby small differences between mean and nodal thermodynamic saturation values, the coefficients

αGZntCEzddtzCE=12ACEAMnαMnsddtzCE   (n = 1 and zCE> zBE)=0(1 < n < NCT)= αCAzddtzCAACAAMnCAs-12αMns)ddtzCA (n = NCT, NCT>1 and zCA< zBA)E47


  • the mean nodal void fraction time-derivative

ddtαMn=  αASntAPntCAzddtzNnCEzddtzNn-1   (n=1,NCTand zCA=zBA) or (n=1, NCT-1, NCT>1 and zCA<zBA), LFTYPE=0E48

with the coefficients

αASntAQntAGnt=1ρMn//hSWMn [qMn- hSWMnGSNnsAMn]=1VMnρMn//hSWMn[QMn(GSNn- GSNn-1) hSWMn]       E49
 αAPnt= CRHPnddtPSYSwith CRHPn=1(ρ//hSW)Mn[(1-a)ρh’P+α(ρh’’P+ρ’’PhSW) - 1]MnE50
αAZntGZnt(see eq.(47))E51
GSNns=ΔGSNnΔzNnGSNnz=GSNnααNnsif αNnαCE= 0 (at n=1)or αNnαCA= 1 (at n= NCT)E52

It can again be expected that at the begin of each (intermediate) time step the mean nodal void fraction values αMn are known. This either from steady state considerations (at the start of the transient calculations) or as a result of the integration procedure. Hence, the additional parameters needed in the relations above can be determined too. From the PAX procedure it follow their nodal boundary void fraction terms αNn together with their slopes α Nn(s)and α Mn(s)resp. gradient α Nn(z)and thus, as shown both in section 2.2.3 but also in the tables given by Hoeld (2001 and 2002a), all the other characteristic two-phase parameters (steam, water or relative velocities, steam qualities etc). Obviously, due to the non-linearity of the basic equations the steady state solution procedure has to be performed in a recursive way.

If, in the transient case, the SC nodal boundary void fraction αNn does (within the entire BC) not reach its limit value (αLIMNn=1 or 0) it follows as total number of SC nodes NCT=NBT–NBCE and zNn (at n=NCT) =zCA=zBT. Otherwise, if this limit is reached (at node n), NCT = n, αNn=1 (or =0) with zCA (< zBT) resulting from the integration. Then, from the procedure above also the time-derivative of the boiling boundary, moving within BC, can be established (as this will be discussed in section 3.9).

Hence, it follows a relation for the steam mass flow gradients

GSNn(α)= ACEvS0ρNn//with vS0= vS (atαNn= 0)    (n=1       and αNnαCE= 0)==ACAvW1ρNn/with vW1= vW(atαNn=1)   (n=NCT  andαNnαCA= 1)E53

If eliminating the term ddtαMn in eq.(45) by inserting from eq.(48) yields a relation between GSNn and GNn

 GNn+ (ρ/ρ//-1)MnGSNn= GXn   (n=1,NCT, LFTYPE=0)E54

The ‘auxiliary’ mass flow term GXn is directly available since it refers only to already known parameters (for example taken from the node before or the power profile)

GXn= GNn-1+ (ρ/ρ//-1)Mn(GSNn-1+QMnhSWMn)) - VMn(ρ/-ρ//)MnAPntGPnt)(n=1,NCT, LFTYPE=0)E55

A similar relation can be established if starting from the drift flux correlation (eq.(10)) and taking advantage of the fact that the needed drift velocity vDNn and the phase distribution parameter C0Nn are independent from the total mass flow GNn (and can thus be determined before knowing GNn). This term (GNn) results then by combining the eqs.(54) and (10)

GNn =GXn-(AαvDρ/CDC)Nn1+(αC0CDC)Nn(n=1,NCT)E56

with the coefficient

CDCNn= (ρ/ρ//-1)Mn (ρ//ρ/CGC)NnE57

From the drift flux correlation package (Hoeld et al., 1992) and (Hoeld, 1994) follow then all the other characteristic two-phase parameters. These are especially the nodal steam mass flow GSNn and, eventually, the slope α Nn(z)resp., according to eq.(53), G SNn(s). Then, finally, from eq.(48) (or eq.(45)) the mean nodal void fraction time-derivative ddtαMn will result, needed for the next integration step.

Obviously, at a mixture flow situation the mean nodal temperature and enthalpy terms are equal to their saturation values

TMn= TSAT(PMn) resp. hMn= h/(PMn) or = h//(PMn)    (n= 1, NCTand LFTYPE= 0)E58

and are thus only dependent on the local resp. system pressure value.

Relations for the steady state case can be derived if setting in the eqs.(45) and (48) (resp. the eq.(49)) the time-derivatives equal to 0. For the total mass flow parameters a similar relation as already given for the single-phase flow (see eq.(48)) is valid, yielding GNn,0 = GBE,0.

Hence, one obtains for the (steady state) nodal steam mass flow

GSNn,0= GSNn-1,0+QMn,0hSWMn,0GNn,0 = GBE,0      (n= 1, NBT–NBCEand LFTYPE= 0)E59

if knowing the term QMCA,,0 = QMn,,0 from eq.(27).

Then the total number (NCT) of SC nodes is given too

NCT= n < NBT- NBCE  and zCA,0(=zNn,0) < zBA       (if GSNn,0= GBE,0 and LFTYPE= 0)       NCT= NBT- NBCE          and zCA,0 = zBA                    (if GSNn,0< GBE,0 and LFTYPE= 0)E60

and also the corresponding steam quality parameter

XNn,0=GSMn,0GNn,0(n = 1, NCTif LFTYPE= 0)E61

The nodal boundary void fraction values αNn,0 can now be determined by applying the inverse drift-flux correlation, the mean nodal void fraction value αMn,0 from the PAX procedure. All of them are needed as starting values for the transient calculation.

3.9. SC boundaries

The SC entrance position zCE (= zNn at n=0) is either (for the first SC within the BC) equal to BC entrance zBE or equal to the SC outlet boundary of the SC before.

In the steady state case the SC outlet boundary (= boiling boundary zBB,0 or mixture level zML,0) can be represented as

zCA,0= zBA  (n = NCTand zCA,0= zBA)E62


zCA,0 = zNn-1 +DzCA,0               (n = NCTand zCA,0< zBA)E63

The (steady state) numbers (NCT) of (single- or two-phase) SC-s are already determined by the eqs.(42) or (60), the nodal power terms QMCA,0 by the eqs.(41) or (59). Hence, after dividing eq.(27) by ΔzBkqLBk-1,0 one gets an algebraic quadratic equation of the form

12qLBk,0-qLBk-1,0qLBk-1,0 (ΔzCA,0ΔzBk)2+ΔzCA,0ΔzBk-QMCA,0zBkqLBk-1,0= 0(n = NCTandk= NBCAif NCT< NBCA)E64

yielding finally as solution

              ΔzCA,0=ΔzBk       (n = NCTand k= NBCA  if  NCT= NBT)= ΔzBkqLBk-1,0qLBk-1,0-qLBk,0[1 -1 - 2(1 - qLBK,0qLBk-1,0 ) QCMA,0ΔzBk  qLBMk-1,0 ](k= NBCA  if  NCT< NBT)QMCA,0qLBk-1,0[1 + (1 - qLBk,0qLBk-1,0)QMCA,0zBkqLBk-1,0]if qLBk,0   qLBk-1,0E65

Then, from the relations in section 3.5, also the other characteristic steady state power terms can be calculated.

In the transient case the SC outlet boundary zCA (=boiling boundary or mixture level) follows (and thus also ΔzCA and NCT ), as already pointed-out, directly from the integration procedure. During the transient this boundary can move along the entire BC (and thereby also cross BC node boundaries). A SC can even shrink to a single node (NCT =1), start to disappear or to be created anew. Then, if ΔzCA -> 0, in the relations above the slope in the vicinity of such a boundary has to be replaced by a gradient (determined in PAX).

The mean nodal coolant temperature or, if LFTYPE=0, void fraction of the last SC node is interrelated by the PAX procedure with the locally varying SC outlet boundary zCA. Hence, in a transient situation the time-derivative of only one of these parameters is demanded. The second one follows then from the PAX procedure after the integration.

If combining (in case of single-phase flow) the eqs.(23) and (31), the wanted relation for the SC boundary time derivative can be expressed by

ddtzCA=ddtzBB=TPAXCAt-TTCAtTTCAz-TPAXCAzor = 0       (n = NCT, zCA< zBAor zCA= zBAif LFTYPE> 0)E66

and if taking for the case of a mixture flow the eqs. (23) and (48) into account

ddtzCA=ddtzML=αPAXCAt-αACAtαCAz-αPAXCAzor = 0      (n = NCT, zCA< zBAor zCA= zBAif LFTYPE= 0)E67

If zCA < zBA, the corresponding time-derivatives ddtTMn or ddtαMn of the last SC node (at n=NCT) follow by inserting the terms above into the eqs.(31) or (48). After the integration procedure then the SC outlet boundary zCA (= boiling boundary zBB or mixture level zML) and thus also the total number NCT of SC nodes are given.

3.10. Pressure profile along a SC (and thus also BC)

After having solved the mass and energy balance equations, separately and not simultaneously with the momentum balance, the now exact nodal SC and BC pressure difference terms (ΔPNn = PNn - PNn-1 and ΔPBNn) can be determined for both single- or two-phase flow situations by discretizing the momentum balance eq.(3) and, if applying a modified ‘finite element method’, integrating the eqs.(4) to (6) over the corresponding SC nodes. The total BC pressure difference ΔPBT = PBA - PBE between BC outlet and entrance follows then from the relation

             ΔPBT =ΔPPBT-ΔPGBT      (withΔPGBT,0= 0 at steady state conditions)  E68

with the part

ΔPPBT=ΔPSBT+ΔPABT+ΔPXBT+ΔPFBT+ΔPDBT(withΔPPBT,0=ΔPBTIN,0 at steady state conditions)E69

comprising, as described in section 2.1.3, terms from static head (ΔPSBT), mass acceleration (ΔPABT), wall friction (ΔPFBT) and external pressure accelerations (ΔPXBT due to pumps or other perturbations from outside) and an (only in the transient situation needed) term ΔPGBT, being represented as

ΔPGBT =0zBTddtGFB(z,t) dz = zBTddtGFBMT    at transient conditions=0at steady stateE70

This last term describes the influence of time-dependent changes in total mass flux along a BC (caused by the direct influence of changing nodal mass fluxes) and can be estimated by introducing a ‘fictive’ mean mass flux term GFBMT (averaged over the entire BC)


Its time derivative can then be represented by

ddtGFBMTGFBMT-GFBMTBΔt(ddtGFBMT)at t=tB  if  Δt=t-tB  0  (Index B = begin of time-step)E72

Looking at the available friction correlations, there arises the problem how to consider correctly contributions from spacers, tube bends, abrupt changes in cross sections etc. as well. The entire friction pressure decrease (ΔPFBT) along a BC can thus never be described in a satisfactory manner solely by analytical expressions. To minimize these uncertainties a further friction term, ΔPDBT, had to be included into these considerations.

The corresponding steady state additive pressure difference term ΔPDBT,0 is (according to the eqs.(69) and (70)) fixed since it can be assumed that the total steady state BC pressure difference ΔPBT,0 = ΔPPBT,0 is known from input. It seems, however, to be reasonable to treat this term as a ‘friction’ (or at least the sum with ΔPFBT) and not as a ‘driving’ force. Thus it must be demanded that these terms should remain negative. Otherwise, input terms such as the entire pressure difference along the BC or corresponding friction factors have to be adjusted in an adequate way.

Assuming the general additive pressure difference term ΔPDBT may have the form


means that ΔPDBT is either supplemented with a direct additive term (index FADD) or the friction part is provided with a multiplicative factor (fFMP,0-1). For the additive part it will be assumed to be the (1-εDPZ)-th part of the total additional pressure difference term and to be proportional to the square of the total coolant mass flow, e.g., at BC entrance

ΔPFADD= - fADDzBT(GF|GF|2ρdH)BE= (1 -εDPZ)ΔPDBT and = 0 (i.e., εDPZ=1) ifΔPDBT,0> 0E74

with the input coefficient εDPZ = εDPZI governing from outside which of them should prevail.

From the now known steady state total ‘additional’ term ΔPDBT,0 the corresponding additive friction factor fADD,0 follows then directly from the equation above, from eq.(73) then the multiplicative one

fFMP,0 =1+ εDPZΔPDBT,0ΔPFBT,0(settingεDPZ=1 ifΔPDBT,0> 0)E75

For the steady state description of a 3D compartment with (identical) parallel channels in a first step a factor should be derived for a ‘mean’ channel. Then this factor can be assumed to be valid for all the other (identical) parallel channels too (See Jewer et al., 2005).

For the transient case there arises the question how the validity of both friction factors could be expanded to this situation too. In the here presented approach this is done by assuming that these factors should remain time-independent, i.e., that fADD = fADD,0 and fFMP = fFMP,0. This allows finally also to determine the wanted nodal pressure decrease term ΔPDBT of eq.(69) for the transient case.

By adding the resulting nodal BC pressure difference terms to the (time-varying) system pressure PSYS(t), given from outside as boundary condition with respect to a certain position (in- or outside of the BC), then finally also the absolute nodal pressure profile PBk along the BC can be established (This term is needed at the begin of the next time step within the constitutive equations).

3.11. BC entrance mass flow (‘Open and closed channel concept’)

As to be seen from the sections above in order to be able to calculate the characteristic nodal and total single- and two-phase parameters along a BC the BC entrance mass flow must be known. This term is for the steady state case given by input (GBE = GBEIN), together with the two pressure entrance and outlet values (PBE = PBEIN, PBA = PBAIN). In the transient situation it can, however, be expected that for the normal case of an ‘open’ channel besides the entrance mass flow only one of these two pressure terms is known by input, either at BC entrance or outlet. The missing one follows then from the calculation of the pressure decrease part.

Such a procedure can not be applied without problems if the channels are part of a complex set of closed loops. Loops consisting of more than one coolant channel (and not driven by an outside source, such as a pump). Then the mass flow values (and especially the entrance term of at least one of the channels) has to be adjusted with respect to the fact that the sum of the entire pressure decrease terms along such a closed circuit must be zero. Usually, in the common thermal-hydraulic codes (see for example the ‘separate-phase’ approaches) this problem is handled by solving the three (or more) fundamental equations for the entire complex system simultaneously, a procedure which affords very often immense computational times and costs. In the here applied module (based on a separate treatment of momentum from mass and energy balance) a more elegant method could be found by introducing an additional aspect into the theory of CCM. It allows, different to other approaches, to take care of this situation by solving this problem by means of a ‘closed channel concept’ (in contrast to the usual ‘open channel’ method).

Choosing for this purpose a characteristic ‘closed’ channel within such a complex loop it can be expected that its pressure difference term ΔPBT = ΔPBA - ΔPBE over this channel is fix (= negative sum of all the other decrease terms of the remaining channels which can be calculated by the usual methods). Thus also the outlet and entrance BC pressure values (PBAIN, ΔPBEIN) are now available as inputs to CCM. Since, according to eq.(69), the term ∆PPBT is known too, then from eq.(68) a ‘closed channel concept criterion’ can be formulated

zBTddtGFBMT = ΔPGBT= ΔPPBT- ΔPBTIN     (at ‘closed channel’ conditions)E76

demanding that in order to fulfil this criterion the total mass flow along a BC (and thus also at its entrance) must be adapted in such a way that the above criterion which is essential for the ‘closed channel’ concept remains valid at each time step, i.e., that the actual time derivative of the mean mass flux GFBMT averaged over the channel must, as recommended by eq.(71), agree in a satisfactory manner with the required one from the equation above.

There exist, obviously, different methods how to deal with this complicated problem. One of them could be to determine the entrance mass flow GBE by changing this value in a recursive way until the resulting term ddtGFBMT agrees with the criterion above.

Another possibility is to find a relation between the time-derivatives of the mean and of special local mass flux values (e.g., at BC entrance), i.e., to establish for example a relation between the terms ddtGFBMT and ddtGFBE. Then the wanted mass flow time-derivative at BC entrance can be determined directly from eq.(76). One practicable method to establish such a relation could follow if considering that a change of the mass flux is propagating along the channel so fast that the time derivative of its mean values could be set (in a first step) almost equal to the time-derivative at its entrance value. This term can, eventually, be provided with a form factor which can be adapted by an adequate recursion procedure until the condition of eq.(75) is fulfilled. The so won entrance mass flow is then governing the mass flow behaviour of the entire loop.

This last and very provisional method has been tested by applying the UTSG-3 code (Hoeld, 2011b) for the simulation of the natural-circulation behaviour of the secondary steam generator loop. Similar considerations have been undertaken for a 3D case where the automatic mass flow distribution into different entrances of a set of parallel channels had to be found (See e.g. Hoeld, 2004a and Jewer et al., 2005). The experience of such calculations should help to decide which of the different possible procedures should finally be given preference.

The ‘open/closed channel concept’ makes sure that measures with regard to the entire closed loop do not need to be taken into account simultaneously but (for each channel) separately. Its application can be restricted to only one ‘characteristic’ channel of a sequence of channels of a complex loop. This additional tool of CCM can in such cases help to handle the variety of closed loops within a complex physical system in a very comfortable way.


4. Code package CCM

Starting from the above presented ‘drift-flux based mixture-fluid theory’ a (1D) thermal-hydraulic coolant channel module, named CCM, could be established. A module which was derived with the intention to provide the authors of different and sometimes very complex thermal-hydraulic codes with a general and easily applicable tool needed for the simulation of the steady state and transient behaviours of the most important single- and two-phase parameters along any type of heated or cooled coolant channel.

For more details about the construction and application of the module CCM see (Hoeld, 2005, 2007a and b and 2011b).


5. Verification and Validation (V & V) procedures

During the course of development of the different versions of the code combination UTSG-3/CCM the module has gone through an appropriate verification and validation (V&V) procedure (with continuous feedbacks being considered in the continual formulation of the theoretical model).

CCM is (similar as done in the separate-phase models) constructed with the objective to be used only as an element within an overall code. Hence, further V&V steps could be performed only in an indirect way, i.e. in combination with such overall codes. This has been done in a very successful way by means of the U-tube steam generator code UTSG-3. Thereby the module CCM could profit from the experiences been gained in decades of years work with the construction of an effective non-linear one-dimensional theoretical model and, based on it, corresponding digital code UTSG-2 for vertical, natural-circulation U-tube steam generators (Hoeld, 1978, 1988b, 1990a and 1990b), (Bencik et al., 1991) and now also the new advanced code version UTSG-3 (Hoeld 2002b, 2004a).

The good agreement of the test calculations with similar calculations of earlier versions applied to the same transient cases demonstrates that despite of the continuous improvements of the code UTSG and the incorporation of CCM into UTSG-3 the newest and advanced version has still preserved its validity.

A more detailed description over these general V&V measures demonstrated on one characteristic test case can be found in (Hoeld, 2011b).


6. Conclusions

The universally applicable coolant channel module CCM allows describing the thermal-hydraulic situation of fluids flowing along up-, horizontal or downwards channels with fluids changing between sub-cooled, saturated and superheated conditions. It must be recognized that CCM represents a complete system in its own right, which requires only BC-related, and thus easily available input values (geometry data, initial and boundary conditions, resulting parameters from integration). The partitioning of a basic channel into SC-s is done automatically within the module, requiring no special actions on the part of the user. At the end of a time-step the characteristic parameters of all SC-s are transferred to the corresponding BC positions, thus yielding the final set of ODE-s together with the parameters following from the constitutive equations of CCM.

In contrast to the currently very dominant separate-phase models, the existing theoretical inconsistencies in describing a two-phase fluid flowing along a coolant channel if changing between single-phase and two-phase conditions and vice versa can be circumvented in a very elegant way in the ‘separate-region’ mixture-fluid approach presented here. A very unique technique has been established built on the concept of subdividing a basic channel (BC) into different sub-channels (SC-s), thus yielding exact solutions of the basic drift-flux supported conservation equations. This type of approach shows, as discussed in (Hoeld, 2004b), distinct advantages vs. ‘separate phase’ codes, especially if taking into account

  • the quality of the fundamental equations (basic conservation equations following directly from physical laws supported by experimentally based constitutive equations vs. split ‘field’ equations with artificial closure terms),

  • the special solution methods due to the detailed interpolation procedure from PAX allowing to calculate the exact movement of boiling boundaries and mixture (or dry-out) levels (different to the ‘donor-cell averaging’ methods yielding mostly only ‘condensed’ levels),

  • the easy replacement of new and improved correlations within the different packages without having to change the basic equations of the theory (for example the complicated exchange terms of a ‘separate-phase’ approach),

  • the possibility to take advantage of the ‘closed-channel concept’ (needed for example for thermal- hydraulic 3D considerations) allowing thus to decouple a characteristic (‘closed’) channel from other parts of a complex system of loops,

  • the speed of the computation,

  • the derivation of the theory in close connection with the establishment of the code by taking advantage of feedbacks coming from both sides,

  • the considerable effort that has been made in verifying and checking CCM (besides an extensive V & V procedure), with respect to the applicability and adjustment and also for very extreme situations,

  • its easy applicability,

  • the maturity of the module, which is continuously enhanced by new application cases,

  • taking advantage of the fact that most of the development work for the coolant channel thermal-hydraulics has already been shifted to the establishment of the here presented module (including the special provisions for extreme situations such as stagnant flow, zero power or zero sub-cooling, test calculations for the verification and validation of the code etc.).

The existence of the resulting widely verified and validated module CCM represents an important basic element for the construction of a variety of other comprehensive thermal-hydraulic models and codes as well. Such models and modules can be needed for the simulation of the steady state and transient behaviour of different types of steam generators, of 3D thermal-hydraulic compartments consisting of a number of parallel channels (reactor cores, VVER steam generators etc.). It shows special advantages in view of the determination of the mass flow distribution into different coolant channels after non-symmetric perturbations see (Hoeld, 2004a) or (Jewer et al., 2005), a problem which is far from being solved in many of the newest 3D studies.

The introduction of varying cross sections along the z axis allows to take care also of thermal non-equilibrium situations by simulating the two separate phases by two with each other interacting basic channels (for example if injecting sub-cooled water rays into a steam dome).

The resulting equations for different channels appearing in a complex physical system can be combined with other sets of algebraic equations and ODE-s coming from additional parts of such a complex model (heat transfer or nuclear kinetics considerations, top plenum, main steam system and downcomer of a steam generator etc.). The final overall set of ODE-s can then be solved by applying an appropriate time-integration routine (Hoeld, 2004a, 2005, 2007a and b, 2011b).

The enormous efforts already made in the verification and validation of the codes UTSG-3, its application in a number of transient calculations at very extreme situations (fast opening of safety valves, dry out of the total channel with SC-s disappearing or created anew) brings the code and thus also CCM to a very mature and (what is important) easily applicable state. However, there is not yet enough experience to judge how the potential of the mixture-fluid models and especially of CCM can be expanded to other extreme cases (e.g., water and steam hammer). Is it justified to prefer separate-phase models versus the drift-flux based (and thus non-homogeneous) mixture fluid models? This depends, among other criteria, also on the quality of the special models and their exact derivation. Considering the arguments presented above it can, however, be stated that in general the here presented module can be judged as a very satisfactory approach.


7. Nomenclature

ABkBC cross sectional area (at BC node boundary k)
ANn, AMnSC cross sectional area (at SC node boundary n, mean value)
ATWMnSurface area of a (single) tube wall along a node n
Dimensionless constant
Phase distribution parameter
dHYmHydraulic diameter
f(z,t), fNn, fMn-General and nodal (boundary and mean) solution functions
fLIMCA-Upper or lower limit of the approx. function f(z,t)
fADD0, fFMP0-Additive and multiplicative friction coefficients
G = GW + GS
GF =GA= v ρ
Mass flow

Mass flux
h, hP, cP = hT
Jkg,m3kg,Jm3kgSpecific enthalpy and its partial derivatives with
respect to pressure and temperature (= specific heat)
hSW = h// – h/, h// ,h/JkgLatent heat, saturation steam and water enthalpy
KEYBC-Characteristic key number of channel BC
LFTYPE= 0, 1 or 2-SC with saturated water/steam mixture, sub-cooled water or superheated steam
LFTBE (=LFTYPE of 1-st NSC)-LFTYPE of 1-st SC within BC
NBT-Total number of BC nodes
NBCA = NCT+NBCE-1, NBCE-BC node numbers containing SC outlet or entrance
NCT=NBCA-NBCE+1-Total number of SC nodes
NSC = NSCE, NSCA-Characteristic number of SC, setting NSCE= 1, 2 or 3
if LFTYPE= 0, 1 or 2 and NSCA=NSCE+NSCT-1
P, ΔPT = PA- PEPa=Jm3=kgms2Pressure and pressure difference (in flow direction)
QBT, QBMkWTotal and nodal power into BC node k
qBMk =UBkABMkqLBMkWm3Mean nodal BC power density into the fluid
(= volumetric heat transfer rate)
qFBk=ABkqBkUBk=qLBkUTWBkWm2Heat flux from (heated) wall to fluid at BC node k
qLBk=QBkΔzBk=ABkqBkWmLinear power along BC node k
qMn =QMnVMn
Wm3Mean nodal SC power density into fluid
(= volumetric0 heat transfer rate)
T, tC, sTemperature, time
UTWmPerimeter of a heated (single) tube wall
VMn=12(ANn+ANn-1)ΔzNnMean nodal SC volume
vS =GFSαρS, vW =GFW(1α)ρWmsSteam and water velocity
X=GSGor =hh/hSW-Steam quality (2- resp. extended to 1-phase flow)
z, ΔzNn=zNn-zNn-1mLocal position, SC node length (zNn-1=zCE at n=0)
zBA-zBE=zBT,zCA-zCE=zCTmBC and SC outlet and entrance positions, total length
zBB, zML, zSPCmBoiling boundary, mixture level resp. supercritical boundary within a BC
α-Void fraction
αTWkWm2CHeat transfer coefficient along a BC wall surface
Δ-Nodal differences
εDPZ-Coefficient controlling the additional friction part
εQTW-Correction factor with respect to QNOM,0
εTWmAbs. roughness of tube wall (εTW/dHW = relative value)
Φ2PF2-Two-phase multiplier
ΦZG-Angle between upwards and flow direction
ρ, ρ P, ρ Tkgm3,kgJ,kgm3CDensity and their partial derivatives with respect to
(system) pressure and temperature
ΘsTime constant
-Partial derivative


0, 0 (=E or BE)Steady state or entrance to SC or BC (n or k =0)
A, E, T (=AE)Outlet, entrance, total (i.e. from outlet to entrance)
B, SBasic channel or sub-channel (=channel region)
A, S, F, D, X
(P=A+S+F+D+X) and
Acceleration, static head, direct and additional
Friction, external pressure differences
(in connection with ΔP) and pressure differences
due to changes in mass flux
Mn, BMkMean values over SC or BC nodes
Nn, BkSC or BC node boundaries (n = 0 or k = 0: Entrance)
S, WSteam, water
P, TDerivative at constant pressure or temperature
TWTube wall surface


/,//Saturated water or steam
P, TPartial derivatives with respect to P or T
(GS), (α), z, sPartial derivatives with respect to
GS, .α. or z (=gradient), slope


Well-known thermal-hydraulic codes (on the basis
a separate-phase approach)
BC, SCBasic(=coolant) channel subdivided
into subchannels (=regions of different flow types)
BWR, PWRBoiling and pressurized reactors
CCFCounter-current flow
CCMCoolant channel model and module (on the basis
a separate-region approach)
GCSMGeneral Control Simulation Module (Part of
GRSGesellschaft für Anlagen- und Reaktorsicherheit
HETRACHeat transfer coefficients package
HTCHeat transfer coefficients
MDSDrift flux code package
MPPWS, MPPETAThermodynamic and transport property package
NPPNuclear power plant
PAX (PAXDRI)Quadratic polygon approximation procedure
(with driver code)
PDE, ODEPartial and ordinary differential equation
UTSGU-tube steam generator code
V & VVerification and validation procedure


  1. 1. AustregesiloHet al2003ATHLET Mod 2.0, Cycle A, Models and methods. GRS-P-1/Vol.4.
  2. 2. BencikVHoeldA1991Experiences in the validation of large advanced modular system codes.1991 SCS Simulation MultiConf., April 1-5New Orleans
  3. 3. BencikVHoeldA1993Steam collector and main steam system in a multi-loop PWR NPP representation.Simulation MultiConf., March 29 April 1, Arlington/Washington
  4. 4. BestionD1990The physical closure laws in CATHARE code. Nucl.Eng.& Design124229245Burwell, J. et al. 1989.The thermohydraulic code ATHLET for analysis of PWR and SWR system, Proc.of 4-th Int. Topical Meeting on Nuclear Reactor Thermal-Hydraulics, NURETH-4, Oct.10-13, Karlsruhe
  5. 5. FabicS1996How good are thermal-hydraulics codes for analyses of plant transients. International ENS/HND Conference,79Oct. 1996, Proc.193-201, Opatija/Croatia
  6. 6. HaarLGallagherJ. SKellG. S1988NBS/NRC Wasserdampftafeln. London/Paris/New York 1988
  7. 7. HannaB. N1998CATENA. A thermalhydraulic code for CANDU analysis.Nucl. Eng. Des.1801998113131
  8. 8. HoeldA1978A theoretical model for the calculation of large transients in nuclear natural circulation U-tube steam generators (Digital code UTSG). Nucl. Eng. Des. (1978), 47, 1.
  9. 9. HoeldA1988aHETRACA heat transfer coefficients package.GRS-A-1446, May.
  10. 10. HoeldA1988bCalculation of the time behavior of PWR NPP during a loss of a feedwater ATWS case. NUCSAFE8827October, Proc. 1037-1047, Avignon, France.
  11. 11. HoeldA1990aUTSG-2. A theoretical model describing the transient behaviour of a PWR natural-circulation U-tube steam generator. Nuclear Techn.9098118April.
  12. 12. HoeldA1990bThe code ATHLET as a valuable tool for the analysis of transients with and without SCRAM, SCS Eastern MultiConf,2326April, Nashville, Tennessee.
  13. 13. HoeldA1996Thermodynamisches Stoffwertpaket MPP für Wasser und Dampf. GRS, Technische Notiz, TN-HOD-1/96, Mai.
  14. 14. HoeldA1999An Advanced Natural-Circulation U-tube Steam Generator Model and Code UTSG-3, EUROTHERM-SEMINAR63September68Genoa, Italy.
  15. 15. HoeldA2000The module CCM for the simulation of the thermal-hydraulic situation within a coolant channel. Intern. NSS/ENS Conference, September1114Bled, Slovenia.
  16. 16. HoeldA2001The drift-flux correlation package MDS.9th International Conference on Nuclear Engineering (ICONE-9), 8-12 April, Nice, France.
  17. 17. HoeldA2002aThe consideration of entrainment in the drift-flux correlation package MDS.10th Int. Conference on Nuclear Engineering (ICONE-10), April 14-18, Arlington, USA.
  18. 18. HoeldA2002bA steam generator code on the basis of the general coolant channel module CCM. PHYSOR 2002, October710Seoul, Korea.
  19. 19. HoeldA2004aA theoretical concept for a thermal-hydraulic 3D parallel channel core model.PHYSOR 2004, April2529Chicago, USA.
  20. 20. HoeldA2004bAre separate-phase thermal-hydraulic models better than mixture-fluid approaches. It depends. Rather not. Int. Conference on Nuclear Engineering for New Europe 2004’. Sept.69Portoroz, Slovenia.
  21. 21. HoeldA2005A thermal-hydraulic drift-flux based mixture-fluid model for the description of single- and two-phase flow along a general coolant channel.The11th International Topical Meeting on Nuclear Reactor Thermal-Hydraulics (NURETH-11). Paper 144. Oct. 2-6, 2005, Avignon, France.
  22. 22. HoeldA2007aCoolant channel module CCM. An universally applicable thermal-hydraulic drift-flux based mixture-fluid 1D model and code.Nucl. Eng. and Desg. 237(2007), 1952-1967
  23. 23. HoeldA2007bApplication of the mixture-fluid channel element CCM within the U-tube steam generator code UTSG-3 for transient thermal-hydraulic calculations.15th Int. Conf. on Nuclear Engineering (ICONE-15). Paper 10206. April 22-26, 2007, Nagoya, Japan
  24. 24. HoeldA2011aCoolant channel module CCM. An universally applicable thermal-hydraulic drift-flux based separate-region mixture-fluid model.Intech Open Access Book ‘Steam Generator Systems: Operational Reliability and Efficiency’’, V. Uchanin (Ed.), March 2011.978-9-53307-303-3InTech, Austria, Available as: u-tube-steam-generator- model-and-code-utsg-3-based-on-the-universally- applicabor
  25. 25. HoeldA2011bThe thermal-hydraulic U-tube steam generator model and code UTSG3 (Based on the universally applicable coolant channel module CCM).See Intech Open Access Book referred above.
  26. 26. HoeldAJakubowskaEMiro, J.E & Sonnenburg,H.G.1992A comparative study on slip and drift- flux correlations. GRS-A-1879, Jan. 1992
  27. 27. IshiiM1990Two-Fluid Model for Two-Phase Flow, Multiphase Sci. and Techn. 5.1, Hemisphere Publ. Corp.
  28. 28. IshiiMandMishimaK1980Study of two-fluid model and interfacial area.NUREG/CR-1873 (ANL-80-111), Dec.
  29. 29. JewerSBeeleyP. AThompsonAHoeldA2005Initial version of an integrated thermal-hydraulic and neutron kinetics 3D code X3D.ICONE13, May1620Beijing. See also NED 236(2006), 1533-1546.
  30. 30. LerchlGet al2009ATHLET Mod 2.2, Cycle B, User’s Manual’. GRS-1Rev. 5, July
  31. 31. LilesD. Ret al1988TRAC-PF1/MOD1 Correlations and Models. NUREG/CR-5069, LA-11208-MS, Dec.
  32. 32. MartinelliR. CNelsonD. B1948Prediction of pressure drop during forced-circulation of boiling water. Trans. ASME 70, 695.
  33. 33. MoodyN. L. F1994Friction factors for pipe flow. Trans. ASME, 66, 671 (See also VDI-Wärmeatlas, 7.Auflage, 1994, VDI- Verlag).
  34. 34. SchmidtEGrigullUed.)1982Properties of water and steam in SI-Units. Springer-Verlag.
  35. 35. SonnenburgH. G1989Full-range drift-flux model based on the combination of drift-flux theory with envelope theory, NURETH-4,10031009Oct.10-13, Karlsruhe.
  36. 36. ShultzR. R2003RELAP5/3-D code manual. Volume V: User’s guidelines. INEEL-EXT-98-00084. Rev.2.2, October.
  37. 37. U. S-N. R. C.2001aTRAC-M FORTRAN90 (Version 3.0). Theory Manual,NUREG/CR-6724, April
  38. 38. U. S-N. R. C.2001bRELAP5/MOD3.3. Code Manual.Volume I, NUREG/CR-5535, December
  39. 39. ZuberNFindlayJ. A1965Average volume concentration in two-phase flow systems, J. Heat Transfer 87, 453

Written By

Alois Hoeld

Submitted: March 6th, 2012 Published: February 13th, 2013