Radionuclide release in severe accidents

The release of radionuclides after a severe accident in a nuclear reactor poses risk to public health, causes serious environmental pollution, and impacts the ecological environment^[9]. Petit et al. used the radionuclide network of a comprehensive nuclear-test-ban treaty organization and collected data from laboratories to detect radionuclides and study the behavior of fission products^[10-11]. In the case of the Fukushima nuclear disaster, radionuclides spread into the environment, prompting a health management investigation and external radiation dose estimates of nearby residents ^[12]. A series of tests on the PHEBUS fission product was led by the Institute for Nuclear Safety and Radiation Protection in France. The fission product behavior of a light-water reactor under severe accidents was studied by simulating the reactor device to test the reactor core and containment and primary cooling systems^[13]. Kontautas et al. studied the aerosol diffusion behavior of the PHEBUS FPT-1 test using the COCOSYS code and analyzed the effects of aerosol density and solubility on aerosol deposition ^[14]. Haste et al. studied the special behavior and form of radionuclide iodine in the PHEBUS PFT-3 test, and the volatile fission products Cs and Mo had an important effect on aerosol iodine ^[15].

Some important phenomena of radionuclide iodine obtained from PHEBUS PFT0 and PFT1 test data were not predicted and calculated; therefore, a collation test was carried out to study the complex phenomena of fission product behavior in pressurized water reactors (PWRs) under severe accident conditions ^[16-17]. Some researchers used the marine PWR as the research object, using the severe accident-integrated analysis program, MELCOR, to establish an accident analysis model, and carried out a simulation study of radionuclide behavior after a severe accident. The behavior of radioactive fission products was studied under severe accident conditions by large-break loss-of-coolant accident (LOCA) combined with ship blackouts ^[18]. The release regularity of radionuclides in small-break LOCAs was studied, and the diffusion behavior of inert gas and volatile fission products was analyzed ^[19]. In a blackout accident, the creep failure of the pressurizer surge tube was confirmed, and the behavior of the source term before and after the break of the surge tube was investigated ^[20]. The combination of LOCA and the blackout accident of a marine nuclear reactor was modeled and calculated, and the accident process and source term were investigated ^[21]. A comparative study on the radioactive effects of different pressurizer safety valve working conditions was carried out under a combined situation of whole ship blackout and safety valve failure ^[22]. The radionuclide release under the condition of steam generator tube ruputure accident coupled with whole ship blackout was investigated ^[23].

Simulation of mathematical model for radionuclide diffusion

Computational fluid dynamics (CFD) combines hydromechanics, numerical mathematics, and computer science^[24]. Physical models, such as steady and unsteady flows, laminar and turbulent flows, incompressible and compressible flows, heat transfer, and chemical reactions, have been widely used to simulate complex phenomena^[25]. In order to simulate indoor environment radionuclide diffusion^[26], the k-ε model has become the most popular method in atmospheric applications among many CFD turbulence models ^[27]. Lyu et al.^[28] analyzed the hydrogen source term and hazard under the condition of a marine reactor with a LOCA and used the GASFLOW^[29-30] CFD code to compute the two-dimensional puff domain and the release and diffusion of hydrogen in the marine environment. GASFLOW Message Passing Interface (MPI), a CFD tool, was used to simulate the turbulent phenomenon of fluid in the reactor containment vessel^[31], and the full-velocity turbulent flowability of GASFLOW-MPI was verified^[32]. The results show that the deterministic sampling method effectively quantifies the uncertainty of CFD calculations ^[33].Ouyang et. al^[34] studied the diffusion of radioactive materials in the atmosphere in the event of a severe accident in a marine reactor using the Lagrange particle tracking method to establish a distribution model of radionuclides in the atmosphere. The diffusion and distribution of radionuclides in offshore waters under the condition of a nuclear leakage accident were studied using the Lagrange and Euler methods ^[35]. In order to evaluate the risk and harm of radionuclides to the public, the diffusion phenomena of radionuclides in the atmosphere during a hypothetical accident at the Haiyang Nuclear Power Plant was studied based on the improved research method of Weather Research and Forecasting Chemistry ^[36]. There is a hypothesis that a severe nuclear accident occurred in the Akkuyu nuclear power plant, and the diffusion analysis of radionuclide ¹³⁷Cs was carried out using a three-dimensional hydrodynamic model coupled with a Lagrange model, and the distribution regularity of the ¹³⁷Cs active concentration was studied ^[37].

