Positive effect of profiling the gas-cooled reactor core within the framework of the GT-MHR project was investigated in (Podgorny and Kuzevanov 2017, Kuzevanov and Podgorny 2017, 2018). The necessity arises to supplement already implemented analysis of equilibrium conditions of core operation with investigation of effects of profiling on the temperature field in transient modes of reactor core operation.

The present paper is dedicated to the investigation of development of transients in gas-cooled nuclear reactor core subject to the implementation of different principles of core profiling.

Investigation of transients in reactor core represents complex problem, solution of which by conducting direct measurements is beyond the resources available to the authors. Besides the above, numerical simulation based on advanced CFD software complexes (ANSYS 2016, 2016a, 2016b, Shaw 1992, Anderson et al. 2009, Petrila and Trif 2005, Mohammadi and Pironneau 1994) is also fairly demanding in terms of required computer resources.

The algorithm for calculating temperature fields using the model where the reactor core is represented as the solid medium with gas voids was developed by the authors and the assumption was made that heat transfer due to molecular heat conductivity can be described by thermal conductivity equation written for continuous medium with thermal physics parameters equivalent to respective parameters of porous object in order to get the possibility of obtaining prompt solutions of this type of problems.

Computer code for calculating temperature field in gas-cooled reactor in transient operation modes was developed based on the suggested algorithm. Proprietary computation code was verified by comparing the results of numerous calculations with results of CFD-modeling of respective transients in the object imitating the core of gas-cooled nuclear reactor. The advantage of the developed computer code is the possibility of real-time calculation of evolution of conditions in complex configurations of gas-cooled reactor cores with different channel diameters. This allows using the computer code in the calculations of transients in loops of reactor facility as a whole, in particular for developing reactor simulators.

Results are provided of calculations of transients for reactor core imitating the core of gas-cooled nuclear reactor within the framework of GT-MHR project performed for different approaches to profiling coolant mass flow. Results of calculations unambiguously indicate the significant difference of temperature regimes during transients in the reactor core with and without profiling and undeniable enhancement of reliability of nuclear reactor (Design of the Reactor Core 2005, International Safeguards 2014, Safety of Nuclear Power Plants 2014) with profiling of coolant mass flow in the reactor core as a whole.

Reactor core profiling, transient processes, temperature field, porous body, thermal conductivity equation, heat sink

Distributions of flowrate and thermal parameters of coolant inside cooling channels of high-temperature helium-cooled nuclear reactor were investigated in details in (Podgorny and Kuzevanov 2017, Kuzevanov and Podgorny 2017, 2018). Analysis was performed of different principles of core profiling intended to maintain constant the following parameters:

• Coolant mass flowrates in parallel channels;
• Coolant heating in the channels;
• Maximum wall temperatures of cooling channels.

Reactor operation on nominal power in steady-state operation mode was addressed.

In the event of perturbation with variation of either neutronics parameters (neutron field) or coolant parameters (flowrate, coolant temperature at the reactor core inlet) transient process develops with establishment of new temperature field in the reactor core. The present study is dedicated to the investigation of effects of coolant mass flowrate profiling on the variation of temperature in the reactor core during transients.

Physical and mathematical models of transients in gas-cooled nuclear reactor are formulated, calculation results are presented and their comparison with results of CFD-modeling is made using the example of GT-MHR nuclear reactor core.

Nuclear реактор within the framework of GT-MHR reactor design project (GT-MHR 2002, Vasyaev et al. 2001, Neylan et al. 1994, Engle 1977, Engle and Johnson 1976) is the pressurized graphite-gas high-temperature reactor with helium coolant and thermal power of the core equal to 600 MW. Reactor core is shaped as a hexagonal collar and consists of fuel assemblies in the form of perforated prismatic hexagonal blocks. Assemblies of the following three types are used including fuel, control and emergency shutdown blocks. Nuclear fuel is manufactured in the form of spherical fuel elements packed in fuel channels with 12.7-mm diameter. Cooling channels are represented by the following two groups: 9984 channels with 15.88-mm diameter and 642 channels with 12.7-mm diameter. Helium mass flowrate through the core is equal to G₀ = 320 kg/s in nominal power operation mode, coolant temperature at the core inlet Т_in = 491^оC.

