5. Neutronic and thermal-hydraulic coupling model (NK-TH)

The description contained in this section is based on a work published by Ceceñas in Ref. [9] about a TH model developed for boiling water reactors. The TH model was modified from a point kinetics approach with an extension of the NK model to 3D and implemented in the development of AZKIND.

The treatment of neutron kinetics in [9] has been improved by coupling a 3D solution of the neutron diffusion equations with an arrangement of TH channels in parallel. Each channel independently contemplates three regions: (1) one phase, (2) subcooled boiling, and (3) bulk boiling. The objective was to implement a detailed model of a nuclear reactor core, which is somehow perturbed to simulate NK-TH coupling. These perturbations are obtained when the power generated in a group of channels changes and thus affecting the TH state of each channel.

The original [9] TH model is based on a generic channel, which is adapted by transferring to it the operational data as flow area, generated power, axial power profile, and subcooling, among other parameters. Each channel is associated with a number of nuclear fuel assemblies and an axial power profile. Although the neutron model is a two-dimensional model for the radial power profile in each z-plane covering all the channels, information related to the axial power distribution is considered for each individual channel. In Ref. [9], it is assumed that this steady-state axial power profile is invariant over time, and it is used to weight the axial averages of macroscopic cross sections and void fractions. To perform the numerical implementation of the model, the arrangement of channels is obtained by grouping the total core assemblies into an appropriate number of thermal-hydraulic channels, which gives a definition of a set of channels per quadrant.

For the implementation in AZKIND of the TH model of Ceceñas, the grouping of fuel assemblies was maintained for generating a reduced number of TH channels; operational data are also used. The main difference is that the NK model recursively computes the axial power profile for each channel, and this thermal power is the updated source of power for TH model. Therefore, a "new" thermal-hydraulic condition is generated, and it is used by the NEMTAB model to update the nuclear data to generate new thermal power profiles with the NK model. The process is iterative, and it stops when the convergence is met. Convergence is achieved when updated conditions do not change in both NK and TH models.

The NK-TH coupling in AZKIND performs core calculations as described above to obtain a steady-state reactor core condition. For transient conditions in a time interval T, the NK-TH coupling process is the same for each time step ΔT in T, that is, a different quasi-steady-state condition for each successive ΔT. Achieving converge for each ΔT with respective reactor core conditions means to produce a time-dependent behavior of the reactor condition over the total time interval T.

The TH model comprises the solution of the mass, momentum, and energy conservation equations in the three regions contemplated by the channel: (1) one phase, (2) subcooled boiling, and (3) bulk boiling. The system receives heat through a non-uniform source whose profile is axially defined plane by plane. This axial use of the power profile allows the inclusion of a wide range of axial profiles, from relatively flat to profiles with their peak value at some axial point in each channel in the core.

In the following subsections, there are several expressions for which the corresponding parameters are defined in Refs. [10, 11].

#### 5.1. Heat transfer in the fuel

In summary, once the NK model is used to generate the neutron flux distribution in the reactor core, expression (12) can be used to calculate the thermal power being generated along all the nodes in a thermal-hydraulic channel of area Δa and height H. This thermal power can be the axial power profile needed by the TH model to produce the thermal-hydraulic state

The description contained in this section is based on a work published by Ceceñas in Ref. [9] about a TH model developed for boiling water reactors. The TH model was modified from a point kinetics approach with an extension of the NK model to 3D and implemented in the

The treatment of neutron kinetics in [9] has been improved by coupling a 3D solution of the neutron diffusion equations with an arrangement of TH channels in parallel. Each channel independently contemplates three regions: (1) one phase, (2) subcooled boiling, and (3) bulk boiling. The objective was to implement a detailed model of a nuclear reactor core, which is somehow perturbed to simulate NK-TH coupling. These perturbations are obtained when the power generated in a group of channels changes and thus affecting the TH state of each channel. The original [9] TH model is based on a generic channel, which is adapted by transferring to it the operational data as flow area, generated power, axial power profile, and subcooling, among other parameters. Each channel is associated with a number of nuclear fuel assemblies and an axial power profile. Although the neutron model is a two-dimensional model for the radial power profile in each z-plane covering all the channels, information related to the axial power distribution is considered for each individual channel. In Ref. [9], it is assumed that this steady-state axial power profile is invariant over time, and it is used to weight the axial averages of macroscopic cross sections and void fractions. To perform the numerical implementation of the model, the arrangement of channels is obtained by grouping the total core assemblies into an appropriate number of thermal-hydraulic channels, which gives a definition

