Migration Simulation of Radioactive Soil Particles in the Atmospheric Environment Using CFD-DEM Coupled with Empirical Formulas

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. 'e 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. 'e coupling model was applied to simulate the radioactive particle distribution and migration on the nonconfined backward facing step (NBFS). 'e 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.


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]. e 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]. erefore, 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. e 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. e empirical model was developed from the data and an updated NCRP 129 model by Loosmore [8] for shorttime 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. e 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. e 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. is 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.

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.

Radioactive Particle Emission.
is 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.

Radioactivity of Single Soil Particle.
It is assumed that the radionuclides uniformly adhere to the spherical soil particles. e radioactivity estimation equation of the radioactive particle with diameter d (μm) is as follows: where Q (Bq) is the radioactivity in the surface area by measurement soil sampled at the observation point, S soil (m 2 ) is the summation of the surface area for all soil particles, and Q par d (Bq) is the radioactivity of a single particle with diameter d in the unit area on the soil surface.

Radioactive Particle Flux.
e solution procedure of the radioactive particles emission consists of two steps: estimating the radioactive particle flux and calculating the emission particle number. e 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. e 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. e radioactive particle flux (μg m − 2 s − 1 ) is estimated by the following equation: e friction velocity u * (m · s − 1 ) can be calculated as follows [14]: 2 Science and Technology of Nuclear Installations where U(m · s − 1 ) is the mean airflow velocity, D h (m) is the hydraulic diameter of the model, v is the aerodynamic viscosity coefficient, and Re is Reynolds number based on U. e number of emission particles per second from the soil surface in the computational domain can be calculated from the flux as follows: where ρ par 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. e incorporation of the resulting forces into DEM will produce information about the motion of individual particles for the next time step [15].

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. e applicability of the proposed model must be considered under nuclear emergency, and the computational efficiency is mainly considered. e 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).
e mass conversation equation is as follows: e momentum conversation equation is as follows: where α air is the volume fraction of the airflow field; α par is the volume fraction of the radioactive particles, and α air + α par � 1; ρ air is the atmosphere density; U air is the velocity of the airflow; τ air is the stress tensor for airflow; V i par is the volume of particle i; U i par is the velocity of particle i; and β air−par 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 C D is the drag coefficient, which is calculated as a function of the particle Reynolds number (Re par ).

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. ree 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. e 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. e 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. e 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).

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 m i par is the mass of particle i; F N,ij par , F T,ij par , and F V,ij par are the normal force between particles i and j, tangential force between particles i and j, and interparticle viscous damping force, respectively; and F D,air−par and F P,air−par are the fluid-particle drag force and pressure force, respectively. e normal force is obtained by a linear spring dashpot model as follows: where K N , δ N , N ij par , D N , and U N,ij par 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. e tangential contact force is obtained by the Coulombtype friction law as follows: where K T , δ T , T ij par , D T , μ, and U T,ij par 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.

Rotational Motion Equation.
e rotational motion of the particle can be described by the following equation: where ω i par is the angular momentum of particle i; I i par is the momentum of inertia of particle i; M T,ij par is the torque of the tangential forces; and M R,ij par is the rolling friction torque. e detailed solution of equations (5)-(12) refers to Deen et al. [20].

Radioactivity of the Migration Radioactive Particles.
e continuous generation and migration of radioactive particles will produce harmful radiation to a clean environment. e notion of the radioactive particle migration flux is proposed, which is applied to estimate the radiation harmfulness due to radioactive particles migration. e 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 T 1/2 is the radioactive half-life; λ is the decay constant of the radionuclide; Q esc par d (Bq) is the radioactivity of the escaped radioactive particles from the computational domain; and N esc par is the number of escaped radioactive particles, which is obtained by simulation.
Radioactive particle migration flux Φ mig par d (Bq·m −2 ) for the unit area of the soil surface can be defined as follows: 3. Validation e distribution and migration of the radioactive soil particles are the focus of this work. e 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. erefore, the coupled model is mainly validated by simulating the motion of the particles. e 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. erefore, the particle distribution simulation from Iqbal and Rauh and particle velocity from the experiment of Test case B3 described in Ref. [21] is used for validation. e 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.5 th 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 0 th second to the 25 th second in the work by Link et al. [21] and this work. e 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 timeaveraged vertical is approximately 20% in both works.

Case Study
e 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. At the NBFS surface layer, 80 Bq·g −1 of 137 Cs 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 137 Cs particles of 60 s on the NBFS is affected by factors such as the location of the emission, near-wall airflow, and wall condition.

Computational Methods.
In this work, a structured mesh for the NBFS model is generated by ANSYS ICEM CFD 14.0. e radioactive particle emission and migration calculation are simulated by OpenFOAM 4.1. e 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. Figure 2: Illustration of the radioactive particle motion from the surface. Science and Technology of Nuclear Installations

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. e density of the soil particle with 100 μm is selected, and its diameter is 2.5 × 10 3 kg·m −3 .
Friction velocity u * 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.

Grid Independence Analysis.
e magnitude of the grids for environmental pollution and particle motion near the surface can be determined to be approximately 10 4 , 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 × 10 4 , 3.5 × 10 4 , and 4.5 × 10 4 . Figure 6 illustrates the airflow velocity and 137 Cs particle velocity at the 30 th 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 × 10 4 grid, 3.5 × 10 4 grid, and 4.5 × 10 4 grid are shown in red, blue, and yellow, respectively. Compared with the 3.5 × 10 4 grid and the 4.5 × 10 4 grid, there is a certain deviation for the airflow field calculated by the 2.5 × 10 4 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 137 Cs particles as shown in Figure 6(b). In addition, the 3.5 × 10 4 grid and 4.5 × 10 4 grid are consistent with each other. Considering the accuracy and efficiency, the 3.5 × 10 4 grid is selected for the following simulations.

Radioactive Particle Simulation.
In this work, the migration mechanism of 137 Cs particles for the NBFS is investigated using the coupled model. We analyzed two aspects: distribution and migration characteristics of 137 Cs particles. Figure 7 shows the airflow field with the center plane from 40 m to 90 m at the steady state for the NBFS. e 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 137 Cs 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 137 Cs 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 137 Cs 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 137 Cs particles from the 0 th second to the 60 th second. Figure 9(a) shows the variation of the 137 Cs 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 137 Cs particles. It can be seen that the migration flux (4.37 × 10 −6 Bq·m −2 ) of 137 Cs at the 60 th second is too small in this work conditions. ere 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 137 Cs 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 137 Cs 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. e transient velocity of particles along the positive Xdirection at the 60 th second on the lower step has been selected to analyze the mechanism of the 137 Cs particle migration because of the similar particle growth curves at different time.
e particle transient velocity along the positive X-direction at the 60 th second on the lower step is shown in Figure 10. Figure 10 Compared with the 137 Cs 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)). e dotted line (straight line) (as shown in Figure 10(b)) consists of many discrete points indicating that the transient velocity of the individual 137 Cs 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 137 Cs 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. us, the particles hardly escape from the separation bubble region and gradually accumulate near the wall.

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. e 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
e data used to support the findings of this study are included within the article. Science and Technology of Nuclear Installations 9