These literature surveys demonstrate the importance of research on radionuclide release and diffusion after a severe accident. Based on the severe accident analysis program, there are separate studies on radionuclide release in a severe accident and radionuclide diffusion and migration in the atmosphere based on CFD. However, no method for studying the diffusion of radionuclides in a limited, confined space has yet been proposed. Therefore, a new approach using MELCOR, which is a severe accident analysis program coupled with fluid simulation, is proposed in this case study. This new method is used to study the diffusion of radionuclides in a confined environment, and it provides technical and data support for decision-making in a nuclear emergency.

MELCOR coupled with scSTREAM

MELCOR is an integrated second-generation system program with a relatively fast-running code. It is a probabilistic risk evaluation and accident analysis instrument developed by Sandia National Laboratories for the NRC. MELCOR can emulate the main phenomena of severe accident processes in pressurized water reactors (PWRs), and the MELCOR code models a wide range of physical phenomena, such as the release and transport of fission products. To accurately simulate the radionuclide diffusion behavior of a marine reactor after a severe accident in a closed environment, the results of the radionuclide release fraction and temperature calculated by MELCOR were coupled with the scSTREAM calculation.

scSTREAM is a general thermal fluid simulation software based on a structured grid (Cartesian or cylindrical coordinate systems). In addition to providing basic fluid/heat problem-solving functions, it also provides special-related functions. Compared with FLUENT, it can calculate fluid, heat transfer and chemical reactions，scSTREAM has the following characteristics and functions: powerful proprietary-building domain-import function, powerful grid-generation function–automatic generation of structured grid, robust and high-speed large-scale computing capability, intuitive graphical user interface, great flexibility as well as user function settings, excellent post-processing and visualization, and better coupling of calculation results using MELCOR.

Analysis process for radionuclide diffusion

The scSTREAM thermofluid analysis system consists of the following three sets of programs, which are independent and do not affect each other in the application and calculation.

Step 1 Preprocessor (Pre)

Step 2 Solver

Step 3 Postprocessor (Post)

To analyze a problem, input data files were first created in Pre. The solver then performs the calculation based on the data setup in Pre. The results were visually evaluated using Post. Several files were used as the interface for each program, as shown in Fig. 1.

Fig. 1Full-screen View

The relationship between files and programs

During radionuclide diffusion studies, the release share and temperature of radionuclides are coupled with the condition set in the preprocessor. When executing the calculation, the transient analysis can be performed only after the steady-state debugging and convergence are achieved. A flowchart that leads to the final result is shown in Fig. 2.

Fig. 2Full-screen View

The flowchart of scSTREAM process for radionuclide diffusion

The computational model for radionuclide diffusion

In numerical simulations to reproduce physical phenomena, there are four types of governing equations: momentum conservation equations (the Navier-Stokes equations), continuity equations, energy conservation equations, and model equations to represent turbulent transport. The finite volume method was used in Solver to solve the above equations on a computer. A finite volume method solves the governing equations by converting them into an integral conservation form that is expressed on each fractional unit of the decomposed elements or control volume.

If the velocity of a fluid is less than the speed of sound, the change in its density is extremely small; therefore, the density can be treated as a constant. Fluids of this type are known to be incompressible fluids. In this study, the diffusion of airborne radionuclides in the cabin is assumed to be incompressible. Mass conservation^[38] for incompressible fluids is simplified to the following Eq. (1). $\frac{\partial u_i}{\partial x_i}=0$ (1)

For incompressible fluids, the momentum conservation equation[39] is given by Equation (2): $\frac{\partial \rho u_i}{\partial t}+\frac{\partial u_j\rho u_i}{\partial x_j}=-\frac{\partial p}{\partial x_i}+\frac{\partial }{\partial x_j}\mu \frac{\partial u_i}{\partial x_j}-\rho g\beta \left(T-T_0\right).$ (2)

For incompressible fluids or low Mach number flows, approximation is applied. The energy conservation equation is given by Eq. (3): $\frac{\partial \rho C_{\text{p}}T}{\partial t}+\frac{\partial u_j\rho C_{\text{p}}T}{\partial x_j}=\frac{\partial }{\partial x_j}K\frac{\partial T}{\partial x_j}+\dot{q}$ (3)

