Advanced pressurized water reactors are the main part of a new generation of nuclear power plant projects under development that provide cost-effective power production for various needs (Yemelyanov et al. 1982, Klimov 2002, Boyko et al. 2005, Baklushin 2011, Bays et al. 2019, Nuclear Technology Review 2019). The innovative technologies are aimed at improving the safety and reliability as well as at reducing the cost of NPPs. At the same time, improvements in design, technological and layout solutions are focused primarily on the reactor core. Assessments of the efficiency of these improvements are preceded by numerical simulations of the processes in the core, in particular heat generation and sink, with account for the difference between the study object and the standard version tested in operational practice.

The authors of the article propose a method for calculating the temperature field in the core of a heterogeneous reactor (using the example of a pressurized water reactor), which makes it possible to quickly assess the level of temperature safety of various changes in the core and has the necessary speed for analyzing transients in real time.

This method is based on the energy equation for an equivalent homogeneous core in the form of a heat equation that takes into account the main features of the simulated heterogeneous structure. The procedure for recovering the temperature field of a heterogeneous reactor uses the analytical relation obtained in this work for the heat sink function, taking into account inter-fuel element heat leakage losses.

Calculations of temperature fields in the model of the PWR type reactor (The Westinghouse Pressurized Water Reactor Nuclear Plant 1984) were carried out in stationary and transient operating modes. The calculation results were compared with the results of CFD simulation. The area of competing use of the temperature field recovery method was indicated.

temperature field, reactor core, thermal conductivity equation, heat sink function, fuel element, real time

Attention to pressurized water reactors is explained not only by their prevalence as an energy source in various types of nuclear power plants (NPPs), both stationary and transportable (including space ones), but also by the constant improvement of existing design, technological and layout solutions, primarily related to reactor core. The obvious unfeasibility of most full-scale tests aimed at confirming the improvement effect has led to the development of simulation bench tests of core element models. However, the experimental study, in any case, is preceded by numerical simulation.

All computing codes of heterogeneous cores are complicated, requiring significant computing power and imposing individual restrictions on the area of use. Another thing is a homogeneous core, which, moreover, is presented as a continuous medium. For such an core, the energy equation can be written in the form of the heat conduction equation (Kuzevanov et al. 2017, Kuzevanov and Podgorny 2019a, b). The numerical solution of this equation can be arranged as an efficient procedure with a small computation time at the computing power available to any user.

In this work, the task is to build an algorithm for recovering the temperature field in a heterogeneous core in stationary and transient processes based on the results of calculating the temperature in an equivalent homogeneous core of a nuclear reactor.

For definiteness and concretization of the relations obtained below, let us consider the core of a PWR type reactor (The Westinghouse Pressurized Water Reactor Nuclear Plant 1984) loaded with fuel assemblies that are square in cross section. The fuel elements are located in the fuel assemblies at the nodes of a square lattice with a step a1. The fuel assemblies have no jackets. The geometric parameters and composition of the fuel elements are presented in Table 1.

• The pressure of the coolant at the core inlet P _in , its temperature T ⁱⁿ₂ and its mass flow rate G are known; the coolant pressure at the core outlet P _out is the same for all the fuel assemblies.
• The coolant exchange between neighboring fuel assemblies occurs on an equal-mass basis. The mass flow rate of the coolant in any fuel assembly does not change along its height z.
• Let us define a virtual channel as an elementary FA channel, in the center of the square cross-section (with side a₁) of which there is a single fuel element with an outer cladding diameter d _sh .
• Volumetric heat releases q _v in each of the fuel elements and mass flow rates of the coolant in the elementary channels of any fuel assembly are the same. The lateral boundaries of the channels are permeable. The connection of hydraulic losses with the mass flow rate of the coolant in any fuel assembly or elementary channel can be represented by the Darcy-Weisbach equation.