In the model representation the core is examined as the solid medium (index 1) with gas voids (index 2). Solid part is multi-component consisting of fuel, graphite and metal; single-component gaseous part consists of helium.

It is assumed that heat transfer due to molecular thermal conductivity can be described as for the continuous medium with thermal physics parameters equivalent to respective parameters of the porous object.

The following is accounted for in the model:

• Thermal conductivity with varying density ρ and varying thermal conductivity coefficient λ of reactor core components;
• Heat generated by nuclear reactions is represented as distributed internal heat sources q _v ;
• Heat dissipated by coolant from the cooled down surface is represented as internal heat sinks q _v.t .
• The following assumptions are made:
• Averaged temperature Т within the calculation section (cell) is determined by the conditions of equivalent continuous medium;
• Temperatures of solid and gaseous components differ from each other but, however, are linked with the calculated average temperature of the equivalent continuous medium;
• Variation of averaged parameters of the continuous medium in tangential direction can be neglected; in this case reactor core is represented in two-dimensional approximation;
• Average porosity (ε) of the equivalent continuous medium within the reactor core volume is equal to the ratio of the volume of gaseous component V₂ to the total volume V of the core;
• Local porosity ε is the function of coordinate;
• Variation of porosity due to temperature effects can be neglected;
• Conditions of gaseous component are described by equation written for ideal gases;
• Pressure differentials P in the system are small – thermodynamic process during gas flow can be regarded as isobaric in the calculation of thermal physics parameters of the gas.
• Basic elements in the model of heat transfer are following:
• Thermal conductivity equation in the system with variable thermal physics parameters is the main energy balance equation;
• Displacement (flow) of the gaseous component is taken into account as the transfer of heat represented in the form of source (sink) member in the thermal conductivity equation;
• Due to the significant thermal inertia of the solid component variation of parameters of the gaseous component is calculated in quasi-stationary approximation.

Thus, the aggregated heat and mass transfer in the core of the nuclear reactor is represented in the model as the superposition of the non-stationary process of thermal conductivity and the quasi-stationary process of transfer (removal) of heat by the variable gas flow.

Let us write the thermal conductivity equation for the reactor core represented as the continuous medium with internal heat sources and sinks and with variable thermal physics properties in the following form (Kuzevanov et al. 2017)

$\genfrac{}{}{0.1ex}{}{\partial }{\partial \tau }\left(\rho h\right)=\nabla ({\lambda }_{ef}\nabla T)+q_v-q_{v,t}$ (1)

where T is the temperature, K; ρh = ρ₁h₁(1 – ε) + ρ₂h₂ ε; h₁, h₂ are the enthalpies of the solid and gaseous components, respectively, J/K; ρ₁ is the density of the multi-component solid fraction, kg/m³; ρ₂ is the helium density, kg/m³; τ is the time, s; q_v, q_v.t are the intensities of internal heat sources and sinks per unit volume of the reactor core, respectively, W/m³; λ_ef is the effective thermal conductivity coefficient for the equivalent medium, W/(m∙Κ) (definition of λ_ef corresponds to the model suggested in (Kuzevanov et al. 2017)).

After transformation of the left side of Equation (1) we obtain

$\genfrac{}{}{0.1ex}{}{\partial }{\partial \tau }\left(\rho h\right)=\Phi \genfrac{}{}{0.1ex}{}{\partial T}{\partial \tau }$ , $\Phi =(1-\epsilon )[c_1{\rho }_1+h_1\genfrac{}{}{0.1ex}{}{\partial {\rho }_1}{\partial T_1}]+[c_{p2}{\rho }_2+h_2\genfrac{}{}{0.1ex}{}{\partial {\rho }_2}{\partial T_2}]$