For the implementation in AZKIND of the TH model of Ceceñas, the grouping of fuel assemblies was maintained for generating a reduced number of TH channels; operational data are also used. The main difference is that the NK model recursively computes the axial power profile for each channel, and this thermal power is the updated source of power for TH model. Therefore, a "new" thermal-hydraulic condition is generated, and it is used by the NEMTAB model to update the nuclear data to generate new thermal power profiles with the NK model. The process is iterative, and it stops when the convergence is met. Convergence is achieved

The NK-TH coupling in AZKIND performs core calculations as described above to obtain a steady-state reactor core condition. For transient conditions in a time interval T, the NK-TH coupling process is the same for each time step ΔT in T, that is, a different quasi-steady-state

when updated conditions do not change in both NK and TH models.

5. Neutronic and thermal-hydraulic coupling model (NK-TH)

corresponding to the generated thermal power.

development of AZKIND.

14 New Trends in Nuclear Science

of a set of channels per quadrant.

The heat transfer and temperature distribution in the fuel and cladding can be calculated by a simple model where the heat diffusion equation is solved in one dimension (radial) for a fuel rod, since the conduction in axial direction is small compared to the radial one, it can be neglected. An energy balance per unit length yields

$$m\_f c\_{pf} \frac{dT\_f}{dt} = q'(t) - \frac{1}{R\_{\text{g}}'} \left[ T\_f(t) - T\_c(t) \right] \tag{15}$$

$$m\_{\varepsilon}c\_{pc}\frac{d\overline{T\_{\varepsilon}}}{dt} = \frac{1}{R\_{\text{g}}'} \left[\overline{T\_{f}}(t) - \overline{T\_{\varepsilon}}(t)\right] - \frac{1}{R\_{\text{c}}'} \left[\overline{T\_{\varepsilon}}(t) - \overline{T\_{m}}(t)\right] \tag{16}$$

where R<sup>0</sup> <sup>g</sup> and R<sup>0</sup> <sup>c</sup> represent thermal resistances per unit length. The coefficient of heat transfer to the refrigerant fluid is calculated by the Dittus-Boelter or Chen correlation, depending on the type of flow, which can be in one or two phases. These equations are used for the radial averaging of the temperatures in the fuel rod.

#### 5.2. Reactor coolant dynamics

The conservation equations of mass, energy, and momentum are applied in this case to a flow of water along a vertical channel, where the dynamics of the fluid heated by the wall of the fuel is modeled. Conservation equations can be expressed as [10]

$$\frac{\partial \rho\_{\rm m}}{\partial t} + \frac{\partial \mathbf{G}\_{\rm m}}{\partial z} = \mathbf{0} \tag{17}$$

$$\frac{\partial G\_m}{\partial t} + \frac{\partial}{\partial z} \left( \frac{G\_m^2}{\rho\_m^+} \right) = -\frac{\partial p}{\partial z} - \frac{\mathbf{f} G\_m |G\_m|}{2 D\_c \rho\_m} - \rho\_m \mathbf{g} \cos \Theta \tag{18}$$

$$\rho\_m \frac{\partial h\_m}{\partial t} + \mathbf{G}\_m \frac{\partial h\_m}{\partial z} = \frac{q'' P\_h}{A\_z} + \frac{\partial p}{\partial t} + \frac{\mathbf{G}\_m}{\rho\_m} \left( \frac{\partial p}{\partial z} + \frac{\mathbf{f} \mathbf{G}\_m |\mathbf{G}\_m|}{2 D\_e \rho\_m} \right) \tag{19}$$

