Corresponding author: Sergey K. Podgorny ( serkonpod@gmail.com ) Academic editor: Boris Balakin
© 2019 Vyacheslav S. Kuzevanov, Sergey K. Podgorny.
This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Citation:
Kuzevanov VS, Podgorny SK (2019) Temperature field in gascooled reactor core in transient conditions under different approaches to mass flow profiling. Nuclear Energy and Technology 5(4): 297303. https://doi.org/10.3897/nucet.5.48392

Positive effect of profiling the gascooled reactor core within the framework of the GTMHR project was investigated in (
The present paper is dedicated to the investigation of development of transients in gascooled 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 (
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 gascooled 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 CFDmodeling of respective transients in the object imitating the core of gascooled nuclear reactor. The advantage of the developed computer code is the possibility of realtime calculation of evolution of conditions in complex configurations of gascooled 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 gascooled nuclear reactor within the framework of GTMHR 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 hightemperature heliumcooled nuclear reactor were investigated in details in (
Reactor operation on nominal power in steadystate 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 gascooled nuclear reactor are formulated, calculation results are presented and their comparison with results of CFDmodeling is made using the example of GTMHR nuclear reactor core.
Nuclear реактор within the framework of GTMHR reactor design project (
In the model representation the core is examined as the solid medium (index 1) with gas voids (index 2). Solid part is multicomponent consisting of fuel, graphite and metal; singlecomponent 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:
Thus, the aggregated heat and mass transfer in the core of the nuclear reactor is represented in the model as the superposition of the nonstationary process of thermal conductivity and the quasistationary 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 (
$\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 = ρ_{1}h_{1}(1 – ε) + ρ_{2}h_{2} ε; h_{1}, h_{2} are the enthalpies of the solid and gaseous components, respectively, J/K; ρ_{1} is the density of the multicomponent solid fraction, kg/m^{3}; ρ_{2} is the helium density, kg/m^{3}; τ 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^{3}; λ_{ef} is the effective thermal conductivity coefficient for the equivalent medium, W/(m∙Κ) (definition of λ_{ef} corresponds to the model suggested in (
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}ρ_{2} + h_{2}∙δρ_{2}/δT_{2} = 0 in the case of zero value of h_{2} at the point T = 0 K; c_{1} is the average specific thermal capacity of the multicomponent 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} = σ_{i}q_{i}_{.d}, (3)
where σ_{i} = F_{i}/V_{i}; V_{i}, F_{i} are the cell volume, m^{3} and the heat exchange surface corresponding to this volume, m^{2}, respectively; q_{i}_{.d} is the heat flux density per unit surface, W/m^{2} (additional index d stands for the effective value of the parameter in the nonstationary 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 nonstationary process.
In steadystate mode q_{d} = q_{n}, where q_{n} refers to steadystate conditions and is easily determined using wellknown expressions, in particular the general energy balance equation:
$\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^{2}; 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_{2}._{d} averaged over the flow crosssection 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 quasistationary 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^{2}∙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 crosssection.
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 _{1} = π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_{2} 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 steadystate 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, quasistationary 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 steplike perturbation.
Let us include q, G, α ανδ T_{2}^{in} in the list of perturbation factors influencing the temperature regime of the reactor core. Let us denominate the ratio of the new steadystate values of these parameters to the initial ones as follows:
k_{q} = q_{1}/q_{0}; k_{α} = α_{1}/α_{0}; k_{G} = G_{1}/G_{0}; k_{T} = T_{2.1}^{in}/T_{2.0}^{in}. (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_{0}(z) for the steadystate mode:
q _{0}(z) = q_{0}^{max} sin (πz/H_{e}), (13)
where q_{0}^{max} is the heat flux density on the wall channel in the middle of the channel along its height in the initial steadystate 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 crosssection heat flux density q_{0} in the steadystate process is the function of two coordinates:
q _{0} = q_{0}(z,r) = q_{0}^{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_{2} from the condition of satisfaction of the following equation at the coolant inlet in the reactor core:
q _{0}(δ,r) = α_{0}[T_{st.0}(δ,r) – T_{2.0}^{in}]. (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 ]\mathrm{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_{0}) / (T_{1} – T_{0}), (19)
where Т is the solution of thermal conductivity equation during the transient; Т_{0} and Т_{1} 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 steadystate 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 – ε)λ_{1}]}_{i},_{j}, (21)
where b_{j} for the channel with diameter d_{j} is the constant:
b_{j} = (d_{e}/4)^{2}[2ln(d_{e}/d_{j}) + (d_{j} /d_{e})^{2} – 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 crosssection, m^{2}; 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 gascooled nuclear reactor under the GTMHR development project.
CFDmodeling 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 CFDmodeling of the nonstationary process developing after exercising the perturbation was performed using the fragment of the reactor core (Fig.
Fragment of reactor core: 1 – channels with 15.88mm diameter; 2 – channels with 12.7mm 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 CFDmodeling for the examined options of profiling the fragment of the reactor core and different types of perturbations are illustrated in Figures
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 CFDmodeling; 3 – values of average temperature obtained using the suggested algorithm; 4 – values of average temperature obtained in CFDmodeling.
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 CFDmodeling; 3 – values of average temperature obtained using the suggested algorithm; 4 – values of average temperature obtained in CFDmodeling.
Advantages of the core with profiling during the transient, essentially emergency, process as compared with the core without profiling is illustrated in Figures
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 steplike 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.
It is evident that the core of nuclear reactor under the development GTMHR 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
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 gascooled reactor during transient processes. Computer code was verified by comparing the results of numerous calculations with results of CFDmodeling of respective transient processes. The code has the advantage associated with the possibility of realtime calculation of evolution of conditions of gascooled 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 hightemperature gascooled nuclear reactor under the GTMHR 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.