where the component c_p.2ρ₂ + h₂∙δρ₂/δT₂ = 0 in the case of zero value of h₂ at the point T = 0 K; c₁ is the average specific thermal capacity of the multi-component solid fraction within the temperature range under study, J/(kg∙Κ); c_p.2 is the isobaric specific thermal capacity of helium (c_p.2 = const), J/(kg∙Κ).

Differential equation

$\Phi \genfrac{}{}{0.1ex}{}{\partial \mathrm{T}}{\partial \tau }=\nabla ({\lambda }_{ef}\nabla \mathrm{T})+{\mathrm{q}}_v-{\mathrm{q}}_{v.t}$ (2)

Describes the heterogenous nuclear reactor core cooled with gas coolant as the equivalent core with distributed solid and gaseous components without internal boundaries, within the volume of which heat sources and sinks are positioned.

Quite naturally the function q_v.t must reflect specific features of the energy balance for the solid and gaseous components of the core. Let us further determine this function by binding it with the heat flow density on the cooling channel wall and with specific features of heterogenous structure of the core.

It is evident that the following condition is valid in the calculation cell i with finite dimensions assigned in the numerical solution of Equation (2):

(q_v.t)_i = σ_iq_i.d, (3)

where σ_i = F_i/V_i; V_i, F_i are the cell volume, m³ and the heat exchange surface corresponding to this volume, m², respectively; q_i.d is the heat flux density per unit surface, W/m² (additional index d stands for the effective value of the parameter in the non-stationary process).

The task of determination of the heat sink function is reduced to the task of determination of effective thermal flux density on the walls of cooling channels in non-stationary process.

In steady-state mode q_d = q_n, where q_n refers to steady-state conditions and is easily determined using well-known expressions, in particular the general energy balance equation:

${\displaystyle \underset{v}{\int }}{q}_{v}dv={\displaystyle \underset{F}{\int }}{q}_{n}dF+{\mathrm{Q}}_{boun}$ (4)

where F is the total heat exchange surface in the reactor core, m²; Q_boun is the total heat flux on the external boundary of the reactor core, W.

During the transient q_d ≠ q_n.

Let us write for arbitrary cooling channel of the reactor core the equation binding the coolant temperature T₂._d averaged over the flow cross-section with normal heat flux density on the channel wall q_d:

${\mathrm{T}}_{\mathrm{2d}}={\mathrm{T}}_2^{in}+\genfrac{}{}{0.1ex}{}{1}{G\cdot c_{p.2}}{\displaystyle \underset{5}{\overset{\epsilon }{\int }}}\pi d{q}_{d}dz$ (5)

where G is the coolant mass flowrate, kg/s; G ≠ G (z) in quasi-stationary approximation; d is the channel diameter, m; index in indicates the value of the parameter at the channel inlet; z is the axial coordinate, m; δ is the effective additive component, m.

Following the above made assumptions the equation below is satisfied

q _d = α(T_st.d – T_2.d), (6)

where α is the average heat transfer coefficient along the channel length for the case when gaseous coolant flows inside round channel with diameter d and coolant mass flowrate G, W/m²∙K; T_st.d is the effective value of the channel wall temperature averaged over the channel perimeter at the point with respective coordinate z of the channel cross-section.

Differential equation for the determination of q_d follows from Equation (5) taking expression (6) into account:

$\genfrac{}{}{0.1ex}{}{d{q}_{\mathrm{d}}}{dz}+a_1\cdot q_{\mathrm{d}}=\alpha \genfrac{}{}{0.1ex}{}{d{T}_{st.d}}{dz}$ (7)

a ₁ = πdα / Gc_p.2 is the parameter constant along the length of the channel in question following the assumptions made.

The function having the following form is the solution of Equation (7):

${\mathrm{q}}_{\mathrm{d}}=e^{-a_{1}z}\{\alpha {\displaystyle \int }\left(\genfrac{}{}{0.1ex}{}{d{T}_{st.\mathrm{d}}}{dz}\right)e^{a_{1}z}dz+a_2\}$, (8)

where a₂ is the integration constant determined by the conditions at the coolant inlet in the reactor core.

