CALCULATION OF VVER-1000 CORE BAFFLE TEMPERATURE DISTRIBUTION FOR IT'S SWELLING ASSESSMENT

В.В. Філонов, Ю.С. Філонова, Я.Р Дубик., А.В. Богдан. Розрахунок температурного поля вигородки реактору ВВЕР-1000 для аналізу її розпухання. В роботі представлено спрощену CFD-модель охолодження вигородки реактору ВВЕР-1000, розроблену для подальшого аналізу величини її об’ємного розпухання. Так як явище розпухання є основнимобмежуючим фактором при продовженні терміну експлуатації енергоблоків з ВВЕР-1000, а моделі розпухання дуже чутливі до температурного поля в металі, його визначенню приділяється особлива увага. В даній роботі пропонується підхід із застосуванням засобів обчислювальної гідродинаміки (CFD), що дозволяє врахувати локальні гідродинамічні особливості потоку теплоносія, а також азимутальні розподіли характерних параметрів. Розроблено аналітичну модель для оцінки характеристичних параметрів спрощеної CFD-моделі, що дозволяє обґрунтовано звузити її межі. Розрахункова модель охолодження вигородки, що обмежена її висотою, має 60-градусну симетрію та включає активну зону, вигородку, шахту, спрощену геометрію з’єднувальних шпильок та омиваючий теплоносій. Активна зона представлена у вигляді еквівалентного гомогенного тіла з урахуванням просторового розподілу об’ємного енерговиділення. Вигородка розглядається як монолітне тіло, в якому враховано об’ємні енерговиділення за рахунок гамма випромінювання. Також, в моделі враховано циркуляцію охолоджуючого теплоносія через проточки в гайках, що дозволяє отримати більш реальне температурне поле шпильок. Отримані коефіцієнти тепловіддачі та температури добре узгоджуються з аналітичною оцінкою та дають більш прийнятні результати в порівнянні з RELAP5. Отримане поле температур використане для оцінки процесу розпухання. Через менш консервативні результати температурного поля величина розпухання вигородки суттєво зменшується. Розроблену модель в подальшому вдосконалено та використано для розрахунку зміни температурного поля при протіканні представницького перехідного процесу порушення умов нормальної експлуатації. Результати моделювання нестаціонарного процесу використано при оцінці необхідності розрахунку на прогресуючу формозміну. Ключові слова: CFD-моделювання, ВВЕР-1000, вигородка, температурне поле, розпухання матеріалу


Introduction
Possibility of lifetime prolongation for the NPPs with VVER-1000 reactors is a critical task for Ukraine. The main limiting factor in the long-term operation of power units is the phenomenon of PVI's (pressure vessel internals) radiation swelling under the action of high-energy neutron flux. Radiation swelling and creep processes models are very sensitive to the temperature distribution in the metal, if the temperature changes within 30 degrees, the swelling value changes 2.5 times.
To estimates the baffle progressive form change, transient calculation of the baffle temperature field during the representative abnormal operation is also necessary.
The determination of the temperature field consists of two sub-tasks: estimation of radiation energy release and solving the thermal conductivity problem. Previously, there was uncertainty in the estimation of radiation energy release, which completely leveled the difference in the formulations of the thermal conductivity problem. The methods for radiation energy release evaluation have improved over the time, and recently do not differ significantly. Due to this fact uncertainty in the temperature field estimation becomes essential. At current stage, three-dimensional formulations of the complete baffle model are considered, but there is uncertainty in the boundary conditions, which are estimated in the RELAP system code.

