- Email: [email protected]
COBRA-IIIP: An improved version of COBRA for full-core light water reactor analysis
Nuclear Engineering and Design 48 (1978) 293-310 © North-Holland Publishing Company
293
COBRA-IIIP: AN IMPROVED VERSION OF COBRA F O R FULL-CORE LIGHT WATER REACTOR ANALYSIS Robert E. MASTERSON * Batelle Pacific Northwest Laboratory, Richland WA 99352, USA and Lothar WOLF Department of Nuclear Engineering, MIT, Cambridge MA 02139, USA Received 10 February 1978
A new type of pressure-crossflow solution method is developed for solving the fuid conservation equations used by the COBRA-IIIC code. This method is extremely efficient from a numerical point of view, and offers sizeable reductions in running time, relative to conventional numerical schemes for predicting the fluid flow and enthalpy distributions in nuclear reactor fuel rod bundles and reactor cores. This numerical method has been integrated into the computational framework of the COBRA-IIIC code, and an improved version, called COBRA-IIIP, has been developed. This code is faster than COBRAIIIC, and is designed specifically for solving larger and more complex problems, such as full-core PWR transients. Several problems are presented in which the results of COBRA-IIIP, and the method on which it is based, are compared to those of COBRA-IIIC. Excellent agreement between the two codes is observed.
1. Introduction In recent years many digital computer programs [ 1 - 5 ] have been written to solve the set o f fluid conservation equations which characterize the steady state and transient thermal hydraulic performance of nuclear systems. Although these computer programs are widely accepted and used for reactor thermal hydraulic analysis, very little effort has been made to increase the efficiency of the numerical methods used by these codes so that extremely large and complex problems can be solved within the limits of existing computational capacity. The purpose of this paper is to present a new numerical method that has been developed for solving the fluid conservation equations used by the COBRA-IIIC code, which is considerably faster and more efficient than previous methods used to analyze operational reactor conditions. Another object o f this paper is to discuss an improved version of the COBRA-IIIC code, called COBRA-IIIP, which has been developed for determing the fluid flow and enthalpy distributions in fuel pin bundles and light water reactors cores during steady state and transient conditions.
2. Mathematical formulation Owing to the widespread acceptance and use o f the COBRA-IIIC code, a considerable amount o f the physical and mathematical formalism employed by the code will also be used as the basis for the numerical method to be developed here. This essentially means that the same basic set of governing equations will be considered and that
* Former address: Department of Nuclear Engineering, MIT, Cambridge, MA 02139, USA.
R.E. Masterson, L. Wolf/COBRA-II1P. improved code Jbr full-core L W R analysis
294
the same type of solution procedure used for the initial value problem will be retained in space and time. The object of the numerical method to be proposed in this paper is to solve the fluid conservation equations as efficiently as possible so that much larger and more complex problems can be handled within the limits of existing computational Capacity. 3. Basic assumptions With these points in mind, the primary assumptions that will be made in the following development are as follows. (1) A lumped parameter approach is valid. (2) Sonic velocity propagation is ignored. (3) The diversion crossflow is usually smaller than the axial mass flow rate. (4) Viscous dissipation is neglected. (5) The liquid and vapor phases during two-phase flow are in thermodynamic equilibrium. (6) Electromagnetic body forces are ignored. (7) A homogeneous flow model can be used to describe the bulk boiling of the coolant. (8) Flow reversals, recirculating flows, and coolant expulsions are not considered.
4. Basic conservation equations Byadopting the control volume approach described in ref. [5], the equations for the conservation of mass, energy, and momentum of the fluid in a lattice composed of a number of computational cells can be written as:
con tinuity : OPi
N ~
Omi_
Ai ~--~-+~ x - - - i = l
(a)
wi] ;
energy." 1 Ohi + a h i _ q ~ t?
u i Ot
Dx
mi
N
/=1
N
N
(ti - ti) "~i+ ~
(hi - h i )wij ~ (hi _ h,) Wij . 1=I mi ]=i mij'
(2)
axial momentum." N
1 ~m i
api + _ _ = _
A-~l a-~- - 2ui ~--t- Ox
+
L 2Di
k,o, +
22A~ ~
O \ A Jil /.~X-'i" -- gPi COS 0 -- f t
.
i:
(.,-.j, 1
~
Wii
=
Ai
N
+ ~ (2u i -- u*) wiJ" j=l Ai '
(3)
transverse momentum."
where the summations indicated are to be performed by summing over all the cells in a cross-sectional plane of the lattice connected to each cell represented by the index i. In this way the differential equations for a typical cell in
R.E. Masterson, L. Wolf / COBRA-IIIP." improved code for full-core L WR analysis
295
RECTANGULARLA~I CE HEXAGONALLA'nI CE
Fig. 1. Cell numbering schemes for generating systems of difference equationS.
the hexagonal lattice shown in fig. 1 (say cell 4) can be written explicitly by simply summing from] = 1-3, and the differential equations for cell 5 in the rectangular lattice shown in the figure can be generated by summing in an analogous manner from ] = 1-4. The equations for other layouts of cells, such as those in round bundles, can be written by simply summing over all the cells in the lattice that interact with one another in such a way that the interchange of mass, energy, and momentum is allowed to occur. As in COBRA-IIIC, this set of equations is dosed by defining an additional equation for the physical state of the coolant of the form
m : p(h~, e*)
(5)
and empirical correlations for turbulent mixing coefficients, friction factors, and heat transfer coefficients. In this equation it is assumed that the pressure throughout the flow field exhibits only small variations about a known mean or system reference pressure P*. This assumption has been used on many occasions by other authors [6,7], and has found applications in the study .of both single phase and two-phase flow [8].
5. Matrix form of the conservation equations For a lattice composed of a large number of computational cells, the set of conservation equations that has been discussed here becomes unwieldy; consequently, a more compact matrix notation is desirable. In ref. 5 a vector form of the conservation equations is presented in which a transformation matrix [S] and its transpose [S] T is used to reduce these equations to a more tractable form. The elements of the matrix [S] are defined through two column vectors and where i and ] represent the indices of adjacent computational cells and K is used to represent the index of the common boundary connecting them. As in COBRA-IIIC, the entries in [S] and [S] T may only have values of 0, 1, and - 1 ; thus, the crossflow is arbitrarily considered to be positive when the dominant direction of the flow is from cell i to cell ], where i is less than ]. Using the definitions of [S] and [S] T presented in ref. 5, the foregoing conservation equations can also be written
i(K) ](K),
as:
continuity; [AI ~
+
energy: ~,
[~_/=k~j{Q}
_[ll[s]T[Ah]{w'}_[1][s]T[At]{c}+[1 ] [[h][S] T ' -
[SIT[h'l] {w};
(7)
R.E. Masterson, L. Wolf/ COBRA-IIIP: iwpoved
296 Table
code for full-core I, WR analJ’sis
I
Symbol Matrices:
Dimension
Meaning
Special properties
IAjl
NP X NP
diagonal
lcjl
NF X NT:
fcjl
NF X NF
[Djl
NF X NF
Pjl
NP x NP
I+1
NF x NF
Pfl
NF X NF
[Qjl
NP x NP
PI
NP x NF
matrix of cross sectional flow areas in the axial direction matrix of thermal conductivities between adjacent computational cells matrix of transverse friction coefficients between adjacent computational cells matrix used to compute the cross flow distribution from the transverse momentum equation matrix of enthalpies from each computational cell matrix of enthalpy differences between adjacent computational cells matrix of the enthalpies carried by the diversion cross flow between adjacent computational cells matrix of volumetric heat generation rates interface or gap connection matrix used as a lateral differencing operator
FIT
NP x NP
interface or gap connection matrix used as a lateral summation operator
NF x NF
matrix of lateral temperature differences between adjacent computation cells matrix containing the axial velocity of the fluid from each computational cell matrix containing the axial velocity of the fluid carried by the diversion crossflow between adjacent computational cells
NP x NP
NF x NF
NP X 1
{bj)
NP x 1
{hj)
NP x 1
vector containing the axial pressure gradients from each channel in the lattice source vector used by the MAT method vector containing the enthalpies from each channel in the lattice
diagonal
diagonal
diagonal diagonal
diagonal contains only entries having l,O, and -1, not diagonal contains only entries having the values of 1, 0, and -1, not diagonal diagonal
diagonal
diagonal
R,E. Masterson, L. Wolf/COBRA.IIIP: improved code for full-core L WR analysis
297
Table 1 (continued) Symbol Matrices:
Dimension
Meaning
Special properties
(mi}
NP X 1
(Pi)
NP X 1
{Pi-1}
NP X 1
{p]}
NP × 1
{u]}
NP X 1
{u]*}
NF X 1
(wj}
NF X 1
vector containing the axial mass flow rates from each channel in the lattice vector containing the pressure field at axial level] vector containing the pressure field at axial level] - 1 vector containing the density of the fluid in each channel in the lattice vector containing the axial velocity of fluid in each channel in the lattice vector containing the axial velocity of the fluid carried by the diversion crossflow vector containing the cross flow between adjacent computational cells
-
-
Miscellaneous constants." s/l
-
x Ax At 0
--
transverse momentum factor axial elevation axial mesh spacing time step size weighting function used by the MAT method
NF = number of boundaries for crossflow between adjacent cells. NP = number of nodes at which the pressure is to be found.
axial momentum: -
2u
/~/-
{a'} +
[[2ul IS] T - [slT[u*I]
(w},
(8)
w h e r e the c o m p o n e n t s o f the pressure d r o p due to frictional forces, gravitational forces, and the e f f e c t s o f turbulent m i x i n g are given b y +--
-l\--A] ~2-D 2&,¢
+A
3x
]
pgcos4~
-ft
[ s ] T [ A u ] (w'}
(9)
and transverse i n t e r a c t i o n s are t a k e n i n t o a c c o u n t b y m e a n s o f the m a t r i x f o r m o f the transverse m o m e n t u m equation d e v e l o p e d for the COBRA-IIIC code:
transverse momentum:
298
R.E. Masterson, L. Wolf/COBRA-IIIP: improved code for full-core L WR analysis
For convenience the same nomenclature and notation used in ref. [5] are also used here. Consequently, although most of this nomenclature is summarized in table 1, ref. [5] should be consulted for more detailed information about the definition and derivation of the individual terms.
6. Conversion of the conservation equations into difference equations The preceding conservation equations must be converted into a system of difference equations before they can also be solved by numerical means. Taking a first order backward differencing scheme for the spatial and temporal derivatives to provide numerical stability allows the equations for the conservation of mass and axial momentum to be written as
continuity:
axial momentum:
For reasons of computational efficiency, the energy equation is approximated using a differencing scheme which is spatially explicit and temporally implicit:
energy: [ l lth/--hJi*[h/-h/-1}
= {mi} -1 { { Q i - t / 2 } - [ S ] [ A h j - t ] ( w ; - t } - [ s ] T t A t j - t ] [ c j - t ] }
+ (mi)--l[[hj_l] [S]x - [S]T[h~_I]I .
(13)
On the other hand, the energy equation can be made fully implicit, and the coolant enthalpy distribution can be determined by a solution scheme similar to that given in ref. [9]. In either case, it is convenient to visualize each control volume of fluid as being bounded by two adjacent control volumes in the axial direction; thus, the amount of mass, energy, and momentum that can be interchanged also depends on the number of computational cells connected to one another in each cross-sectional plane of the core.
7. The transverse momentum equation The transverse momentum equation to be considered in this analysis is similar to the transverse momentum equation used by the COBRA-IIIC code [5]. However, it differs from previous versions of the equation because it is based on the concept of splitting the pressure field (P) and the transverse friction terms [C] (w} which drive the crossflow distribution into a sum of spatial implicit and explicit parts. The algorithm which governs the form of the terms is
(P) = 0 (Pi} + (1 - O) (Pi_ l) ,
(14a)
[C] {w} = 0 [Cj] {wi} + (1 - O)[Ci_ 1] (wi- ,},
(14b)
where 0 is a weighting function having an arbitrary value between 0 and 1.0. As a result of this algorithm the trans-
R.E. Masterson, L. Wolf/COBRA-IIIP: improved code for full-core L WR analysis
299
verse momentum equation, eq, (10), can be written in finite difference form as
---XT--] +
ax
+
(O[©l(wj)+(l-O)[Q_ll(wj_~}}= [sl(o{Pj)+(l-O)(pj_l)}, (15)
where 0 ~<0 ~< 1.0. Again, it is worth mentioning that there are no numerical restrictions on the time step size, At, since eq. (15) and all the other equations which govern the behavior of the fluid are temporally implicit.
8. Solution of the multidimensional COBRA-IIIC conservation equations with the MAT method In order to solve the conservation equations presented in this paper, it is necessary to develop a numerical method for computing the fluid flow and enthalpy fields in an operational reactor core. Since the conservation equations have already been cast into a form suitable for digital computer analysis, the development to follow will consist of identifying techniques for solving these equations which can bring about sizeable reductions in computing time, relative to conventional numerical methods. In ref. [10] a pressure-velocity method is presented in which enormous amounts of computer time can be saved by developing an iteration scheme in which the pressure fields in the lattice are found as a function of the axial mass flow rate. Unfortunately, problems with numerical stability are encountered with this scheme when it is used for a flow blockage analysis because the frictional and gravitational terms in the axial momentum equation, eq. (8), are approximated by spatially explicit relationships of the form {a') = - { K j } ( m / _ l } 2 - {fl'-l),
(16)
where {Kj} is the coefficient of the flow-squared terms and all the other terms (including the gravitational terms) are lumped into (t)_ l }. In this paper a modified version of this numerical method will be presented in which this shortcoming is overcome by replacing eq. (16) by a differencing scheme that is both spatially and temporally implicit: {a') = -{Kj}(mj) 2 - {fj} .
(17)
The numerical method which results from this approximation is a new type of pressure-velocity method called the MAT method (Modified and Advanced Theta method). In what follows it will be shown that the MAT method combines the computational efficiency of the method presented in ref. [10] with the stability of the crossflow solution scheme used by the COBRA-IIIC code. Finally, it will be demonstrated that this combination of computational characteristics makes it possible to save a remarkable amount of computer time when a method of this type is used for solving large problems in which there are either single-phase or two-phase flows.
9. The MAT method: derivation of the difference equations for the pressure field If the continuity equation, eq. (11), is written as (18) it is possible to express the mass flow rate at axial level ] in terms of the mass flow rate at axial level ] - 1 by means of a relationship of the form
(m/) = { m / _ l ) + (Am),
(19)
300
R.k'. Masterson, L. Wolf~ COBRA-II1P: improved code jor fidl-core L WR ana(vsis
where
{~xm} = --~xx [S] t ( w j ) - ~ ¢ [A/] l~J ~?P~J} k
(20)
l
is the change in the axial mass flow rate between tile top and the bottom of each plane of computational cells. Using this relationship, the terms in the axial momentum equation which account for the pressure drop due to frictional and gravitational effects can be written as (2 l a)
(a;} = - {k/}(m~} - ~t).} = -(ki}(mi-1
(21b)
+ A m } 2 -- (J)}
= - { k i } ( r n ~ _ , + 2rnj_ ,Av'n + ~ n 2}
~.},
(21c)
where (Am} is generally smaller than {rni_ 1}. In many cases the stability of this difference approximation can be improved by applying the continuity equation {~n} = {m/} - {m/_ l}
(22)
to write the flow-squared terms as a linear function of { ~ z ) : (rn/} 2 = (my_ ,}2 + {mj + rni_ 1}{&m) .
(23)
The terms in eq. (2l) which account for the effects of frictional and gravitational forces on the axial pressure gradient are then given by {a'}
= -{kj}{mj_ 1} 2 - {kj}{rttj + rrtj_ 1}{~11} -- {f/-} .
(24)
Applying the same procedure to the temporal acceleration term in the axial momentum allows this term to be written as [Aft -1 At {mj
-
~i} -
[Aj]-I At
{mj -- m i _ 1 + m / _ 1 - in/} - [ Aat j ] - 1 { m i - mi-1} +
[Aft-'
at - { m / - 1 - mfl '
(2s)
and since the continuity equation requires that {ZXrn) = {mj} - {my_ l ) ,
(26)
the temporal acceleration term can also be expressed as a linear function of {Am): [ AAt I ] - I { m i _ ~ i } = [A/]_ 1 {~mt ~-}+[A/]
-1
{m J - l~- Ftn J
}.
(27)
Substituting these linear relationships for the frictional and temporal acceleration terms into the axial momentum equation, eq. (12), allows the conservation of momentum in the axial direction to be written as ax
at
J
[Ail
at
+ [Aj] -I [[2uj] IS] T - [s]T[uTII {w/).
- {ki}{mj_ 1} = - 0')} -
+ {/9}{,ni + m i -
(28)
The terms containing {am} can then be eliminated from this system of equations by combining eq. (28) with eq. (20). The result of this operation is a system of difference equations which simultaneously satisfies the conservation of mass and momentum at each axial level of the core at each instant of time:
R.E. Masterson, L. Wolf/COBRA-IIIP." improved code for full-core L WR analysis
+,A,,,
301
,s,w,u.,1
The effects of lateral m o m e n t u m transfer can be incorporated into eq. (29) by combining it with the transverse m o m e n t u m equation, eq. (10), to eliminate the crossflow-dependent terms. Since the transverse m o m e n t u m equation can be written as 1 +--u; + s 0[Ci ] {w/} = At Ax
+ ui-~xi-1 +
(0-1)[Cy_l]{Wi_l}+
[S]{O{P/}+(1-O){Pi_l}}
,(30)
the crossflow distribution {wj} at axial level j can be expressed as a function of the pressure fields at levels/" and j - 1 by means of a relationship of the form
{wi} = [O3 -'{0 i} + (~)[Oj] -1 [Sl{O{Pj} + (1 - O)(Pj_ ,}},
[Dj] -1 is the
where
(31)
diagonal matrix defined by
[Oyl-' = [ l ' l - U ] * + LAt 2v¢
(32)
O[Cjl
and
[[atJ[
Ax
J
is a source vector consisting of the other forcing terms in the transverse m o m e n t u m equation. If eq. (31) is then combined with eq. (29) to eliminate the crossflow-dependent terms, it is possible to derive the following system of difference equations for the pressure field as a function of the axial mass flow rate:
{Pj - PI-1 } = ~
+ ~[A/I-'
[Ail-'{ff9 - m j _ l } -- ~ {kj } {my2_~} -- A,v. ~.} + - ~ {pj -- ~j} 2uj + - - ~ + ,5x kjA j(rnj + m j _ l )
2 u / + - - + Ax[Aj]kj(mj + mj_0] [S] T - [s]T[u; At
[Dyl-'{Oy}+
[Dy]-I[S]{0{Py}
x
+ (1
- o) {Py_,}}}.
(34)
This system of equations for the pressure field can also be written in matrix notation as
[I+ (1 - O)Myl{Py_ 1} = [ I - OMyl{Py} + {by},
(35)
where I is the identity matrix,
[My] = A x ( / )
[[Aj]-'
[By][sIT[D/]-'
iS] -
[A/l-'
[S] T[U]*] [D/l-' [Sl]
(36)
is a matrix used for determining the pressure field,
{by} =--~ [ A j ] - ' { ~ y - m y _ ,} - 2~x{kj}{ml 2._,} - Ax ~.} +Ax at {py - -~y} 2uy + + Ax [A/] -1 [By] [ s ] T [ D j ] - ' { O j }
+ axkjA1{m j + my_,)
- A x [ A i ] - 1 [S] T[u~] [1)/] -l{Oj}
(37)
is a source vector consisting of terms that contribute to the axial pressure gradient, and
[By] = [2uy] + ~ - + Ax[A/]ky(m/+ my_,)
(38)
R.E. Masterson, L. Wolf/COBRA-IIIP." improved code for full-core L WR analysis
302
is a matrix which is used in the construction of both the source vector and [Mj]. Although the source vector that drives the pressure field consists of only one entry from each computational cell, it is worth noting that this vector contains all the physical forces that can act upon an elemental volume of fluid. These forces can be due to frictional effects, gravitational effects, and the spatial and temporal acceleration of the flow. Finally, although the axial mass flow rates used in eq. (35) are not actually known, it is possible to develop simple procedures for estimating the values of these flow rates and updating them through iteration. Procedures of this kind are now used by the COBRA-IIIC and COBRA-IV-I codes [5,9]. 10. Discussion of the difference equations governing the pressure field It is interesting to observe that the use of the weighting function 0 in the transverse momentum equation enables a system of difference equations to be derived in which the crossflow distribution can be driven by any combination of the pressure fields that exist at the top and the bottom of a plane of computational cells. This can be seen more clearly by writing eq. (35) as
{P/- 1} = [I + (1 --0 )1~]]-1 [I -- OMj] {Pj} + [I + (1 -O)Mj] -l{b/},
(39)
and by examining the structure of the matrix equations governing the pressure field {Pi_ 1} for various values of 0. Setting 0 = 0 allows eq. (39) to be written as: 0=0:
(P]- 1) = [I +/111/]--1 [I - Mi] (Pi) + [I + Mj] -l(b/) .
(40)
If 0 is set to ½ and 1.0, respectively, eq. (39) then becomes 0=1.
2'
{Pj_ ,} = [I + ~,Mj]-1 [I - ~ j ] {ej} + [I + ~,Mj]-l{bj}
;
(41)
0 = 1.0:
(P/- 1} = [I - Mj] ~pj} + {bj} .
(42)
These equations demonstrate that a direct relationship exists between the form of the difference equations for the pressure field and the forces used to drive the crossflow distribution between computational cells in the transverse plane. For example, setting 0 = 1 allows the crossflow distribution to be driven by the pressure field that exists at the top of each plane of computational cells. From eq. (42) it can be seen that this results in a differencing scheme that is spatially explicit but temporally implicit. Setting 0 = 0 means that the crossflow distribution is to be driven by the pressure field that exists at the bottom of each plane of computational cells. In this case the difference equations governing the behavior of the fluid are both spatially and temporally implicit. Finally, when 0 is set equal to ½, eq. (41) shows that it is possible to generate a system of equations where the crossflow distribution is driven by the average of the pressure fields that exist at the top and the bottom of each plane of computational cells. In this case the difference equations are temporally implicit but have a spatial component whose structure is analogous to that of the Crank-Nicholson method. Of course, it should be evident that other spatial differencing schemes can be constructed by simply choosing other values for 0.
11. Solution schemes for the pressure distribution If the difference equations governing the computation of the pressure field are written as
(Pj-I} = [I+ (1 - O)Mj]-l[I - OM/]~p/} + [I+ (1 - O)Mj]-Ifbj}
(43)
R.~: Masterson, L. Wolf/COBRA-IIIP: improved code for full-core L WR analysis
303
it is possible to invert the coefficient matrix [K + (1 - 0)
Mj]
(44)
and solve these equations for {Pj_ l) as a function of (Pi). Since the pressure field at axial level j - 1 is upstream of the pressure field at axial level j, this allows disturbances in the flow field to be felt at upstream locations. One of the primary computational advantages of the pressure-oriented solution scheme outlined here is that the coefficient matrix governing the computation of the pressure field is a Stieltjes matrix. In other words, it is possible to show that the coefficient matrix in eq. (43) is a diagonally dominant, irreducible, positive definite matrix with a simple and predictable band structure for a variety of subchannel numbering schemes. Since this enables the pressure distribution to be found quite easily by both direct and iterative solution schemes, the computation of the pressure field requires considerably less effort than in other approaches where the crossflow distribution is found directly as a function of the axial mass flow rate [11].
12. Derivation of the crossflow distribution from the pressure distribution Whereas the numerical methods used to solve the momentum equations in previous versions of COBRA have been based upon the concept of solving directly for the crossflow distribution as a function of the axial mass flow rate [5,9,11 ], the pressure-velocity method discussed in this paper is based upon the alternative numerical approach of first solving for the pressure field as a function of the axial mass flow rate and then using this pressure field to determine the crossflow distribution from the transverse momentum equation. In other words, this approach assumes that substantial transverse pressure gradients can exist in a reactor lattice and that these gradiants are the primary physical mechanism for driving the crossflow distribution. After the inlet flow and enthalpy distributions have been specified, eq. (35) can be solved for {Pj_ 1) and a pressure field can be constructed for driving the crossflow distribution between the cells in each computational plane. This pressure field can be found by blending together the pressure fields from the top and the bottom of each plane of cells to form the composite pressure field given by eq. (14a) as the numerical solution scheme sweeps downstream. The crossflow distribution can then be found from the transverse momentum equation by writing it in the form
(w/} = [DjI-I(oj) + (~)[Ojl-I[s](o(Pj}
+(1 - o ) { e j _
1)),
(45)
where [Dj]-I is the designed matrix defined by eq. (32) and (Oj) is a vector consisting of the other forcing terms in eq. (30). Since the inverse of [Dj] can be found by inspection, this additional step in the solution scheme requires very little computational effort to determine the flow and enthalpy fields once the pressure distribution is known.
13. Solution schemes for the energy equation After the equations for the conservation of mass and momentum have been solved, it is necessary to interface the thermal response of the fluid with the hydraulics through heat transfer coefficients. This can be done by soNing either an explicit or implicit energy equation for the enthalpy of the fluid at axial level j as a function of the enthalpy of the fluid at the previous axial level and the amount of heat produced in each rod in the lattice. An energy balance can be performed using eq. (13) and the axial mass flow rates (mj} found from a previous solution of the momentum equations. Thus, the solution of the energy equation can be looked upon as the 'outer iteration' in this two-step solution procedure. For problems in which there are substantial departures from operational reactor conditions, several additional iterations between the energy and the momentum equations may have to be performed.
304
R.E. Masterson, L. Wo(f / COBRA-IIIP. improved code jbr fidl-core L WR ana(vsis
14. Overall iterative strategies At this point it may be helpful to summarize how the MAT method can be used to find tile flow and enthalpy fields in an operational reactor where there can be strong coupling between the energy and momentum equations. Referring to the flow chart shown in fig. 2, a typical calculation consists of sweeping downstream several times in succession between the inlet and the outlet of each channel in the core. Once the inlet flow and enthalpy distribution have been specified, the enthalpy can be advanced from axial level / - I to axial level / by solving either an implicit or explicit version of eq. (7). For the first axial iteration, the flow rate (mi} for axial level / is set equal to the flow rate from level j - 1; otherwise the value of {rn/} from the previous iteration is used. After eq. (35) is solved for the pressure distribution {Pj_ 1} the pressure fields from adjacent axial levels are blended together using eq. (14a) and the transverse m o m e n t u m equation is solved for {w/}. Since the value of {Pj} in eq. (35) is not initially known, it is set equal to a constant (such as zero to avoid the necessity of computing the source term [1 0m]
ESTABLISHINITIAL CONDITIONS ONh, w, ANDm
SOLVETHEENERGYEQUATION,
EQ. (/)TO GET(hi)
I FORFIRSTAXIAL ITERATION ASSUME: (rni) (rnj.l) •
i
SOLVEEQ. (35)FORTHE PRESSUREFIELD(Pi_lI
I BLENDTOGETHERTHEPRESSURE FIELDSUSINGEQ. ll4t
I I
SOLVEEQ. 145)FORTHE ~ ] CROSSFLOWDI STRI BUTION
I
] FINOA NEWVALUEOFIm ) FROM I THECONTINUITYEQUATON, EQ. i18 ] i
l
l
SWEEPDOWNSTREAMTOTHE NEXTAXIALLEVELAND REPEATTHISLOOPAGAIN
I
DOANOTHERAX,AL
1.0[---WHENTHEITERATIONREACHES ~
ITERATIONWITH NEWVALUESFOR(mjl I
TES'FORCO VERGENCE
l
THEENDOFTHECHANNEL
~, YES PRINTRESULTS CONDITIONSAND CONTINUETHECALCULATION
DOANOTHERTIMESTEP?
Fig. 2. Flow chart for the MAT method.
R.E. Masterson, L. Wolf/COBRA-IIIP: improved code for full-core L WR analysis
305
(Pj} for the first axial iteration) and then updated so that its new values can be used for successive iterations. After the axial mass flow rate {mj} is updated using the continuity equation, eq. (18), the numerical scheme is allowed to sweep downstream to the next axial level, where this entire process is repeated again. When the calculation reaches the exit of each channel, a check is made to see if the flow distribution has converged at every mesh point in the lattice to the tolerance desired. If the solution has converged, the computation stops; otherwise, the entire iterative process is repeated again, starting at the first axial level (e.g. the core inlet) and sweeping downstream to the last axial level (e.g. the core outlet). For all iterations after the first, the most recent values of the axial mass flow rates are used for computing the flow and enthalpy fields. This eliminates any errors that may have been introduced into the solution by estimating a value of (m/) for the first axial iteration. For all practical problems, it has been found that this iterative procedure yields a unique, stable, and physically realistic solution to the fluid conservation equations in a nuclear reactor core.
15. The COBRA-IIIP code In order to take advantage of the computational capabilities of the MAT method, a new version of the COBRAIIIC code, called COBRA-IIIP *, has been written for light water reactor analysis. This version of the code is based upon the same physical formulation as previous versions of COBRA, but due to the use of the MAT method and substantial improvements in computer coding, the new code has the capability to solve much larger and more complex problems than COBRA-IIIC with much greater speed. For example, the code has been used to model an entire PWR fuel pin bundle during a 5-min loss-of-flow transient. It has also considered problems having enough computational cells to represent over 600 fuel assemblies in a reactor core explicitly. To more clearly understand the capabilities of the new code, it may be helpful to summarize the types of problems to which the code can be applied. These problems are generally the same as those that can be solved by COBRAIV-I. However, since the code does use a 'MAC' methodology for solving the momentum equations, it cannot be used for solving problems where there are substantial flow reversals and recirculations. Problems of this type can only be solved by numerical schemes such as those developed for the TRAC code [ 12] and the explicit part of COBRA-IV-I [91.
16. Comparisons of computation time To demonstrate the computational advantages of the MAT method several problems were run in which the computation times of the COBRA-IIIP code was compared to that of COBRA-IIIC. In these problems the number of cells in each plane of a pressurized water reactor core was varied between 16 and 400 and only rectangular mesh grids were used. In all cases the flow fields were converged to a tolerance of 1%; that is, the relative change in the axial mass flow rate, dm/m, at each point in the core was constrained to change by less than 1% between successive axial iterations. Moreover, the momentum equations in COBRA-IIIP were solved by both direct and iterative solution techniques to demonstrate the advantages of different methods of matrix inversion as a function of problem size. The results of the timing run comparison are summarized in table 2. In all cases it was found that an enormous reduction in computation time could be achieved by using the COBRA-IIIP code and the MAT method in place of the crossflow solution technique in COBRA-IIIC. Whereas the crossflow solution time of the COBRA-IIIC code was found to increase as the number of channels cubed, the crossftow solution times of COBRA-IIIP increased as only the 1.5 power of the number of channels for a direct inversion scheme, and as almost a linear function of the number of channels when the method of successive over-relaxation was used. For large problems, such as full* The letter P is used to designate the fact that the code is based on a pressure solution scheme rather than the crossflow solution schemes used by COBRA-IIICand COBRA-IV-I.
306
R. bZ Masterson, L. Wolf / COBRA-IIIP. improved code for full.core L WR analysis
Table 2 A comparison of the crossflow solution times a) of the COBRA-IIIC and COBRA-IIIP codes Code
COBRA-II1C
COBRA-IIIP
COBRA-IIIP
Method of solution
Gauss elimination with no compression of crossflow coefficient matrix Gauss elimination with full compression of pressure coefficient matrix successive overrelaxation with optimized relaxation factor
Number of channels per plane 16
64
128
200
300
400
1.7
128.49 b)
987.36 b)
3243.61 b)
10 250 b)
31 980 b)
0.17
0.76
2.16
4.14
7.95 b)
16.53 b)
0.15
0.44
1.23
1.82
2.65 b)
3.62 b)
Notes: All results are for ten axial levels on an IBM 370/168 computer with the H compiler. Note that the results of these timing runs may vary by -+10% during the course of a day due to changes in the work load on the system. a) In seconds of central processing unit time b) For economic reasons these results are estimated by extrapolation.
scale simulations o f PWR fuel assemblies, the COBRA-IIIP code and the MAT method were found to be over three orders of magnitude faster than COBRA-IIIC. It should be emphasized that these comparisons were made by using exactly the same types o f iteration schemes in both codes. In other words, the flow and enthalpy fields were recomputed over the entire axial height of the core whenever a single point was found where the solutions were not fully converged. For problems where the effects of natural convection are much more important than those o f forced convection, iteration schemes o f this kind are desirable because they allow details o f the downstream flow distribution to be felt at upstream locations. However, for many other practical problems these schemes have the drawback that they consume an enormous amount of computation time because they require the flow fields to be recomputed at every mesh point in the core - whether or not this recomputation is actually needed. In order to overcome this difficulty, an alternative iteration scheme has been developed for COBRA-IIIP in which the effects o f upstream pressure propagation are suppressed by performing additional iterations at only those points in the lattice where they are needed. This iteration scheme consists o f solving eq. (35) for the pressure field and then solving eq. (45) for the crossflow distribution to complete the first radial itejration. If the change in the axial mass flow rate, {Am) = {mj - m / _ 1), between the top and the b o t t o m o f each node in less than some prescribed convergence criterion, the flow field is assumed to converge, and the numerical method is allowed to sweep downstream to the next axial level, where this process is repeated again. If, however, the flow field has not converged, the axial mass flow rates are updated via the continuity equation, eq. (I 8), and a new pressure field is found b y solving eq. (35). A new crossflow distribution is then computed from the transverse m o m e n t u m equation, and the computational cycle is repeated again - with feedback from the energy equation, if necessary - until a converged flow and enthalpy distribution is found. Since the flow and enthalpy fields are only updated at points where additional iterations are needed, it is often possible to reduce the computation time o f COBRA-IIIP b y another factor o f 3 or 4 compared to COBRA-IIIC. In other words, a significant reduction in computing time can be achieved by simply developing a more efficient numerical method and a more efficient iteration scheme for solving the same set o f m o m e n t u m equations that is
R.E. Masterson, L. Wolf / COBRA-IIIP: improved code for full-core L WR analysis
307
used by the COBRA-IIIC code. For extremely large problems, such as full core loss-of-flov~ transients with 600 channels, this saving in computer time may be as great as three and one half orders of magnitude.
17. Results To demonstrate that the predictions of the COBRA-IIIP code agree with those of COBRA-IIIC, a number of sample problems were constructed in which the flow and enthalpy fields were computed by both codes in an operational reactor core. Since the COBRA-IIIC code was not capable of solving problems as large as those which could be solved by COBRA-IIIP, the comparisons were restricted to problems in which reasonable numbers of computational cells could be used. The first sample problem to be considered is shown in fig. 3. In this problem, the departure from nuclate boiling ratio (DNBR) was computed by both codes for the hottest rod in the lattice of the PWR whose operating parameters are shown in table 3. The results of these calculations are summarized in fig. 4. Inspection of these curves reveals that the COBRA-IIIP code gives a departure from nucleate boiling ratio that is almost exactly the same as that given by COBRA-IIIC. Moreover, fig. 4 shows that the Westinghouse (W-3) correlation apparently gives more conservative estimates of the departure from nucleate boiling ratio than the Babcock and Wilcox (B&W-2) correlation over the entire axial length of the core. The predictions of the two correlations can be seen to deviate from one another by as much as 20% in the hottest region of the channel for the operating conditions given in table 3. Fig. 5 shows a cross-sectional view of a second sample problem that was used to compare the predictions of the COBRA-IIIC and COBRA-IIIP code. This problem consists of two identical fuel pin bundles of a PWR with different radial power peaking factors whose size and operational characteristics are described in table 4. In this problem, a grid with a loss coefficient of 2.5 is located midway between the inlet and the outlet of each assembly and both assemblies are subjected to a power transient in which the power level is doubled uniformly in 1 s. To follow the behavior of the assemblies during the transient, it was assumed that each of the assemblies could be homogenized (i.e. represented by a single computational cell) and that an axial mesh spacing of 2 in. could be used. Table 5 shows the coolant temperatures, enthalpies, densities, and axial mass flow rates predicted by both codes at the exit of the hot channel as a function of time. It can be seen that the predictions of both codes are in excellent agreement with one another for this particular case. The crossflow distribution between the assemblies is also plotted at the end of the transient in fig. 6. Inspection of these curves reveals that the magnitudes and the shapes of the crossflow distributions are essentially the same along the entire length of the core. However, some small
Table 3 Input parameters for sample problem 1 Type of reactor - pressurized water Axial power distribution Radial power peaking factors Nominal operating conditions: System outlet pressure: Inlet temperatures (uniform for all channels): Average inlet mass flux (uniform for all channels): Average reactor heat flux: Axial mesh spacing: Channel length:
skewed cosine as shown in fig. 3 2100 psia 541.0°F 2.51 X 106 lbm](h ft 2) 0.21 X 106 Btu/(h ft 2) 6 in. 180 in.
R. l'#. Masterson, L. Wolf~ COBRA-IIIt'. improved code for/idl-core L WR aptalvs'is
308
DISTANCE FROM INLET, cm
0 9.O I
100
200
D0
300
T
----COBRA IIIP
I
8.0~
--COBRA
IIIC
i! ~
x, . ~ . / H O T ASSEMBLY (FOR DETAILS SEE BELOW] ........
B&W-2 ~
_o ~ 6.0
,
.194
..... L .,,__.~'__~._
o.9~
g
:-.__~___~ q---1---'---~--+.._, "--, ,
1.351
7.o i
3.0
w ~3
CORRELATION
1.5,~...~m 11.5,~I a.5~ ~...~l.ssl ~.
:~
"
N
2..0
HOT ROD
1.0
00 ~ 4 9 ~ •
~
G CTORS I
1
1
40
80
120
160
DI STANCE FROM INLET, in.
Fig. 3. A one-eighth full-core symmetry section of a PWR (powers normalized to 1.0). Fig. 4. Axial distribution o f the DNBR predicted by COBRA-IIIC and COBRA-IIIP for the hottest of the lattice shown in fig. 3.
Table 4 Input parameters for sample problem 2 Type of reactor
pressurized water
Axial power distribution Radial power peaking factors Nominal operating conditions: System outlet pressure: inlet enthalpies (uniform for all changes): Average inlet mass flux (uniform for all channels): Average reactor heat flux: Channel length: Axial mesh spacing: Flow area for each cell: Number of rods per cell:
uniform as shown in fig. 5 2 100 psia 538.0 Btu/lb m 2.48 X 106 lbm/(h ft 2) 0.2 × 106 Btu/(h ft 2) 120 in. 2 in. 0.266 ft 2 225
R.E. Masterson, L. Wolf / COBRA-IIIP: improved code for full-core L WR analysis
N
8"
.i./ / / /
"
309
RADIALPOWER . PEAKINGFACTORS DISTANCEFROMINLET,cm 100 200
8.0 7.0
300
COBRA-IIIP COBRA-IIIC
---- -
6.0 "COLD"~
/
ASSEMBLY~ ; E A C T O R O U T L E T /
'I-lOT"
5.0
ASSEMBLY ASS
4.0
3.0
2.0
|,
,, ]
GR,D 1.0
//
L/
/x
,,
GRIDLOCATION.
O.O
/
.
.
ONSETOF . ~ ~
I
I
i
3O
6O
9O
120
DISTANCEFROMINLET, in.
REACTORINLET Fig. 5. A cross-sectional view of two homogenized PWR fuel assemblies.
Fig. 6. Axial crossflow distributions predicted by COBRA-IIIC and COBRA-IIIP for the problem shown in fig. 5.
Table 5 Values of the fluid properties predicted by COBRA-IIIC and COBRA-IIIP at the outlet of the hot channel as a function of time Time (s)
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
Enthalpies(Btu/lbm)
Densities ( l b m / f t 3)
Mass fluxes (Mlbm/ h ft 2)
Temperatures (°F)
COBRAIIIC
COBRAIIIP
COBRAIIIC
COBRAIIIP
COBRAIIIC
COBRAIIIP
COBRAIIIC
COBRAIIIP
627.70 629.08 631.84 635.90 641.21 647.56 654.78 662.67 671.03 679.78 688.73
627.72 629.10 631.86 635.90 641.20 647.55 654.76 662.64 670.99 679.75 688.69
41.65 41.57 41.40 41.16 40.84 40.44 39.99 39.48 38.93 38.34 35.94
41.65 41.57 41.40 41.16 40.84 40.45 39.99 39.49 38.94 38.36 35.96
2.464 2.479 2.490 2.499 2.505 2.509 2.511 2.511 2.510 2.508 2.501
2.462 2.480 2.493 2.501 2.506 2.510 2.513 2.513 2.512 2.510 2.509
607.30 608.24 610.08 612.80 616.29 620.42 625.01 629.96 635.08 640.29 64.2.81
607.33 608.26 610.09 612.79 616.29 620.40 624.98 629.93 635.05 640.26 642.78
310
R.E. Masterson, L. Wolf/COBRA.IIIP. improl,ed code for full.core L WR analysis
differences between the predictions of the codes can be seen at the location of the onset of bulk boiling and in the vicinity of the grid. These differences are due primarily to the differences in the axial iteration schemes discussed in section 16. Nevertheless, in view of the fact that they are considerably smaller than the uncertainties in some of the empirical correlations used by the codes, it is felt that these differences are quite small a price to pay for such an enormous increase in computational efficiency.
18. Conclusion In conclusion, a new numerical method has been developed for solving the set of fluid conservation equations used by the COBRA-IIIC code. This method is considerably faster and more efficient than previous methods that have been used for reactor thermal-hydraulic analysis. The new method of solution has been integrated into the computational framework of COBRA-IIIC, and a computer code called COBRA-IIIP has been written for light water reactor analysis. The predictions of this code have been compared to those of COBRA-IIIC and it has been found that the results of the codes agree extremely well with one another for a variety of three-dimensional problems. However, numerical tests have also revealed that the COBRA-IIIP code is considerably faster and more powerful than COBRA-IIIC, particularly for problems with very large numbers of computational cells. In view of the enormous amounts of computer time that can be saved by using the COBRA-IIIP code to its full potential, it is hoped that this code will be extensively used as an alternative to COBRA-IIIC for full-core light water reactor analysis.
AcknOwledgements This research was partly supported by the Electric Power Research Institute, Palo Alto, CA, USA. The authors wish to express their gratitude to EPRI for this support. The advice given by members of the faculty of the Department of Nuclear Engineering at MIT and the staff of Batelle Pacific Northwest Laboratory during the course of this research is gratefully acknowledged.
References [1] R.W. Bowring, AEEW-R524, Atomic Energy Establishment, UK (1967).
[21 P.T. Chu, H. Chelemer and L.E. Hochreiter, WCAP-7956, Westinghouse Electric Corporation (1973). [3] [4] [5] [6] [7] [8] [9] [10] [ 11 ] [12]
A.D. Grossman et al., AEEW-R905, UK Atomic Energy Authority, Winfrith (1973). W.T. Sha and R.C. Sehmitt, ANL-8112, Argonne National Laboratory (1975). D.S. Rowe, BNWL-B82 (1971). T.A. Porsching, Nucl. Sci. Eng. 64 (1977) 177. J.E. Meyer, Nucl. Sci. Eng. 10 (1961) 269. C.W. Stewart et al., BNWL-2214, July (1977). C.L. Wheeler et al., BNWL-1962, Batelle, Pacific Northwest Laboratories, Mar. (1976). R.E. Masterson and L. Wolf, Nuel. Sci. Eng. 64 (1977) 222. R.E. Masterson, Sc. D. Thesis, Department of Nuclear Engineering, Massachusetts Institute of Technology, June (1977). W.H. Reed et al., Proc. ANS Thermal Reactor Safety Meeting, Sun Valley, Idaho, July-Aug. (1977).
293
COBRA-IIIP: AN IMPROVED VERSION OF COBRA F O R FULL-CORE LIGHT WATER REACTOR ANALYSIS Robert E. MASTERSON * Batelle Pacific Northwest Laboratory, Richland WA 99352, USA and Lothar WOLF Department of Nuclear Engineering, MIT, Cambridge MA 02139, USA Received 10 February 1978
A new type of pressure-crossflow solution method is developed for solving the fuid conservation equations used by the COBRA-IIIC code. This method is extremely efficient from a numerical point of view, and offers sizeable reductions in running time, relative to conventional numerical schemes for predicting the fluid flow and enthalpy distributions in nuclear reactor fuel rod bundles and reactor cores. This numerical method has been integrated into the computational framework of the COBRA-IIIC code, and an improved version, called COBRA-IIIP, has been developed. This code is faster than COBRAIIIC, and is designed specifically for solving larger and more complex problems, such as full-core PWR transients. Several problems are presented in which the results of COBRA-IIIP, and the method on which it is based, are compared to those of COBRA-IIIC. Excellent agreement between the two codes is observed.
1. Introduction In recent years many digital computer programs [ 1 - 5 ] have been written to solve the set o f fluid conservation equations which characterize the steady state and transient thermal hydraulic performance of nuclear systems. Although these computer programs are widely accepted and used for reactor thermal hydraulic analysis, very little effort has been made to increase the efficiency of the numerical methods used by these codes so that extremely large and complex problems can be solved within the limits of existing computational capacity. The purpose of this paper is to present a new numerical method that has been developed for solving the fluid conservation equations used by the COBRA-IIIC code, which is considerably faster and more efficient than previous methods used to analyze operational reactor conditions. Another object o f this paper is to discuss an improved version of the COBRA-IIIC code, called COBRA-IIIP, which has been developed for determing the fluid flow and enthalpy distributions in fuel pin bundles and light water reactors cores during steady state and transient conditions.
2. Mathematical formulation Owing to the widespread acceptance and use o f the COBRA-IIIC code, a considerable amount o f the physical and mathematical formalism employed by the code will also be used as the basis for the numerical method to be developed here. This essentially means that the same basic set of governing equations will be considered and that
* Former address: Department of Nuclear Engineering, MIT, Cambridge, MA 02139, USA.
R.E. Masterson, L. Wolf/COBRA-II1P. improved code Jbr full-core L W R analysis
294
the same type of solution procedure used for the initial value problem will be retained in space and time. The object of the numerical method to be proposed in this paper is to solve the fluid conservation equations as efficiently as possible so that much larger and more complex problems can be handled within the limits of existing computational Capacity. 3. Basic assumptions With these points in mind, the primary assumptions that will be made in the following development are as follows. (1) A lumped parameter approach is valid. (2) Sonic velocity propagation is ignored. (3) The diversion crossflow is usually smaller than the axial mass flow rate. (4) Viscous dissipation is neglected. (5) The liquid and vapor phases during two-phase flow are in thermodynamic equilibrium. (6) Electromagnetic body forces are ignored. (7) A homogeneous flow model can be used to describe the bulk boiling of the coolant. (8) Flow reversals, recirculating flows, and coolant expulsions are not considered.
4. Basic conservation equations Byadopting the control volume approach described in ref. [5], the equations for the conservation of mass, energy, and momentum of the fluid in a lattice composed of a number of computational cells can be written as:
con tinuity : OPi
N ~
Omi_
Ai ~--~-+~ x - - - i = l
(a)
wi] ;
energy." 1 Ohi + a h i _ q ~ t?
u i Ot
Dx
mi
N
/=1
N
N
(ti - ti) "~i+ ~
(hi - h i )wij ~ (hi _ h,) Wij . 1=I mi ]=i mij'
(2)
axial momentum." N
1 ~m i
api + _ _ = _
A-~l a-~- - 2ui ~--t- Ox
+
L 2Di
k,o, +
22A~ ~
O \ A Jil /.~X-'i" -- gPi COS 0 -- f t
.
i:
(.,-.j, 1
~
Wii
=
Ai
N
+ ~ (2u i -- u*) wiJ" j=l Ai '
(3)
transverse momentum."
where the summations indicated are to be performed by summing over all the cells in a cross-sectional plane of the lattice connected to each cell represented by the index i. In this way the differential equations for a typical cell in
R.E. Masterson, L. Wolf / COBRA-IIIP." improved code for full-core L WR analysis
295
RECTANGULARLA~I CE HEXAGONALLA'nI CE
Fig. 1. Cell numbering schemes for generating systems of difference equationS.
the hexagonal lattice shown in fig. 1 (say cell 4) can be written explicitly by simply summing from] = 1-3, and the differential equations for cell 5 in the rectangular lattice shown in the figure can be generated by summing in an analogous manner from ] = 1-4. The equations for other layouts of cells, such as those in round bundles, can be written by simply summing over all the cells in the lattice that interact with one another in such a way that the interchange of mass, energy, and momentum is allowed to occur. As in COBRA-IIIC, this set of equations is dosed by defining an additional equation for the physical state of the coolant of the form
m : p(h~, e*)
(5)
and empirical correlations for turbulent mixing coefficients, friction factors, and heat transfer coefficients. In this equation it is assumed that the pressure throughout the flow field exhibits only small variations about a known mean or system reference pressure P*. This assumption has been used on many occasions by other authors [6,7], and has found applications in the study .of both single phase and two-phase flow [8].
5. Matrix form of the conservation equations For a lattice composed of a large number of computational cells, the set of conservation equations that has been discussed here becomes unwieldy; consequently, a more compact matrix notation is desirable. In ref. 5 a vector form of the conservation equations is presented in which a transformation matrix [S] and its transpose [S] T is used to reduce these equations to a more tractable form. The elements of the matrix [S] are defined through two column vectors and where i and ] represent the indices of adjacent computational cells and K is used to represent the index of the common boundary connecting them. As in COBRA-IIIC, the entries in [S] and [S] T may only have values of 0, 1, and - 1 ; thus, the crossflow is arbitrarily considered to be positive when the dominant direction of the flow is from cell i to cell ], where i is less than ]. Using the definitions of [S] and [S] T presented in ref. 5, the foregoing conservation equations can also be written
i(K) ](K),
as:
continuity; [AI ~
+
energy: ~,
[~_/=k~j{Q}
_[ll[s]T[Ah]{w'}_[1][s]T[At]{c}+[1 ] [[h][S] T ' -
[SIT[h'l] {w};
(7)
R.E. Masterson, L. Wolf/ COBRA-IIIP: iwpoved
296 Table
code for full-core I, WR analJ’sis
I
Symbol Matrices:
Dimension
Meaning
Special properties
IAjl
NP X NP
diagonal
lcjl
NF X NT:
fcjl
NF X NF
[Djl
NF X NF
Pjl
NP x NP
I+1
NF x NF
Pfl
NF X NF
[Qjl
NP x NP
PI
NP x NF
matrix of cross sectional flow areas in the axial direction matrix of thermal conductivities between adjacent computational cells matrix of transverse friction coefficients between adjacent computational cells matrix used to compute the cross flow distribution from the transverse momentum equation matrix of enthalpies from each computational cell matrix of enthalpy differences between adjacent computational cells matrix of the enthalpies carried by the diversion cross flow between adjacent computational cells matrix of volumetric heat generation rates interface or gap connection matrix used as a lateral differencing operator
FIT
NP x NP
interface or gap connection matrix used as a lateral summation operator
NF x NF
matrix of lateral temperature differences between adjacent computation cells matrix containing the axial velocity of the fluid from each computational cell matrix containing the axial velocity of the fluid carried by the diversion crossflow between adjacent computational cells
NP x NP
NF x NF
NP X 1
{bj)
NP x 1
{hj)
NP x 1
vector containing the axial pressure gradients from each channel in the lattice source vector used by the MAT method vector containing the enthalpies from each channel in the lattice
diagonal
diagonal
diagonal diagonal
diagonal contains only entries having l,O, and -1, not diagonal contains only entries having the values of 1, 0, and -1, not diagonal diagonal
diagonal
diagonal
R,E. Masterson, L. Wolf/COBRA.IIIP: improved code for full-core L WR analysis
297
Table 1 (continued) Symbol Matrices:
Dimension
Meaning
Special properties
(mi}
NP X 1
(Pi)
NP X 1
{Pi-1}
NP X 1
{p]}
NP × 1
{u]}
NP X 1
{u]*}
NF X 1
(wj}
NF X 1
vector containing the axial mass flow rates from each channel in the lattice vector containing the pressure field at axial level] vector containing the pressure field at axial level] - 1 vector containing the density of the fluid in each channel in the lattice vector containing the axial velocity of fluid in each channel in the lattice vector containing the axial velocity of the fluid carried by the diversion crossflow vector containing the cross flow between adjacent computational cells
-
-
Miscellaneous constants." s/l
-
x Ax At 0
--
transverse momentum factor axial elevation axial mesh spacing time step size weighting function used by the MAT method
NF = number of boundaries for crossflow between adjacent cells. NP = number of nodes at which the pressure is to be found.
axial momentum: -
2u
/~/-
{a'} +
[[2ul IS] T - [slT[u*I]
(w},
(8)
w h e r e the c o m p o n e n t s o f the pressure d r o p due to frictional forces, gravitational forces, and the e f f e c t s o f turbulent m i x i n g are given b y +--
-l\--A] ~2-D 2&,¢
+A
3x
]
pgcos4~
-ft
[ s ] T [ A u ] (w'}
(9)
and transverse i n t e r a c t i o n s are t a k e n i n t o a c c o u n t b y m e a n s o f the m a t r i x f o r m o f the transverse m o m e n t u m equation d e v e l o p e d for the COBRA-IIIC code:
transverse momentum:
298
R.E. Masterson, L. Wolf/COBRA-IIIP: improved code for full-core L WR analysis
For convenience the same nomenclature and notation used in ref. [5] are also used here. Consequently, although most of this nomenclature is summarized in table 1, ref. [5] should be consulted for more detailed information about the definition and derivation of the individual terms.
6. Conversion of the conservation equations into difference equations The preceding conservation equations must be converted into a system of difference equations before they can also be solved by numerical means. Taking a first order backward differencing scheme for the spatial and temporal derivatives to provide numerical stability allows the equations for the conservation of mass and axial momentum to be written as
continuity:
axial momentum:
For reasons of computational efficiency, the energy equation is approximated using a differencing scheme which is spatially explicit and temporally implicit:
energy: [ l lth/--hJi*[h/-h/-1}
= {mi} -1 { { Q i - t / 2 } - [ S ] [ A h j - t ] ( w ; - t } - [ s ] T t A t j - t ] [ c j - t ] }
+ (mi)--l[[hj_l] [S]x - [S]T[h~_I]I .
(13)
On the other hand, the energy equation can be made fully implicit, and the coolant enthalpy distribution can be determined by a solution scheme similar to that given in ref. [9]. In either case, it is convenient to visualize each control volume of fluid as being bounded by two adjacent control volumes in the axial direction; thus, the amount of mass, energy, and momentum that can be interchanged also depends on the number of computational cells connected to one another in each cross-sectional plane of the core.
7. The transverse momentum equation The transverse momentum equation to be considered in this analysis is similar to the transverse momentum equation used by the COBRA-IIIC code [5]. However, it differs from previous versions of the equation because it is based on the concept of splitting the pressure field (P) and the transverse friction terms [C] (w} which drive the crossflow distribution into a sum of spatial implicit and explicit parts. The algorithm which governs the form of the terms is
(P) = 0 (Pi} + (1 - O) (Pi_ l) ,
(14a)
[C] {w} = 0 [Cj] {wi} + (1 - O)[Ci_ 1] (wi- ,},
(14b)
where 0 is a weighting function having an arbitrary value between 0 and 1.0. As a result of this algorithm the trans-
R.E. Masterson, L. Wolf/COBRA-IIIP: improved code for full-core L WR analysis
299
verse momentum equation, eq, (10), can be written in finite difference form as
---XT--] +
ax
+
(O[©l(wj)+(l-O)[Q_ll(wj_~}}= [sl(o{Pj)+(l-O)(pj_l)}, (15)
where 0 ~<0 ~< 1.0. Again, it is worth mentioning that there are no numerical restrictions on the time step size, At, since eq. (15) and all the other equations which govern the behavior of the fluid are temporally implicit.
8. Solution of the multidimensional COBRA-IIIC conservation equations with the MAT method In order to solve the conservation equations presented in this paper, it is necessary to develop a numerical method for computing the fluid flow and enthalpy fields in an operational reactor core. Since the conservation equations have already been cast into a form suitable for digital computer analysis, the development to follow will consist of identifying techniques for solving these equations which can bring about sizeable reductions in computing time, relative to conventional numerical methods. In ref. [10] a pressure-velocity method is presented in which enormous amounts of computer time can be saved by developing an iteration scheme in which the pressure fields in the lattice are found as a function of the axial mass flow rate. Unfortunately, problems with numerical stability are encountered with this scheme when it is used for a flow blockage analysis because the frictional and gravitational terms in the axial momentum equation, eq. (8), are approximated by spatially explicit relationships of the form {a') = - { K j } ( m / _ l } 2 - {fl'-l),
(16)
where {Kj} is the coefficient of the flow-squared terms and all the other terms (including the gravitational terms) are lumped into (t)_ l }. In this paper a modified version of this numerical method will be presented in which this shortcoming is overcome by replacing eq. (16) by a differencing scheme that is both spatially and temporally implicit: {a') = -{Kj}(mj) 2 - {fj} .
(17)
The numerical method which results from this approximation is a new type of pressure-velocity method called the MAT method (Modified and Advanced Theta method). In what follows it will be shown that the MAT method combines the computational efficiency of the method presented in ref. [10] with the stability of the crossflow solution scheme used by the COBRA-IIIC code. Finally, it will be demonstrated that this combination of computational characteristics makes it possible to save a remarkable amount of computer time when a method of this type is used for solving large problems in which there are either single-phase or two-phase flows.
9. The MAT method: derivation of the difference equations for the pressure field If the continuity equation, eq. (11), is written as (18) it is possible to express the mass flow rate at axial level ] in terms of the mass flow rate at axial level ] - 1 by means of a relationship of the form
(m/) = { m / _ l ) + (Am),
(19)
300
R.k'. Masterson, L. Wolf~ COBRA-II1P: improved code jor fidl-core L WR ana(vsis
where
{~xm} = --~xx [S] t ( w j ) - ~ ¢ [A/] l~J ~?P~J} k
(20)
l
is the change in the axial mass flow rate between tile top and the bottom of each plane of computational cells. Using this relationship, the terms in the axial momentum equation which account for the pressure drop due to frictional and gravitational effects can be written as (2 l a)
(a;} = - {k/}(m~} - ~t).} = -(ki}(mi-1
(21b)
+ A m } 2 -- (J)}
= - { k i } ( r n ~ _ , + 2rnj_ ,Av'n + ~ n 2}
~.},
(21c)
where (Am} is generally smaller than {rni_ 1}. In many cases the stability of this difference approximation can be improved by applying the continuity equation {~n} = {m/} - {m/_ l}
(22)
to write the flow-squared terms as a linear function of { ~ z ) : (rn/} 2 = (my_ ,}2 + {mj + rni_ 1}{&m) .
(23)
The terms in eq. (2l) which account for the effects of frictional and gravitational forces on the axial pressure gradient are then given by {a'}
= -{kj}{mj_ 1} 2 - {kj}{rttj + rrtj_ 1}{~11} -- {f/-} .
(24)
Applying the same procedure to the temporal acceleration term in the axial momentum allows this term to be written as [Aft -1 At {mj
-
~i} -
[Aj]-I At
{mj -- m i _ 1 + m / _ 1 - in/} - [ Aat j ] - 1 { m i - mi-1} +
[Aft-'
at - { m / - 1 - mfl '
(2s)
and since the continuity equation requires that {ZXrn) = {mj} - {my_ l ) ,
(26)
the temporal acceleration term can also be expressed as a linear function of {Am): [ AAt I ] - I { m i _ ~ i } = [A/]_ 1 {~mt ~-}+[A/]
-1
{m J - l~- Ftn J
}.
(27)
Substituting these linear relationships for the frictional and temporal acceleration terms into the axial momentum equation, eq. (12), allows the conservation of momentum in the axial direction to be written as ax
at
J
[Ail
at
+ [Aj] -I [[2uj] IS] T - [s]T[uTII {w/).
- {ki}{mj_ 1} = - 0')} -
+ {/9}{,ni + m i -
(28)
The terms containing {am} can then be eliminated from this system of equations by combining eq. (28) with eq. (20). The result of this operation is a system of difference equations which simultaneously satisfies the conservation of mass and momentum at each axial level of the core at each instant of time:
R.E. Masterson, L. Wolf/COBRA-IIIP." improved code for full-core L WR analysis
+,A,,,
301
,s,w,u.,1
The effects of lateral m o m e n t u m transfer can be incorporated into eq. (29) by combining it with the transverse m o m e n t u m equation, eq. (10), to eliminate the crossflow-dependent terms. Since the transverse m o m e n t u m equation can be written as 1 +--u; + s 0[Ci ] {w/} = At Ax
+ ui-~xi-1 +
(0-1)[Cy_l]{Wi_l}+
[S]{O{P/}+(1-O){Pi_l}}
,(30)
the crossflow distribution {wj} at axial level j can be expressed as a function of the pressure fields at levels/" and j - 1 by means of a relationship of the form
{wi} = [O3 -'{0 i} + (~)[Oj] -1 [Sl{O{Pj} + (1 - O)(Pj_ ,}},
[Dj] -1 is the
where
(31)
diagonal matrix defined by
[Oyl-' = [ l ' l - U ] * + LAt 2v¢
(32)
O[Cjl
and
[[atJ[
Ax
J
is a source vector consisting of the other forcing terms in the transverse m o m e n t u m equation. If eq. (31) is then combined with eq. (29) to eliminate the crossflow-dependent terms, it is possible to derive the following system of difference equations for the pressure field as a function of the axial mass flow rate:
{Pj - PI-1 } = ~
+ ~[A/I-'
[Ail-'{ff9 - m j _ l } -- ~ {kj } {my2_~} -- A,v. ~.} + - ~ {pj -- ~j} 2uj + - - ~ + ,5x kjA j(rnj + m j _ l )
2 u / + - - + Ax[Aj]kj(mj + mj_0] [S] T - [s]T[u; At
[Dyl-'{Oy}+
[Dy]-I[S]{0{Py}
x
+ (1
- o) {Py_,}}}.
(34)
This system of equations for the pressure field can also be written in matrix notation as
[I+ (1 - O)Myl{Py_ 1} = [ I - OMyl{Py} + {by},
(35)
where I is the identity matrix,
[My] = A x ( / )
[[Aj]-'
[By][sIT[D/]-'
iS] -
[A/l-'
[S] T[U]*] [D/l-' [Sl]
(36)
is a matrix used for determining the pressure field,
{by} =--~ [ A j ] - ' { ~ y - m y _ ,} - 2~x{kj}{ml 2._,} - Ax ~.} +Ax at {py - -~y} 2uy + + Ax [A/] -1 [By] [ s ] T [ D j ] - ' { O j }
+ axkjA1{m j + my_,)
- A x [ A i ] - 1 [S] T[u~] [1)/] -l{Oj}
(37)
is a source vector consisting of terms that contribute to the axial pressure gradient, and
[By] = [2uy] + ~ - + Ax[A/]ky(m/+ my_,)
(38)
R.E. Masterson, L. Wolf/COBRA-IIIP." improved code for full-core L WR analysis
302
is a matrix which is used in the construction of both the source vector and [Mj]. Although the source vector that drives the pressure field consists of only one entry from each computational cell, it is worth noting that this vector contains all the physical forces that can act upon an elemental volume of fluid. These forces can be due to frictional effects, gravitational effects, and the spatial and temporal acceleration of the flow. Finally, although the axial mass flow rates used in eq. (35) are not actually known, it is possible to develop simple procedures for estimating the values of these flow rates and updating them through iteration. Procedures of this kind are now used by the COBRA-IIIC and COBRA-IV-I codes [5,9]. 10. Discussion of the difference equations governing the pressure field It is interesting to observe that the use of the weighting function 0 in the transverse momentum equation enables a system of difference equations to be derived in which the crossflow distribution can be driven by any combination of the pressure fields that exist at the top and the bottom of a plane of computational cells. This can be seen more clearly by writing eq. (35) as
{P/- 1} = [I + (1 --0 )1~]]-1 [I -- OMj] {Pj} + [I + (1 -O)Mj] -l{b/},
(39)
and by examining the structure of the matrix equations governing the pressure field {Pi_ 1} for various values of 0. Setting 0 = 0 allows eq. (39) to be written as: 0=0:
(P]- 1) = [I +/111/]--1 [I - Mi] (Pi) + [I + Mj] -l(b/) .
(40)
If 0 is set to ½ and 1.0, respectively, eq. (39) then becomes 0=1.
2'
{Pj_ ,} = [I + ~,Mj]-1 [I - ~ j ] {ej} + [I + ~,Mj]-l{bj}
;
(41)
0 = 1.0:
(P/- 1} = [I - Mj] ~pj} + {bj} .
(42)
These equations demonstrate that a direct relationship exists between the form of the difference equations for the pressure field and the forces used to drive the crossflow distribution between computational cells in the transverse plane. For example, setting 0 = 1 allows the crossflow distribution to be driven by the pressure field that exists at the top of each plane of computational cells. From eq. (42) it can be seen that this results in a differencing scheme that is spatially explicit but temporally implicit. Setting 0 = 0 means that the crossflow distribution is to be driven by the pressure field that exists at the bottom of each plane of computational cells. In this case the difference equations governing the behavior of the fluid are both spatially and temporally implicit. Finally, when 0 is set equal to ½, eq. (41) shows that it is possible to generate a system of equations where the crossflow distribution is driven by the average of the pressure fields that exist at the top and the bottom of each plane of computational cells. In this case the difference equations are temporally implicit but have a spatial component whose structure is analogous to that of the Crank-Nicholson method. Of course, it should be evident that other spatial differencing schemes can be constructed by simply choosing other values for 0.
11. Solution schemes for the pressure distribution If the difference equations governing the computation of the pressure field are written as
(Pj-I} = [I+ (1 - O)Mj]-l[I - OM/]~p/} + [I+ (1 - O)Mj]-Ifbj}
(43)
R.~: Masterson, L. Wolf/COBRA-IIIP: improved code for full-core L WR analysis
303
it is possible to invert the coefficient matrix [K + (1 - 0)
Mj]
(44)
and solve these equations for {Pj_ l) as a function of (Pi). Since the pressure field at axial level j - 1 is upstream of the pressure field at axial level j, this allows disturbances in the flow field to be felt at upstream locations. One of the primary computational advantages of the pressure-oriented solution scheme outlined here is that the coefficient matrix governing the computation of the pressure field is a Stieltjes matrix. In other words, it is possible to show that the coefficient matrix in eq. (43) is a diagonally dominant, irreducible, positive definite matrix with a simple and predictable band structure for a variety of subchannel numbering schemes. Since this enables the pressure distribution to be found quite easily by both direct and iterative solution schemes, the computation of the pressure field requires considerably less effort than in other approaches where the crossflow distribution is found directly as a function of the axial mass flow rate [11].
12. Derivation of the crossflow distribution from the pressure distribution Whereas the numerical methods used to solve the momentum equations in previous versions of COBRA have been based upon the concept of solving directly for the crossflow distribution as a function of the axial mass flow rate [5,9,11 ], the pressure-velocity method discussed in this paper is based upon the alternative numerical approach of first solving for the pressure field as a function of the axial mass flow rate and then using this pressure field to determine the crossflow distribution from the transverse momentum equation. In other words, this approach assumes that substantial transverse pressure gradients can exist in a reactor lattice and that these gradiants are the primary physical mechanism for driving the crossflow distribution. After the inlet flow and enthalpy distributions have been specified, eq. (35) can be solved for {Pj_ 1) and a pressure field can be constructed for driving the crossflow distribution between the cells in each computational plane. This pressure field can be found by blending together the pressure fields from the top and the bottom of each plane of cells to form the composite pressure field given by eq. (14a) as the numerical solution scheme sweeps downstream. The crossflow distribution can then be found from the transverse momentum equation by writing it in the form
(w/} = [DjI-I(oj) + (~)[Ojl-I[s](o(Pj}
+(1 - o ) { e j _
1)),
(45)
where [Dj]-I is the designed matrix defined by eq. (32) and (Oj) is a vector consisting of the other forcing terms in eq. (30). Since the inverse of [Dj] can be found by inspection, this additional step in the solution scheme requires very little computational effort to determine the flow and enthalpy fields once the pressure distribution is known.
13. Solution schemes for the energy equation After the equations for the conservation of mass and momentum have been solved, it is necessary to interface the thermal response of the fluid with the hydraulics through heat transfer coefficients. This can be done by soNing either an explicit or implicit energy equation for the enthalpy of the fluid at axial level j as a function of the enthalpy of the fluid at the previous axial level and the amount of heat produced in each rod in the lattice. An energy balance can be performed using eq. (13) and the axial mass flow rates (mj} found from a previous solution of the momentum equations. Thus, the solution of the energy equation can be looked upon as the 'outer iteration' in this two-step solution procedure. For problems in which there are substantial departures from operational reactor conditions, several additional iterations between the energy and the momentum equations may have to be performed.
304
R.E. Masterson, L. Wo(f / COBRA-IIIP. improved code jbr fidl-core L WR ana(vsis
14. Overall iterative strategies At this point it may be helpful to summarize how the MAT method can be used to find tile flow and enthalpy fields in an operational reactor where there can be strong coupling between the energy and momentum equations. Referring to the flow chart shown in fig. 2, a typical calculation consists of sweeping downstream several times in succession between the inlet and the outlet of each channel in the core. Once the inlet flow and enthalpy distribution have been specified, the enthalpy can be advanced from axial level / - I to axial level / by solving either an implicit or explicit version of eq. (7). For the first axial iteration, the flow rate (mi} for axial level / is set equal to the flow rate from level j - 1; otherwise the value of {rn/} from the previous iteration is used. After eq. (35) is solved for the pressure distribution {Pj_ 1} the pressure fields from adjacent axial levels are blended together using eq. (14a) and the transverse m o m e n t u m equation is solved for {w/}. Since the value of {Pj} in eq. (35) is not initially known, it is set equal to a constant (such as zero to avoid the necessity of computing the source term [1 0m]
ESTABLISHINITIAL CONDITIONS ONh, w, ANDm
SOLVETHEENERGYEQUATION,
EQ. (/)TO GET(hi)
I FORFIRSTAXIAL ITERATION ASSUME: (rni) (rnj.l) •
i
SOLVEEQ. (35)FORTHE PRESSUREFIELD(Pi_lI
I BLENDTOGETHERTHEPRESSURE FIELDSUSINGEQ. ll4t
I I
SOLVEEQ. 145)FORTHE ~ ] CROSSFLOWDI STRI BUTION
I
] FINOA NEWVALUEOFIm ) FROM I THECONTINUITYEQUATON, EQ. i18 ] i
l
l
SWEEPDOWNSTREAMTOTHE NEXTAXIALLEVELAND REPEATTHISLOOPAGAIN
I
DOANOTHERAX,AL
1.0[---WHENTHEITERATIONREACHES ~
ITERATIONWITH NEWVALUESFOR(mjl I
TES'FORCO VERGENCE
l
THEENDOFTHECHANNEL
~, YES PRINTRESULTS CONDITIONSAND CONTINUETHECALCULATION
DOANOTHERTIMESTEP?
Fig. 2. Flow chart for the MAT method.
R.E. Masterson, L. Wolf/COBRA-IIIP: improved code for full-core L WR analysis
305
(Pj} for the first axial iteration) and then updated so that its new values can be used for successive iterations. After the axial mass flow rate {mj} is updated using the continuity equation, eq. (18), the numerical scheme is allowed to sweep downstream to the next axial level, where this entire process is repeated again. When the calculation reaches the exit of each channel, a check is made to see if the flow distribution has converged at every mesh point in the lattice to the tolerance desired. If the solution has converged, the computation stops; otherwise, the entire iterative process is repeated again, starting at the first axial level (e.g. the core inlet) and sweeping downstream to the last axial level (e.g. the core outlet). For all iterations after the first, the most recent values of the axial mass flow rates are used for computing the flow and enthalpy fields. This eliminates any errors that may have been introduced into the solution by estimating a value of (m/) for the first axial iteration. For all practical problems, it has been found that this iterative procedure yields a unique, stable, and physically realistic solution to the fluid conservation equations in a nuclear reactor core.
15. The COBRA-IIIP code In order to take advantage of the computational capabilities of the MAT method, a new version of the COBRAIIIC code, called COBRA-IIIP *, has been written for light water reactor analysis. This version of the code is based upon the same physical formulation as previous versions of COBRA, but due to the use of the MAT method and substantial improvements in computer coding, the new code has the capability to solve much larger and more complex problems than COBRA-IIIC with much greater speed. For example, the code has been used to model an entire PWR fuel pin bundle during a 5-min loss-of-flow transient. It has also considered problems having enough computational cells to represent over 600 fuel assemblies in a reactor core explicitly. To more clearly understand the capabilities of the new code, it may be helpful to summarize the types of problems to which the code can be applied. These problems are generally the same as those that can be solved by COBRAIV-I. However, since the code does use a 'MAC' methodology for solving the momentum equations, it cannot be used for solving problems where there are substantial flow reversals and recirculations. Problems of this type can only be solved by numerical schemes such as those developed for the TRAC code [ 12] and the explicit part of COBRA-IV-I [91.
16. Comparisons of computation time To demonstrate the computational advantages of the MAT method several problems were run in which the computation times of the COBRA-IIIP code was compared to that of COBRA-IIIC. In these problems the number of cells in each plane of a pressurized water reactor core was varied between 16 and 400 and only rectangular mesh grids were used. In all cases the flow fields were converged to a tolerance of 1%; that is, the relative change in the axial mass flow rate, dm/m, at each point in the core was constrained to change by less than 1% between successive axial iterations. Moreover, the momentum equations in COBRA-IIIP were solved by both direct and iterative solution techniques to demonstrate the advantages of different methods of matrix inversion as a function of problem size. The results of the timing run comparison are summarized in table 2. In all cases it was found that an enormous reduction in computation time could be achieved by using the COBRA-IIIP code and the MAT method in place of the crossflow solution technique in COBRA-IIIC. Whereas the crossflow solution time of the COBRA-IIIC code was found to increase as the number of channels cubed, the crossftow solution times of COBRA-IIIP increased as only the 1.5 power of the number of channels for a direct inversion scheme, and as almost a linear function of the number of channels when the method of successive over-relaxation was used. For large problems, such as full* The letter P is used to designate the fact that the code is based on a pressure solution scheme rather than the crossflow solution schemes used by COBRA-IIICand COBRA-IV-I.
306
R. bZ Masterson, L. Wolf / COBRA-IIIP. improved code for full.core L WR analysis
Table 2 A comparison of the crossflow solution times a) of the COBRA-IIIC and COBRA-IIIP codes Code
COBRA-II1C
COBRA-IIIP
COBRA-IIIP
Method of solution
Gauss elimination with no compression of crossflow coefficient matrix Gauss elimination with full compression of pressure coefficient matrix successive overrelaxation with optimized relaxation factor
Number of channels per plane 16
64
128
200
300
400
1.7
128.49 b)
987.36 b)
3243.61 b)
10 250 b)
31 980 b)
0.17
0.76
2.16
4.14
7.95 b)
16.53 b)
0.15
0.44
1.23
1.82
2.65 b)
3.62 b)
Notes: All results are for ten axial levels on an IBM 370/168 computer with the H compiler. Note that the results of these timing runs may vary by -+10% during the course of a day due to changes in the work load on the system. a) In seconds of central processing unit time b) For economic reasons these results are estimated by extrapolation.
scale simulations o f PWR fuel assemblies, the COBRA-IIIP code and the MAT method were found to be over three orders of magnitude faster than COBRA-IIIC. It should be emphasized that these comparisons were made by using exactly the same types o f iteration schemes in both codes. In other words, the flow and enthalpy fields were recomputed over the entire axial height of the core whenever a single point was found where the solutions were not fully converged. For problems where the effects of natural convection are much more important than those o f forced convection, iteration schemes o f this kind are desirable because they allow details o f the downstream flow distribution to be felt at upstream locations. However, for many other practical problems these schemes have the drawback that they consume an enormous amount of computation time because they require the flow fields to be recomputed at every mesh point in the core - whether or not this recomputation is actually needed. In order to overcome this difficulty, an alternative iteration scheme has been developed for COBRA-IIIP in which the effects o f upstream pressure propagation are suppressed by performing additional iterations at only those points in the lattice where they are needed. This iteration scheme consists o f solving eq. (35) for the pressure field and then solving eq. (45) for the crossflow distribution to complete the first radial itejration. If the change in the axial mass flow rate, {Am) = {mj - m / _ 1), between the top and the b o t t o m o f each node in less than some prescribed convergence criterion, the flow field is assumed to converge, and the numerical method is allowed to sweep downstream to the next axial level, where this process is repeated again. If, however, the flow field has not converged, the axial mass flow rates are updated via the continuity equation, eq. (I 8), and a new pressure field is found b y solving eq. (35). A new crossflow distribution is then computed from the transverse m o m e n t u m equation, and the computational cycle is repeated again - with feedback from the energy equation, if necessary - until a converged flow and enthalpy distribution is found. Since the flow and enthalpy fields are only updated at points where additional iterations are needed, it is often possible to reduce the computation time o f COBRA-IIIP b y another factor o f 3 or 4 compared to COBRA-IIIC. In other words, a significant reduction in computing time can be achieved by simply developing a more efficient numerical method and a more efficient iteration scheme for solving the same set o f m o m e n t u m equations that is
R.E. Masterson, L. Wolf / COBRA-IIIP: improved code for full-core L WR analysis
307
used by the COBRA-IIIC code. For extremely large problems, such as full core loss-of-flov~ transients with 600 channels, this saving in computer time may be as great as three and one half orders of magnitude.
17. Results To demonstrate that the predictions of the COBRA-IIIP code agree with those of COBRA-IIIC, a number of sample problems were constructed in which the flow and enthalpy fields were computed by both codes in an operational reactor core. Since the COBRA-IIIC code was not capable of solving problems as large as those which could be solved by COBRA-IIIP, the comparisons were restricted to problems in which reasonable numbers of computational cells could be used. The first sample problem to be considered is shown in fig. 3. In this problem, the departure from nuclate boiling ratio (DNBR) was computed by both codes for the hottest rod in the lattice of the PWR whose operating parameters are shown in table 3. The results of these calculations are summarized in fig. 4. Inspection of these curves reveals that the COBRA-IIIP code gives a departure from nucleate boiling ratio that is almost exactly the same as that given by COBRA-IIIC. Moreover, fig. 4 shows that the Westinghouse (W-3) correlation apparently gives more conservative estimates of the departure from nucleate boiling ratio than the Babcock and Wilcox (B&W-2) correlation over the entire axial length of the core. The predictions of the two correlations can be seen to deviate from one another by as much as 20% in the hottest region of the channel for the operating conditions given in table 3. Fig. 5 shows a cross-sectional view of a second sample problem that was used to compare the predictions of the COBRA-IIIC and COBRA-IIIP code. This problem consists of two identical fuel pin bundles of a PWR with different radial power peaking factors whose size and operational characteristics are described in table 4. In this problem, a grid with a loss coefficient of 2.5 is located midway between the inlet and the outlet of each assembly and both assemblies are subjected to a power transient in which the power level is doubled uniformly in 1 s. To follow the behavior of the assemblies during the transient, it was assumed that each of the assemblies could be homogenized (i.e. represented by a single computational cell) and that an axial mesh spacing of 2 in. could be used. Table 5 shows the coolant temperatures, enthalpies, densities, and axial mass flow rates predicted by both codes at the exit of the hot channel as a function of time. It can be seen that the predictions of both codes are in excellent agreement with one another for this particular case. The crossflow distribution between the assemblies is also plotted at the end of the transient in fig. 6. Inspection of these curves reveals that the magnitudes and the shapes of the crossflow distributions are essentially the same along the entire length of the core. However, some small
Table 3 Input parameters for sample problem 1 Type of reactor - pressurized water Axial power distribution Radial power peaking factors Nominal operating conditions: System outlet pressure: Inlet temperatures (uniform for all channels): Average inlet mass flux (uniform for all channels): Average reactor heat flux: Axial mesh spacing: Channel length:
skewed cosine as shown in fig. 3 2100 psia 541.0°F 2.51 X 106 lbm](h ft 2) 0.21 X 106 Btu/(h ft 2) 6 in. 180 in.
R. l'#. Masterson, L. Wolf~ COBRA-IIIt'. improved code for/idl-core L WR aptalvs'is
308
DISTANCE FROM INLET, cm
0 9.O I
100
200
D0
300
T
----COBRA IIIP
I
8.0~
--COBRA
IIIC
i! ~
x, . ~ . / H O T ASSEMBLY (FOR DETAILS SEE BELOW] ........
B&W-2 ~
_o ~ 6.0
,
.194
..... L .,,__.~'__~._
o.9~
g
:-.__~___~ q---1---'---~--+.._, "--, ,
1.351
7.o i
3.0
w ~3
CORRELATION
1.5,~...~m 11.5,~I a.5~ ~...~l.ssl ~.
:~
"
N
2..0
HOT ROD
1.0
00 ~ 4 9 ~ •
~
G CTORS I
1
1
40
80
120
160
DI STANCE FROM INLET, in.
Fig. 3. A one-eighth full-core symmetry section of a PWR (powers normalized to 1.0). Fig. 4. Axial distribution o f the DNBR predicted by COBRA-IIIC and COBRA-IIIP for the hottest of the lattice shown in fig. 3.
Table 4 Input parameters for sample problem 2 Type of reactor
pressurized water
Axial power distribution Radial power peaking factors Nominal operating conditions: System outlet pressure: inlet enthalpies (uniform for all changes): Average inlet mass flux (uniform for all channels): Average reactor heat flux: Channel length: Axial mesh spacing: Flow area for each cell: Number of rods per cell:
uniform as shown in fig. 5 2 100 psia 538.0 Btu/lb m 2.48 X 106 lbm/(h ft 2) 0.2 × 106 Btu/(h ft 2) 120 in. 2 in. 0.266 ft 2 225
R.E. Masterson, L. Wolf / COBRA-IIIP: improved code for full-core L WR analysis
N
8"
.i./ / / /
"
309
RADIALPOWER . PEAKINGFACTORS DISTANCEFROMINLET,cm 100 200
8.0 7.0
300
COBRA-IIIP COBRA-IIIC
---- -
6.0 "COLD"~
/
ASSEMBLY~ ; E A C T O R O U T L E T /
'I-lOT"
5.0
ASSEMBLY ASS
4.0
3.0
2.0
|,
,, ]
GR,D 1.0
//
L/
/x
,,
GRIDLOCATION.
O.O
/
.
.
ONSETOF . ~ ~
I
I
i
3O
6O
9O
120
DISTANCEFROMINLET, in.
REACTORINLET Fig. 5. A cross-sectional view of two homogenized PWR fuel assemblies.
Fig. 6. Axial crossflow distributions predicted by COBRA-IIIC and COBRA-IIIP for the problem shown in fig. 5.
Table 5 Values of the fluid properties predicted by COBRA-IIIC and COBRA-IIIP at the outlet of the hot channel as a function of time Time (s)
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
Enthalpies(Btu/lbm)
Densities ( l b m / f t 3)
Mass fluxes (Mlbm/ h ft 2)
Temperatures (°F)
COBRAIIIC
COBRAIIIP
COBRAIIIC
COBRAIIIP
COBRAIIIC
COBRAIIIP
COBRAIIIC
COBRAIIIP
627.70 629.08 631.84 635.90 641.21 647.56 654.78 662.67 671.03 679.78 688.73
627.72 629.10 631.86 635.90 641.20 647.55 654.76 662.64 670.99 679.75 688.69
41.65 41.57 41.40 41.16 40.84 40.44 39.99 39.48 38.93 38.34 35.94
41.65 41.57 41.40 41.16 40.84 40.45 39.99 39.49 38.94 38.36 35.96
2.464 2.479 2.490 2.499 2.505 2.509 2.511 2.511 2.510 2.508 2.501
2.462 2.480 2.493 2.501 2.506 2.510 2.513 2.513 2.512 2.510 2.509
607.30 608.24 610.08 612.80 616.29 620.42 625.01 629.96 635.08 640.29 64.2.81
607.33 608.26 610.09 612.79 616.29 620.40 624.98 629.93 635.05 640.26 642.78
310
R.E. Masterson, L. Wolf/COBRA.IIIP. improl,ed code for full.core L WR analysis
differences between the predictions of the codes can be seen at the location of the onset of bulk boiling and in the vicinity of the grid. These differences are due primarily to the differences in the axial iteration schemes discussed in section 16. Nevertheless, in view of the fact that they are considerably smaller than the uncertainties in some of the empirical correlations used by the codes, it is felt that these differences are quite small a price to pay for such an enormous increase in computational efficiency.
18. Conclusion In conclusion, a new numerical method has been developed for solving the set of fluid conservation equations used by the COBRA-IIIC code. This method is considerably faster and more efficient than previous methods that have been used for reactor thermal-hydraulic analysis. The new method of solution has been integrated into the computational framework of COBRA-IIIC, and a computer code called COBRA-IIIP has been written for light water reactor analysis. The predictions of this code have been compared to those of COBRA-IIIC and it has been found that the results of the codes agree extremely well with one another for a variety of three-dimensional problems. However, numerical tests have also revealed that the COBRA-IIIP code is considerably faster and more powerful than COBRA-IIIC, particularly for problems with very large numbers of computational cells. In view of the enormous amounts of computer time that can be saved by using the COBRA-IIIP code to its full potential, it is hoped that this code will be extensively used as an alternative to COBRA-IIIC for full-core light water reactor analysis.
AcknOwledgements This research was partly supported by the Electric Power Research Institute, Palo Alto, CA, USA. The authors wish to express their gratitude to EPRI for this support. The advice given by members of the faculty of the Department of Nuclear Engineering at MIT and the staff of Batelle Pacific Northwest Laboratory during the course of this research is gratefully acknowledged.
References [1] R.W. Bowring, AEEW-R524, Atomic Energy Establishment, UK (1967).
[21 P.T. Chu, H. Chelemer and L.E. Hochreiter, WCAP-7956, Westinghouse Electric Corporation (1973). [3] [4] [5] [6] [7] [8] [9] [10] [ 11 ] [12]
A.D. Grossman et al., AEEW-R905, UK Atomic Energy Authority, Winfrith (1973). W.T. Sha and R.C. Sehmitt, ANL-8112, Argonne National Laboratory (1975). D.S. Rowe, BNWL-B82 (1971). T.A. Porsching, Nucl. Sci. Eng. 64 (1977) 177. J.E. Meyer, Nucl. Sci. Eng. 10 (1961) 269. C.W. Stewart et al., BNWL-2214, July (1977). C.L. Wheeler et al., BNWL-1962, Batelle, Pacific Northwest Laboratories, Mar. (1976). R.E. Masterson and L. Wolf, Nuel. Sci. Eng. 64 (1977) 222. R.E. Masterson, Sc. D. Thesis, Department of Nuclear Engineering, Massachusetts Institute of Technology, June (1977). W.H. Reed et al., Proc. ANS Thermal Reactor Safety Meeting, Sun Valley, Idaho, July-Aug. (1977).