The diffusive species conservation equation is shown in Eq. (4). $\frac{\partial C}{\partial t}+\frac{\partial u_{j}C}{\partial x_j}=\frac{\partial }{\partial x_j}{D}_{\text{m}}\frac{\partial C}{\partial x_j}+d$ (4) where $u_i$ is velocity of flow in $x_i$ direction (m/s); $x_i$ is coordinates; $i$=1,2,3 correspond to x, y, and z directions, respectively (m); $\rho$ is density of a fluid or a solid (kg/m³); $t$ is time (s); $p$ is pressure of a fluid (N/m²); $g$ is gravity (m/s²); $\beta$ is thermal expansion coefficient (1/K); $T$ is temperature of a fluid or a solid (K); $T_0$ is reference temperature of a fluid (K); $\mu$ is viscosity ($\mu ={\mu }_1+{\mu }_2$, ${\mu }_1$ is molecular viscosity and ${\mu }_2$ is eddy viscosity; kg/(m·s)); $C_{\text{p}}$ is specific heat at constant pressure (J/(kg·K)); $K$ is thermal conductivity (J/(m·s·K)); $\dot{q}$ is heat source (J/(m³·s)); $C$ is concentration of a diffusive species; $D_{\text{m}}$ is diffusion coefficient (m²/s); and $d$ is source terms of diffusive species (1/s).

The improvements are known as the RNG k–ε model, which stands for renormalization group analysis. The model constants that appear in the standard k–ε model have been established experimentally; however, the RNG k–ε model obtains the constants theoretically using Fourier analysis. For this reason, the first step is to modify the momentum equation and perform a Fourier transformation to obtain the k–ε model in the wavenumber space. However, after the Fourier transformation, the equation becomes one, in which the low wavenumber range is subjected to the influence of the higher wavenumber range, and the high wavenumber range is, in turn, affected by the even higher wavenumber range. In contrast, the k–ε model, in the first place, only deals with eddies of a large scale; in other words, the impact of the high wavenumber range is incorporated into the low wavenumber range. The equations of turbulent kinetic energy and its dissipation rate (k–ε equations) are as follows: $\frac{\partial \rho k}{\partial t}+\frac{\partial u_i\rho k}{\partial x_i}=\frac{\partial }{\partial x_i}\left(\frac{{\mu }_t}{{\sigma }_k}\frac{\partial k}{\partial x_i}\right)+G_{\text{S}}+G_{\text{T}}-\rho \epsilon$ (5) $\begin{array}{c}\frac{\partial \rho \epsilon }{\partial t}+\frac{\partial u_i\rho \epsilon }{\partial x_i}=\frac{\partial }{\partial x_i}\left(\frac{{\mu }_t}{{\sigma }_{\epsilon }}\frac{\partial \epsilon }{\partial x_i}\right)+C_1\frac{\epsilon }{k}(G_{\text{S}}+G_{\text{T}})\left(1+C_3{R}_{\text{f}}\right)\\ -C_2\frac{\rho {\epsilon }^2}{k}\end{array}$ (6) where, $G_{\text{S}}={\mu }_{\text{t}}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)\frac{\partial u_i}{\partial x_j}$ (7) $G_{\text{T}}=g\beta \frac{{\mu }_{\text{t}}}{{\sigma }_{\text{t}}}\frac{\partial T}{\partial x_i}$ (8) $R_{\text{f}}=-\frac{G_{\text{T}}}{G_{\text{S}}+G_{\text{T}}}$ (9) ${\mu }_{\text{t}}=C_{\mu }\rho \frac{k^2}{\epsilon }$ (10) where $k$ is the turbulent kinetic energy (m²/s²) and $\epsilon$ is the turbulent dissipation rate (m²/s³).

These equations include many empirical constants such as ${\sigma }_{\text{k}}$, ${\sigma }_{\epsilon }$, and $C_1$. However, since virtually the same values are used for these empirical constants in a large number of reports, it would be fair to conclude that the use of empirical constants yields greater accuracy than empirically determining the eddy viscosity for each problem. The model constants of the standard k–ε and RNG k–ε models are listed in Table 1.

where, $C_1\left(\eta \right)=1.42-\frac{\eta \left(1-\eta /4.38\right)}{1+0.012{\eta }^3}$ (11) $\eta =\frac{k}{\epsilon }S$ (12) $\text{S}={\left\{\frac{1}{2}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)\right\}}^{\frac{1}{2}}$ (13)

Establishment of a three-dimensional model of a marine reactor

Based on the practical operational experience and design characteristics of a different marine reactor types, a pressurized water reactor was selected as the marine reactor in this study, and a simple three-dimensional model of the marine reactor was established. Radionuclide will be released from the reactor after the LOCA of the marine reactor. It is assumed that the containment of the third protective barrier of the reactor has a breach failure and that the radionuclide diffuses to other spaces from the breach. Four specific areas were selected to study the influence of radionuclides on staff and the pollution of the air environment in the cabin. The distribution of the cabins in the simplified model of the marine reactor is shown in Fig. 3. Radionuclides will have different levels of radiation effects on workers or equipment. Based on the consideration of radiation safety of staff and radiation damage of instruments and equipment, the staff area (chamber 1), equipment area of the steam generator (chambers 2 and 3), and radiation monitoring area (chamber 4) are regarded as four specific areas for radionuclide diffusion impact analysis.