• In the calculations of changes in the temperature field T in time τ, the core as a whole is represented as a continuous medium. At the same time, in the numerical calculation, the elementary volume coincides with the elementary channel as its part with the size Δ z along the axial coordinate z. The parameters of the continuous medium are weighted average. Each calculated volume has heat sources and sinks.
• In the analysis of coolant temperature changes and determination of the sink term power, the core is considered as heterogeneous, composed of elementary channels, including discrete components “ i”. As components of the heterogeneous calculated volume, we consider the coolant (water, any coolant parameter index is “2”), fuel rod cladding (metal, the parameter index is “ sh”), gas gap between the cladding and fuel in the fuel rod (gas, the parameter index is “g”), and fuel (uranium dioxide, the parameter index is “ f”).
• A single elementary channel is a model of the corresponding fuel assembly in the analysis of the temperature field of the reactor core.
• Heat and mass transfer between neighboring equivalent channels is described in the quasi-stationary approximation.

The core is described as a continuous medium with weighted average temperatures in the calculated cells. The heat conduction equation in the form (Kuzevanov et al. 2017, Kuzevanov and Podgorny 2019a, b) is used as the energy conservation equation:

$\Phi \frac{\partial T}{\partial \tau }=div\left({\lambda }_{ef}gradT\right)+q_v-q_{v,st}$ (1)

where Φ = Σε_i; λ_ef is the effective thermal conductivity of the simulated system W/(m∙°C); q_v is the specific volumetric energy release, W/m³; q_v._st is the heat sink function reflecting the heat removal from the fuel rod surface per unit volume, W/m³; weighted average temperature in the calculated cell of the calculated volume T = ΣT_i∙ε^*_i, ε^*_i = ε_i/Φ; ε_i = (ρ∙c∙φ)_i, where ρ, c and φ are the are the density, specific heat capacity and volume fraction of the heterogeneous core component in the calculated volume, respectively.

The analysis showed that for the core of a pressurized water reactor, it is possible to use the relation for determining the heat removal from the surface of a fuel element in a non-stationary process, obtained in (Kuzevanov and Podgorny 2019) for the core of a high-temperature gas-cooled reactor:

q_v.st = σ{q⁰(z, r)f ^~ + P^~}, (2)

where σ = F/V; F is the heat exchange surface of a single fuel element, m²; V is the volume of the calculation cell, m³; q⁰ is the heat flux density per unit area of the fuel element, W/m² (index “0” is the stationary (initial) state); f^~ is the a function of time, spatial coordinates (z, r) and disturbing effects; P^~ is the function reflecting the influence of boundary conditions.

One of the features of using relation (2) is the need to determine the true heating of the coolant in each elementary channel, taking into account thermal leakage between them, which is possible only if the design and hydrodynamic features of the heterogeneous core are considered.

With the adopted physical model of the core taken into account, the heat balance equation for the channel “j” of the core with square fuel assemblies in the stationary mode looks like this:

(3)

where Δh_j = c₂,_jΔT_j, h is the specific enthalpy of the coolant, J/kg; 〈q_n,j〉 is the average density of the heat flux from other fuel assemblies, W/m²; a₁ is the cross-size of a cell containing one fuel element, m; H is the core height, m; ΔТ is the coolant heating in the channel, °C.

For an arbitrary cross section z of the channel j, we define the heat flux density q_n,j as q_n = q_n,j.1 – q_n,j.2. Assume that the components of the heat flux q_n,j.1 and q_n,j.2 through the virtual side surface of the channel can be represented as:

$q_{n,j.1}={\lambda }_{ef.j-1}\frac{T_{2,j-1}-T_{2,j}}{a_2}$ , $q_{n,j.2}={\lambda }_{ef.j}\frac{T_{2,j}-T_{2,j+1}}{a_2}$ (4)

where λ_ef is the effective heat transfer coefficient averaged over the channel height, taking into account the molecular and turbulent components of heat transfer, W/(m∙°C); а₂ is the effective distance between adjacent fuel assemblies, m.

And now we shall determine the effective thermal conductivity coefficient λ_ef based on the following considerations. Let us assume that the change in the intensity of heat transfer on the heating surface during the transition from the laminar flow regime to the turbulent one is directly related to the general change in the heat-conducting properties of the medium. Then we obtain:

