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

Bu-Er Wang ? Shi-Chao Zhang ? Zhen Wang ID ? Jiang-Tao Jia ?Zhi-Bin Chen

Abstract 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,were 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.

Keywords Supersonic jet · Radioactive dust · Loss of vacuum accident · Mach disk · Friction velocity

1 Introduction

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].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]. 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,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]. 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], 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–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–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]. The pressurization rates, local air velocity,temperature, and resuspended 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]. 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–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, 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

2.1 LOVA description

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, resuspends 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

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[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 m2.

3 Numerical model and validation

3.1 Governing equations

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]. 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):

Table 1 Parameters for the simplified geometry of the VV

The momentum equations is defined as Eq. (2):

The energy equation is defined as Eq. (3):÷

The ideal gas equation of state is defined as Eq. (4):

The stress tensor τ is defined as Eq. (5).

where ρf isthe density ofthe fluid,is thevelocity vector,τis thestresstensor,p is the static pressure,ρfis the gravitational body force, htotis the total enthalpy, λ is the thermal conductivity coefficient, μ 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.

3.2 Numerical mesh

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. 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. 2a. 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], 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.

Fig. 1 (Color online) Mesh sensitivity analysis. a Mass flow rate; b pressure transient at monitoring point

Fig. 2 (Color online) Mesh used in the simulation

3.3 Numerical methods

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, 12, 14]. First-order implicit unsteady transient solver was chosen with an extremely small time step, such as 10-8s, 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.

3.4 Boundary and initial conditions

As a supersonic flow is expected, the static and total pressure should be prescribed at the pressure inlet [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 simulation, 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 (ps) can be obtained using Eq. (7) follows:where ptis 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. 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 noslip boundary condition was applied to all walls.

Table 2 Numerical methods used in the calculation

The initial conditions are listed in Table 3. In actual conditions, the temperature distribution within the VV is not well known, while the initial pressure is as low as 10-5Pa 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 Knis sufficiently low(6.16 × 10-6?1) based on Eq. (8) and the continuumassumption still holds.The initial temperature was assumed to be uniformly distributed at 573 K in VV.

Table 3 Boundary and initial conditions

where kBis 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

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:

where A is the flow area, ptis the total pressure, Ttis the total temperature,R is the gas constant,and γ is the ratio of the specific heat of air.

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,includingthe total pressure ptand total temperature Tt. 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.

Table 4 Critical mass flow rate

4 Simulation results and discussion

4.1 Supersonic jet flow characteristics

4.1.1 Flow field evolution

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. 3a. According to the physical nature of a highly under-expanded jet (η ＞4) [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 Fig.3b–f.Due to the large distance of the inboard 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 impinges the walls around the center of the axis (Fig. 3d) and flows further along the walls of the VV,(Fig. 3e).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. 3f.

4.1.2 Mass flow rate and pressure transient in VV

The mass flow rate remained approximately constant during the simulation, as shown in Fig. 4a 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. Figure 4b shows the pressure transient in VV. The pressure at monitoring point P1near the inlet sharply increases to ～40,000 Pa, while the pressure at the other monitoring points (P2, P3) slowly increases to less than 550 Pa. The volume average pressure in VV increases from 500 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

Figure 5a 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. Figure 5b shows that the variations in 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 LMD=0.67×Deηe0.62(‘‘LMD’’ is the position of the Mach disk, Deis the exit diameter of the nozzle, ηeis the ratio between the exit pressure and the ambient pressure) was noted to be 1.76 [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. 3 (Color online) Evolution of the flow field during the LOVA: a 1.7 ms; b 5 ms; c 10.1 ms; d 16.1 ms; e 427 ms; f 1.152 s

Fig. 4 (Color online) a Mass flow rate; b pressure transient

Fig. 5 a Velocity distribution on the symmetry planes close to the inlet region at 427 ms; b variations of the Mach number along the jet axis

Fig. 6 (Color online) Streamlines at the VV at 427 ms. a Vertical perspective; b horizontal perspective

Figure 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. Figure 7a 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 7b 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.

Fig. 7 (Color online) a Friction velocity fields at 427 ms; b variations of the friction velocity at the bottom surface with time

4.1.4 Pressure and temperature field

Figure 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. 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.

Fig. 8 (Color online) Static pressure fields at 427 ms: a vertical cross section; b horizontal cross section

Fig. 9 (Color online) Temperature fields in the VV at 427 ms: a vertical cross section; b horizontal cross section

4.2 Dust transport behavior

4.2.1 Main assumptions

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].

where CDis the drag coefficient, Re is the relative Reynolds number for the continuous and dispersed phases, ρpis the density of the particle,and dpis the particle diameter.

The gravity force and buoyancy force can then be written as Eq. (11):

4.2.2 Particle trajectories

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]. 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. 7b,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.Figure 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(St=ρpD2p/（18μD0U0)),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 that 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.

Fig. 10 (Color online) Particle trajectories in the ‘‘frozen’’ flow field at 427 ms. a Particle velocity magnitude; b particle temperature

5 Summary and outlook

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 densitybased 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 behaves like an extremely free underexpanded jet. The volume average pressure in the VV slowly increases from 500 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,are 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.

Author contributionsAll 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.