Numerical analysis of supersonic jet flow and dust transport induced by air ingress in a fusion reactor

Numerical analysis of supersonic jet flow and dust transport induced by air ingress in a fusion reactor 1 2 first-author 1 corresp Corresponding author, shichao.zhang@inest.cas.cn Corresponding author, shichao.zhang@inest.cas.cn shichao.zhang@inest.cas.cn 1 corresp zhen.wang@inest.cas.cn zhen.wang@inest.cas.cn zhen.wang@inest.cas.cn https://orcid.org/0000-0003-4968-9139 1 1 Institute of Nuclear Energy Safety Technology, Hefei Institutes of Physical Science, Chinese Academy of Sciences, Hefei 230031, China Institute of Nuclear Energy Safety Technology, Hefei Institutes of Physical Science, Chinese Academy of Sciences, Hefei 230031, China University of Science and Technology of China, Hefei 230026, China University of Science and Technology of China, Hefei 230026, China During a loss of vacuum accident (LOVA), the air ingress into a vacuum vessel (VV) may lead to radioactive dust resuspension, migration, and even explosion, thereby posing a great threat to the safe operation of future fusion reactors; thus, it is crucial to understand the flow characteristics and radioactive dust transport behavior induced by LOVA. However, only a few studies have identified the characteristics of the highly under-expanded jet flow at a scale of milliseconds during LOVA. Particularly, the occurrence and behavior of a Mach disk is yet to be captured in existing studies. In this study, we used a more advanced model with a finer mesh and adaptive mesh strategies to capture the Mach disk in a VV during LOVA. In detail, a computational fluid dynamics–discrete phase model one-way coupled multiphase approach was established using the computational fluid dynamics code ANSYS FLUENT and applied to the analysis during the first seconds of LOVA. The results showed that air ingress into the VV behaved like a highly free under-expanded jet at the initial stage and Mach disk was formed at ~6 ms. Moreover, the flow field dramatically changed at the position of the Mach disk. The jet core before the Mach disk had a maximum velocity of ~8 Mach with the corresponding lowest static pressure (~100 Pa) and temperature (few tens of K). The friction velocities in the lower part of the VV, which is an area of concern due to dust deposition, was generally larger than 15 m/s near the inlet region. Lastly, the crude prediction of the particle trajectories demonstrated the spiral trajectories of the dust following the air motion. Therefore, this study provided a basis for further safety analysis and accident prevention related to dust transport and explosion in future fusion reactors. Supersonic jet Radioactive dust Loss of vacuum accident Mach disk Friction velocity 1. Introduction 1 Nuclear fusion is widely used as the ultimate solution to address global energy requirements and the associated environmental issues. To achieve this, tokamak-type fusion reactors fueled by deuteron–tritium (D–T) have been proven to be the most reliable method [ 1 1 ]. In tokamaks, dust is mainly produced by the energetic plasma–wall interactions and many mechanisms are responsible for dust generation, including edge localized modes and breaking of the deposited layer [ 2 2 ]. This resulting dust particles have compositions similar to that of the plasma-facing walls, such as tungsten and beryllium, and various shapes and sizes ranging from nanometers to hundreds of microns. In addition, this dust is generally tritiated and activated, and has strong chemical reactivity and toxicity [ 2 2 , 3 3 ]. Although dust inventory in existing tokamaks is relatively small in the range of grams (e.g., JET, JT-60) or kilograms (e.g., TFTR), large amounts of dust are expected in International Thermonuclear Experimental Reactor (ITER) and fusion DEMO reactors, such as CFETR and EU DEMO, which are expected to produce approximately several hundreds of kilograms of dust [ 4 4 ]. In the case of a loss of vacuum accident (LOVA), in which air ingress into the vacuum vessel (VV) due to a failure of the windows/valves in a large penetration line of the VV [ 5 5 ], resuspension, migration, and even explosion of the dust pose great threat to the safe operation of future fusion reactors. Therefore, it is highly important to understand the flow characteristics and radioactive dust transport behavior induced by LOVA [ 6 6 - 8 8 6-8 ]. To date, there have been several studies on LOVA in fusion reactors, which can be classified into two categories: studies carried out in scaled experiments, and numerical studies for ITER and fusion DEMO reactors. For the studies carried out in scaled experiments, most of the representative experimental trials have been performed on STARDUST and STARDUST-Upgrade facilities to investigate the dust mobilization phenomena [ 9 9 - 11 11 9-11 ]. These analyses focused on the velocity of the continuous phase as the key parameter for dust mobilization using different precision grids and turbulence models, such as ZE, k-ω , and shear stress transport (SST), to analyze the behavior of the supersonic flow and verify the turbulence by comparing the simulated Mach disk position with the empirical formula [ 12 12 ]. The pressurization rates, local air velocity, temperature, and re-suspended dust fraction in different accident conditions were detected, and a 2D thermofluid dynamic model was developed using the computational fluid dynamics (CFD) commercial code FLUENT, in which the renormalization group-based k–ε model was adopted to analyze the effect of LOVA on the generated dust [ 13 13 ]. The comparison results showed a substantial agreement between the numerical and experimental results in terms of the pressurization rate and local air velocity; however, further model development is still needed to improve the numerical accuracy. In the numerical studies of ITER and fusion DEMO reactors [ 14 14 - 17 17 14-17 ], CFD code is typically used to simulate the airflow characteristics during LOVA, focusing on the wall friction velocity field, which directly influences the subsequent dust transport. Furthermore, resuspension and particle transport models based on Euler–Euler two-phase flow were incorporated to simulate dust distribution [ 18 18 , 19 19 ]. However, existing studies have mainly focused on the long-term changes in the pressure and temperature in the vacuum chamber at the scale of hundreds of seconds, and only a few studies have identified the characteristics of the highly under-expanded jet flow at the scale of milliseconds during LOVA. In addition, Mach disk formation is expected in the air jet flow field during a LOVA event; however, existing studies have failed to capture this phenomenon, resulting in the possible inaccuracy of the simulated dust transport. In this study, to understand the highly under-expanded jet flow induced by LOVA and the consequent dust migration behavior, CFD–discrete phase model (DPM) one-way coupling model based on a density-based solver was established using ANSYS FLUENT, and used in the LOVA analysis of a fusion reactor. In Sect. 2, LOVA is described and a simplified geometry model of the fusion reactor is presented. In Sect. 3, a detailed numerical model, including the governing equations, numerical methods, and solver, and boundary and initial conditions are provided. In addition, the model validation was specified. In Sect. 4, the results of the implementation of the validated model in the analysis of the supersonic jet characteristics, dust transport behavior, and formation and development of Mach disk are presented. Finally, Sect. 5 presents a summary of the study and outlook for future research. This study provides a basis for further safety analysis and accident prevention in future fusion reactors. 2. LOVA description and geometry model 1 2.1 LOVA description 2 LOVA is induced by the failure of one or more VV/cryostat penetration lines. Air ingress through the penetration line pressurizes VV, resulting in the immediate termination of the fusion reaction. This airflow behaves like a highly under-expanded jet owing to its high pressure ratio (η ≈ 100) between the interior and exterior of the VV, and a supersonic flow is expected. Consequently, even with an explosive dust mixture, the dust deposited on underlying surfaces, such as the divertor surface, re-suspends and redistributes in the VV. If the VV is overpressurized, the burst disk is triggered to relieve pressure, resulting in the dispersion of dust into the VV pressure suppression system. In addition, the dust may directly leak out through the port cell through the fracture, resulting in the potential release of these radioactive dust particles into the environment. In this study, we only analyzed on the transient phase before the burst disk was opened, focusing on the supersonic jet flows induced by LOVA and their influence on dust transport behavior. 2.2 Geometry model 2 The VV of the fusion reactor is generally designed to be a torus with a D-shaped cross section. In this study, a simplified geometry for the VV based on the European fusion power plant conceptual model C (PPCS C) was established. The complex structure of the divertor at the bottom was deemed negligible. The specific parameters of the VV are listed in Table 1 Table 1 [ 20 20 ]. The penetration line was assumed to fail in the equatorial plane of the VV and the bypass size was assumed to be 0.02 m 2 . 1 1 Items 1 1 Values 1 1 Major radius (m) 1 1 7.5 1 1 Minor radius (m) 1 1 2.5 1 1 Aspect ratio 1 1 3 1 1 Triangularity 1 1 0.47 1 1 Total volume (m 3 ) 1 1 1722 3. Numerical model and validation 1 3.1 Governing Equations 2 The numerical analysis was performed using the CFD code ANSYS FLUENT. The simulations were performed based on the DPM model Eulerian–Lagrangian method. The air phase was considered as the continuous phase, which was simulated using the Eulerian framework, while the tungsten dust was regarded as the discrete phase, which was traced by the Lagrangian method. The volume fraction of the tungsten dust was within the range of 0% to 10%; thus, the Lagrangian method was deemed suitable to track dust in air [ 21 21 ]. The discrete phase exhibits a relatively dispersed distribution, indicating the limited influence of tungsten dust on the air phase; thus, only the influence of air flow on the entrained dust was considered. Based on these ideas, the CFD–DPM one-way coupled multiphase approach was established, in which the unsteady airflow (continuum phase) was first simulated by the CFD model and the steady dust particle tracking (dispersed phase) was subsequently calculated using the DPM model to obtain the particle trajectories in the "frozen" flow field at a certain time. In the highly under-expanded air jet induced by LOVA, supersonic flow is expected. Therefore, air must be treated as a compressible fluid, whose dynamic behaviors satisfy the conservation equations of mass, momentum, and energy, respectively, as shown below. To keep the equation closed, the ideal gas equation of state was used to calculate the pressure in the computational domain. The continuity equation is defined as Eq. (1): <math display="block" id="M1"><mrow><mfrac><mrow><mo>∂</mo><msub><mi>ρ</mi><mtext>f</mtext></msub></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>+</mo><mo>∇</mo><mo>⋅</mo><mrow><mo>(</mo><mrow><msub><mi>ρ</mi><mi>f</mi></msub><mover accent="true"><mi>v</mi><mo>→</mo></mover></mrow><mo>)</mo></mrow><mo>=</mo><mn>0.</mn></mrow></math> The momentum equations is defined as Eq. (2): <math display="block" id="M2"><mrow><mfrac><mrow><mo>∂</mo><msub><mi>ρ</mi><mtext>f</mtext></msub><mover accent="true"><mi>v</mi><mo>→</mo></mover></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>+</mo><mo>∇</mo><mo>⋅</mo><mrow><mo>(</mo><mrow><msub><mi>ρ</mi><mtext>f</mtext></msub><mover accent="true"><mi>v</mi><mo>→</mo></mover><mover accent="true"><mi>v</mi><mo>→</mo></mover></mrow><mo>)</mo></mrow><mo>=</mo><mo>−</mo><mo>∇</mo><mi>p</mi><mo>+</mo><mo>∇</mo><mo>⋅</mo><mrow><mo>(</mo><mrow><mover accent="true"><mover accent="true"><mi>τ</mi><mo stretchy="true">¯</mo></mover><mo stretchy="true">¯</mo></mover></mrow><mo>)</mo></mrow><mo>+</mo><msub><mi>ρ</mi><mtext>f</mtext></msub><mover accent="true"><mi>g</mi><mo>→</mo></mover><mo>.</mo></mrow></math> The energy equation is defined as Eq. (3): <math display="block" id="M3"><mrow><mo> </mo><mfrac><mrow><mo>∂</mo><msub><mi>ρ</mi><mtext>f</mtext></msub><msub><mi>h</mi><mrow><mtext>tot</mtext></mrow></msub></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>−</mo><mfrac><mrow><mo>∂</mo><mi>p</mi></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>+</mo><mo>∇</mo><mo>⋅</mo><mrow><mo>(</mo><mrow><msub><mi>ρ</mi><mtext>f</mtext></msub><mover accent="true"><mi>v</mi><mo>→</mo></mover><msub><mi>h</mi><mrow><mtext>tot</mtext></mrow></msub></mrow><mo>)</mo></mrow><mo>=</mo><mo>∇</mo><mo>⋅</mo><mrow><mo>(</mo><mrow><mi>λ</mi><mo>∇</mo><mi>T</mi></mrow><mo>)</mo></mrow><mo>+</mo><mo>∇</mo><mo>⋅</mo><mrow><mo>(</mo><mrow><mover accent="true"><mi>v</mi><mo>→</mo></mover><mo>⋅</mo><mover accent="true"><mover accent="true"><mi>τ</mi><mo stretchy="true">¯</mo></mover><mo stretchy="true">¯</mo></mover></mrow><mo>)</mo></mrow><mo>+</mo><msub><mi>ρ</mi><mtext>f</mtext></msub><mover accent="true"><mi>g</mi><mo>→</mo></mover><mo>⋅</mo><mover accent="true"><mi>v</mi><mo>→</mo></mover><mo>.</mo></mrow></math> The ideal gas equation of state is defined as Eq. (4): <math display="block" id="M4"><mrow><mi>P</mi><mo>=</mo><mi>ρ</mi><mi>R</mi><mi>T</mi><mo>.</mo></mrow></math> The stress tensor <math id="M5"><mrow><mover><mi>τ</mi><mo>=</mo></mover></mrow></math>   is defined as Eq. (5). <math display="block" id="M6"><mrow><mover accent="true"><mover accent="true"><mi>τ</mi><mo stretchy="true">¯</mo></mover><mo stretchy="true">¯</mo></mover><mo>=</mo><mi>μ</mi><mrow><mo>[</mo><mrow><mrow><mo>(</mo><mrow><mo>∇</mo><mover accent="true"><mi>v</mi><mo>→</mo></mover><mo>+</mo><mo>∇</mo><msup><mover accent="true"><mi>v</mi><mo>→</mo></mover><mstyle mathvariant="bold" mathsize="normal"><mi>T</mi></mstyle></msup></mrow><mo>)</mo></mrow><mo>−</mo><mfrac><mn>2</mn><mn>3</mn></mfrac><mo>∇</mo><mo>⋅</mo><mover accent="true"><mi>v</mi><mo>→</mo></mover><mi>I</mi></mrow><mo>]</mo></mrow></mrow></math> where <math id="M7"><mrow><msub><mi>ρ</mi><mtext>f</mtext></msub></mrow></math> is the density of the fluid, <math id="M8"><mover accent="true"><mi>v</mi><mo>→</mo></mover></math> is the velocity vector, <math id="M9"><mrow><mover><mi>τ</mi><mo>=</mo></mover></mrow></math> is the stress tensor, p is the static pressure, <math id="M10"><mrow><msub><mi>ρ</mi><mtext>f</mtext></msub><mover accent="true"><mi>g</mi><mo>→</mo></mover></mrow></math> is the gravitational body force, <math id="M11"><mrow><msub><mi>h</mi><mrow><mtext>tot</mtext></mrow></msub></mrow></math> is the total enthalpy, μ is the thermal conductivity coefficient, <math id="M12"><mrow><mi>μ</mi><mo> </mo></mrow></math> is the molecular viscosity, and I is the unit tensor. After the simulation of the continuum phase, the trajectory of the discrete phase particles can be computed using Eq. (6). The detailed DPM model is explained in Sect. 4.2. <math display="block" id="M13"><mrow><mfrac><mrow><mtext>d</mtext><msub><mover accent="true"><mi>v</mi><mo>→</mo></mover><mtext>p</mtext></msub></mrow><mrow><mtext>d</mtext><mi>t</mi></mrow></mfrac><mo>=</mo><msub><mover accent="true"><mi>F</mi><mo>→</mo></mover><mrow><mtext>buoyancy</mtext></mrow></msub><mo>+</mo><msub><mover accent="true"><mi>F</mi><mo>→</mo></mover><mrow><mtext>drag</mtext></mrow></msub><mo>+</mo><msub><mover accent="true"><mi>F</mi><mo>→</mo></mover><mrow><mtext>gravity</mtext></mrow></msub></mrow></math> 3.2 Numerical mesh 2 The mesh sensitivity analysis was performed on three hexahedral structured mesh resolutions created by ICEM (Grid 1, Grid 2, and Grid 3), as shown in Fig. 1 Fig. 1 . There is minimal difference on the mass flow rates of the three mesh scales, as confirmed by overlap of the curves, while the pressure transient curves at the monitoring point near the inlet show the overlap of Grid 2 and Grid 3. Therefore, the results demonstrated that two million elements are acceptable for the analysis, which afforded computational accuracy and efficiency. Thus, the subsequent simulations were conducted based on this mesh size. The selected mesh with a total of 2 242 589 elements is shown in Fig. 2(a) Fig. 2(a) . As the inlet is a crucial location throughout the transient of the LOVA, the local mesh at the inlet was refined with constant elements. In addition, owing to the nature of the LOVA phenomenon, a dedicated mesh adaptivity strategy based on the gradient approach [ 22 22 ], which is recommended for problems with supersonic flows, was adopted to refine the mesh in areas with large Mach number gradient changes to capture the Mach disk. Particularly, the specific mesh adaptation is based on the normalized Mach number gradient with a refined threshold of 0.5 and coarsen threshold of 0.2, indicating that meshes with a Mach number gradient higher than 0.5 will be marked and refined, while meshes below 0.2 will be marked and coarsened. 3.3 Numerical methods 2 The density-based solver was selected for the calculation as the pressure–velocity–temperature fields are strongly coupled in supersonic flow transients. The k-ω SST turbulence model was applied to the simulation because this model was confirmed to be the most suitable model for this type of air jet [ 8 8 , 12 12 , 14 14 ]. First-order implicit unsteady transient solver was chosen with an extremely small time step, such as 10 −8 s, used at the beginning. As the calculation proceeded, the time step was progressively increased manually. Convergence was ensured during the computation. In the calculation, the residual standard deviation was set to 10 −5 . Once the residual standard deviation was reached, the iteration stopped or the calculation continued until the maximum iteration per time step, which was manually adjusted to 30–300, according to the residual error. For the spatial discretization, the second-order upwind scheme was adopted for the flow equations and least squares cell-based discretization of the gradients, while the other equations used the first-order upwind scheme by default. The numerical methods used in the calculations are summarized in Table 2 Table 2 . 1 1 Items 1 1 Solution methods 1 1 Solver type 1 1 Density-based 1 1 Turbulence model 1 1 k-ω SST 1 1 Time discretization 1 1 First-order implicit unsteady transient solver Second-order upwind scheme for the flow equations 1 1 Spatial discretization 1 1 Least squares cell-based discretization for the gradients First-order upwind scheme for the turbulent kinetic energy First-order upwind scheme for the specific dissipation rate 3.4 Boundary and initial conditions 2 As a supersonic flow is expected, the static and total pressure should be prescribed at the pressure inlet [ 23 23 ]. This aspect is highly important in the calculation of compressible flow; however, only a few studies have incorporated this in their LOVA analysis. In this simulati on, the inlet region is a typical converging–diverging (CD) nozzle structure, where the outside pressure is significantly larger than the inside pressure. The ideal flow through the inlet can be considered as isentropic. Therefore, the static pressure ( <math id="M14"><mrow><msub><mi>p</mi><mtext>s</mtext></msub></mrow></math> ) can be obtained using Eq. (7) follows: <math display="block" id="M15"><mrow><msub><mi>p</mi><mtext>t</mtext></msub><mo>=</mo><msub><mi>p</mi><mtext>s</mtext></msub><msup><mrow><mo stretchy="false">(</mo><mn>1</mn><mo>+</mo><mfrac><mrow><mi>γ</mi><mo>−</mo><mn>1</mn></mrow><mn>2</mn></mfrac><msup><mrow><mo stretchy="false">(</mo><mi>M</mi><mi>a</mi><mo stretchy="false">)</mo></mrow><mn>2</mn></msup><mo stretchy="false">)</mo></mrow><mrow><mfrac><mi>γ</mi><mrow><mi>γ</mi><mo>−</mo><mn>1</mn></mrow></mfrac></mrow></msup><mo>,</mo></mrow></math> where <math id="M16"><mrow><msub><mi>p</mi><mtext>t</mtext></msub></mrow></math> is the total pressure, γ is the ratio of the specific heat of air, and Ma is the Mach number, which is the ratio of the fluid speed at the boundary to the speed of sound in the medium. As the inlet, such as the throat of a CD nozzle, has the smallest area, sonic flow is expected, thereby Ma can be considered equal to 1. The boundary conditions are listed in Table 3 Table 3 . It should be noted that there is no outlet boundary because the computation domain is closed only when the pressure in the VV exceeds the opening limit threshold values (150 kPa for rupture disks and 94 kPa for bleed line). Particularly, as the analyzed time span for this simulation is only limited to the first second of the transient, which is extremely short, it was assumed that the internal pressure cannot reach 94 kPa or 150 kPa. The no-slip boundary condition was applied to all walls. 1 1 Items 1 1 Values 1 1 Air pressure at inlet (Pa) 1 1 101 325 1 1 Air temperature at inlet (K) 1 1 298 1 1 Initial pressure in VV (Pa) 1 1 500 1 1 Initial temperature in VV (K) 1 1 573 1 1 Initial wall temperature (K) 1 1 573 1 1 Ratio of the specific heat of air 1 1 1.4 1 1 Boltzmann constant (J/K) 1 1 1.38×10 −23 1 1 Mean molecule diameter (m) 1 1 3.40×10 −10 [ 18 18 ] 1 1 Characteristic length (m) 1 1 5 1 1 Knudsen number 1 1 6.16×10 −6 The initial conditions are listed in Table 3 Table 3 . In actual conditions, the temperature distribution within the VV is not well known, while the initial pressure is as low as 10 −5 Pa during the steady-state operation of a tokamak. However, due to the limitations of the numerical model, the initial pressure in VV was assumed to be 500 Pa, in which the Knudsen number K n is sufficiently low (6.16×10 −6 <math id="M17"><mrow><mtext> </mtext><mo>≪</mo><mn>1</mn></mrow></math> ) based on Eq. (8) and the continuum assumption still holds. The initial temperature was assumed to be uniformly distributed at 573 K in VV. <math display="block" id="M18"><mrow><msub><mi>K</mi><mi>n</mi></msub><mo>=</mo><mfrac><mrow><msub><mi>k</mi><mtext>B</mtext></msub><mi>T</mi></mrow><mrow><msqrt><mn>2</mn></msqrt><mi>π</mi><msup><mi>σ</mi><mn>2</mn></msup><mi>p</mi><mi>L</mi></mrow></mfrac><mo>,</mo></mrow></math> where <math id="M19"><mrow><msub><mi>k</mi><mtext>B</mtext></msub></mrow></math> is the Boltzmann constant, T is the temperature, σ is the mean molecule diameter, p is the pressure, and L is the characteristic length. 3.5 Model validation 2 To ensure the accuracy of the model, we validated the model by comparing the results of the theoretical calculation and model simulation. In the initial stages of LOVA, the flow is choked and the mass flow rate is approximately constant as the pressure is below the critical value. Therefore, theoretically, the critical mass flow rate can be obtained when the flow through the inlet is isentropic and the fluid is an ideal gas, as follows: <math display="block" id="M20"><mrow><mtext>Critical mass flow rate</mtext><mo>=</mo><mi>A</mi><msub><mi>p</mi><mtext>t</mtext></msub><msqrt><mrow><mfrac><mi>γ</mi><mrow><mi>R</mi><msub><mi>T</mi><mtext>t</mtext></msub></mrow></mfrac><msup><mrow><mrow><mo>(</mo><mrow><mfrac><mn>2</mn><mrow><mi>γ</mi><mo>+</mo><mn>1</mn></mrow></mfrac></mrow><mo>)</mo></mrow></mrow><mrow><mfrac><mrow><mi>γ</mi><mo>+</mo><mn>1</mn></mrow><mrow><mi>γ</mi><mo>−</mo><mn>1</mn></mrow></mfrac></mrow></msup><mo>,</mo></mrow></msqrt></mrow></math> where A is the flow area, <math id="M21"><mrow><msub><mi>p</mi><mtext>t</mtext></msub></mrow></math> is the total pressure, <math id="M22"><mrow><msub><mi>T</mi><mtext>t</mtext></msub></mrow></math> is the total temperature, R is the gas constant, and <math id="M23"><mrow><mo> </mo><mi>γ</mi></mrow></math> is the ratio of the specific heat of air. Table 4 Table 4 shows that the critical mass flow rate simulated by the numerical model is slightly lower than the theoretical value with a relative error of approximately 0.8%, which is attributed to the marginal effects due to the presence of the boundary layer. In addition, the theoretical formula is derived from the ideal state, in which the mass flow rate mainly depends on the total conditions, including the total pressure <math id="M24"><mrow><msub><mi>p</mi><mtext>t</mtext></msub></mrow></math> and total temperature <math id="M25"><mrow><msub><mi>T</mi><mtext>t</mtext></msub></mrow></math> . Meanwhile, in the model simulation, no-slip conditions were applied to the wall and the effect of the boundary layer due to viscous phenomena was considered. Therefore, these results confirmed the validity of the model. 1 1 Items 1 1 Values 1 1 Theoretical calculation 1 1 4.7448 1 1 Model simulation 1 1 4.7067 1 1 Relative error 1 1 −0.8% 4. Simulation results and discussion 1 4.1 Supersonic jet flow characteristics 2 4.1.1 Flow field evolution 3 Figure 3 Figure 3 shows the evolution of the flow field during the LOVA, and the structure of the under-expanded jet and its transient development. The high-pressure air, in which the pressure ratio η between the air at the inlet and in the VV is approximately 100, expands after exiting the nozzle, and a supersonic jet core with a spherical shape leading front appears in the region from the inflow boundary, as shown in Fig. 3(a) Fig. 3(a) . According to the physical nature of a highly under-expanded jet (η > 4) [ 24 24 ], supersonic flow experiences the Prandtl–Meyer expansion to recover the ambient pressure value. The expansion waves are then reflected into compression waves at the jet boundary, which coalesce to form an intercepting shock. Once the intercepting shock wave is reflected at a large angle, a normal shock (also called a Mach disk) and reflected shock are generated, which was observed in this work. From the interception point of the Mach disk, intercepting shock, and reflected shock, which is referred to as the “triple point”, a slipstream is generated, as shown in Figs. 3(b)–(f) Figs. 3(b)–(f) . Due to the large distance of the in-board walls from the nozzle exit ( x / De ≈30), the jet behaves like a free under-expanded jet at the initial stage. Subsequently, the slipstream first impinge the walls around the center of the axis ( Fig. 3(d) Fig. 3(d) ) and flows further along the walls of the VV, ( Fig. 3(e) Fig. 3(e) ). As the average pressure in the VV increases, the airflow velocity decreases and the jet shrinks, which is reflected in the receding position and decreasing diameter of the Mach disk, as shown in Fig. 3(f) Fig. 3(f) . 4.1.2 Mass flow rate and pressure transient in VV 3 The mass flow rate remained approximately constant during the simulation, as shown in Fig. 4(a) Fig. 4(a) due to the average pressure in the VV, which is maintained below the critical value. The initial slight decrease in the mass flow rate can be attributed to the effect of the boundary layer, as discussed in Sect. 3.5. Fig. 4(b) Fig. 4(b) shows the pressure transient in VV. The pressure at monitoring point P 1 near the inlet sharply increases to ~40000 Pa, while the pressure at the other monitoring points ( P 2 , P 3 ) slowly increases to less than 550 Pa. The volume average pressure in VV increases from 500 Pa to 565.8 Pa in 1.152 s. Although this value is far from the pressure equilibrium, there is no need for further calculations as it would be time-consuming; thus, the current simulation time used in this study is sufficient. 4.1.3 Velocity field 3 Fig. 5(a) Fig. 5(a) shows the velocity distribution on the symmetry planes close to the inlet region at 427 ms, in which the maximum velocity is noted at the core of the jet. Fig. 5(b) Fig. 5(b) shows that the variations of the Mach number along the jet axis. The Mach disk is observed at 1.65 m from the inlet. In a previous study, the position of the Mach disk at steady-state conditions calculated by <math id="M26"><mrow><mtext>LMD=</mtext><mn>0.67</mn><mo>×</mo><mi>D</mi><mi>e</mi><msub><mi>η</mi><mtext>e</mtext></msub><msup><mrow/><mrow><mn>0.62</mn></mrow></msup></mrow></math> was noted to be 1.76 [ 24 24 ]. Therefore, the Mach disk positions obtained in the current and previous study agree well with a difference of less than 7%, which can be attributed to the transient effects (increasing pressure in the VV) and influence of the inboard VV wall in this study. Fig. 6 Fig. 6 (Color online) shows the streamlines at 427 ms. The initially static flow field is largely agitated by the air jet, resulting in the formation of two swirling flows in the upper and lower parts of the VV, respectively. The friction velocity of the airflow directly influences the aerodynamic forces, which is important for the subsequent dust transport. Fig. 7(a) Fig. 7(a) shows the friction velocity fields at 427 ms. The maximum friction velocity is obtained at the wall during the direct impact of the jet directly. Meanwhile, the friction velocities in the lower part of the VV, which should be considered due to dust deposition, are larger than 24 m/s near the inlet region. Figure 7(b) Figure 7(b) shows that the friction velocity at the bottom surface varies with time. The average friction velocity is greater than 10 m/s after a certain period, which is expected to have a significant effect on the dust resuspension in the first second. 4.1.4 Pressure and temperature field 3 Fig. 8 Fig. 8 shows the static pressure fields at 427 ms. The air jet with a pressure of 1 bar rapidly expands, resulting to a significant pressure drop to an extremely low pressure. Meanwhile, the pressure in the jet core is approximately 100 Pa, which is consistent with the corresponding high speed of the jet. In addition, there is a sharp increase in the pressure behind the Mach disk and the position where the jet impinges the wall has a relatively high pressure of approximately 850 Pa. The temperature distribution in the VV at 427 m is shown in Fig. 9 Fig. 9 . The pressure and temperature exhibit similar trends due to the energy conversion from internal (high pressure, high temperature) to kinetic (high speed) through jet expansion. The air jet core cooled to a few tens of Kelvin. Meanwhile, there is also sharp increase in the temperature behind the Mach disk. Because the air temperature (298 K) is lower than the initial temperature in the VV (573 K), and the air is also relatively cold where the air flows at a high velocity. 4.2 Dust transport behavior 2 4.2.1 Main assumptions 3 The DPM model implemented in Fluent uses a Lagrangian approach, in which the dispersed phase is solved by tracking a large number of particles to obtain the individual theoretical particle trajectories. In the DPM calculation, the dust particles were assumed to be spherical in shape with equal sizes of approximately 1 μm. Moreover, the chemical reactions, and breakup and coalescence of particles during LOVA were not considered. Drag, gravity, and buoyancy forces were considered to be the main factors of the dust transport behavior, while other forces, including virtual mass and pressure gradient forces, were not considered because of the large density difference between the tungsten dust particles and air. The drag force per unit particle mass can be expressed as Eq. (10), which adopted the Schiller–Naumann drag correlations [ 25 25 ]. <math display="block" id="M27"><mrow><msub><mover accent="true"><mi>F</mi><mo>→</mo></mover><mrow><mtext>drag</mtext></mrow></msub><mo>=</mo><mfrac><mrow><mn>3</mn><mi>μ</mi><msub><mi>C</mi><mtext>D</mtext></msub><mi>R</mi><mi>e</mi></mrow><mrow><mn>4</mn><msub><mi>ρ</mi><mtext>p</mtext></msub><msubsup><mi>d</mi><mtext>p</mtext><mn>2</mn></msubsup></mrow></mfrac><mrow><mo>(</mo><mrow><mover accent="true"><mi>v</mi><mo>→</mo></mover><mo>−</mo><msub><mover accent="true"><mi>v</mi><mo>→</mo></mover><mtext>p</mtext></msub></mrow><mo>)</mo></mrow></mrow></math> where, <math id="M28"><mrow><msub><mi>C</mi><mtext>D</mtext></msub><mo> </mo></mrow></math> is the drag coefficient, Re is the relative Reynolds number for the continuous and dispersed phases, <math id="M29"><mrow><msub><mi>ρ</mi><mtext>p</mtext></msub></mrow></math> is the density of the particle, and <math id="M30"><mrow><msub><mi>d</mi><mtext>p</mtext></msub></mrow></math> is the particle diameter. The gravity force and buoyancy force can then be written as Eq. (11): <math display="block" id="M31"><mrow><msub><mover accent="true"><mi>F</mi><mo>→</mo></mover><mrow><mtext>buoyancy</mtext></mrow></msub><mo>+</mo><msub><mover accent="true"><mi>F</mi><mo>→</mo></mover><mrow><mtext>gravity</mtext></mrow></msub><mo>=</mo><mfrac><mrow><mrow><mo>(</mo><mrow><msub><mi>ρ</mi><mtext>p</mtext></msub><mo>−</mo><mi>ρ</mi></mrow><mo>)</mo></mrow></mrow><mrow><msub><mi>ρ</mi><mtext>p</mtext></msub></mrow></mfrac><mover accent="true"><mi>g</mi><mo>→</mo></mover><mo>.</mo></mrow></math> 4.2.2 Particle trajectories 3 In this study, detailed simulations of actual dust transport behavior in a fusion reactor were not carried out owing to the limitations of the current model and the need for computational resources. Meanwhile, it is unnecessary to conduct the said simulation at the current stage owing to the extensive uncertainty on the amount and properties of dust in future tokamaks [ 26 26 ]. To determine whether the dust can be blown up by the set airflow, the dust particle was tracked in a "frozen" flow field at a steady state of 427 ms with a large friction velocity, as shown in Fig. 7(b) Fig. 7(b) , using the DPM model to obtain a crude prediction. In the "frozen" flow field at 427 ms calculated by compressible CFD model, it was assumed that the particles were injected at 0 m/s from the bottom of the VV with a total mass flow rate 100 kg/s. Fig. 10 Fig. 10 shows the particle trajectories in the VV in the "frozen" flow field at 427 ms. The dust particles follow the airflow, showing spiral trajectories, which are highly similar to the streamlines of the continuum phase. This can be explained by the Stokes number ( <math id="M32"><mrow><mi>S</mi><mi>t</mi><mo>=</mo><msub><mi>ρ</mi><mtext>p</mtext></msub><msubsup><mi>D</mi><mtext>p</mtext><mn>2</mn></msubsup><mo>/</mo><mo stretchy="false">(</mo><mn>18</mn><mtext>μ</mtext><msub><mi>D</mi><mn>0</mn></msub><msub><mi>U</mi><mn>0</mn></msub></mrow></math> )), which is a dimensionless number that characterizes the behavior of particles suspended in a fluid flow. For the 1 μm tungsten dust in under LOVA, St is approximately 0.05, which is less than 1, indicating the the dust particles mostly followed the fluid motion, acting as a tracer. Further, velocity and temperature of the particles are also consistent with those of the flow fields. 5. Summary and outlook 1 In this work, to understand the air jet flow induced by LOVA and the subsequent dust migration behavior, a CFD–DPM one-way coupling model with on a density-based solver was established by using ANSYS FLUENT. The airflow characteristics were simulated by a compressible CFD model based on the unsteady density-based solver and steady dust particle tracking was calculated by the DPM model to determine the particle trajectories in the "frozen" flow field at a specific time of 4 ms. The airflow characteristic simulations showed that air ingress into the VV behave like an extremely free under-expanded jet. The volume average pressure in the VV slowly increases from 500 Pa to 565.8 Pa in 1.152 s. The flow at the inlet was choked and the mass flow rate was kept relatively constant during the simulation due to the average pressure in the VV, which was below the critical value. The supersonic jet core with a spherical leading front appeared at the region from the inflow boundary at ~5 ms and underwent Prandtl–Meyer expansion. A Mach disk is formed under multiple reflections of the expansion waves, where the flow field dramatically changed. As the LOVA evolved, the air jet shrunk, as indicated by the reduction in the Mach disk position and diameter. The friction velocities in the lower part of the VV, which is an area of concern due to the dust deposition in this region, is almost always larger than 24 m/s near the inlet region during the simulations. The crude prediction of the particle trajectories demonstrated that the dust particles followed spiral trajectories that approximately followed the fluid motion and acting as a tracer. Therefore, this study provided a basis for the LOVA analysis model and improved an understanding of the LOVA accident process. Future work will focus on obtaining the results for longer transients and improving computational models, such as incorporating dust resuspension models, and adopting the two-way coupling model of the fluid and particle. In addition, the dust and hydrogen concentration distributions need to be further calculated using a reliable CFD model during LOVA to provide a guide for accident prevention by determining the locations for igniters or inert gas injection points. References J. Ongena , R. Koch , R. Wolf et al. , Magnetic-confinement fusion . Nature Physics 12 , 398 - 410 ( 2016 ). doi: 10.1038/nphys3745 http://doi.org/10.1038/nphys3745 Magnetic-confinement fusion J.P. Sharpe , D.A. Petti , H.W. Bartels , A review of dust in fusion devices: Implications for safety and operational performance . Fusion Engineering Design , 63-64 : 153 - 163 ( 2002 ). doi: 10.1016/S0920-3796(02)00191-6 http://doi.org/10.1016/S0920-3796(02)00191-6 A review of dust in fusion devices: Implications for safety and operational performance A. Rondeau , S. Peillon , A. Roynette et al. , Characterization of dust particles produced in an all-tungsten wall tokamak and potentially mobilized by airflow . Journal of Nuclear Materials 463 : 873 - 876 ( 2015 ). doi: 10.1016/j.jnucmat.2014.12.051 http://doi.org/10.1016/j.jnucmat.2014.12.051 Characterization of dust particles produced in an all-tungsten wall tokamak and potentially mobilized by airflow G. Mazzini , T. Kaliatka , M.T. Porfiri , Tritium and dust source term inventory evaluation issues in the European DEMO reactor concepts . Fusion Engineering and Design 146 : 510 - 513 ( 2019 ). doi: 10.1016/j.fusengdes.2019.01.008 http://doi.org/10.1016/j.fusengdes.2019.01.008 Tritium and dust source term inventory evaluation issues in the European DEMO reactor concepts G. Caruso , M. Nobili , L. Ferroni , Modelling of dust resuspension in tokamak devices during an air inflow event . Journal of Fusion Energy 34 ( 5 ): 1039 - 1050 ( 2015 ). doi: 10.1007/s10894-015-9921-8 http://doi.org/10.1007/s10894-015-9921-8 Modelling of dust resuspension in tokamak devices during an air inflow event Y. Wu , Z. Chen , L. Hu et al. , Identification of safety gaps for fusion demonstration reactors . Nature Energy 1 : 16154 ( 2016 ). doi: 10.1038/nenergy.2016.154 http://doi.org/10.1038/nenergy.2016.154 Identification of safety gaps for fusion demonstration reactors M. Lukacs , L.G. Williams , Nuclear safety issues for fusion power plants . Fusion Engineering and Design 150 : 111377 ( 2020 ). doi: 10.1016/j.fusengdes.2019.111377 http://doi.org/10.1016/j.fusengdes.2019.111377 Nuclear safety issues for fusion power plants Y. Xu , S. Liu , X. Ma et al. , Numerical simulation of airflow characteristics during the loss of vacuum accident of CFETR . International Journal of Hydrogen Energy 43 ( 24 ): 11160 - 11172 ( 2018 ). doi: 10.1016/j.ijhydene.2018.04.228 http://doi.org/10.1016/j.ijhydene.2018.04.228 Numerical simulation of airflow characteristics during the loss of vacuum accident of CFETR S. Paci , N. Forgione , F. Parozzi et al. , Bases for dust mobilization modelling in the light of STARDUST experiments . Nuclear Engineering & Design 235 ( 10-12 ): 1129 - 1138 ( 2005 ). doi: 10.1016/j.nucengdes.2005.01.015 http://doi.org/10.1016/j.nucengdes.2005.01.015 Bases for dust mobilization modelling in the light of STARDUST experiments L.A. Poggi , A. Malizia , J.F. Ciparisse et al. , First experimental campaign to demonstrate STARDUST-Upgrade facility diagnostics capability to investigate lova conditions . Journal of Fusion Energy , 34 ( 6 ): 1320 - 1330 ( 2015 ). doi: 10.1007/s10894-015-9964-x http://doi.org/10.1007/s10894-015-9964-x First experimental campaign to demonstrate STARDUST-Upgrade facility diagnostics capability to investigate lova conditions L.A. Poggi , A. Malizia , J.F. Ciparisse et al. , STARDUST-U experiments on fluid-dynamic conditions affecting dust mobilization during LOVAs . Journal of Instrumentation 11 ( 07 ): C07012 ( 2016 ). doi: 10.1088/1748-0221/11/07/C07012 http://doi.org/10.1088/1748-0221/11/07/C07012 STARDUST-U experiments on fluid-dynamic conditions affecting dust mobilization during LOVAs I. Lupelli , P. Gaudio , M. Gelfusa et al . Numerical study of air jet flow field during a loss of vacuum . Fusion Engineering & Design 89 ( 9-10 ): 2048 - 2052 ( 2014 ). doi: 10.1016/j.fusengdes.2014.03.064 http://doi.org/10.1016/j.fusengdes.2014.03.064 Numerical study of air jet flow field during a loss of vacuum P. Gaudio , A. Malizia , L. Lupelli , RNG k-ε modelling and mobilization experiments of loss of vacuum in small tanks for nuclear fusion safety applications . International Journal of Systems Engineering, Applications and Development . 5 ( 3 ): 287 - 305 ( 2011 ). RNG k-ε modelling and mobilization experiments of loss of vacuum in small tanks for nuclear fusion safety applications T. Gélain , A. Rondeau , S. Peillon et al. , CFD modelling of the wall friction velocity field in the ITER tokamak resulting from airflow during a loss of vacuum accident—Consequences for particle resuspension . Fusion Engineering Design 100 : 87 - 99 ( 2015 ) doi: 10.1016/j.fusengdes.2015.04.043 http://doi.org/10.1016/j.fusengdes.2015.04.043 CFD modelling of the wall friction velocity field in the ITER tokamak resulting from airflow during a loss of vacuum accident—Consequences for particle resuspension J.F. Ciparisse , R. Rossi , A. Malizia et al. , 3D Simulation of a Loss of Vacuum Accident (LOVA) in ITER (International Thermonuclear Experimental Reactor): Evaluation of static pressure, mach number, and friction velocity . Energies , 11 ( 4 ): 856 ( 2018 ). doi: 10.3390/en11040856 http://doi.org/10.3390/en11040856 3D Simulation of a Loss of Vacuum Accident (LOVA) in ITER (International Thermonuclear Experimental Reactor): Evaluation of static pressure, mach number, and friction velocity F. Feuillebois , F. Gensdarmes , T. Gelain , Particle behavior close to the lower wall of a tokamak type geometry during a loss of vacuum accident . Fusion Engineering and Design 153 , 111500 ( 2020 ). doi: 10.1016/j.fusengdes.2020.111500 http://doi.org/10.1016/j.fusengdes.2020.111500 Particle behavior close to the lower wall of a tokamak type geometry during a loss of vacuum accident Y.Y. Xu , S.L. Liu , X.M. Cheng et al. , Numerical analysis of loss of vacuum accident (LOVA) and preliminary discussion about dust resuspension for CFETR . Fusion Engineering and Design 143 : 82 - 90 ( 2019 ). doi: 10.1016/j.fusengdes.2019.03.159 http://doi.org/10.1016/j.fusengdes.2019.03.159 Numerical analysis of loss of vacuum accident (LOVA) and preliminary discussion about dust resuspension for CFETR L.L. Tong , Numerical analysis of the dust distribution during LOVA . Annals of Nuclear Energy , 87 : 454 - 461 ( 2016 ). doi: 10.1016/j.anucene.2015.09.025 http://doi.org/10.1016/j.anucene.2015.09.025 Numerical analysis of the dust distribution during LOVA T. Gelain , F. Gensdarmes , S. Peillon et al. , CFD modelling of particle resuspension in a toroidal geometry resulting from airflows during a loss of vacuum accident (LOVA) . Fusion Engineering and Design 151 , 111386 ( 2020 ). doi: 10.1016/j.fusengdes.2019.111386 http://doi.org/10.1016/j.fusengdes.2019.111386 CFD modelling of particle resuspension in a toroidal geometry resulting from airflows during a loss of vacuum accident (LOVA) J. F. Daunenhofer and J. R. Baron , Grid Adaption for the 2D Euler equations. Technical Report AIAA-85-0484 . American Institute of Aeronautics and Astronautics , 1985 . Grid Adaption for the 2D Euler equations. Technical Report AIAA-85-0484 D. Maisonnier , I. Cook , S. Pierre et al. , The European power plant conceptual study . Fusion Engineering and Design 75 , 1173 - 1179 ( 2005 ). doi: 10.1016/j.fusengdes.2005.06.095 http://doi.org/10.1016/j.fusengdes.2005.06.095 The European power plant conceptual study D. Fang , L.F. Li , J. Li et al. , Full-scale numerical study on the thermal-hydraulic characteristics of steam-water separation system in an advanced PWR UTSG. Part two: Droplets separation process . Progress in Nuclear Energy , 118 : 103139 ( 2020 ). doi: 10.1016/j.pnucene.2019.103139 http://doi.org/10.1016/j.pnucene.2019.103139 Full-scale numerical study on the thermal-hydraulic characteristics of steam-water separation system in an advanced PWR UTSG. Part two: Droplets separation process A. Zappatore , A. Froio , G.A. Spagnuolo et al. , 3D transient CFD simulation of an in-vessel loss-of-coolant accident in the EU DEMO fusion reactor . Nuclear Fusion , 60 ( 12 ): 126001 ( 2020 ). doi: 10.1088/1741-4326/abac6b http://doi.org/10.1088/1741-4326/abac6b 3D transient CFD simulation of an in-vessel loss-of-coolant accident in the EU DEMO fusion reactor E. Franquet , V. Perrier , S. Gibout et al . Free underexpanded jets in a quiescent medium: A review . Progress in Aerospace Sciences , 77 : 25 - 53 ( 2015 ). doi: 10.1016/j.paerosci.2015.06.006 http://doi.org/10.1016/j.paerosci.2015.06.006 Free underexpanded jets in a quiescent medium: A review L. Schiller , A. Naumann , A drag coefficient correlation . Z. Ver. Deutsch. Ing 77 , 318 - 320 ( 1935 ). A drag coefficient correlation R. Rossi , P. Gaudio , L. Martellucci et al . Numerical simulations of radioactive dust particle releases during a Loss Of Vacuum Accident in a nuclear fusion reactor . Fusion Engineering and Design , 163 , 112161 ( 2021 ). doi: 10.1016/j.fusengdes.2020.112161 http://doi.org/10.1016/j.fusengdes.2020.112161 Numerical simulations of radioactive dust particle releases during a Loss Of Vacuum Accident in a nuclear fusion reactor All authors contributed to the study conception and design. Material preparation, data collection and analysis were performed by Zhen Wang, Bu-Er Wang and Shi-Chao Zhang. The first draft of the manuscript was written by Zhen Wang and Bu-Er Wang, and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript. 10.1007/s41365-021-00912-z This work was supported by the National Natural Science Foundation of China (No. 51906249) , National Key R&D Program of China (No. 2019YFE0191600) and President Foundation of Hefei Institute of Physical Science, Chinese Academy of Sciences (No. YZJJ2019QN21) . 2021-01-26 2021-05-03 2021-05-22 2021-07-15 2021-07-01 2021-08-05 © China Science Publishing & Media Ltd. (Science Press), Shanghai Institute of Applied Physics, the Chinese Academy of Sciences, ChineseNuclear Society and Springer Nature Singapore Pte Ltd. 2021 This is a post-peer-review, pre-copyedit version of an article published in Nuclear Science and Techniques. The final authenticated version is available online at: http://dx.doi.org/10.1007/s41365-021-00912-z http://dx.doi.org/10.1007/s41365-021-00912-z . 2021 Nuclear Science and Techniques 07 32