α_l = α_t∙λ/λ_ef, (5)

where α and λ are the are the coefficients of heat transfer and molecular thermal conductivity, respectively (the indexes “l” and “t” refer to the laminar and turbulent coolant flow regimes, respectively).

For the laminar regime on the stabilized section in the round pipe, the solution of the integral Lyon relation for a laminar fluid flow leads to the equality Nu_l = A = const.

Extending the relation Nu_l = A to channels of arbitrary shape with an equivalent diameter d_e, we obtain from (5) the following expression for the effective thermal conductivity:

λ_ef = α_t∙d_e/A. (6)

The transformation of equation (3) into a system of algebraic equations for the connection of flow, hydraulic and thermodynamic parameters of channels, convenient for analysis, was carried out using B. Petukhov’s formula for calculating Nu_t (Petukhov and Kirillov 1958). If we combine as identical (i.e. having the same heat load and coolant flow rate) fuel assemblies under the same channel group number “j”, designating the number of such groups as “m”, then, in the case of numbering from the center of the core, the transformed system of equations will be look like this:

$G_j{c}_{p.j}\Delta T_j\left(1+0.5C_1{\xi }_1\left({\delta }_1+{\delta }_2{Y}_j\right)\right)=$

$=v_j\frac{Q_p}{N·n^2}+0.5G_j{c}_{p.j}{C}_1{\xi }_j\left[{\delta }_2{Y}_j\Delta T_{j-1}+{\delta }_1\Delta T_{j+1}\right]$ ; (7)

1 ≥ j ≥ m

Here ξ are friction resistance coefficients;