Fig. 3Full-screen View

Distribution of cabins in marine reactor simplified model

Severe accident calculations based on MELCOR

The MELCOR program was used to simulate the process of a serious accident, and the changes in the release share of radionuclide ¹³⁷Cs and temperature with time were obtained, as shown in Fig. 4. The temperature under the initial conditions represents the temperature trend at the containment break. These conditions were coupled to SCSTREAM as boundary conditions for the study and analysis of radionuclide diffusion in a closed environment.

Fig. 4Full-screen View

Release share of radionuclide 137Cs (a) and temperature trends over time (b)

Boundary conditions for radionuclide diffusion

The boundary conditions include physical and initial conditions when using scSTREAM to simulate radionuclide diffusion. The containment had a breach in an area of 1 m × 1 m and a distance of 12.5 m from the cabin floor. The ventilation system in the cabin was ventilated via internal circulation. In addition to the above two physical conditions, the boundary conditions also includes the leakage velocity and wind speed of the ventilation system. The setting of the initial conditions for scSTREAM includes an assumed incompressible fluid, selection of an improved k-ε diffusion model, and transient operation analysis.

Results analysis and discussion

The main research work was a comparative analysis of single variables in the study of radionuclide diffusion. A comparative analysis of the influence of radionuclide diffusion under different ventilation velocities was conducted. During cabin design, different ventilation velocities were set according to different working conditions. The most commonly used velocity was 3 m/s; the minimum ventilation velocity was 1 m/s, and the maximum was 11 m/s. The middle ventilation velocity was adjusted at any time, in combination with the actual situation. Therefore, this study is arranged according to the equal-difference sequence to study the diffusion rule of radionuclides in the cabin under different ventilation velocities. A comparative analysis of the effects of radionuclide diffusion at different leakage speeds is conducted.

At the initial stage, the diffusion of radionuclides covers the entire cabin space. In the subsequent process, the diffusion of radionuclides increases in a positively proportional manner to the release share. Therefore, only the previous time period was selected for discussion and research. In Fig. 5, 6, 7 and Fig. 9, 10, the units are represented by dimensionless mass fractions. The depth of the color represents the concentration of radionuclides that diffuse in the region.

Fig. 5Full-screen View

(Color online) Cloud chart of radionuclide diffusion: The ventilation velocity is 1m/s, and the diffusion time is 90 s (a), 150 s (b) and 210 s (c) respectively. The diffusion time is 90 s, and the ventilation velocity is 3 m/s (d), 5 m/s (e) and 7 m/s (f).

Fig. 6Full-screen View

(Color online) Cloud chart of radionuclide diffusing into chamber 1: The ventilation velocity is 1 m/s (a), 3 m/s (b), 5 m/s (c), 7 m/s (d), 9 m/s (e) and 11 m/s (f) respectively.

Fig. 7Full-screen View

(Color online) Cloud chart of radionuclide diffusing into chamber 2: The ventilation velocity is 1 m/s (a), 3 m/s (b), 5 m/s (c), 7 m/s (d), 9 m/s (e) and 11m/s (f) respectively.

Fig. 9Full-screen View

(Color online) Cloud chart of radionuclide diffusing into chamber 2: The diffusion speed is 0.25 m/s (a), 0.5 m/s (b), 0.75 m/s (c) and 1 m/s (d) respectively.

Fig. 10Full-screen View

(Color online) Cloud chart of radionuclide diffusing into chamber 4: The diffusion speed is 0.25 m/s (a), 0.5 m/s (b), 0.75 m/s (c) and 1 m/s (d) respectively.

When the ventilation velocity was 1 m/s, the time of radionuclide ¹³⁷Cs diffusion from the containment break was 90 s, 150 s, and 210 s as shown in Fig. 5(a) and (c), respectively, and this diffusion area was the outer end of the containment space. When the radionuclide diffusion time is 90 s, radionuclide diffusion cloud charts at 3 m/s, 5 m/s, and 7 m/s, as shown in Fig. 5(d), (e) and (f), respectively. Figure 5(a)-(c) indicate that the diffusion process of radionuclides is slow with an increase in diffusion time, and the diffusion phenomenon in a closed environment is not obvious and not ideal. A comparative analysis of Fig. 5(a), (d)-(f) shows that with increasing wind velocity of the ventilation system, the diffusion and migration of radionuclide ¹³⁷Cs becomes faster, and the distribution of radionuclide ¹³⁷Cs in the cabin shows a sedimentation trend caused by turbulent diffusion. The main reason for this diffusion process of the radionuclides released from the containment is the tardiness due to the low wind velocity in the ventilation system. Most radionuclides are trapped in space through steam condensation and gravity-induced deposition.