In this work, the conservation equations are solved by the Integral Moment method [11], according to which it is assumed that the refrigerant is incompressible but thermally expandable, and the density is a function of enthalpy at a constant pressure

$$\frac{\partial \rho\_m}{\partial t} = \frac{\partial \rho\_m}{\partial h\_m}\bigg|\_p \frac{\partial h\_m}{\partial t} + \frac{\partial \rho\_m}{\partial p}\bigg|\_{h\_m} \frac{\partial p}{\partial t} = \mathbf{R}\_l \frac{\partial h\_m}{\partial t} + \mathbf{R}\_p \frac{\partial p}{\partial t} \tag{20}$$

Neglecting terms related to pressure changes and wall friction forces, the energy equation is simplified as

$$
\rho\_{\rm m} \frac{\partial \mathbf{h\_{m}}}{\partial t} + \mathbf{G\_{m}} \frac{\partial \mathbf{h\_{m}}}{\partial z} = \frac{\mathbf{q'} \mathbf{P\_{h}}}{\mathbf{A\_{z}}} \tag{21}
$$

where Gi is the flow rate for channel i, the index k represents the number of the iteration, w is an arbitrary weight to control the convergence, and P is the average pressure drop of all

> 1¼1 Pk

It is observed that even though the pressures are equaled, the value of the pressure drop in the core is not imposed as a boundary condition. Convergence is achieved when the following

iteration, the enthalpy and void fraction profiles are affected. It is necessary to recalculate the TH solution at each iteration for all channels, achieving convergence when every parameter

Although reference [12] has important issues to be considered in the development of an NK-TH-coupled model, those issues are not repeated here, but taken into account. The most direct way of coupling NK module and TH module, as implemented in AZKIND, consists simply in that axially both NK mesh and TH mesh have the same partition, making possible to assign an NK node at position z to the TH node in the same position. This relationship is a one-to-one

As it can be seen in Figure 4, before initiating the NK-TH feedback process, the initial nuclear parameters and kinetics parameter (XS) are loaded from files constructed in NEMTAB format, previously generated by means of a lattice code. Then, following the reading of the nuclear reactor burn-up state and thermo-physics initial conditions, the XS parameters are obtained

The process continues as follows. The corresponding neutron flux is calculated in the NK module with the mgcs numerical solver, and this power (initial neutron flux) is the heat source to be assigned to the TH model. The axial power profile can be that of each fuel assembly assigned to a unique TH channel or the power profile of a set of fuel assemblies assigned to a TH channel. The axial power profile is the heat source for each node in the z-direction. Once the axial power profiles have been constructed in the TH module, an initial thermal-hydraulic state of the reactor system is calculated. The thermal-hydraulic state is calculated for each node

The important variables sent to the NK module are the fuel temperature (Tf), moderator temperature (Tm), and moderator density (Dens). The XS parameters are updated using these 3D variables for interpolation in the NEMTAB tables. The next step is to calculate new 3D power profiles to be sent to the TH module. This cyclic NK-TH calculation continues and stops when the TH criterion and neutron-flux criterion are met. Stopping the cyclic calculation means that the reactor power and thermal-hydraulic conditions have reached a steady state.

from the Nemtab multi-dimensional tables by means of interpolation calculations.

in the TH channels from the bottom to the top of the reactor core.

<sup>i</sup> (25)

http://dx.doi.org/10.5772/intechopen.79723

Nuclear Reactor Simulation

17

� <sup>&</sup>lt; <sup>ε</sup> . By changing the flow rate of the channel for each

P k ¼ 1 N X N

channels at iteration k, obtained as

<sup>i</sup>¼<sup>1</sup> <sup>P</sup> k � Pk i

� � � � �

involved in the thermal-hydraulic calculation remains unchanged.

5.3. Neutron kinetics: thermal-hydraulics (NK-TH) coupling model

relationship is met: P<sup>N</sup>

node correspondence.

where the axial flow variation can be obtained by

$$\frac{\partial \mathbf{G\_m}}{\partial \mathbf{z}} = -\frac{\mathbf{R\_h}}{\rho\_\mathbf{m}} \left( \frac{\mathbf{q'} \mathbf{P\_h}}{\mathbf{A\_z}} - \mathbf{G\_m} \frac{\partial \mathbf{h\_m}}{\partial \mathbf{z}} \right) \tag{22}$$

This equation provides the flow variations with respect to an average value imposed as a boundary value or provided by the dynamics of the coolant recirculation system. Three regions are defined by which the coolant circulates as it ascends into the channel: a one-phase region, a subcooled boiling region, and a bulk boiling region. The first region begins at the bottom of the channel, where the coolant enters with known enthalpy and ends at the point of separation of the bubbles Zsc. The bulk temperature at this point is obtained by the Saha and Zuber correlation. The subcooled boiling region ends when the bulk temperature reaches the saturation temperature, and its axial location is determined by an energy balance. The enthalpy distribution allows the calculation of the thermodynamic equilibrium quality, used to calculate the flow quality. The axial distribution of the void fractions is calculated by iteratively solving the equation for void fraction α and the Bankoff correlation slip (S):

$$\alpha = \frac{\mathbf{x}}{S\left(\frac{\rho\_{\mathbf{s}}}{\rho\_{\uparrow}}\right) + \mathbf{x}\left(1 - S\left(\frac{\rho\_{\mathbf{s}}}{\rho\_{\uparrow}}\right)\right)}, S = \frac{1 - \alpha}{\mathbf{k}\_{\text{s}} - \alpha + (1 - \mathbf{k}\_{\text{s}})\alpha^{\mathbf{r}}},\tag{23}$$

where, in this case, the parameters ks and r are functions of the system pressure:

$$\mathbf{k}\_s = 0.71 + 1.2865 \times 10^{-3}p,\text{ and, } r = 3.33 - 2.56021 \times 10^{-3}p + 9.306 \times 10^{-5}p^2.$$

The total pressure drop in the channel is made up of the contributions of each region. Every term in each region includes the contribution by acceleration, gravity, and friction. For the channel arrangement, the steady state is obtained by iterating the coolant flow rate of each channel to obtain the same pressure drop for all of them. This iteration consists of a correction to the flow defined by the deviation of the pressure drop of the channel with respect to the average of all the channels:

$$\mathbf{G}\_{\rm i}^{\rm k+1} = \mathbf{G}\_{\rm i}^{\rm k} + w \mathbf{G}\_{\rm i}^{\rm k} \left( \frac{\overline{P}^{\rm k} - P\_{\rm i}^{\rm k}}{P\_{\rm i}^{\rm k}} \right) \tag{24}$$

where Gi is the flow rate for channel i, the index k represents the number of the iteration, w is an arbitrary weight to control the convergence, and P is the average pressure drop of all channels at iteration k, obtained as

$$\overline{P}^k = \frac{1}{N} \sum\_{1=1}^N P\_i^k \tag{25}$$

It is observed that even though the pressures are equaled, the value of the pressure drop in the core is not imposed as a boundary condition. Convergence is achieved when the following relationship is met: P<sup>N</sup> <sup>i</sup>¼<sup>1</sup> <sup>P</sup> k � Pk i � � � � � � <sup>&</sup>lt; <sup>ε</sup> . By changing the flow rate of the channel for each iteration, the enthalpy and void fraction profiles are affected. It is necessary to recalculate the TH solution at each iteration for all channels, achieving convergence when every parameter involved in the thermal-hydraulic calculation remains unchanged.

#### 5.3. Neutron kinetics: thermal-hydraulics (NK-TH) coupling model

In this work, the conservation equations are solved by the Integral Moment method [11], according to which it is assumed that the refrigerant is incompressible but thermally expand-

> ∂r<sup>m</sup> ∂p � � � � hm

Neglecting terms related to pressure changes and wall friction forces, the energy equation is

q 00 Ph Az

This equation provides the flow variations with respect to an average value imposed as a boundary value or provided by the dynamics of the coolant recirculation system. Three regions are defined by which the coolant circulates as it ascends into the channel: a one-phase region, a subcooled boiling region, and a bulk boiling region. The first region begins at the bottom of the channel, where the coolant enters with known enthalpy and ends at the point of separation of the bubbles Zsc. The bulk temperature at this point is obtained by the Saha and Zuber correlation. The subcooled boiling region ends when the bulk temperature reaches the saturation temperature, and its axial location is determined by an energy balance. The enthalpy distribution allows the calculation of the thermodynamic equilibrium quality, used to calculate the flow quality. The axial distribution of the void fractions is calculated by

iteratively solving the equation for void fraction α and the Bankoff correlation slip (S):

rf

The total pressure drop in the channel is made up of the contributions of each region. Every term in each region includes the contribution by acceleration, gravity, and friction. For the channel arrangement, the steady state is obtained by iterating the coolant flow rate of each channel to obtain the same pressure drop for all of them. This iteration consists of a correction to the flow defined by the deviation of the pressure drop of the channel with respect to the

> <sup>i</sup> <sup>þ</sup> wG<sup>k</sup> i P k � <sup>P</sup><sup>k</sup> i Pk i

!

� � � � , S <sup>¼</sup> <sup>1</sup> � <sup>α</sup>

<sup>þ</sup> <sup>x</sup> <sup>1</sup> � <sup>S</sup> <sup>r</sup><sup>g</sup>

G<sup>k</sup>þ<sup>1</sup> <sup>i</sup> <sup>¼</sup> <sup>G</sup><sup>k</sup>

where, in this case, the parameters ks and r are functions of the system pressure:

p, and, <sup>r</sup> = 3.33 � 2.56021 � <sup>10</sup>�<sup>3</sup>

∂hm <sup>∂</sup><sup>z</sup> <sup>¼</sup> <sup>q</sup> 00 Ph Az

� Gm

� �

∂hm ∂z

∂p <sup>∂</sup><sup>t</sup> <sup>¼</sup> <sup>R</sup><sup>h</sup> ∂hm <sup>∂</sup><sup>t</sup> <sup>þ</sup> <sup>R</sup><sup>p</sup> ∂p

<sup>k</sup><sup>s</sup> � <sup>α</sup> <sup>þ</sup> ð Þ <sup>1</sup> � <sup>k</sup><sup>s</sup> <sup>α</sup><sup>r</sup> , (23)

p2 .

<sup>p</sup> + 9.306 � <sup>10</sup>�<sup>5</sup>

<sup>∂</sup><sup>t</sup> (20)

(21)

(22)

(24)

able, and the density is a function of enthalpy at a constant pressure

rm ∂hm <sup>∂</sup><sup>t</sup> <sup>þ</sup> Gm

∂Gm <sup>∂</sup><sup>z</sup> ¼ � Rh rm

<sup>α</sup> <sup>¼</sup> <sup>x</sup> S <sup>r</sup><sup>g</sup> rf � �

ks = 0.71 + 1.2865 � <sup>10</sup>�<sup>3</sup>

average of all the channels:

∂hm ∂t þ

∂r<sup>m</sup> <sup>∂</sup><sup>t</sup> <sup>¼</sup> <sup>∂</sup>r<sup>m</sup> ∂hm � � � � p

where the axial flow variation can be obtained by

simplified as

16 New Trends in Nuclear Science

Although reference [12] has important issues to be considered in the development of an NK-TH-coupled model, those issues are not repeated here, but taken into account. The most direct way of coupling NK module and TH module, as implemented in AZKIND, consists simply in that axially both NK mesh and TH mesh have the same partition, making possible to assign an NK node at position z to the TH node in the same position. This relationship is a one-to-one node correspondence.

As it can be seen in Figure 4, before initiating the NK-TH feedback process, the initial nuclear parameters and kinetics parameter (XS) are loaded from files constructed in NEMTAB format, previously generated by means of a lattice code. Then, following the reading of the nuclear reactor burn-up state and thermo-physics initial conditions, the XS parameters are obtained from the Nemtab multi-dimensional tables by means of interpolation calculations.

The process continues as follows. The corresponding neutron flux is calculated in the NK module with the mgcs numerical solver, and this power (initial neutron flux) is the heat source to be assigned to the TH model. The axial power profile can be that of each fuel assembly assigned to a unique TH channel or the power profile of a set of fuel assemblies assigned to a TH channel. The axial power profile is the heat source for each node in the z-direction. Once the axial power profiles have been constructed in the TH module, an initial thermal-hydraulic state of the reactor system is calculated. The thermal-hydraulic state is calculated for each node in the TH channels from the bottom to the top of the reactor core.

The important variables sent to the NK module are the fuel temperature (Tf), moderator temperature (Tm), and moderator density (Dens). The XS parameters are updated using these 3D variables for interpolation in the NEMTAB tables. The next step is to calculate new 3D power profiles to be sent to the TH module. This cyclic NK-TH calculation continues and stops when the TH criterion and neutron-flux criterion are met. Stopping the cyclic calculation means that the reactor power and thermal-hydraulic conditions have reached a steady state.

assemblies of an LWR. Fine discretization means that each fuel assembly was subdivided in a mesh of size 10 10. As an example, an arrangement of 6 6 fuel assemblies consists of a square with 36 fine-discretized fuel assemblies. The corresponding algebraic system for each fuel arrangement was solved with parallel processing performed by the bicgstab solver mentioned earlier. In Tables 2 and 3, the speedup of the different cases is shown [1] with a remarkable performance. Despite the speedup for small matrices that is comparable for the three computer architectures used, it is also important to notice that the speedup values listed in Table 3 do not present a linear behavior, and the reason is because although more GPU processor cores are used with massive data transference to and from the GPU, a data traffic delay is present in the communication bus between the GPU and the CPU. For the analysis of the computing acceleration or "speedup," a definition of speedup is used in [15], known as relative speedup or speedup ratio: S = T1/Tn, where T1 is the computing time using a single processor (serial calculation) and Tn is the computing time using n processor cores. The "no memory" insert listed in Table 2 is because for those large matrix dimensions, there is not

Figure 5 [1] shows the distribution of nuclear fuel assemblies in the core of a boiling water reactor. Excepting the blue-shaded zone, colors are for different types of fuel assemblies. In the plane xy, the mesh is 24 24, according to each fuel zone, and axially, there are 25 nodes. The matrix for this coarse mesh (1,274,304 nnz) is comparable to the matrix of the fine mesh created

As described in [1], a reactor power transient was simulated as the capability to remove neutrons was highly increased in the perturbed assembly shown in Figure 5. An increase as

> 4 4 1,947,200 21,174,400

6 6 4,363,200 47,606,400 10 10 12,080,000 132,160,000

Nuclear Reactor Simulation

19

http://dx.doi.org/10.5772/intechopen.79723

2 2 492,800 5,305,600

Serial 24 124 372 994 2471 GTX 860 M 2.1 7.9 31.3 No memory No memory Tesla K20c 1.3 4.0 16.6 40.1 No memory

GTX TITAN X 1.0 2.6 10.4 26.7 95.4

Assemblies 1 1 2 2 4 4 6 6 10 10 GTX 860 M 11 16 12 – – Tesla K20c 18 31 22 24 – GTX TITAN X 24 48 36 37 26

enough memory to load the matrix and solvers.

Assemblies array: Matrix dimension (n): nnz elements:

Table 3. Speedup comparison (S) [1].

for the case of a unique assembly (case 1 1 listed in Table 2).

1 1 126,200 1,332,400

Table 2. Parallel processing time (seconds) in different architectures [1].

Figure 4. The NK-TH feedback process in AZKIND.
