Abstract
CIPS is a shift in the axial power towards the bottom half of the core, also known as axial offset anomaly (AOA), which results from the deposited of corrosion products during an operation. The main reason of CIPS is the solute particles especially boron compounds concentrated inside the porous deposit. The impact of CIPS is that the axial power distribution control may be more difficult and the shutdown margin can be decreased simultaneously. Besides, it also requires estimated critical condition (ECC) calculations to account for the effects of AOA. In this article, thermal-hydraulic subchannel code and boron deposit model have been combined to analyze the CIPS risk. The neutronics codes deal with the generation of homogenized neutron cross section as well as the calculation of local power factor. A simple rod assembly is analyzed with this combined method and simulation results are presented. Simulation results provide the boron hideout amount inside crud deposits and power shapes. The obtained results clearly show the power shape suppression in regions where crud deposits exist, which is a clear indication of CIPS phenomenon. And the CIPS effects on CHF have also been investigated. Result shows a margin of DNBR decrease in the crud case.
1. Introduction
It is universally acknowledged that axial offset anomaly (AOA), also known as crud induced power shift (CIPS), can suppress the neutron flux at the upper half of the core. The most famous AOA occurred in the late 1990s when Callaway experienced a severe AOA during Cycle 9 [1]. In order to recover shutdown margin, Callaway had to reduce its power to 70%. It must be noted that AOA results from a combination of multiple physical phenomena. One of the main reasons is crud, which is referred to be the deposition of corrosion products onto fuel cladding releasing from reactor coolant system surfaces. The porous crud is formed as a consequence of sub-cooled nucleate boiling (SNB) in the upper part of the assemblies. As the deposits are formed, boron concentrates within the crud while SNB is proceeding. Because a large amount of boron accumulates in the crud layer, the axial power distribution of nuclear core shift due to boron-10 has a high absorption cross section. Hence, neutron flux and the burnup peaking factor shift from top half to bottom half. Interestingly, near the end of cycle (EOC), excess burnup in the bottom of the core and reduced boron concentration in the primary system cause the power shape shift to reverse, partially restoring the burnup distribution. The appearance of CIPS makes the control of power distribution, especially axial, more difficult. Therefore, under the background of the objective of prolonged core life as well as the improvement of NPP’s safety and economy, it is essential to assess the CIPS risk of the core.
The CIPS phenomenon involves the thermal hydraulic, the formation and deposition of crud, boron advection and diffusion inside crud deposits, and neutron transport. Both the single mechanism and the coupling effects should be understood to perform a reasonable CIPS prediction. The relationship between evaporation and deposition has been studied in the past years, presenting the experimental evidences of deposits formed during boiling [2, 3]. The process of heat and mass transfer within crud deposits is a self-multi-physics mechanism, including a special capillary-derived boiling phenomenon, which is referred to as wick-boiling [4]. Water is absorbed into the deposit by capillary suction, evaporates into steam at the “mouth” of the chimney, and finally disappears in the coolant as steam bubbles. Upon the recent years, researchers have built upon this concept and developed it into different crud calculation models [5–11]. In the experimental research, scientists observed the fractalline properties of PWR fuel crud [12] and then found that a systematic theory of fractal porous media [13] can be used to refine some parameters and assumptions in the crud model.
Simulations between crud deposition models and thermal-hydraulic codes have also been investigated, such as the coupled solution between STAR-CCM+ and the crud model [14], the coupled simulation between COBRA-IV and the crud model [15], and the coupled calculation between CTF and the crud model, which models deposits growth, erosion, and thermal feedback [16]. In nuclear simulated circumstances, the neutronics and thermal-hydraulics coupling are the most classical combination to resolve core physical phenomena. By adding the crud prediction model, multi-coupling the simulation between STAR-CCM+, crud model, and neutronics code DeCART is also investigated [17–19]. In fact, although STAR-CCM+ is a very powerful CFD tool, it is still time-consuming, whereas subchannel codes are much more economical and acceptable for CIPS risk assessment and prediction. BOA 3.0 [20], a tool for CIPS risk assessment [21] which is codeveloped by Westinghouse Electric Company and the Electric Power Research Institute (EPRI), adopted the VIPRE-01 [22] subchannel code as the thermal hydraulic code.
The effort of this paper is to implement the crud and thermal hydraulic subchannel code together to simulate boron-accumulated amount along the assembly. The boron inventory can act as an input to the neutronics code, which could offer the power shape distribution along the assembly. In the following section, the combined methodology is presented in detail. The methodology consists of three parts: a boron deposit model based on Cohen’s model [5], the subchannel code COBRA-IV which is a sophisticated code for reactor core simulation, and the neutronics calculation codes DRAGON [23] and SKETCH-N [24]. Section 3 demonstrates the simulation results and shows the CIPS effect on a 33 fuel pin assembly. Additionally, the effect of various crud thicknesses is discussed. And the effects on CHF are simulated and discussed. Section 4 concludes this paper.
2. Methodology
2.1. Overview of the Framework
The present methodology in this paper provides a combined framework for modeling the CIPS phenomenon. This framework includes three modules (TH, crud, and neutron) and the simulation is processed in sequence. Firstly, an initial power distribution along axial fuel pin is calculated in the SKETCH-N code. In this stage, the fuel pellets are arranged uniformly along the axial height. So a cosine-shape power distribution can be solved and the power factor at each axial node is known. This power distribution shape is up-down symmetric and is also used for valuing the power shift at the end of stage. Secondly, coupled calculation between the thermal-hydraulic code and the crud model begins. While COBRA-IV code sweeping is at the axial level, channel sweeping and rod sweeping are the subcycles of the axial sweeping. Fluid temperature is solved in the channel sweeping cycle and the deposit location is determined by the SNB criterion in the rod sweeping cycle. If the criterion indicates the sign of deposit, the crud module solves the boron deposit amount on the local nodal. Lastly, the neutron absorption cross section is solved in the DRAGON code and then axial power distribution shape can be solved in the SKETCH-N code. The power shape diagram can be plotted, which could be regarded as an indication of CIPS risk. This calculation sequence is shown in Figure 1.
Due to the difference in solution meshes, the thermal hydraulic results are channel-centered while the crud solution domain is rod-centered. Because the rod sweeping cycle activates after the channel sweeping cycle, it is important to approximate the spatial transform from COBRA-IV code to the crud model. Here, volume-averaged density and mass-weighted temperature of the surrounding channels (northwest; northeast; southwest; southeast) are used as coolant information for crud model simulation. The schematic diagram is shown in Figure 2.
2.2. Module Descriptions
In the present work, the subchannel analysis code COBRA-IV is chosen for the thermal-hydraulic analysis. The main objective of this module is to calculate the fluid temperature and the cladding temperature. The fluid temperature serves as the simulating boundary of the crud model whereas the cladding temperature is calculated to identify the subcooled nucleate boiling region. This paper proposes a conservative assumption that crud will only occur in the SNB locations. Corrosion products source from corrosion surfaces, or from previously deposited crud back to the bulk coolant. The deposition rate of the crud caused by SNB is much higher than the nonboiling surfaces and it is assumed that the deposits have uniform thickness in the SNB region because this framework does not include deposited dynamics for now.
The subchannel code COBRA-IV uses Levy model [25] to calculate the subcooled void fraction and then predicts the outer cladding temperature and the local heat transfer coefficient by Thom correlation [26] (see (1)). The subcooled nucleate boiling occurs when the predicted outer cladding temperature is just above the saturated temperature. The Thom correlation is suitable for up to about 20MPa under conditions where the nucleate boiling contribution domains over forced convection. The correlation is used for estimate temperature difference when the local heat flux is given.where is the wall temperature elevation above the saturation temperature, K, q isthe heat flux, MW/m2, and P is the pressure of water, MPa.
It is assumed that the soluble metal ions only exist in liquid phase, so if bubbles evaporate from the surface of fuel pin, solute particles will precipitate and then coat on the boiling surface. The deposit has special morphology. The area, always filled with vapor during wick boiling phenomena, is called chimney. Microscopic exanimation of the deposits shows the existence of steam chimneys. The chimneys are across the deposit thickness orientation and point toward the coolant. Crud can grow to 10-100μm thick with the average pore size ranging within 0.1-1μm, porosity range 40-50% [27]. The chimney population represents the number of chimneys per mm2, 2000-5000 typically and each with a diameter of 2-6 μm.
The model used in this study calculates temperature and boron species concentration distributions in the crud deposit. The simulation unit cell consists of one steam chimney and its surrounding porous shell. Heat transfer is assumed to occur in the porous shell media. Conduction and convection take place in the medium while evaporation and condensation occur at the chimney boundary. Flow of liquid is modeled as diffusive flow while velocities are related to pressure gradients according to Darcy’s law.where is the Darcy velocity, m/s, is permeability the heat flux, MW/m2, μ is dynamic viscosity, Pa·s, and P is the pressure gradient vector of water, Pa/m.
And solute species are modeled by equilibrium between advective current and diffusive current.where is the net solute flux, is diffusive coefficient, m2/s, and C is boron concentration, mol/dm3.
A schematic diagram is shown in Figure 3, displaying a two-dimensional cylindrical simulation unit and its boundaries. The above equations are discreted by upwind differencing, finite-volume approach and iteratively solved until convergence. Temperature and pressure fields are solved separately and the results affect the concentration solutions. The solution process is couple-solved. This iteratively solved model is introduced by Haq et al. [6] and Short et al. [11] and has been reproduced by Li and Liu [15].
Obviously, deposit layer is established between fuel pin cladding and coolant. Therefore, there are two nature boundaries of a single simulation unit, which are referred to as coolant boundary and cladding boundary. The boundary conditions vary among different rods and various coolant conditions, such as rod powers at different rods and axial levels, or coolant temperatures and concentrations in different channels and various rod-coolant heat transfer coefficient. It should be noted that the origin COBRA-IV missed the boron transport calculation ability. In our previous research [28], the boron transport simulation capacity has been implemented in the COBRA-IV code. After the boundaries are all determined, boron hideout amount of any rod/assembly could be calculated. The boron hideout amount is proportional to deposit thickness, local heat flux, bulk concentration, and some porous characters due to some sensitive researches [9, 15, 29].
The precipitation of lithium metaborate in the crud layer is considered to be a root cause of CIPS. The equilibrium equation of this reaction is [7]
The neutron-kinetics codes DRAGON and SKETCH-N are selected as the neutron solver for the neutron-physical calculation. In order to calculate neutron cross sections, it is necessary to identify the compositions in the deposit. The majority compositions in the crud solid structure are Fe3O4, NiO, NiFe2O4, and ZrO2 [27]. Both water and concentrated solutes are in the pore region. After providing the composited element, the DRAGON code would interpolate the absorption cross section which is supplied by standard libraries. The cross sections of the crud layer are divided into seven groups and each group correspond to a range of energy. Finally, the SKETCH-N module would use these cross sections to print the shifted power shape diagram.
3. Numerical Results
In this section, a 3×3 fuel rods assembly with square lattice is solved under this combined method and the results are presented. A top view of the assembly geometry is schemed in Figure 4. Typical geometry parameters are listed in Table 1. Due to the symmetric geometry of the assembly, one-eighth of the assembly can represent the flow and heat transfer field of the whole assembly. Channel can be identified as corner, side, and center channel, respectively.
As the COBRA-IV input, the 3657.6mm-height fuel rod is divided into 36 equal axial levels, each with 101.6mm height. Channel exit pressure is set equal to 15.5Mpa, and inlet mass velocity is 3600kg/(m2s) with a temperature of 270°C (543.15 K). A cosine shape axial power distribution with the average heat flux value of 1000kw/m2 is used. In addition, the coolant boron concentration is set equal to 1200ppm.
According to the preceding part of the content, some parameters should be determined as the input of the crud calculation module. First, uniform deposit thickness is assumed at each location where SNB takes place. Then, the statistical majority number [27], the porosity is set equal to 0.6, chimney population is set to 2000 per square millimeter, chimney diameter is set to 4 μm, and the other porous medium character values are calculated based on fractalline porous medium concept. In order to compare different deposit thickness effects on boron hideout inventory, two cases are calculated in this work: one case with the thickness of 30 μm representing thin deposit and the other case with the thickness of 50 μm representing a thick one.
The corrosion source in primary system mainly comes from SG tubes. Crud solid almost consists of nickel, iron, and there oxide. Table 2 shows the volume fractions of crud solid which are assumed in this study. Crud is mainly composed of nickel ferrate, whose density varies slightly with temperature. At 600k, the density is 3.9 g/cm3 [31]. Each compound element in the crud and the boron amount is transferred into mass fraction. A pin geometry with crud layer, which is also a schematic drawing for SKETCH-N, is shown in Figure 5. The geometry is the same as the thermal hydraulic input. Because of the radial symmetry of the component, the fuel assembly is treated as a fuel pin with crud and coolant layer around it. The axial active zone is of 3657.6mm-height and it is divided into 18 levels. Reflect boundary layers are arranged at both the top and the bottom of the fuel pin.
3.1. SNB Region and Boron Hideout Amount
After completely sweeping the axial levels along the fuel assembly, the calculation of combined TH-crud module is finished. Boron hideout mass along the assembly is figured out. It is found that boron inventory distributes relative uniform in the radial direction due to the high symmetric of fuel pin; the axial distribution is more uniform. The axially boron hideout mass result is drawn in Figure 6. This figure also indicates the locations where SNB takes place. SNB takes place from the middle upper part of the assembly and ends before the channel exit. The SNB could not be maintained at the end of channel because of the decreased local heat flux.
3.2. CIPS Prediction and AO Value
The axial power distribution, shown in Figure 7, indicates a signal of CIPS phenomenon. Local power reduces as a result of neutron flux suppression. Also, because of the power constant in the whole reactor, the suppression of power in the crud pin region leads to an increased power at the nondeposit region. Therefore, it results in a downward shift in the axial power distribution. Moreover, it can be found that thick crud will bring more axial power shift than thin crud.
The axial offset anomaly value is defined by Westinghouse Company [1]. AO difference indicates the power difference between top half and bottom half of core:where is power integration for the top half of the core and is power integration for the bottom half of the core.
If AO is negative, power shift will occur in the core. The AO differences can indicate AOA levels of the power plant. These levels can be divided into “moderate AOA” with a difference of -5% and “severe AOA” with a difference of -10%. In this simulated case, AO value is equal to -10.4% for thin crud and -14.8% for thick crud. Both of them achieve the “severe AOA” difference. Typically, core design has sufficient excess shutdown margin to come out against moderate AOA; however severe AOA may treat the full power operation of the plant.
3.3. CIPS Effects on CHF
The shifted axial power shape will not only change the heat source along the rod, but also affect the thermal properties of the heat transfer mechanism. Subchannel codes are always adopted to calculate the critical heat flux (CHF) in the reactor design. The Westinghouse-3 (W-3) CHF correlation [32] is widely used for PWR conditions, where departure from nucleate boiling (DNB) is the dominant CHF mechanism. The departure from nucleate boiling ratio (DNBR) is a local value, which varies along the rod height. The minimum value of the DNBR (MDNBR) is a significant parameter to identify the worst location and its dangerous level. In this section, the DNBR are calculated both in the normal cosine power shape and in the shifted cosine power shape. Different thickness effects are also considered.
The results of the DNBR distribution along the assembly are presented in Figures 8 and 9. Curves clearly show the MDNBR location is changed and the value is also decreased in the shifted power cases. In the “no crud” case, MDNBR occurs at relative location of about 0.7 with the value of 3.32. In the “thin crud” case, MDNBR occurs at relative location of about 0.64 with the value of 3.11. And in the “thick crud” case, MDNBR occurs at relative location of about 0.61 with the value of 3.0. The distribution trend of MDBNR is the same as that of CIPS: towards the bottom of the core. And the thicker deposit worsens this situation.
4. Conclusion
In this study, a combination of thermal-hydraulic, crud, and neutronics modules is described. The aim of this framework is to predict the CIPS phenomena. A 33 fuel pin assembly is solved by this combined work and the axial power shape is given. The results clearly show the power is suppressed in the crud-existing region and rises in the normal region, which is a clear sign of axial power shift. Thick crud induces more power shift than thin crud. According to the AO value definition, the calculated case in this study achieves the severe AOA level. The CHF calculation results show the CIPS effects on thermal-hydraulic design. CIPS will cause the location of MDNBR to change and reduces the margin. However, this combined work is based on several assumptions, such as equal thickness deposit in SNB region and neglecting the redistribution effect caused by spacer grid and mixing vane. Advanced CHF correlation is also important because the surface of crud is rough and hydrophilic. And more efforts are needed to improve the simulation ability, especially validation and verification.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors have declared that no conflicts of interest exist.