The spatial and chemical composition of the particles are represented by the size distribution function n ( v , t ) , where n ( v , t ) d v is the number of particles in the particle range [ v , v + d v ] and the number of particles in the volume range v to v + d v . The collide and coagulate are described by the general particulate conservation (Gelbard and Seinfeld, 1978). ∂ n ( v , t ) ∂ t = − ∂ ∂ v [ I ( v , t ) n ( v , t ) ] + 1 2 ∫ 0 v β ( u , v − u ) n ( u , t ) n ( v − u , t ) d u − n ( v , t ) ∫ 0 ∞ β ( v , u ) n ( u , t ) d u + S ( v , t ) (1) where I ( v , t ) is the rate of change of the particle volume due to mass exchange with the liquid phase. β ( v , u ) is the coagulation coefficient for particles of volumes v and u . S is the net increase rate of particles flowing into the system.

By discretizing the aerosol particles into m sections and s components, the mass of the section l at a given time is defined as Q l : Q l ( t ) = ∑ k = 1 s Q l , k ( t ) = ∫ v l − 1 v l v n ( v , t ) d v (2) where Q l , k ( t ) is the mass of the component k in the section l . v l and v l − 1 are the upper and lower boundaries of the section l , respectively. n ( v , t ) is the distribution function of the section l ; the mass concentration change rate for the component k of the aerosol particle in the section l can be written as the following equation (Gelbard and Seinfeld, 1980). d Q l , k ( t ) d t = 1 2 ∑ i = 1 l − 1 ∑ j = 1 l − 1 [ β ¯ 1 a i , j , l Q j , k Q i + β ¯ 1 b i , j , l Q i , k Q j ] − ∑ i = 1 l − 1 [ β ¯ 2 a i , l Q i Q l , k − β ¯ 2 b i , l Q l Q i , k ] − 1 2 β ¯ 3 l , l Q l Q l , k − Q l , k ∑ i = l + 1 m β ¯ 4 i , l Q i + S ¯ l , k − Ψ ¯ l , k Q l , k (3) where d Q l , k ( t ) / d t is the aerosol mass change rate of the component k in the section l at a time t , β ¯ 1 a i , j , l , β ¯ 1 b i , j , l is agglomeration rata of particles in section i and section j to form a particles in section l , β ¯ 2 a i , l is agglomeration rata of particles in section i and section l to form a particle larger than section l , β ¯ 2 b i , l is agglomeration rata of particles in section i and section l to form a particle in section l , β ¯ 3 l , l is agglomeration rata of particles both in section l to form a particle in section l , β ¯ 4 i , l is agglomeration rata of particles in section l and section i ( i > l ). The coagulation coefficient of the traditional sectional method is shown in Table 1. S ¯ l , k is the source of the component k of the aerosol particle in section l . Ψ ¯ l , k is the aerosol removal rate of component k of the aerosol particle in section l .

Moreover, the first and second terms in the right of Eq. 3 indicate the flux of the radioactive aerosol component k into the size section l from the lower size sections by coagulating with each section. The third term represents the flux of the component k out of the section l due to the coagulation of particles within the section l . For l < m in the fourth term, the flux of the component k leaves section l by coagulation of particles within the section l and those of higher sections. The fifth term represents the source term of the component k in the section l entering the containment. The sixth term is the removal flux of the component k in the section l by a natural settlement mechanism. The detailed theory of the natural settlement mechanism is as follows.

For the calculation of the deposition term in this study, gravity, Brownian diffusion, thermophoresis, and diffusiophoresis are considered for each heat structure surface of the containment, and the removal rate of all heat structure surfaces is added up to the total removal rate (Murata et al., 1997; Humphries et al., 2017). Ψ ¯ l , k = ∑ j = 1 N s t r λ j , l Q l , k (4) where N s t r is the total number of heat structure surfaces for aerosol deposition in the control volume, and λ j , l is the deposition rate for the heat structure j for the aerosol section l . The natural deposition coefficients can be calculated by using the deposition velocity. λ j , l = A j V ( v s + v d i f f + v t h e r m + v d i f f u s i o ) (5) where A j is the area of the heat structure surface j , and V is the control volume atmosphere volume. The velocity of aerosol deposition on a heat structure surface is defined as follows.