Analysis of recent publications and problem description
In recent years the use of CFD analysis has become widespread in calculations to clarify the distribution of thermal and hydraulic parameters in the reactor [1] and to obtain profiles of the coolant temperature distribution during the transient processes [2,3]. It must be noted that the core is represented in the form of a porous body [1,4]. In order to predict temperature distribution in the core baffle metal, CFD analysis is used only for PWR reactors of a different type (not for VVER's) in order to predict the temperature embrittlement processes of structural elements [5]. However, with regard to the baffle swelling problem, which is currently an important industry issue, there are uncertainties in preliminary solution of the temperature problem and, as a consequence, there is a need in development of improved approaches.
The purpose and objectives of the study Therefore, the purpose of this work is to develop an improved approach for the baffle metal temperature evaluation during operation at the nominal power level for further calculations of swelling, as well as during transient process of abnormal operation for estimation of the progressive form change and cyclic loadings. The use of CFD to solve this problem allows us to take into account the local hydrodynamic features of the coolant flow, as well as the azimuthal distributions of characteristic parameters.

Preliminary analytical estimation
At the stage of problem formulation and defining the CFD-model boundaries, the approach of preliminary analytical estimation of necessary boundary conditions was developed. It allowed us to narrow the boundaries of the model and reduce the need for computing resources. The distribution of coolant flow through PVI, as well as the upper limit of HTC values from the baffle inner surface (in the assumption of equivalence to annular channel, taking into account the direction of heat flux) were determined analytically. This paper provides general information on the approaches used, a more detailed description is provided in [6].
Analytical model and evaluation procedure description. The evaluation of thermal hydraulical properties was carried out by the method of determination the average and maximum loaded core in one-dimensional approximation (determining axial coordinate). The maximum loaded fuel assembly has characteristics (mass flow rate and energy release amplitude) that are proportional to the energy release unevenness coefficients (for heat flux -K V , for mass flow rate -K r ). Unevenness coefficients are assumed to be equal to the maximum values that can be observed during normal operation.
The model takes into account the core volume with a 5 mm gap between the peripheral fuel assemblies and the baffle inner surface. It is assumed that each fuel assembly forms an isolated prismatic volume, thus mass transfer of the coolant between fuel assemblies is not considered. The main assumption of the analytical model is the equality of pressure in the outlet cross section. The pressure distribution along the height of the core is assumed to be linear.
In order to validate the model, quantitative parameters are monitored, it is a necessary condition for the correct calculation, especially: the fuel core temperature (during normal operation it should not exceed 1600 °С), reserve to the boiling crisis of the first kind in the maximum loaded fuel assembly and coolant heating value.
The minimum value of the reserve to the heat transfer crisis in the maximum loaded cell is 1.56, the maximum temperature of the fuel core is 1512 °С. Such results indicate the adequacy of analytical model. Coolant heating is equal to 29.5 °С, which does not exceed the TVSA project limits (also valid for TVS-WR). The results of the estimation are applicable for the normal operation regime, without accounting for the measurement errors.
Determination of the coolant mass flow rate distribution in the PVI channels Estimates of flow rates were performed using empirical and semi empirical expressions, accounting for coolant flow specific. An equality of static pressure drop across all hydraulic elements in the reactor was made as a basic simplifying assumption. The core, the baffle and the annulus between baffle and core barrel were represented by a simplified equivalent hydraulic scheme. Each channel of a specific group, that takes into account the adjacent channel in the core basket belt and the annular space between the baffle and the barrel, can be represented by the geometric similarity in a simplified form of their real geometry. This approach allows to distinguish between different groups of pressure drops, namely loss on entrance, friction loss, expansion and contraction of the flow. The resulting mass flow rates are used in the CFD-model as inlet boundary conditions.

Estimation of HTC values from the baffle internal surface.
To determine the distribution of HTC from the baffle inner surface, the core (in equivalence to a simple system) can be treated as: -smooth-walled circular channel, which is cooled by coolant with cross-section averaged velocity; -equivalent circular channel with considering of core friction surface and with cross-section averaged coolant velocity; -annular channel, taking into account the direction of the heat flux and with cross-section averaged coolant velocity.
It is interesting to note that similar estimations of the HTC from the baffle inner surface to the coolant, carried out in RELAP5 gave an average value -14500 W/(m 2 K). This value was compared with the results obtained during the analytical evaluation (assuming the equivalences listed above). The closest values are obtained in predicting the equivalence of a smooth-walled circular channel (using an equation similar to Dittus-Boelter's formula in RELAP5 [7]). This indicates that the flow characteristics in the core are not taken into account, which leads to underestimated values of HTC in RELAP5 results.
For further calculations the upper limit of HTC value was obtained from estimation, where the gap between the peripheral fuel assemblies and baffle inner surface is schematized as an annular channel considering the direction of heat flux (received minimum/maximum value -29382/31898 W/(m 2 K)).
Reactor baffle cooling three-dimensional CFD-model CFD analysis is necessary to clarify the temperature distribution in the baffle. It permits to: -take into account the influence of the heat generation (gamma+neutron) in the reactor internals on the distribution of thermal hydraulic characteristics at the baffle boundary; -obtaine a temperature distribution in the connecting studs of the baffle rings, and the threaded rod; -consider the energy field distribution influence on the most heated baffle cross-section position. The resulting edge effects decrease the maximum temperature due to the presence of thermal conductivity in the axial direction and also the axial profile of the heat generation distribution in the baffle metal.
Model description. The CFD model has a two-plane 60-degree symmetry. Geometry is limited by the baffle height in the axial direction. Our three-dimensional CFD modeling has some geometric simplifications and physical assumptions: -The core is represented as a porous body. Porous body integral characteristics and intensities of heat transfer are estimated by the analytical model.
-The baffle is presented as monolithic body, in which there are no cooling fins. Thermal resistance in the contacts is neglected.
-The cooling channels cross-sectional area changes and, consequently, the actual vortex structure of the flow is disregarded, as this will unreasonably increase the calculation mesh density. Mass flow rate of the coolant through the channels determines the coolant speed, as previously estimated analytically (with local pressure losses).
-The coolant is considered to be a singlephase mixture, the content of a gas phase, which can occur during water radiolysis processes, as well as the vapor phase (up to 0.5…1.1 %) in the presence of hypothetical abnormal thermalhydraulic processes [8] are disregarded.
-Changes in the thermophysical properties of the internals steel under influence of radiation degradation are also ignored.
It should also be noted that there are physical assumptions that are impossible to account due to the absence of extensive and intensive process characteristics (for example, consideration of the actual thermal resistance between the baffle rings). The general view of the CFD model is presented in Fig. 1.
Boundary conditions. The boundary conditions for the calculated model are of Dirichlet type for the pressure and mass BC and also of the Neumann type for thermal conditions. Energy release in the core is equal to 104 % of nominal power level.
Boundary conditions specifications for the external surfaces of the model are shown in Fig. 2.

Fig. 2. Boundary conditions of simplified baffle cooling CFD-model
Mass boundary conditions at inlet and pressure boundary conditions at the outlet were applied to all groups of cooling channels. For the core, in addition to the mass and pressure BC, the porous body parameters were: the porosity coefficient -0.5815, the specific pressure loss coefficient (hydraulic resistance coefficient, which corresponds to the core height unit) -4.28 1/m, the specific heat ex- change area (heat transfer area divided by the volume of the porous core domain) -183.9 1/m and the distribution of the heat transfer coefficients. In the baffle, internal sources of heat generation having an amplitude of 8.5…9 W/cm 3 were specified. Additionally, heat generation of about 0.9…1.1 W/cm 3 was set in the core barrel, to get a more realistic temperature distribution. Outside of the core barrel, conditions of a third type were used, namely -the temperature of the coolant and the heat transfer coefficient, calculated in RELAP5.
Computational mesh. The boundary layer defines the processes of energy dissipation and heat transfer intensity, so special attention was paid to its discretization. To estimate the discretization, we used semi empirical dependencies, neglecting the edge effects in the channels. The thickness of the first boundary layer element was determined from the following expression: The coefficient k, that characterizes the height growth degree of each subsequent boundary layer element, must not exceed 1.5. Then, the number of the boundary layer divisions can be determined based on the sum of the geometric progression, as well as the boundary layer total thickness у.
For a two-parameter turbulence model of the k-ε family which has analytic wall functions, the appropriate value of y + should not exceed 300 [9].
In the model, functions of energy sources (radiation energy release in a metal due to attenuation of a gamma-ray flux) are spatially dependent, therefore insufficient to determine the quality of the boundary layer. It is also necessary to investigate the sensitivity of the model on different meshes [10].
Calculations for four mesh discretization variants were carried out. Independent results from mesh discretization were achieved in calculations for mesh#3 and mesh#4 [6]. Therefore, Mesh#4 was chosen for further calculation because it generally has less elements, but better discretization in baffle and studs. A detailed description of the calculation meshes characteristics and the comparison of obtained results are given in [6].

Results of baffle temperature field calculation at nominal power level (considering 4 % error)
According to the results, the maximum loaded cross section is at the level of 3.325 m from the bottom end of the baffle. Temperature distributions in the baffle cross sections are shown in Fig. 3. The maximum temperature in the baffle metal reaches a value of 368.9 °С.

. Temperature distributions in the most heated (а) and other baffle cross-sections (b)
General temperature ranges for other characteristic elements of the model, such as the outer surfaces of the baffle and connecting studs are shown in Fig. 4.

Fig. 4. Temperature distribution on the outer surfaces of characteristic elements
The results show the presence of reliable heat transfer from the studs by circulating coolant (Fig. 5), the temperature of which reaches the maximum value of 322 °С in contrast to the ÚJV Řež results [11]. This difference is explained by the presence of grooves in the nuts of studs in our model, whereas in [11] a stagnant zone around the studs is postulated and the boiling of the coolant is achieved.
The results of CFD analysis correspond with analytical estimation, which used the boundary conditions from RELAP. Comparison of obtained average values of heat transfer coefficients and coolant temperatures in channels is performed in Table. The main differences are related to the consideration of spatial features of heat transfer processes. When comparing the heat transfer coefficients, the main difference is caused by different reference temperatures and the influence of the heat flux on the temperature function profile in boundary layer. Similar estimations of the HTC from the baffle inner surface to the coolant, which were carried out for Ukrainian NPPs in RELAP5, are insufficient. It can be explained by comparison of the average HTC value obtained in RELAP5 (14500 W/(m 2 K)) with estimations based on empirical formulas. Results of CFD-analysis, were used for baffle swelling calculating. Strength assessment. Calculated temperature fields of the baffle are less conservative compared to the RELAP data and also are more detailed, in particular for the studs. However, the final results of the baffle calculation are the values of volumetric swelling, equivalent creep deformation and equivalent stresses, which determine the baffle material degradation and affect load capacity of the PVI elements. Therefore, a comparison of the listed target functions will be representative. The swelling calculation is based on the formula [Ошибка! Источник ссылки не найден., 13]:  (2) is the simple law of swelling, also called the law of free swelling. In general, swelling also depends on stress 1 ( ) eff f σ and equivalent plastic deformation 2 (ae ) p f : . Thus, the full calculation on the Eq. (3) is quite complex and requires simultaneous consideration of the irradiation-induced swelling, creep and plasticity effects. Analysis of (2) shows, that temperature is one of the determining factors because the dependence between swelling and temperature is exponential and, even a slight change in it, results in significant differences in swelling value.
Comparison of calculations based on RELAP data (see Table) and obtained from CFD calculations is shown in Fig. 6. It should be noted that the neutron irradiation dose is assumed to be the same. Two-dimensional model of the most loaded baffle cross-section using hypotheses of a plain strain state is created. The problem is substantially nonlinear, so load was increased slowly (with the step -1/1000 of the fuel campaign time) to achieve better convergence of results. Fig. 6 compares the results obtained for the predicted lifetime. Obtained values of volumetric swelling and creep differ 1.5…1.6 times even with a slight difference in temperatures for a 2D-model (the difference in maximum temperatures does not exceed 10 °C). Difference in equivalent stresses is less. The correct values of the studs temperatures are very important as their strength is one of the limiting factors for the long term operation. Such a difference in the target parameters is obviously significant for the lifetime extension of the Ukrainian NPP units.

Simulation of the baffle cooling during transient process
To perform calculations of cyclic strength and estimate the possible baffle progressive form change, transient calculation of the baffle temperature field during the abnormal operation is carried out. A regime with maximum changes in pressure and temperature gradients across the baffle thickness is selected for analysis.
In order to minimize computational time, we simulated a part of the transient process. The upper time limit of the calculation is determined based on the following requirements: ISSN 2223-3814 (online) ENERGETICS 43 -the average static pressure in the primary circuit is stabilizes (while the pressure drop remains unchanged); -the temperatures at the core inlet and outlet are constant (coolant flow is persistent); -the time limits of calculation regime include all characteristic events that affect the energy release and heat transfer from the baffle elements (emergency protection -130 s, shutdown of the MCP's -755…950 s); -the independence of average and maximum baffle metal temperatures on time is achieved. The execution time of the first three criteria is determined on the basis of RELAP input data (temperature and pressure changes at the boundaries of the core and inlet flow rate). The last criterion is monitored during calculations.
Boundary conditions estimation. Developed calculation model is geometrically and physically simplified. The time-dependent BC's (values of the coolant mass flow rates through PVI, as well as the heat transfer coefficients from the inner baffle surface) were formulated using analytical estimation. The necessary initial data for preliminary analytical estimation were taken from the RELAP5 calculation data. It is assumed that each time step from RELAP data characterizes the pseudo-steady state of the general transient process. The inputs to the analytical estimation are: static pressure at the core outlet, pressure drop in the core, coolant temperature at the inlet, and mass flow rate through the core. Coolant heating value in the core selected as an analytical solution adequacy criterion.
After emergency protection operation, the power generation in the core decreases to the level of residual energy and in an hour ranges from 1 % to 1.5 % of the nominal power level. In order to account for the emergency protection action in the calculation model, the energy release in the core is multiplied by the time-factor. This time-factor characterizes the decrease in the energy release in the porous body model.
Since the pressure drop as well as the average velocity of the coolant in the core vary over time, changes in the pressure loss factor (PLF) are also taken into account. The result of the analytical BC's estimation is the time-dependent values of the mass flow rates through the baffle cooling channels and the HTC from the fuel rods and from the baffle to the coolant from the core side (Fig. 7). Transient calculation results. CFD calculations are made using the mesh#4. The initial distribution of thermal hydraulic parameters for all domains corresponds to the results obtained for the reactor nominal power level (considering 4 % of measurement error). Radiation energy release in PVI's metal was multiplied by the decrease factor after the emergency protection operation (with/without considering of 25 % error). Figures 8 and 9 show the obtained changes in the baffle temperature, with the reduction factor of radiation energy release after the emergency protection operation. The steady-state and transient calculations were performed using ANSYS.  Fig. 9 shows the maximum stresses during the abnormal operation regime. The obtained results allowed to prove (in accordance with [14]) the absence of the need to calculate the progressive formchange.
Conclusions This paper presents a simplified CFD-model of VVER-1000 baffle cooling during operation at the nominal power level to address swelling problem. Developments of approach with preliminary analytical estimation of necessary boundary conditions allowed to simplify and narrow the model boundaries.
The calculation model has a 60-degree symmetry, which includes: a simplified core in the form of a porous body with regard to the volumetric energy release, the metal of internals (baffle, threaded rods, studs and barrel) taking into account gamma heating and the coolant. Model geometry is limited by the height of the baffle. Analyses of the mesh discretization in the boundary layer and calculations on four mesh variants with different spatial discretization were made. Since independent results from mesh discretization were achieved in calculations for mesh#3 and mesh#4, mesh#4 was chosen for further calculation because it generally has less elements, but better discretization in baffle and studs. According to the results, the maximum loaded cross section is at the level of 3.325 m from the bottom end of the baffle. The maximum temperature in the baffle metal reaches a value of 368.9 °С. The reliability of heat transfer from the stud's metal was also investigated. The results show the presence of reliable heat transfer from the studs by circulating coolant with estimated flow rate of 0.1…0.15 m 3 /h. The temperature of coolant around the studs reaches the maximum value of 322 °С, so the boiling of coolant is not reached due to the absence of a "stagnant" zone around the studs.
The obtained temperature distribution in the baffle is less conservative, compared with the previously used results, so it was used for baffle swelling calculating. Due to the reduced conservatism in temperature distribution calculation, the estimated values of volumetric swelling and creep deformation decreased 1.5…1.6 times for the two-dimensional model. The calculations for the full three-dimensional model will be even less conservative, especially in terms of assessing the strength of the studs.
To perform calculations of cyclic strength and to clarify the need of the baffle progressive form change estimates, the developed CFD-model was enhanced to estimate the temperature field during transients. The most representative regime of abnormal operation with maximum changes in pressure and temperature gradients across the baffle thickness was selected for analysis. The procedure for estimating the boundary conditions was adapted to obtain the time dependent values of mass flow rates through the cooling channels and HTC from the baffle inner surface. The reduction of volumetric energy release in the core due to the emergency protection operation in the chosen transient regime was also accounted. The obtained changes in the temperature distribution over time allowed us to estimate the maximum stresses during the representative abnormal operation regime. The obtained results showed (in accordance with [14]) the absence of the need to calculate the progressive form-change.