Research Article  Open Access
A Combined Method for Predicting the Boron Deposited Mass and the CIPS Risk
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, thermalhydraulic 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 subcooled 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 boron10 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 selfmultiphysics mechanism, including a special capillaryderived boiling phenomenon, which is referred to as wickboiling [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 thermalhydraulic codes have also been investigated, such as the coupled solution between STARCCM+ and the crud model [14], the coupled simulation between COBRAIV 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 thermalhydraulics coupling are the most classical combination to resolve core physical phenomena. By adding the crud prediction model, multicoupling the simulation between STARCCM+, crud model, and neutronics code DeCART is also investigated [17–19]. In fact, although STARCCM+ is a very powerful CFD tool, it is still timeconsuming, 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 VIPRE01 [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 boronaccumulated 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 COBRAIV which is a sophisticated code for reactor core simulation, and the neutronics calculation codes DRAGON [23] and SKETCHN [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 SKETCHN code. In this stage, the fuel pellets are arranged uniformly along the axial height. So a cosineshape power distribution can be solved and the power factor at each axial node is known. This power distribution shape is updown symmetric and is also used for valuing the power shift at the end of stage. Secondly, coupled calculation between the thermalhydraulic code and the crud model begins. While COBRAIV 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 SKETCHN 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 channelcentered while the crud solution domain is rodcentered. Because the rod sweeping cycle activates after the channel sweeping cycle, it is important to approximate the spatial transform from COBRAIV code to the crud model. Here, volumeaveraged density and massweighted 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 COBRAIV is chosen for the thermalhydraulic 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 COBRAIV 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/m^{2}, 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 10100μm thick with the average pore size ranging within 0.11μm, porosity range 4050% [27]. The chimney population represents the number of chimneys per mm^{2}, 20005000 typically and each with a diameter of 26 μ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/m^{2}, μ 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, m^{2}/s, and C is boron concentration, mol/dm^{3}.
A schematic diagram is shown in Figure 3, displaying a twodimensional cylindrical simulation unit and its boundaries. The above equations are discreted by upwind differencing, finitevolume approach and iteratively solved until convergence. Temperature and pressure fields are solved separately and the results affect the concentration solutions. The solution process is couplesolved. 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 rodcoolant heat transfer coefficient. It should be noted that the origin COBRAIV missed the boron transport calculation ability. In our previous research [28], the boron transport simulation capacity has been implemented in the COBRAIV 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 neutronkinetics codes DRAGON and SKETCHN are selected as the neutron solver for the neutronphysical 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 Fe_{3}O_{4}, NiO, NiFe_{2}O_{4}, and ZrO_{2} [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 SKETCHN 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, oneeighth 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 COBRAIV input, the 3657.6mmheight 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/(m^{2}s) with a temperature of 270°C (543.15 K). A cosine shape axial power distribution with the average heat flux value of 1000kw/m^{2} 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/cm^{3} [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 SKETCHN, 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.6mmheight 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 THcrud 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 Westinghouse3 (W3) 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 thermalhydraulic, 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 crudexisting 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 thermalhydraulic 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.
References
 J. Deshon, PWR Axial Offset Anomaly (AOA) Guidelines, Revision 1, EPRI, Palo Alto, Calif,USA, 2004.
 N. B. Hospeti and R. B. Mesler, “Deposits formed beneath bubbles during nucleate boiling of radioactive calcium sulfate solutions,” AIChE Journal, vol. 11, no. 4, pp. 662–665, 1965. View at: Publisher Site  Google Scholar
 E. P. Partridge and A. H. White, “Mechanism of formation of calcium sulfate boiler scale,” Industrial & Engineering Chemistry, vol. 21, no. 9, pp. 834–838, 1929. View at: Publisher Site  Google Scholar
 R. V. Macbeth, “Boiling on surfaces overlayed with a porous deposit: heat transfer rates obtainable by capillary action,” Atomic Energy Establishment AEEWR711, Winfrith, England, 1971. View at: Google Scholar
 P. Cohen, “Heat and mass transfer for boiling in porous deposits with chimneys,” in Proceedings of the in AICHE Symposium series 70, pp. 71–80, American Institute of Chemical Engineers, 1974. View at: Google Scholar
 I. U. Haq, N. Cinosi, M. Bluck, G. Hewitt, and S. Walker, “Modelling heat transfer and dissolved species concentrations within PWR crud,” Nuclear Engineering and Design, vol. 241, no. 1, pp. 155–162, 2011. View at: Publisher Site  Google Scholar
 J. Henshaw, J. C. McGurk, H. E. Sims, A. Tuson, S. Dickinson, and J. Deshon, “A model of chemistry and thermal hydraulics in PWR fuel crud deposits,” Journal of Nuclear Materials, vol. 353, no. 12, pp. 1–11, 2006. View at: Publisher Site  Google Scholar
 M. Jin and M. Short, “Multiphysics modeling of twophase film boiling within porous corrosion deposits,” Journal of Computational Physics, vol. 316, pp. 504–518, 2016. View at: Publisher Site  Google Scholar  MathSciNet
 C. Pan, B. G. Jones, and A. J. Machiels, “Concentration levels of solutes in porous deposits with chimneys under wick boiling conditions,” Nuclear Engineering and Design, vol. 99, no. C, pp. 317–327, 1987. View at: Publisher Site  Google Scholar
 C. Pan, B. G. Jones, and A. J. Machiels, “Wick boiling performance in porous deposits with chimneys,” in Proceedings of the ASME, National Heat Transfer Conference Symposium on Multiphase flow and Heat Transfer, 1985. View at: Google Scholar
 M. P. Short, D. Hussey, B. K. Kendrick, T. M. Besmann, C. R. Stanek, and S. Yip, “Multiphysics modeling of porous CRUD deposits in nuclear reactors,” Journal of Nuclear Materials, vol. 443, no. 13, pp. 579–587, 2013. View at: Publisher Site  Google Scholar
 I. Dumnernchanvanit, V. K. Mishra, N. Q. Zhang et al., “The fractalline properties of experimentally simulated PWR fuel crud,” Journal of Nuclear Materials, vol. 499, pp. 294–300, 2018. View at: Publisher Site  Google Scholar
 B. M. Yu, “Analysis of flow in fractal porous media,” Applied Mechanics Reviews, vol. 61, no. 5, Article ID 050801, 2008. View at: Publisher Site  Google Scholar
 V. Petrov, B. K. Kendrick, D. Walter, A. Manera, and J. Secker, “Prediction of CRUD deposition on PWR fuel using a stateoftheart CFDbased multiphysics computational tool,” Nuclear Engineering and Design, vol. 299, pp. 95–104, 2016. View at: Publisher Site  Google Scholar
 S. Li and X. Liu, “Development of boron tracking and boron hideout (CRUD) model based on subchannel approach,” Nuclear Engineering and Design, vol. 338, pp. 166–175, 2018. View at: Publisher Site  Google Scholar
 R. Salko et al., “Development of COBRATF for modeling fullcore, reactor operating cycles,” Advances in Nuclear Fuel Management V (ANFM 2015), 2015. View at: Google Scholar
 D. J. Walter, B. K. Kendrick, V. Petrov, A. Manera, B. Collins, and T. Downar, “Proofofprinciple of highfidelity coupled CRUD deposition and cycle depletion simulation,” Annals of Nuclear Energy, vol. 85, pp. 1152–1166, 2015. View at: Publisher Site  Google Scholar
 D. J. Walter and A. Manera, “CRUD, boron, and burnable absorber layer 2D modeling requirements using MOC neutron transport,” Annals of Nuclear Energy, vol. 87, pp. 388–399, 2016. View at: Publisher Site  Google Scholar
 L. Zou, H. Zhang, J. Gehin, and B. Kochunas, “Coupled thermalhydraulic/neutronics/crud framework in prediction of crudinduced power shift phenomenon,” Nuclear Technology, vol. 183, no. 3, pp. 535–542, 2013. View at: Publisher Site  Google Scholar
 BoronInduced Offset Anomaly Risk Assessment Tool Version 3.0, EPRI, Palo Alto, Calif, USA, 2010.
 J. Deshon, D. Hussey, B. Kendrick, J. McGurk, J. Secker, and M. Short, “Pressurized water reactor fuel crud and corrosion modeling,” JOM: The Journal of The Minerals, Metals & Materials Society (TMS), vol. 63, no. 8, pp. 64–72, 2011. View at: Publisher Site  Google Scholar
 C. W. Stewart, “VIPRE01: a thermalhydraulic code for reactor cores,” Tech. Rep., Electric Power Research Institute, Pacific Northwest Laboratory, Palo Alto, Calif, USA, Richland, Washington, USA. View at: Google Scholar
 G. Marleau, A. Hébert, and R. Roy, “A user guide for DRAGON 3.06,” Rep. IGE174 Rev 7, vol. 7, 2008. View at: Google Scholar
 V. G. Zimin, SKETCHN: A Nodal Neutron Diffusion Code for Solving SteadyState and Kinetics Problems, JAERI Report, Tokyo, Japan, 2011.
 S. Levy, “Forced convection subcooled boilingprediction of vapor volumetric fraction,” International Journal of Heat and Mass Transfer, vol. 10, no. 7, pp. 951–965, 1967. View at: Publisher Site  Google Scholar
 W. M. Rohsenow, J. P. Hartnett, and Y. I. Cho, Handbook of Heat Transfer, McGrawHill, New York, NY, USA, 1998.
 J. Buongiorno, “Can corrosion and CRUD actually improve safety margins in LWRs?” Annals of Nuclear Energy, vol. 63, pp. 9–21, 2014. View at: Publisher Site  Google Scholar
 S. Z. Li, X. J. Liu, and X. Cheng, “Development of Boron transport model based on subchannel approach,” in Proceedings of the 11th International Topical Meeting on Nuclear Reactor ThermalHydraulics, Operation and Safety (NUTHOS11), Gyeongju, Korea, 9 October 2016. View at: Google Scholar
 D. J. Walter, A High Fidelity Multiphysics Framework for Modeling CRUD Deposition on PWR Fuel Rods, University of Michigan, 2016.
 A. T. Godfrey, “VERA Core Physics Benchmark Progression Problem Specifications,” CASLU20120131004, p. 189, 2014. View at: Google Scholar
 A. T. Nelson, J. T. White, D. A. Andersson et al., “Thermal expansion, heat capacity, and thermal conductivity of nickel ferrite (NiFe2O4),” Journal of the American Ceramic Society, vol. 97, no. 5, pp. 1559–1565, 2014. View at: Publisher Site  Google Scholar
 L. S. Tong, “Prediction of departure from nucleate boiling for an axially nonuniform heat flux distribution,” Journal of Nuclear Energy, vol. 21, no. 3, pp. 241–248, 1967. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2019 Shengzhe Li et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.