The agglomeration kernels of the four mechanisms are proposed in the study by Humphries et al. (2017). β B = 2 π ( D i + D j ) ( γ i d i + γ j d j ) F (11) β S = ε g π 4 C s ( γ i d i + γ j d j ) 2 | v s , i − v s , j | (12)   β T 1 = π 2 ε t 120 v g C s ( γ i d i + γ j d j ) 3 (13) β T 2 = 0.04029 ρ g 1 / 4 ε t 3 / 4 μ g 5 / 4 C s ( γ i d i + γ j d j ) 2 | ρ i C i d i 2 χ i − ρ j C j d j 2 χ j | (14) where d is the aerosol particle diameter, γ is the agglomeration shape factor, F is the correction coefficient, C s is the aerosol particle sticking coefficient, and ε g is the collision efficiency, which is defined below. ε g = 1.5 [ min ( d i , d j ) ( d i + d j ) ] 2 (15) where v s is the settling velocity of the gravitational mechanism, which is defined as v s = ρ i g d i 2 C i 18 μ g χ i (16) where ρ is the particle density, g is the acceleration of gravity, μ g is the gas viscosity, χ i is the dynamic shape factor, and C is the Cunningham slip correction factor in Eq. 7. D is the particle diffusivity, which is defined as D = k T g 3 π d μ g χ C (17)

According to the traditional section method, if the coagulation coefficient is resolved at each time step, it will waste the computation time. Therefore, it is necessary to simplify the calculation process of the coagulation coefficient.

Such as β ¯ 1 a i , j , l in Table 1, it can be written as β ¯ 1 a i , j , l = ∑ n ∫ x i − 1 x j ∫ x j − 1 x j θ ( v l − 1 < u + v < v l ) u β n ( u , v ) u v ( x i − x i − 1 ) ( x j − x j − 1 ) d y d x (18) where the subscript n represents the agglomeration mechanisms. According to Eq. 11, the variable related to the aerosol particle density is separated from the double integral sign in Eq. 18. The Brownian agglomeration can be further written as β ¯ 1 a B , i , j , l = ρ j ∗ ρ j ∫ x i − 1 x j ∫ x j − 1 x j θ ( v l − 1 < u + v < v l ) u β B ∗ ( u , v ) u v ( x i − x i − 1 ) ( x j − x j − 1 ) d y d x = k B ∫ x i − 1 x j ∫ x j − 1 x j θ ( v l − 1 < u + v < v l ) u β B ∗ ( u , v ) u v ( x i − x i − 1 ) ( x j − x j − 1 ) d y d x (19) where the superscript ∗ represents the standard agglomeration coefficients and density. By analogy, the gravitational agglomeration can be written as β ¯ 1 a S , i , j , l = | ρ i ρ j ∗ ρ i ∗ ρ j ∫ x i − 1 x j ∫ x j − 1 x j θ ( v l − 1 < u + v < v l ) u β S ( 1 ) ∗ ( u , v ) u v ( x i − x i − 1 ) ( x j − x j − 1 ) d y d x + ρ j ρ j ∗ ρ j ∗ ρ j ∫ x i − 1 x j ∫ x j − 1 x j θ ( v l − 1 < u + v < v l ) u β S ( 2 ) ∗ ( u , v ) u v ( x i − x i − 1 ) ( x j − x j − 1 ) d y d x | = | k s 1 ∫ x i − 1 x j ∫ x j − 1 x j θ ( v l − 1 < u + v < v l ) u β S ( 1 ) ∗ ( u , v ) u v ( x i − x i − 1 ) ( x j − x j − 1 ) d y d x + k s 2 ∫ x i − 1 x j ∫ x j − 1 x j θ ( v l − 1 < u + v < v l ) u β S ( 1 ) ∗ ( u , v ) u v ( x i − x i − 1 ) ( x j − x j − 1 ) d y d x | (20)