${\delta }_1=\left\{\begin{array}{l}0,\text{ if }j=m\\ 1,\text{ if }j\ne m\\ 2,\text{ if }j=1\end{array}\right.$ ${\delta }_2=\left\{\begin{array}{l}0,\text{ if }j=1\\ 1,\text{ if }j\ne 1\end{array}\right.$ (8)

Y_j = ξ_m–1G_m–1/ξ_mG_m; С₁ = (Hb^*)/L; L = Aa^*₂/a₁ is the heat exchange constant between fuel assemblies; a^*₂ = a₂n; b^* = C_t /2πd[K+ε(Pr)]; K and ε(Pr) are the temperature correction and the coefficients of Petukhov’s formula for Nut (Petukhov and Kirillov 1958); n² is the number of fuel elements in a square fuel assembly; N is the number of fuel assemblies.

The system of equations (7) is supplemented by a system of equations for the pressure drop in a group of identical fuel assemblies in the form of the Darcy-Weisbach equations (Kuzevanov and Podgorny 2018).

For the system of interconnected channels, it is proposed to determine q⁰_v.st as follows:

q ⁰ _v.st = k (T⁰ – T⁰₂), (9)

where for a square lattice of fuel elements arrangement:

$k=\frac{\pi }{a_1^2}{\left\{\left(\frac{R_{l1}}{2}+R_{l2}\right){\epsilon }_{sh}^{*}+\left[R_{l1}+R_{l2}+\frac{1}{8{\lambda }_f}\right]{\epsilon }_f^{*}\right\}}^{-1}$ (10)

Note that in relation (10) R_l1 and R_l2 are the linear thermal resistance of the cladding (including the gas gap) and heat transfer, respectively.

Using expression (9) in the equation

div(λ_ef grad T) + q⁰_v – k (T – T⁰₂) = 0 (11)

together with the system of equations (7) makes it easy to determine the stationary temperature distribution in the reactor core.

The dimensionless time function f_τ is included in the relations for determining f^~ and P^~ (2) (Kuzevanov and Podgorny 2019, 2019a). In the case of high thermal conductivity of a heterogeneous system (k/Φ >> 1) with known stationary fields T⁰(τ = 0), T^∞(τ → ∞ after exposure to disturbing factors), f_τ = f^*_τ = (T – T⁰)/(T^∞– T⁰), which is confirmed by calculations of the temperature fields of a gas-cooled reactor (Kuzevanov and Podgorny 2019, 2019a).

Analytical and computational studies conducted by the authors have shown the possibility of using the following dependencies when calculating the function f_τ for cores of PWR reactors:

$f_{\tau }=\left\{\begin{array}{l}{f}_{\tau }^{*}\text{, if }{k}_T=1;k_G=1\\ A_1^{-1}{f}_{\tau }^{*}+1-A_1^{-1}-\left(1-A_1^{-1}-\delta f_{\tau }\right)\exp \left(-A_2\tau /B\right),\text{ if }{k}_T\ne 1\text{ and }{k}_G\ne 1\end{array}\right.$ (12)

The following notations are used here:

$\left.\begin{array}{c}\delta f_{\tau }=\left(1-{\epsilon }_f^{*}\right)\left(\delta T_2-{\Delta }_2\right)/\delta T;\delta T_2=T_2^{\infty }-T_2^0;\delta T=T^{\infty }-T^0;\\ B=\delta T\left[1-\left(1-{\epsilon }_f^{*}\right){\Delta }_2/\delta T\right];A_2=\left(\delta T_{sh}-{\Delta }_2\right)k_{\alpha }\sigma {\alpha }_0/\Phi ;\\ \delta T_{sh}={\left(T_{sh}^{’}\right)}^{\infty }-{\left(T_{sh}^{’}\right)}^0;{\Delta }_2=\delta T_2-\left(k_T-\frac{1}{k_G}\right){\left(T_2^0\right)}^{in}+\left(1-\frac{1}{k_G}\right)T_2^0;\\ A_1=\left[\left(1-{\epsilon }_f^{*}\right){\Delta }_2+A_2\Phi k^{-1}\right]/\delta T;\end{array}\right\}$ (13)

k_i is the ratio of the new stationary values of the disturbing parameters to the initial ones; T ′_sh is the temperature of the outer surface of the fuel element cladding.

The results of calculating the weighted average temperature when the core is represented as an equivalent homogeneous medium formed the basis for the procedure for recovering the temperature field in the elements of any calculation cell, i.e., the coolant, cladding and fuel. In the coolant, the temperature field was not detailed; only its average temperature in the cross section of the elementary channel and the equality of coolant and cladding temperatures on the outer surface of the fuel element were taken into account. It was assumed that the temperature profile in the fuel cladding remains logarithmic, while in the fuel it was described by a power function during the entire transient process. Within these model approximations, the procedure for recovering the temperature field in any calculated cell of the core looks quite simple.

Indeed, at the time τ after the start of the transient process in the core, caused by an abrupt change in any of the parameters or several parameters that affect the temperature distribution in the core, the following fields are directly known as a result of calculating the equivalent homogeneous core:

• average temperatures Т₂;
• values of the heat sink function q _v.st ; and
• weighted average temperature Т.

Neglecting the thermal inertia of the thin cladding, we additionally calculate the temperature on the outer T ′_sh and inner T ″_sh surfaces of the cladding:

T ′_sh = T₂ + q_v.st∙a₁²∙R_l2/π; T ″_sh = T ′_sh + q_v.st∙a₁²∙R_l1/π (14)

and find the average value T_sh. The average value of the fuel temperature Тf is found from the determination of the weighted average temperature with known Т₂, Т_sh and Т. We consider the quadratic function as approximating the temperature profile in the fuel.

The maximum temperature value in the fuel of the calculation cell T_f max is determined according to the dependence:

T_f ^max = 2T_f – T ″_sh. (15)

We considered a calculation version of the core model, which consists of m groups of elementary cells, which were identical, square in cross section, with a size of а₁ (Podgorny and Kuzevanov 2020).

In terms of thermal, structural, flow and temperature characteristics, the calculated core corresponds to the PWR core (The Westinghouse Pressurized Water Reactor Nuclear Plant 1984), except for the radial dimensions: one elementary channel in the model fuel assembly (n² = 1) of n² = 289 identical elementary channels in a real fuel assembly.

Figs 1, 2 show the results calculating changes in the weighted average temperature fields after step-like disturbances were introduced into the stationary operation of the reactor. We studied the complex disturbance and the disturbance in one parameter, i.e., in terms of the reactor thermal power. The complex disturbance manifested itself as a simultaneous change in the reactor power, core inlet temperature and coolant flow rate, as well as the heat transfer coefficient, consistent with the change in the named parameters. The distribution of the weighted average temperature in the core model was calculated by the method of establishing the implicit solution to equation (1) using the system of equations (7) to determine the coolant heating in the elementary channels (Podgorny and Kuzevanov 2020).

Figure 1. 

Comparison of the average core temperature with thermal power surges by 50%, coolant mass flow rate by 20%, and coolant temperature at the core inlet by 20%: 1 – CFD simulation; 2 – developed algorithm

Figure 2. 

Comparison of the average temperature in the core with a thermal power surge by 50%: 1 – CFD simulation; 2 – developed algorithm

Figs 3, 4 show the recovered temperature fields for the diametrical cross section of fuel elements in the calculated cells with the indicated relative coordinates of the cell centers for different times of the transient process. For comparison, the same figures show the temperature values obtained as a result of CFD simulation of the transient process (Shaw 1992, Mohammadi and Pironneau 1994, Petrila and Trif 2005, ANSYS Fluent 2016a, b, c, Anderson et al. 2009). Since the maximum difference between the compared temperature values at similar points of the cross-section of the elementary cells does not exceed 50 °C, the proposed method for recovering the temperature field of a heterogeneous reactor can be recognized as a completely satisfactory approach to describing temperature changes in the elements of the core of a nuclear reactor in transient processes.

Figure 3. 

Comparison of the recovered temperature fields of the fuel elements with the results of CFD simulation for a time of 5 seconds from the beginning of the transient process shown in Fig. 1: 1 – fuel temperature of the central fuel element (in the central plane of the core); 2 – temperature of the central fuel element cladding (in the central plane of the core); 4 – fuel temperature of the fuel element of Channel 4 at a distance of 0.976 m below the central plane of the core; 5 – temperature of the fuel element cladding of Channel 4 at a distance of 0.976 m below the central plane of the core; 1, 2, 4, 5 – developed algorithm; 3, 6 – CFD simulation

Figure 4. 

Comparison of the recovered temperature fields of the fuel elements with the results of CFD simulation for a time of 10 seconds from the beginning of the transient process shown in Fig. 1: 1 – fuel temperature of the central fuel element (in the central plane of the core); 2 – temperature of the central fuel element cladding (in the central plane of the core); 4 – fuel temperature of the fuel element of Channel 4 at a distance of 1.22 m above the central plane of the core; 5 – temperature of the fuel element cladding of Channel 4 at a distance of 1.22 m above the central plane of the core; 1, 2, 4, 5 – developed algorithm; 3, 6 – CFD simulation

The proposed method for recovering the temperature field of a heterogeneous reactor does not claim to increase the level of detail of the temperature distribution in the core components in comparison with the resulting description of the temperature field using CFD simulation. However, in some cases, the authors’ approach described in this paper can be useful, since it has the following advantages:

• mobile availability of the computing power required for the calculations;
• short time for the complete calculation of the temperature distribution in the local region of the core of interest to the researcher, i.e., approximately two orders of magnitude less than when the basic CFD simulation algorithms are used; and
• operational preliminary calculation of a set of options for structurally different cores to select a limited number of them for the purpose of subsequent refining analysis.

Note that the time of calculating the temperature field by the recovery method is less than the time of the transient process. In this case, such a computational procedure can be an element of a complex program that describes the dynamics of the reactor circuit, e.g., in the software package of a nuclear power plant simulator. In addition, it is possible to use the recovery algorithm in the control systems of NPPs to correct the control based on the forecast of changes in the temperature field.