Abstract
Radioactive particle migration from the soil surface is an unignorable factor for the radioactive material distribution prediction after a nuclear accident, especially for the decision support of radioactive disposal. Considering the continuous emission, collision, and reattachment of radioactive particles, a creative simulation method with a coupled model was proposed, which combines an empirical model and the CFD-DEM method, and was established to simulate the secondary emission and motion of radioactive particles. The source term of the radioactive particles is estimated by an empirical model as the input of the CFD-DEM. Regarding the characteristics of the particle motion, the spout-fluidized bed simulation by the coupled model is consistent with the referred simulation results and experimental data, which validates the correctness of this model. The coupling model was applied to simulate the radioactive particle distribution and migration on the nonconfined backward facing step (NBFS). The simulation reveals that the distribution features and migration flux of the radioactive particles can be estimated in detail by the proposed model, which can help to provide more actual information for radioactive disposal after nuclear accidents.
1. Introduction
Over the past almost 10 years since the Fukushima Daiichi Nuclear Power Plant (FDNPP) accident, the secondary emission of deposited particles, due to factors such as aerodynamic entrainment after impaction and breakage of soil particles, is the main source of radiocesium in the atmosphere [1, 2]. The migration of radioactive particles from the soil surface can cause persistent recontamination of the cleaned ground surface. It has harmful radiation for the environment and increases the decontamination cost [3]. One of the lessons learnt from Fukushima accident is that the technologies of reducing decontamination expense, including temporary storage of decontamination waste, removed soil in containers, transport of waste, volume reduction of waste, and temporary, interim, and final storage of waste, should be developed for the enormous cost of radioactive decontamination after the nuclear accident [4]. Therefore, the impact of the migration of the soil radioactive particles for the environment after the nuclear accident must be assessed.
Numerical simulation can intuitively assess the harmful radiation under the nuclear accident condition and quickly predict the region of radioactive particle migration [5]. In this work, two categories of models, empirical formulas and particle motion model [6], are coupled to calculate the source term from the surface and motion of radioactive particles, respectively, which are two key points to predict the migration radioactive particles from the ground.
The previous method as the empirical formulas estimates the radioactivity concentration of the particles in the atmosphere after the nuclear accident, which aims to predict the emission using formulas or correlations. For example, the resuspension factor proposed by Hatano and Hatano [7] is used to accurately estimate the deposit date of contaminant. The empirical model was developed from the data and an updated NCRP 129 model by Loosmore [8] for short-time prediction of the resuspension. A new resuspension model, which is named the size-resolved, one-dimensional resuspension scheme, was proposed by Ishizuka et al [9] to estimate the radioactivity concentration of the particles in the atmosphere near FDNPP because of the secondary emission after the nuclear accident. Based on the description of macroscopic properties, the empirical model can quickly evaluate the resuspension factors of radioactive particles; however, it cannot simulate the evolution process of the resuspension radioactive particles, which is necessary for the calculation of the particle distribution [6].
Recently, increasing attention is paid to the development of a numerical method for the particle motion to predict the particle distribution. The radioactive soil particles are transported in the atmospheric environment after their break-up of the particle-surface contact by airflow, and this is the phenomenon of gas-solid coupling. The coupled computational fluid dynamics (CFD) with discrete element method (DEM) is the most used model of particle motion, which can capture the complex physics of gas-solid processes. This technique has been widely applied to many fields such as fluidization, pipeline flow, and near-wall flow. Zhao and Shan [10] solved two classic geomechanics problems using the CFD-DEM model and analyzed the sand heap formed characteristics. Park et al. [11] used the CFD-DEM to simulate the liquid-gas-particle mixture in dam break and compared the results of the numerical simulation and experiments to validate the application to multiphase flow. Iqbal and Rauh [12] investigated the details of the particle motion using the CFD-DEM, which prove consistency with experimental data.
In this work, the proposed model is considered under the nuclear emergency as following: (1) the proposed model would be used to evaluate migration process of radioactive soil particles from soil surface and (2) the decontaminant area would be differentiated accurately based on evaluation by the proposed model during the middle and last stage of nuclear accident. Considering the continuous emission, collision, and reattachment of radioactive particles, the coupled model, which combines the CFD-DEM and empirical model, was proposed to simulate the radioactive particle distribution and migration in bench terrace near the nuclear power plant in China. In the proposed model, the source term of radioactive particles, which is estimated by the empirical model as the input of the CFD-DEM, can provide more actual information of the particles and solve the limitations of the empirical model. In addition, considering nuclear emergency and computational efficiency, the progress of the continuous emission of radioactive particles from the soil surface has been simplified, and the soil surface is treated as a wall instead of a multilayered bed.
2. Methodology
Radioactive particles are transported after their break-up of the particle-surface contact by the airflow. In this model (as shown in Figure 1), the emission of radioactive soil particles from the surface is accurately estimated by the radioactive particle emission model and used as the input of the CFD-DEM simulation; the particle-surface, particle-fluid, and particle-particle interactions are calculated by CFD-DEM during the process of particle migration. Because the boundary condition of emission of the radioactive soil particles is difficult to determine, coupling an empirical formula for migration with a particle tracking algorithm can quickly estimate the flux of radioactive particle emission on the soil surface under nuclear accident condition.
2.1. Radioactive Particle Emission
This model aims to describe the continuous radioactive particle emission from the ground surface due to aeolian processes, and it composes of two parts: radioactivity estimation of single soil particle and radioactive particle flux calculation.
2.1.1. Radioactivity of Single Soil Particle
It is assumed that the radionuclides uniformly adhere to the spherical soil particles. The radioactivity estimation equation of the radioactive particle with diameter (μm) is as follows:where (Bq) is the radioactivity in the surface area by measurement soil sampled at the observation point, (m2) is the summation of the surface area for all soil particles, and (Bq) is the radioactivity of a single particle with diameter in the unit area on the soil surface.
2.1.2. Radioactive Particle Flux
The solution procedure of the radioactive particles emission consists of two steps: estimating the radioactive particle flux and calculating the emission particle number.
The radioactive particle flux is obtained by a long-term steady dust flux model, which describes the relationship between the friction velocity and the dust flux [13]. For sandy loam and silty loam, this model has been successfully applied to estimate the radioactive particle concentration in the atmosphere near FDNPP [9] although the empirical formula has a certain error. Quantization of the error of this formula is not the focus of this work. The empirical formula aims to estimate the number of radioactive soil particles generated on the surface, which is only used as the input for the CFD-DEM simulation.
The radioactive particle flux () is estimated by the following equation:
The friction velocity can be calculated as follows [14]:where is the mean airflow velocity, is the hydraulic diameter of the model, is the aerodynamic viscosity coefficient, and is Reynolds number based on .
The number of emission particles per second from the soil surface in the computational domain can be calculated from the flux as follows:where is the density of soil particle.
Two-way coupling, which indicates the interaction between fluid phase and particulate phase, is considered in the following CFD-DEM method. In this method, the motion of the particulate phase is obtained by the DEM method, while the fluid phase is determined by the CFD on a computational cell scale. At each time step, DEM provides information to evaluate the porosity and volumetric fluid drag force in a computational cell, such as the positions and velocities of the radioactive soil particle. For CFD, the fluid field will be determined by using these data, which yields the fluid drag forces that act on individual particles. The incorporation of the resulting forces into DEM will produce information about the motion of individual particles for the next time step [15].
2.2. Solution of the Airflow Field
To illustrate the effect of the particulate phase on the gas phase, the notion of the gas phase volume fraction is applied to solve the airflow field, which is calculated in the computational domain. The applicability of the proposed model must be considered under nuclear emergency, and the computational efficiency is mainly considered. The advantage of the standard k-ε model from a numerical viewpoint is that it is robust and reliable, especially for external flows [16]. Based on the standard k-ε model, the airflow field in the presence of the solid phase is solved by averaged Navier–Stokes equations (5) and (6).
The mass conversation equation is as follows:
The momentum conversation equation is as follows:where is the volume fraction of the airflow field; is the volume fraction of the radioactive particles, and ; is the atmosphere density; is the velocity of the airflow; is the stress tensor for airflow; is the volume of particle ; is the velocity of particle ; and is the momentum exchange coefficient between particles and airflow.
In this study, the momentum exchange coefficient is described by the following equations [17, 18]:where is the drag coefficient, which is calculated as a function of the particle Reynolds number ().
2.3. Radioactive Particle Tracking
During the radioactive particle motion, the migration particles from the surface are governed by the airflow field and reattached on the surface or collided with other particles. Three types of forces are considered (as shown in Figure 2) in this work: fluid force (drag force and pressure force), gravity force, and contact force (normal force and tangential force). Particles smaller than approximately 100 μm in diameter are mainly considered for application of the coupling model because they are fully resuspended in the atmosphere at normal wind speeds. The larger millimeter-sized particles are not considered because they are restricted to rolling or creeping across the surface [19], and the lift force is not used in this model. The effect of all forces on the particle results in the translation and rotation of the particles, which are governed by Newton’s second law and the balance of angular momentum. The forces and torques that act on each particle during the process of particle migration are solved by the translation motion equation and rotational motion equation, respectively. In addition, the contact force is solved based on the soft-sphere model [20] (as shown on the right of Figure 2).
2.3.1. Translation Motion Equation
According to Newton’s second law of motion, the translation motion of the particles can be solved by the following equation:where is the mass of particle ; , , and are the normal force between particles and , tangential force between particles and , and interparticle viscous damping force, respectively; and and are the fluid-particle drag force and pressure force, respectively.
The normal force is obtained by a linear spring dashpot model as follows:where are the normal spring stiffness, particle overlap with the normal direction at the contact point, normal unit vector, normal damping, and normal velocity at the contact point, respectively.
The tangential contact force is obtained by the Coulomb-type friction law as follows:where are the tangential spring stiffness, particle displacement in the tangential direction, tangential unit vector, tangential damping, friction coefficient, and tangential velocity at the contact point, respectively.
2.3.2. Rotational Motion Equation
The rotational motion of the particle can be described by the following equation:where is the angular momentum of particle ; is the momentum of inertia of particle ; is the torque of the tangential forces; and is the rolling friction torque.
The detailed solution of equations (5)–(12) refers to Deen et al. [20].
2.4. Radioactivity of the Migration Radioactive Particles
The continuous generation and migration of radioactive particles will produce harmful radiation to a clean environment. The notion of the radioactive particle migration flux is proposed, which is applied to estimate the radiation harmfulness due to radioactive particles migration. The prediction of radioactivity of the migration radioactive particles can help emergency decision makers take action to reduce loss as soon as possible. Considering the radionuclide decay, the radioactivity model of the migration radioactive particles can be described by the following equation:where is the radioactive half-life; is the decay constant of the radionuclide; (Bq) is the radioactivity of the escaped radioactive particles from the computational domain; and is the number of escaped radioactive particles, which is obtained by simulation.
Radioactive particle migration flux (Bq·m−2) for the unit area of the soil surface can be defined as follows:
3. Validation
The distribution and migration of the radioactive soil particles are the focus of this work. The radioactivity of a single soil particle and the radioactive particle flux can be estimated through equations (1)–(4), and the radioactivity variability of the radioactive particles from the soil surface is mainly affected by the motion of the soil particles. Therefore, the coupled model is mainly validated by simulating the motion of the particles. The spout-fluidized bed with the same geometry and boundary conditions is used in both studies by Iqbal and Rauh [12] and Link et al. [21]. Iqbal and Rauh described the migration distribution of the particles, which is validated by the experiment of Link et al. Therefore, the particle distribution simulation from Iqbal and Rauh and particle velocity from the experiment of Link et al. are applied to validate the correctness of the proposed model.
From Figure 3, the spout-fluidized bed in Iqbal and Rauh has a volume of 0.154 m × 0.084 m, a height of 1.0 m, and a spout area of 0.022 m × 0.012 m. The domain of the model has 21 × 14 × 100 cells with a structured mesh; 44800 particles, which have a diameter of 0.004 m and a density of 2526 kg m−3, are uniformly set in the computational domain.
Test case B3 described in Ref. [21] is used for validation. The spout gas velocity is 65 m·s−1, and the background gas velocity is 3.5 m·s−1. Figure 4(a) shows the particle distribution at the 0.5th second from the work by Iqbal and Rauh (left side of the graph) and this work (right side of the graph). Figure 4(b) illustrates the time-averaged vertical velocity of the particles at z = 0.15 m from the 0th second to the 25th second in the work by Link et al. [21] and this work.
(a)
(b)
The result of two works in Figure 4(a) indicates that there are fewer particles in the upper part of the spout-fluidized bed than in the lower part, and there are few particles at the center of the spout-fluidized bed lower part because of the effect of the spout flow. Figure 4(b) illustrates that the time-averaged vertical velocity of the particles calculated by this work is consistent with that of Link et al. Although there is a slight difference due to the solver scheme and partial undeclared parameter setting (such as the air density and aerodynamic viscosity coefficient), the particle distribution and time-averaged vertical velocity of the particles are consistent, and the fractional error of the time-averaged vertical is approximately 20% in both works.
4. Case Study
The nonconfined backward facing step (NBFS) [22], which has the characteristics of the bench terrace airflow, is used as the classical test model to simulate the migration of radioactive particles from the soil surface to provide more actual information for nuclear emergency makers.
4.1. Geometry Model of the NBFS
Figure 5(a) shows the geometry of the NBFS, which includes the upper step (50 m × 17 m) and lower step (45 m × 17 m); the NBFS has an inlet (1.8 m × 17 m), an outlet (4.3 m × 17 m), and a step (2.5 m × 17 m). There are three main regions for the flow wake in NBFS: separation bubble, reattachment, and vortex shedding, as shown in Figure 5(b).
(a)
(b)
At the NBFS surface layer, 80 Bq·g−1 of 137Cs particles is measured, which corresponds to 0.000105 Bq per single particle with 100 μm, which is emitted from the ground on the NBFS. A continuous emission of 137Cs particles of 60 s on the NBFS is affected by factors such as the location of the emission, near-wall airflow, and wall condition.
4.2. Computational Methods
In this work, a structured mesh for the NBFS model is generated by ANSYS ICEM CFD 14.0. The radioactive particle emission and migration calculation are simulated by OpenFOAM 4.1. The airflow field in the NBFS and motion of radioactive particles are simulated by DPMFoam with the time step of 0.0001 s to simulate the transient migration of radioactive particles based on the standard k-ε model.
4.3. Boundary and Initial Conditions
In this work, the isothermal state of the ambient temperature (approximately 20°C) is assumed for the 1.205 kg·m−3 atmosphere density, the airflow viscosity is 1.8 × 10−5 kg m−1·s−1, and the initial inlet velocity is 2 m·s−1. A no-slip condition at the soil surface is applied. The density of the soil particle with 100 μm is selected, and its diameter is 2.5 × 103 kg·m−3.
Friction velocity can be calculated by equation (3), which is approximately 0.23 ms−1, and the emission number (∼57) of radioactive particles per second from the surface of the NBFS is estimated. In addition, 5 cm depth contaminated soil is assumed, which contains sufficiently many radioactive particles to ensure the radioactive particle emission.
4.4. Results and Discussion
4.4.1. Grid Independence Analysis
The magnitude of the grids for environmental pollution and particle motion near the surface can be determined to be approximately 104, and the nonuniform structured grids with 1.2 times expansion near the wall are used to mesh to guarantee the required grid resolution [12, 23]. To guarantee that the simulation is independent from the grids, three groups of uniform structured grids are calculated: 2.5 × 104, 3.5 × 104, and 4.5 × 104. Figure 6 illustrates the airflow velocity and 137Cs particle velocity at the 30th second, respectively. Figure 6(a) shows the airflow velocity at z = 0.5 m with the y = 8.5 m plane, and Figure 6(b) shows the particle velocity in the domain from the y = 8.5 m plane to the y = 9 m plane, where the simulations of 2.5 × 104 grid, 3.5 × 104 grid, and 4.5 × 104 grid are shown in red, blue, and yellow, respectively. Compared with the 3.5 × 104 grid and the 4.5 × 104 grid, there is a certain deviation for the airflow field calculated by the 2.5 × 104 grid because of the effect caused by the coarse grid as shown in Figure 6(a), which would cause an obvious disparities for the velocity of 137Cs particles as shown in Figure 6(b). In addition, the 3.5 × 104 grid and 4.5 × 104 grid are consistent with each other. Considering the accuracy and efficiency, the 3.5 × 104 grid is selected for the following simulations.
(a)
(b)
4.4.2. Radioactive Particle Simulation
In this work, the migration mechanism of 137Cs particles for the NBFS is investigated using the coupled model. We analyzed two aspects: distribution and migration characteristics of 137Cs particles.
Figure 7 shows the airflow field with the center plane from 40 m to 90 m at the steady state for the NBFS. The simulation indicates that the airflow velocity gradually decreases behind the NBFS because of the expansion of the geometry cross-section width. In addition, a large-scale turbulent vortex in the separation bubble (from 50 m to 58 m) near the NBFS step and the airflow velocity in the separation bubble are very small, which would cause some of the 137Cs particles trapped by the turbulent vortex. To analyze the migration characteristics of the particles, the lower step is uniformly divided into 10 segments from 50 m to 90 m to investigate the effects of turbulent vortices and airflow fields on the particle distribution and migration.
Figure 8 shows the distribution evolution of 137Cs particles in different segments at different time in the lower step, and the red points express the number of particles generated from 50 m to 54 m. Increasingly many 137Cs particles are emitted from the surface; then, the increase in number of particles generated from 50 m to 54 m is faster than that in the other segments because most particles generated at the upper step are trapped by the turbulent vortex. Finally, the particles gradually accumulate near the step wall in the separation bubble region because of the effect of the airflow field and wall. In contrast, there are fewer particles in the other segments.
Figure 9 illustrates the change in radioactivity of 137Cs particles from the 0th second to the 60th second. Figure 9(a) shows the variation of the 137Cs particle migration flux in the atmosphere due to particle emission, and Figure 9(b) shows the ratio between migration flux and total emission flux of the 137Cs particles. It can be seen that the migration flux (4.37 × 10−6 Bq·m−2) of 137Cs at the 60th second is too small in this work conditions. There are three reasons for these results: (1) the emission number of radioactive particles per second assumed from the surface of the NBFS is less; (2) most of radioactive particles are trapped by the turbulent vortex in the separation bubble, which would cause the particles hardly escape from the separation bubble region, so as to cause the separation bubble a rapid radioactive particles growth corresponding to the results of Figure 8; and (3) there is the short simulation time. However, the migration flux of the 137Cs particles in the atmosphere linearly increases with time because the number of radioactive particles increases, and in contrast Figure 9(b) indicates that the ratio between migration flux and total emission flux of 137Cs particles first quickly increases and tends to a certain value (about 0.02) as times go on in this work conditions because of the radioactive particles trapped by the turbulent vortex near the step wall.
(a)
(b)
The transient velocity of particles along the positive X-direction at the 60th second on the lower step has been selected to analyze the mechanism of the 137Cs particle migration because of the similar particle growth curves at different time. The particle transient velocity along the positive X-direction at the 60th second on the lower step is shown in Figure 10. Figure 10(a) shows the average transient velocity of 137Cs particles in different segments; Figure 10(b) shows the transient velocity of individual 137Cs particles from 50 m to 52 m.
(a)
(b)
Compared with the 137Cs particle from 50 m to 54 m, the particles from 54 m to 58 m have much faster average transient velocity, and the particles migrate along the negative X-direction due to the turbulent vortex effect. From 58 m to 78 m, the particles have larger average transient velocity along the positive X-direction further from the step, and the particles have lower average transient velocity when the airflow weakens from 78 m to 90 m (as shown in Figure 10(a)). The dotted line (straight line) (as shown in Figure 10(b)) consists of many discrete points indicating that the transient velocity of the individual 137Cs particles corresponding to different locations from 50 m to 52 m near the NBFS step in the separation bubble. Due to a large-scale turbulent vortex and small airflow velocity, most of the 137Cs particles in the separation bubble region fluctuate in velocity at approximately zero, and the particle motion is mainly a Brownian movement due to the effect of the wall and turbulent vortex. Thus, the particles hardly escape from the separation bubble region and gradually accumulate near the wall.
5. Conclusion
In this work, considering the continuous emission and migration of radioactive particles, the distribution features and migration flux of the radioactive particles can be estimated in detail by the proposed model. One of the important factors from the distribution and migration of radioactive particles is that the turbulent vortex causes accumulation near the step in the separation bubble region. In addition, there are wide variations caused by different factors such as the parameters, steps with different heights, and wind speed. However, the method to study the radioactive particle migration after a nuclear accident is the focus of this work instead of the mechanism of the coupled model. The distribution and migration of radioactive soil particles under different boundary conditions and parameters will be investigated in the next work.
Considering the economy and timeliness, the decontamination area in the complex soil surface should be accurately differentiated for various meteorological conditions, which can help nuclear decision makers to evaluate the cost and effectiveness of the radiation decontamination to provide more actual information for radioactive disposal after a nuclear accident.
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This work was supported by the Natural Science Foundation of the Anhui Higher Education Institutions of China (grant no. KJ2020A0110), the Natural Science Foundation of Anhui Province (grant no. 2008085MA23), and Informatization Project of Chinese Academy of Sciences (grant no. XXH13506-104).