For determining q_d using the relation (8) reflecting the connection with channel wall temperature T_st.d it is necessary to know the functional dependence of this temperature on the coordinate and time.

The assumption of validity of the following relation was made within the framework of the suggested model

(T_st. – T_st.0) / (T_st.1 – T_st.0) = f_τ. (9)

Here index 0 refers to the stationary (initial) mode of operation of the reactor core before the development of the perturbation (variation) of parameters influencing the temperature regime; T_st.1 is the channel wall temperature averaged over the perimeter in newly established steady-state mode with modified parameters which affected the temperature regime; f_τ is the function of time and radial coordinate not dependent on the coordinate z.

Assumption (9) is made based on the following considerations. Firstly, the reactor core structure does not change along any cooling channel (along the axial coordinate z). Secondly, quasi-stationary approximation is used in cases when there exists lagging of the perturbation factor as refers to flowrate and thermal conductivity coefficient along the direction of coolant flow.

Solution of Equation (8) taking into account the expression (9) allows finding the form of the function q_d = q_d(z,τ) for any abrupt step-like perturbation.

Let us include q, G, α ανδ T₂ⁱⁿ in the list of perturbation factors influencing the temperature regime of the reactor core. Let us denominate the ratio of the new steady-state values of these parameters to the initial ones as follows:

k_q = q₁/q₀; k_α = α₁/α₀; k_G = G₁/G₀; k_T = T_2.1ⁱⁿ/T_2.0ⁱⁿ. (10)

It follows from expression (9) that

T _st.d = T_st.0 + (T_st.1 – T_st.0) f_τ. (11)

Let us take into consideration that stationary values T_st.0 and T_st.1 are determined as

$T_{st.0}=T_{2.0}^{in}+\genfrac{}{}{0.1ex}{}{q_0}{{\alpha }_0}+{\Phi }_0\left(z\right)$,

$T_{stl}=k_T{T}_{20}^{in}+\genfrac{}{}{0.1ex}{}{k_q}{k_n}\genfrac{}{}{0.1ex}{}{q_0}{{\alpha }_0}+\genfrac{}{}{0.1ex}{}{k_q}{k_G}{\Phi }_0\left(z\right)$, (12)

${\Phi }_0\left(z\right)=\genfrac{}{}{0.1ex}{}{\pi d}{G_0{c}_{p.2}}{\displaystyle \underset{\delta }{\overset{z}{\int }}}{q}_{0}dz$

Let us accept the “cosine” energy release dependence along the height of the core. If the coordinate origin is defined at the distance δ from the reactor core inlet we will have the following form of function q₀(z) for the steady-state mode:

q ₀(z) = q₀^max sin (πz/H_e), (13)

where q₀^max is the heat flux density on the wall channel in the middle of the channel along its height in the initial steady-state process; H_e = H + 2δ; H is the reactor core height, m.

Let us note that with heat release density varying over the reactor core cross-section heat flux density q₀ in the steady-state process is the function of two coordinates:

q ₀ = q₀(z,r) = q₀^max(r) sin (πz/H_e).

We integrate the first term of summation braced in figure brackets in Equation (8) and determine the integration constant a₂ from the condition of satisfaction of the following equation at the coolant inlet in the reactor core:

q ₀(δ,r) = α₀[T_st.0(δ,r) – T_2.0ⁱⁿ]. (14)

After not complicated transformations we obtain the following form of function q_d(z,r,τ):

$q_d(z,r;\tau )=q_0(z,r)\cdot \tilde{f}+\tilde{P}$(15)

where

$\tilde{f}=A(1-f_{\tau })(k_{\alpha }-k_G)[\genfrac{}{}{0.1ex}{}{\pi }{H_{\mathrm{e}}}\mathrm{ctg}\left(\pi \genfrac{}{}{0.1ex}{}{z}{H_{\mathrm{e}}}\right)-a_1]+(1-f_{\tau })k_{\alpha }+f_{\tau }{k}_q$,

$\tilde{P}=[k_{\alpha }{\alpha }_0(1-f_{\tau })(1-k_T)T_{\mathrm{2,0}}^{in}-\psi ]\exp [-a_1(z-\delta )]$. (16)