The turbulent agglomeration can be written as β ¯ 1 a T 1 , i , j , l = ( ρ g ρ g ∗ ) 0.5 ρ j ∗ ρ j ∫ x i − 1 x i ∫ x j − 1 x j θ ( v l − 1 < u + v < v l ) u β T 1 ( 1 ) ∗ ( u , v ) u v ( x i − x i − 1 ) ( x j − x j − 1 ) d y d x = k t 1 ∫ x i − 1 x i ∫ x j − 1 x j θ ( v l − 1 < u + v < v l ) u β S ( 1 ) ∗ ( u , v ) u v ( x i − x i − 1 ) ( x j − x j − 1 ) d y d x (21) β ¯ 1 a T 2 , i , j , l = | ( ρ g ρ g ∗ ) 0.25 ρ i ρ j ∗ ρ i ∗ ρ j ∫ x i − 1 x i ∫ x j − 1 x j θ ( v l − 1 < u + v < v l ) u β T 2 ( 1 ) ∗ ( u , v ) u v ( x i − x i − 1 ) ( x j − x j − 1 ) d y d x + ( ρ g ρ g ∗ ) 0.25 ρ j ρ j ∗ ρ j ∗ ρ j ∫ x i − 1 x i ∫ x j − 1 x j θ ( v l − 1 < u + v < v l ) u β T 2 ( 1 ) ∗ ( u , v ) u v ( x i − x i − 1 ) ( x j − x j − 1 ) d y d x | = | k t 2 ( 1 ) ∫ x i − 1 x i ∫ x j − 1 x j θ ( v l − 1 < u + v < v l ) u β S ( 1 ) ∗ ( u , v ) u v ( x i − x i − 1 ) ( x j − x j − 1 ) d y d x + k t 2 ( 2 ) ∫ x i − 1 x i ∫ x j − 1 x j θ ( v l − 1 < u + v < v l ) u β S ( 1 ) ∗ ( u , v ) u v ( x i − x i − 1 ) ( x j − x j − 1 ) d y d x | (22)

According to the above equation, the standard agglomeration coefficients are only calculated before the iterative process. The coefficients of k B , k s 1 , k s 2 , k t 1 , k t 2 ( 1 ) , and k t 2 ( 2 ) need an update at every time step. Therefore, the aerosol model is developed based on the aforementioned work. The flow chart of the radioactive aerosol analysis is shown in Figure 1, where t is time, Δ t is the time step, l is the section number, m is the total section number, ρ p is the radioactive aerosol particle density, v l is the lower boundary of the section l , and v l + 1 is the upper boundary of the section l . Detailed information for calculating the coagulation coefficient is discussed in the study by Gelbard (1982). The radioactive aerosol density is updated at each time step. The Runge–Kutta integration methods are adopted to numerically solve the equation and update the radioactive aerosol conditions.

The STORM-SR11 test from the International Standard Problem 40 (ISP-40) is selected for aerosol model validation. The experiment consisted of two phases: the first focusing on deposition due to natural deposition mechanisms and the second on the resuspension process of aerosols under conditions of increased airflow. This study focuses on the first phase for validation.

Figure 5 shows a schematic diagram of the STORM test facility. The test section is a 5.0055 m long straight tube with an internal diameter of 63 mm made of stainless steel. In the deposition test, the supplied carrier gas and aerosol are mixed in a mixing tank and flow into the test section, where the exhaust from the test section is connected to a cleaning and filtration system. The aerosol used is tin oxide (SnO₂), and the carrier gas is a mixture of nitrogen, steam, argon, helium, and air.

The constant mass flow rate of the aerosol in the test is 3.83 × 10–4 kg/s. Assuming that its initial distribution follows a log-normal distribution, the number distribution particle size f(d) is defined as follows. f ( d ) = 1 d ⁡ ln ⁡ σ g 2 π e x p [ − 1 2 ( ln ⁡ d − ln ⁡ d g ln ⁡ σ g ) 2 ] (34) where d g = ln ⁡ d ¯ is the geometric mean diameter. ln ⁡ σ g = [ [ ln ⁡ d − ln ⁡ d g ] 2 ¯ ] 1 / 2 is the standard geometric deviation. d is the particle diameter. The aerosol distribution in the initial state is shown in Figure 6.

The time of the numerical simulation is the same as the experiment time of 9,000 s. Figure 7 shows the comparison between the simulated results and the experiment values. There are three distinct peaks, which are located at 1.02, 3.27, and 4.36 m. The maximum measured peak value is 0.292 kg/m². The peak of the measured values in the experiment is caused by the uneven surface of the test pipe joints. The simulated results and the measured values of the smooth pipes agree very well, and they all show a constant decrease in the aerosol mass along the pipe. Table 2 depicts the deposition mass of the simulated and measured values. The deposition mass changes with the number of nodes. As the number of nodes increases, the error between the simulated value and the measured value also decreases, which gave the minimum error 6.19%.