A cloud chart of radionuclide diffusion into chamber 1 under the wind velocities of different ventilation systems is shown in Fig. 6. As shown in Fig. 6, with an increase in the wind velocity of the ventilation system, the time when 137Cs began to diffuse to chamber 1 gradually decreased. This is because chamber 1 is above the radionuclide release outlet, close to the outlet, and has no barrier walls.

Figure 7 shows a cloud chart of radionuclides diffusing into chamber 2 under the wind velocity conditions of different ventilation systems. It can be seen from Fig.7 that time of radionuclide diffusion from containment external environment and diffusion into chamber 2 was 122 s, 75 s, 60 s, and 59 s, when the wind velocity of the ventilation system was 1 m/s, 3 m/s, 5 m/s, and 7 m/s, respectively, showing a trend of gradual decline. However, with further increase in the wind velocity of the ventilation system, the time of radionuclide diffusion to Chamber 2 increased.

Figure 8 (a) and (b) indicate that the time of radionuclide diffusion to chamber 1 and Chamber 4 gradually decreases with the increase in ventilation velocity, while the time of radionuclide diffusion to chamber 2 and chamber 3 decreases initially and then increases. Chambers 1 and 4 are located on the side of the diffusion outlet due to the difference in location of the four chambers, and the radionuclides are directly diffused into the corresponding environment. The position of chambers 2 and 3 is, however, different from the diffusion outlet. Eddy currents form when the acute airflow meets the wall, lengthening the time it takes for radionuclides to diffuse to chambers 2 and 3. It can be seen from Fig. 8 (c) and (d) that when the ventilation wind velocity is 1 m/s, 3 m/s, and 5 m/s, the time trend of radionuclide diffusion to chamber 1-4 is similar and shows an upward trend for the different positions of the four chambers. When the ventilation rate was 7 m/s, 9 m/s and 11 m/s, the time variation trend of radionuclide diffusion to chamber 1-4 was similar, with rising initially and then falling, which was caused by the common reason of different positions of the four chambers and the increase in ventilation velocity. Figure 8 (a) takes ventilation velocity of 3 m/s as an example for analysis. The diffusion time of radionuclides from the core to chambers 1-4 is 36 s, 75 s, 97 s, and 197 s, respectively, showing a gradual upward trend. At 2700 s, the release of radionuclides is 0, and the diffusion of radionuclides in the cabin reaches equilibrium.

Fig. 8Full-screen View

(Color online) Time-trend graph of radionuclide diffusion into different chambers at different ventilation velocity

A comparative experimental investigation of different diffusion speeds was performed under a constant wind speed of 5 m/s in the ventilation system. Figure 9 reveals that the time for the radionuclide to spread out of chamber 2 gradually decreased with an increase in the diffusion speed. The radionuclide diffusion behaviors in Fig. 9(b), (c), and (d) are similar; however the diffusion cloud distribution is quite different than that in Fig. 9(a). This is because the diffusion velocity was sprayed horizontally from the containment breach, resulting in a short gas flow stroke at the initial stage of radionuclide diffusion when the diffusion speed was low.

Figure10(b), (c) and (d) show that the phenomena of the radionuclide diffusion is similar, and in all cases, the radionuclides diffuse from the lower end of the diffusion port to chamber 4, with an increase in diffusion speed, while the diffusion behavior in Fig. 10(a) is much different and radionuclide diffuses from the upper end of the diffusion port to chamber 4. Because of the increase in diffusion speed and the internal space structure of the marine reactor, the diffusion time to chamber 4 is different, and radionuclides diffuses in the form of different vortices.

In a series of experiments under different ventilation velocity and diffusion speed conditions, it was found that the two boundary conditions caused opposite results for radionuclide diffusion behavior (Fig. 8(c), (d), and Fig. 11(b)). When the ventilation velocity was low or the diffusion speed was high, the time change of radionuclide diffusion to chamber 1-4 shows an upward trend (Fig. 8(c) and Fig. 11 (b)). However, when the ventilation speed was high or the diffusion speed was low, the time change of radionuclide diffusion to chamber 1-4 first increases and then decreases (Fig. 8(d) and Fig. 11 (b)).

Fig. 11Full-screen View

(Color online) Cloud chart of radionuclide diffusing into chamber 4 at different diffusion speeds