In their turn,

$a_1=\genfrac{}{}{0.1ex}{}{k_{\alpha }}{k_G}\genfrac{}{}{0.1ex}{}{\pi d}{G_0{c}_{p.2}}{\alpha }_0$;

$\begin{array}{c}\psi =\psi (\delta ,r,\tau )=q_0(\delta ,r)\left\{A(1-f_{\tau })(k_{\alpha }-k_G)\right\}\times \hfill \\ \times [\genfrac{}{}{0.1ex}{}{\pi }{H_{\mathrm{e}}}\mathrm{ctg}\left(\pi \genfrac{}{}{0.1ex}{}{\delta }{H_{\mathrm{e}}}\right)-a_1]\hfill \end{array}$ (17)

Now it is easy to write down the expression for the determination of the heat sink function during the transient for calculation cell i containing arbitrary number of channels of different diameter:

${\left(q_{v,t}\right)}_i={\displaystyle \sum _{j=1}^n}{\sigma }_{i,j}{\left({\mathrm{q}}_{j,\mathrm{d}}\right)}_i$ (18)

where n is the number of channels included in the calculation cell of the reactor core; σ_i,_j = F_i,_j /V_i; F_i,_j is the heat exchange surface of individual channel j belonging to calculation cell i; V_i is the total volume of the calculation cell.

According to its physical meaning the function f_τ reflects the local thermal inertia of the reactor core.

We assume that function f_τ cannot be expressed analytically in the general case and determine the function in arbitrary point as

f _τ = (T – T₀) / (T₁ – T₀), (19)

where Т is the solution of thermal conductivity equation during the transient; Т₀ and Т₁ are the respective solutions of the stationary equation

$\nabla ({\lambda }_{ef}\nabla T)+q_{\nu }-q_{v,t}=0$ (20)

With boundary conditions identical to the initial (index 0) and newly established (index 1) modes.

Numerical solution of Equation (2) gives the values of average temperatures of the equivalent continuous medium in the calculation cells. Coolant temperatures averaged over the volume of gaseous component in all cooling channels or in part thereof incorporated in the cells, as well as wall temperatures averaged over the perimeter of the channels are determined in each of the cells in the course of solution. Values of local maximum temperatures T ^max in each of the cells are additionally determined for the purpose of obtaining the whole picture of the temperature regime of the reactor core both in steady-state and in transient processes. Calculation expressions for cell i containing channel j have the following form:

T_i,_j^max = (T_st)_i,_j + ΔT_i,_j^max, ΔT_i,_j^max = b_i{q_v/[(1 – ε)λ₁]}_i,_j, (21)

where b_j for the channel with diameter d_j is the constant:

b_j = (d_e/4)²[2ln(d_e/d_j) + (d_j /d_e)² – 1]. (22)

Here, the value d_э as the diameter of conventional isolated area with cooling channel d_j in the center is determined by the dependence

d _e = [4S/(πN)]^1/2, (23)

where S is the total area of the reactor core cross-section, m²; N is the total number of all cooling channels in the reactor core.

The following numerical experiments were performed for the object imitating the core of gas-cooled nuclear reactor under the GT-MHR development project.

CFD-modeling of the reactor core requires significant computational capacity not available for the authors. Because of this reason comparison of results of calculation of transient processes using the methodology developed by the authors with those obtained using the algorithm with numerical CFD-modeling of the non-stationary process developing after exercising the perturbation was performed using the fragment of the reactor core (Fig. 1). Neumann’s boundary conditions corresponding to thermally insulated fragment of the reactor core were used both in the calculation by the authors and in CFD-modeling.

Figure 1. 

Fragment of reactor core: 1 – channels with 15.88-mm diameter; 2 – channels with 12.7-mm diameter. Dimensions are given in mm; only a half of the fragment is shown in the figure because of mirror symmetry; total length of the fragment is equal to 1074.34 mm; height of the fragment is the same as the height of the reactor core equal to 7930 mm.