In this research, an LBLOCA accident is selected to study the particle distribution and behavior of aerosols, with a breach diameter of 50 cm and a breach location in the hot leg. To maximize the release rate of aerosols under accident conditions, all safe injection systems in the primary system are shut down. The sequence of events is shown in Table 3, the pressure vessel failing at 5,489 s.

According to the source term evaluation results after the Fukushima accident, when a serious accident occurs, the monitoring data for I¹³¹, Cs¹³⁴, and Cs¹³⁷ are the most nuclides, and they are also the widely important nuclides in the radiation evaluation. The air source term of I¹³¹ is 60–390 PBq, the air source term of Cs¹³⁴ is 15–20 PBq, and the air source term of Cs¹³⁷ is 5–50 PBq, after the accident (Bois et al., 2014; Cervone and Franzese, 2014; Koo et al., 2014). Therefore, this study focuses on the fission product CsI. Figure 8 shows the variation of the mass of CsI aerosols with time for the accident sequence. The aerosol starts to increase at 8 s due to core exposure; when the temperature exceeds the melting temperature at 1,664 s, the core starts to melt and the aerosol release rate starts to increase dramatically and reaches a peak, with a small increase at 5,489 s due to pressure vessel failure.

As the accident proceeds, the aerosol mass will peak at 2,000 and 6,710 s. Figure 9 shows the distribution of the aerosol mass at 2,000, 6,710, and 10,000 s. According to the mass distribution, the CsI aerosol particles mainly concentrated at 0.41 um. The mass of the CsI aerosol mass at 2,000 s is maximum in the containment, as shown in Figure 8, and the max aerosol mass in Figure 9 is 2 kg. When the pressure vessel failure is at 5,489 s, the CsI aerosol mass rises again and reaches a maximum at 6,710 s, in which the max aerosol particle mass is 1.43 kg. At the end of the calculation, the CsI aerosol mass drops to 1.22 kg. It can also be seen that according to Figure 9, the distribution of the CsI aerosol mass is declining with the accident process.

The natural deposition is the primary removal mechanism of aerosols. Thus, the effect of different natural deposition on aerosols is important. As shown in Figure 10, at the end of the simulation, gravity accounts for 54% of the total deposition mass, thermophoresis 28%, diffusion 3%, and Brownian diffusion 15%. Comparing the natural deposition mass of aerosols under the accident conditions, gravitational effect is more than the other three deposition mechanisms.

Different deposition mechanisms can also have an impact on aerosol particle distribution. Figure 11 shows the effect of the different mechanisms on the aerosol particle. As can be seen from the figure, the natural deposition mechanism has a certain impact on the particles between 0.01  μ m and 0.03  μ m and has an obvious effect on the particle gravity between 2  μ m and 20  μ m , especially because the excessive gravity will increase, leading to the removal of large-diameter particles.

In this section, the influence of the radioactive aerosol particle density of a seven-component problem is discussed. The detailed density of the seven components is listed in Table 4. The initial components’ mass concentration and the initial radioactive aerosol particle numbers are shown in Figure 12. The different line types represent different component mass distribution and the radioactive aerosol particle number distribution. CsI is defined as component 1, UO₂ as component 2, H₂O as component 3, Te as component 4, B₂O₃ as component 5, Cd as component 6, and Pb as component 7. Number distribution is the axis on the left. Other curves are the axis on the right.

The particle number and component mass distribution in the presence of the aerosol coagulation process have been tabulated in Figures 13, 14. The difference in the radioactive aerosol particle number distribution can be seen in Figure 13. The aerosol particle number of the fixed density peaks at 0.45  μ m , but the aerosol particle number of the changed density peaks at 0.65  μ m . Therefore, the distribution of the changed density moves to larger particles, where the aerosol particle number decreases 16% relative to the fixed density. Figure 14 depicts the effect of different densities on the seven components’^, mass distribution. The mass distribution of seven components increases as the radioactive aerosol particles coagulate into bigger aerosol particles. The aerosol particle mass of the fixed density peaks at 0.63  μ m , and the aerosol particle mass of the changed density peaks at 0.80  μ m . So we can conclude that the density has a great influence on the multicomponent aerosol particle distribution.