Qualitative and quantitative correspondence of calculated values of temperatures and those obtained as the result of CFD-modeling for the examined options of profiling the fragment of the reactor core and different types of perturbations are illustrated in Figures 2 and 3. Discrepancy between the temperature values calculated using the suggested model and those obtained by CDF-modeling does not exceed 23 K for any type of transient.

Figure 2. 

Comparison of maximum and average temperatures in the reactor core with profiling under the condition of similar heating in case of power surge by 50%: 1 – values of maximum temperature obtained using the suggested algorithm; 2 – values of maximum temperature obtained in CFD-modeling; 3 – values of average temperature obtained using the suggested algorithm; 4 – values of average temperature obtained in CFD-modeling.

Figure 3. 

Comparison of maximum and average temperatures in the reactor core with profiling under the condition of similar wall temperature in case of drop of coolant mass flowrate by 50%: 1 – values of maximum temperature obtained using the suggested algorithm; 2 – values of maximum temperature obtained in CFD-modeling; 3 – values of average temperature obtained using the suggested algorithm; 4 – values of average temperature obtained in CFD-modeling.

Advantages of the core with profiling during the transient, essentially emergency, process as compared with the core without profiling is illustrated in Figures 4 and 5. Indeed, temperature in the reactor core without profiling coolant mass flowrate increases faster than in the core with mass flowrate profiling. For example, value of maximum temperature equal to 1700 K following the power surge of 50% with respect to the nominal power level will be reached after 140 s in the core without mass flowrate profiling, while in the core with profiling following the condition of either similar coolant heating or the condition of similar maximum cooling channel wall temperatures it is reached after 252 s (see Fig. 4).

Figure 4. 

Comparison of maximum temperature in the reactor core with different profiling conditions for the case of power surge by 50%: 1 – maximum temperature in the core without profiling; 2 – maximum temperature in the core with profiling under the condition of similar heating; 3 – maximum temperature in the core with profiling under the condition of similar maximum wall temperatures.

Similarly, in the case of sharp step-like drop of coolant mass flowrate by 50% with respect to the nominal flowrate value the difference in time for reaching maximum temperature equal to 1800 K amounts to about 150 s (see Fig. 5).

It is evident that the core of nuclear reactor under the development GT-MHR project with coolant mass flowrate profiling is more reliable against the core without profiling in the situations requiring certain time for initiating response measures for accident localization.

Let us note (see Figs 4, 5) that transient processes practically coincide in reactor cores with profiling to satisfy different conditions (maintenance of similar coolant heating or similar channel wall temperatures). However, it has to be taken into consideration, as it was noted in (Kuzevanov and Podgorny 2017), that pressure drop in the core with mass flowrate profiling under the condition of equality of maximum temperature in the channels is smaller than in the core with conventional profiling for maintaining similar coolant heating in the cooling channels of the reactor core. In our case the difference is by 11%.

Figure 5. 

Comparison of maximum temperature in the reactor core with different profiling conditions for the case of drop of mass flowrate by 50%: 1 – maximum temperature in the core without profiling; 2 – maximum temperature in the core with profiling under the condition of similar heating; 3 – maximum temperature in the core with profiling under the condition of similar maximum wall temperatures.

Algorithm was developed on the basis of which computer code was written for calculating temperature field in gas-cooled reactor during transient processes. Computer code was verified by comparing the results of numerous calculations with results of CFD-modeling of respective transient processes. The code has the advantage associated with the possibility of real-time calculation of evolution of conditions of gas-cooled cores with complex configuration containing channels with different diameters, which allows using the code in the calculations of transient processes in the cooling loops of nuclear facilities as a whole, in particular, in the development of reactor simulators.

Results of calculations of transient processes are presented for the reactor core imitating the core of high-temperature gas-cooled nuclear reactor under the GT-MHR development project for different conditions of profiling mass flowrate. Results of calculations unambiguously demonstrate significant differences between temperature regimes in the reactor cores with and without profiling during transient processes and indisputable enhancement of reliability of nuclear reactors with profiling coolant mass flowrate in the reactor core as a whole.