# Numerical Simulations as Tool to Predict Chemical and Radiological Hazardous Diffusion in Case of Nonconventional Events

Authors: J.-F. Ciparisse; A. Malizia; L. A. Poggi; O. Cenciarelli; M. Gelfusa; M. C. Carestia; D. Di Giovanni; S. Mancinelli; L. Palombi; C. Bellecci; P. Gaudio
Journal: Modelling and Simulation in Engineering
 (2016)

---

## Abstract

CFD (Computational Fluid Dynamics) simulations are widely used nowadays to predict the behaviour of fluids in pure research and in industrial applications. This approach makes it possible to get quantitatively meaningful results, often in good agreement with the experimental ones. The aim of this paper is to show how CFD calculations can help to understand the time evolution of two possible CBRNe (Chemical-Biological-Radiological-Nuclear-explosive) events: (1) hazardous dust mobilization due to the interaction between a jet of air and a metallic powder in case of a LOVA (Loss Of Vacuum Accidents) that is one of the possible accidents that can occur in experimental nuclear fusion plants; (2) toxic gas release in atmosphere. The scenario analysed in the paper has consequences similar to those expected in case of a release of dangerous substances (chemical or radioactive) in enclosed or open environment during nonconventional events (like accidents or man-made or natural disasters).


---

## Body

## 1. Introduction

Nowadays the perception of risks is completely changed if compared to that in the past. The geographical boundaries are thick such as the geopolitical ones. The social boundaries have totally disappeared leaving space to a world with different roles compared to the past. The globalization and the simplicity to access information from one side have improved technologies and quality of life but from the other side have created a higher facility to gather information to offend or create terror. The cases of the dispersion of toxic gas in the past and in the recent history [1–5], the DAESH progression in the Arabic Countries [6–10], the Ebola diffusion through the aerial transportation [11–17], and the radiological dispersion due to failures in fission or fusion devices [18–24] are all examples on how the new merge between the well-known risks and the new “dangerous ideas” has exponentially increased the uncertainties of the human health all around the world. The University of Rome Tor Vergata, counting on the synergic work between the Faculty of Engineering and the Faculty of Medicine and Surgery, is facing these new problems with innovative working approaches voted to use the expertise developed in classic scientific disciplines in an unconventional way. In this paper the authors will demonstrate how the numerical simulations, in particular CFD techniques [25–39], can be a powerful tool to predict the consequences of particular accidents giving to the decision makers the chances to better manage the phases of accidents or terroristic events, reducing the risks associated. The authors will start to face the problem of toxic dust resuspension analysing first the problem itself that one of the possible causes of this phenomenon can be the LOVA and will show how these accidents can be approached with the numerical simulations. After that the authors will demonstrate how powerful are numerical simulations to simulate a toxic gas release as a consequence of a terroristic attack. The three cases analysed in this paper will be used to demonstrate that through the numerical simulations it is possible to predict accidental releases scenario (both radiological and chemical) in enclosed environment like a tank (that can be considered as a section of an industrial plant or like an office or a meeting room with an air conditioning plant) or in larger environment like chemical plant or a subway. The paper wants to demonstrate that every single accident, like the one studied in [40–52], can be faced not only from a technical point of view but also from an emergency management point of view. It is important to merge the different point of view and knowledge in order to use the numerical simulation as predicting tool for nonconventional event (like a dangerous release of a substances both man-made or caused by an accident or a natural event) in order to facilitate the phases of emergency management and restoration of normality.The tool used to run the numerical simulation was COMSOL Multiphysics® version 5.1; the results have been analysed and are here discussed by the authors. For details about theory and settings, please refer to the COMSOL Multiphysics 5.1 User’s Guide [53].


## 2. Hazardous Dust Mobilization

### 2.1. Nuclear Fusion Devices Security

Vacuum technologies are routinely used in many industrial and research applications. In some applications the macroscopic erosion of hazardous materials (radioactive, toxic, or explosive) and consequent accumulation of dust represent a risk for the workers and the population, for example, if the dust is mobilized during a loss of vacuum inside the vacuum vessel due to an accident in the system of a nuclear fusion power device. Improvement of accident management guidelines and mitigation measures of experimental nuclear fusion plants will significantly enhance the level of safety of the fusion facilities currently designed and the foreseen ITER and DEMO concepts [19].Possible hazardous consequences from accidental dust mobilization prompted greater attention in the safety analyses of future high energy density fusion machines, like ITER or DEMO. One of the priorities that have been identified by the scientific community among the safety issues in this filed was the simulation of dust transport caused by a continuum phase ingress, water or air, into the vacuum vessel (VV) [54].The Quantum Electronics and Plasma Physics (QEP) Research Group of the University of Rome Tor Vergata realized an experimental facility, STARDUST (Small Tank for Aerosol Removal and Dust) after modified in STARDUST-U (Upgrade) [55], in order to replicate experimental conditions comparable to particular thermofluid dynamic scenario expected in plants like ITER in case of LOVA (Loss Of Vacuum Accidents).The experimental campaigns have been set at the beginning to collect data about the thermofluid dynamic behaviour in case of LOVA replications inside STARDUST and STARDUST-U. Together with these experimental campaigns the QEP research group carried out experiments to analyse the quantity of dust mobilized [56, 57] and determine the velocities direction and modules by the use of optical techniques [58, 59].The numeric tool that will be developed will be a powerful instrument to predict the resuspension of dust (or particles) in enclosed environment. It will be possible to simulate several typologies of accidents like the ones expected in the nuclear fusion plants (like ITER). ITER has guidelines to determine how severe is an accident; these are on the Generic Site Safety Report [60].The simulations, and the relative models, developed and presented in this section are essential tool to support the final aim of the work on safety and security of nuclear fusion plants that is finding a numerical multiphase model to predict the dust resuspension in case of accidents. For all simulations, an unstructured grid with wall refinement to capture the boundary layer was used. The convergence was considered to have been achieved when the maximum value of the residuals is below1
e
-
4 and show an asymptotic trend.


### 2.2. Numerical Simulations to Support Radiological Risks Prevision

The failure of a window or a breach in the metallic vessel of a subatmospheric tank is one of the major concerns about safety in future Tokamaks-based power plants. The understanding of the expansion into vacuum of the air surrounding the vessel is crucial in order to predict dust mobilization. Here, air is injected into a two-dimensional vessel, in which the initial pressure is 300 Pa. As can be seen in Figure1, a 1 cm diameter pipe allows the air to flow downwards into the vacuum chamber, and the specific mass flow rate of it varies over time according to this law:(1)
m
=
a
·
1
-
exp
⁡
-
b
·
t
,where a
=
0.00056 kg/s and b
=
0.962 s−1.Figure 1
STARDUST facility picture showing the section (“sez.A”) that is implemented in the numerical simulation software, as shown in Figures2, 3, and 4.Figure 2
Mach number with streamlines after 1 s.Figure 3
Mach number and streamlines after 5 s.Figure 4
Mach number and streamlines after 10 s.The static temperature of the air is 293.15 K and the duration of the simulation is 10 s. A fully compressible, turbulent flow was considered, and the governing equations are as follows.ConsiderContinuity Equation
(2)
∂
ρ
∂
t
+
∇
→
·
ρ
V
→
=
0
.Navier-Stokes Equation
(3)
ρ
·
∂
V
→
∂
t
+
V
→
∇
→
V
→
=
∇
→
·
-
p
-
2
3
μ
+
μ
T
∇
→
·
V
→
-
2
3
ρ
k
·
I
+
μ
+
μ
T
∇
→
V
→
+
∇
→
V
→
T
.Energy Equation
(4)
ρ
·
C
p
·
∂
T
∂
t
+
V
→
·
∇
→
T
=
∇
→
·
λ
∇
→
T
+
μ
·
∇
→
V
→
+
∇
→
V
→
T
-
2
3
·
∇
→
·
V
→
·
I
:
∇
→
V
→
+
V
→
·
∇
→
p
.Turbulent Kinetic Energy Conservation Equation
(5)
ρ
·
∂
k
∂
t
+
V
→
∇
→
k
=
∇
→
·
μ
+
μ
T
σ
k
∇
→
k
+
P
k
-
ρ
ϵ
.Turbulent Kinetic Energy Dissipation Rate Conservation Equation
(6)
ρ
·
∂
ϵ
∂
t
+
V
→
∇
→
ε
=
∇
→
·
μ
+
μ
T
σ
ϵ
∇
→
ϵ
+
C
ϵ
1
ϵ
k
P
k
-
C
ϵ
2
ρ
ϵ
2
k
.Turbulent Kinetic Energy Production Rate
(7)
P
k
=
μ
T
·
∇
→
V
→
:
∇
→
V
→
+
∇
→
V
→
T
-
2
3
·
∇
→
·
V
→
2
-
2
3
ρ
k
·
∇
→
·
V
→
.Turbulent Viscosity Equation
(8)
μ
T
=
ρ
C
μ
k
2
ϵ
.Perfect Gases State Equation
(9)
p
=
ρ
R
T
.In order to simulate this kind of flow, the High Mach Number Flow module of COMSOL was used. For further information, please see the CFD Module, High Mach Number Flow Interface, pages 282–304.Figures2, 3, and 4 show the Mach number fields with streamlines, respectively, after 1 s, 5 s, and 10 s from the air inlet in the chamber.As can be seen, a subsonic expansion of the gas occurs. The Mach number initially rises because of the mass flow rate increase as the valve opens. Then, while the vessel internal pressure rises along with the gas density, at a constant and fixed mass flow rate, the velocity decreases and therefore also the Mach number decreases. A great vortex, entrained by the jet of air, completely fills the vessel. This structure of the flow is manifestly very effective in mobilizing dust contained in the cavity.Figures5 and 6 present a comparison between numerical results obtained with the model presented and experimental results measured in STARDUST facility. The graph in Figure 5 shows the pressurization rate inside the chamber for the first 10 seconds. Please note that the experimental results are presented with error bars. The graph in Figure 6 shows the air velocity over time, calculated from the Ma number. The results are still not perfectly equal due to insufficiently fine grid and numerical viscosity effects but represent a first step towards the validation of the numerical model with the experimental data collected in STARDUST facility.Figure 5
Pressurization of the chamber in the first 10 s.Figure 6
Local air velocity inside the STARDUST facility in the first 10 s.The second case study is about tungsten powder removal by an air jet flowing inside a cavity containing the dust. Figure7 shows the physical two-dimensional domain (all lengths in mm).Figure 7
Physical domain.The air enters at 100 m/s from the left upper side and exits from the upper right side that is imposed to be at atmospheric pressure. Initially, the air is at rest at 1 atm and a layer of tungsten powder, 3 cm thick, lies at the bottom of the cavity. The diameter of the solid particles is 0.4µm and the tungsten density is 19500 kg/m3. The equations describing the behaviour of this two-phase flow are as follows.ConsiderAverage Density Equation
(10)
ρ
=
ρ
c
·
1
-
φ
d
+
ρ
d
·
φ
d
.Continuity Equation
(11)
∂
ρ
∂
t
+
∇
→
·
ρ
V
→
=
0
.Navier-Stokes Equation
(12)
ρ
·
∂
V
→
∂
t
+
V
→
∇
→
V
→
=
∇
→
·
-
p
-
2
3
μ
+
μ
T
∇
→
·
V
→
-
2
3
ρ
k
·
I
+
μ
+
μ
T
∇
→
V
→
+
∇
→
V
→
T
-
∇
→
·
ρ
·
c
d
·
1
-
c
d
·
u
→
s
l
i
p
-
D
m
d
1
-
c
d
·
∇
→
φ
d
φ
d
·
u
→
s
l
i
p
-
D
m
d
1
-
c
d
·
∇
→
φ
d
φ
d
T
+
ρ
·
g
→
.Mass Balance between the Two Phases Equation
(13)
ρ
c
-
ρ
d
·
∇
→
·
φ
d
·
1
-
c
d
·
u
→
s
l
i
p
-
D
m
d
·
∇
→
φ
d
-
∇
→
φ
d
·
u
→
d
+
ρ
c
·
∇
→
·
V
→
=
0
.Momentum Balance between the Two Phases Equation
(14)
u
→
d
=
V
→
+
1
-
c
d
·
u
→
s
l
i
p
-
D
m
d
·
∇
→
φ
d
φ
d
,
c
d
=
ρ
d
·
φ
d
ρ
,
D
m
d
=
μ
T
ρ
·
σ
T
,
3
·
C
d
·
ρ
c
4
d
d
·
u
→
s
l
i
p
·
u
→
s
l
i
p
=
-
ρ
-
ρ
d
·
-
∂
V
→
∂
t
-
V
→
·
∇
→
·
V
→
+
g
→
.Effective Viscosity Equation
(15)
μ
=
μ
c
1
-
φ
d
/
φ
m
a
x
2.5
·
φ
m
a
x
.Drag Coefficient Calculation (Schiller-Naumann Model)
(16)
C
d
=
24
R
e
p
·
1
+
0.15
·
R
e
p
0.687
,
for
R
e
p
<
1000
,
0.44
,
for
R
e
p
>
1000
.Particle-Based Reynolds Number
(17)
R
e
p
=
d
d
·
ρ
c
·
u
→
s
l
i
p
μ
.Energy Equation
(18)
ρ
·
C
p
·
∂
T
∂
t
+
V
→
·
∇
→
T
=
∇
→
·
k
∇
→
T
+
μ
·
∇
→
V
→
+
∇
→
V
→
T
-
2
3
·
∇
→
·
V
→
·
I
:
∇
→
V
→
+
V
→
·
∇
→
p
.Turbulent Kinetic Energy Conservation Equation
(19)
ρ
·
∂
k
∂
t
+
V
→
∇
→
k
=
∇
→
·
μ
+
μ
T
σ
k
∇
→
k
+
P
k
-
ρ
ϵ
.Turbulent Kinetic Energy Dissipation Rate Conservation Equation
(20)
ρ
·
∂
ϵ
∂
t
+
V
→
∇
→
ε
=
∇
→
·
μ
+
μ
T
σ
ϵ
∇
→
ϵ
+
C
ϵ
1
ϵ
k
P
k
-
C
ϵ
2
ρ
ϵ
2
k
.Turbulent Kinetic Energy Production Rate
(21)
P
k
=
μ
T
·
∇
→
V
→
:
∇
→
V
→
+
∇
→
V
→
T
-
2
3
·
∇
→
·
V
→
2
-
2
3
ρ
k
·
∇
→
·
V
→
.Turbulent Viscosity Calculation
(22)
μ
T
=
ρ
C
μ
k
2
ϵ
.Perfect Gases State Equation
(23)
p
=
ρ
R
T
.The full model takes into account the drift (or slip) speed between the air and the particle, and the Schiller-Naumann drag model was used. To simulate this flow, the multiphase module using the turbulent Euler/Euler transport model was used. For further details, see the CFD Module, Multiphase Flow, the Euler/Euler Model Interface, pages 360–371.Figures8, 9, and 10 show the dust volume fraction and the streamlines at 0 s, 0.5 s, and 1 s.Figure 8
Dust initial volume fraction.Figure 9
Dust volume fraction with streamlines after 0.5 s.Figure 10
Dust volume fraction with streamlines after 1 s.As can be seen, the vortex that forms in this driven cavity rapidly entrains the dust, so its fallout and the contamination of the surroundings of the release point of the powder have to be expected. The high mobilization speed is due to the velocity of the inlet air jet. The order of magnitude of the speed (100 m/s) is common when dust is blown by an air jet during an accident/terroristic attack. Figure11 shows the evolution in time of the mass content of dust in the cavity, normalised with respect to its initial value.Figure 11
Normalised dust content versus time.The mobilization speed has been confirmed by experiments carried out in the STARDUST-U facility of the Industrial Engineering Department of the University of Rome Tor Vergata. This kind of accident can occur in Tokamaks, when a LOVA happens or when hazardous dust is accidentally or voluntarily exposed to the wind. It is clear from Figure11 that almost the total amount of dust is resuspended in the first second after the air inlet in accordance with the peak of air velocities measured (see [19, 54, 57]) and the dust mobilized during the experimental campaigns (see [56, 58, 59]).Once the multiphase model has been developed and validated through a comparison with the experimental data, changing the boundary conditions (geometry, temperatures, flow rates, etc.), it can be used to simulate the release or mobilization of other dangerous compounds (like chemical agents or biological agents).


## 2.1. Nuclear Fusion Devices Security

Vacuum technologies are routinely used in many industrial and research applications. In some applications the macroscopic erosion of hazardous materials (radioactive, toxic, or explosive) and consequent accumulation of dust represent a risk for the workers and the population, for example, if the dust is mobilized during a loss of vacuum inside the vacuum vessel due to an accident in the system of a nuclear fusion power device. Improvement of accident management guidelines and mitigation measures of experimental nuclear fusion plants will significantly enhance the level of safety of the fusion facilities currently designed and the foreseen ITER and DEMO concepts [19].Possible hazardous consequences from accidental dust mobilization prompted greater attention in the safety analyses of future high energy density fusion machines, like ITER or DEMO. One of the priorities that have been identified by the scientific community among the safety issues in this filed was the simulation of dust transport caused by a continuum phase ingress, water or air, into the vacuum vessel (VV) [54].The Quantum Electronics and Plasma Physics (QEP) Research Group of the University of Rome Tor Vergata realized an experimental facility, STARDUST (Small Tank for Aerosol Removal and Dust) after modified in STARDUST-U (Upgrade) [55], in order to replicate experimental conditions comparable to particular thermofluid dynamic scenario expected in plants like ITER in case of LOVA (Loss Of Vacuum Accidents).The experimental campaigns have been set at the beginning to collect data about the thermofluid dynamic behaviour in case of LOVA replications inside STARDUST and STARDUST-U. Together with these experimental campaigns the QEP research group carried out experiments to analyse the quantity of dust mobilized [56, 57] and determine the velocities direction and modules by the use of optical techniques [58, 59].The numeric tool that will be developed will be a powerful instrument to predict the resuspension of dust (or particles) in enclosed environment. It will be possible to simulate several typologies of accidents like the ones expected in the nuclear fusion plants (like ITER). ITER has guidelines to determine how severe is an accident; these are on the Generic Site Safety Report [60].The simulations, and the relative models, developed and presented in this section are essential tool to support the final aim of the work on safety and security of nuclear fusion plants that is finding a numerical multiphase model to predict the dust resuspension in case of accidents. For all simulations, an unstructured grid with wall refinement to capture the boundary layer was used. The convergence was considered to have been achieved when the maximum value of the residuals is below1
e
-
4 and show an asymptotic trend.


## 2.2. Numerical Simulations to Support Radiological Risks Prevision

The failure of a window or a breach in the metallic vessel of a subatmospheric tank is one of the major concerns about safety in future Tokamaks-based power plants. The understanding of the expansion into vacuum of the air surrounding the vessel is crucial in order to predict dust mobilization. Here, air is injected into a two-dimensional vessel, in which the initial pressure is 300 Pa. As can be seen in Figure1, a 1 cm diameter pipe allows the air to flow downwards into the vacuum chamber, and the specific mass flow rate of it varies over time according to this law:(1)
m
=
a
·
1
-
exp
⁡
-
b
·
t
,where a
=
0.00056 kg/s and b
=
0.962 s−1.Figure 1
STARDUST facility picture showing the section (“sez.A”) that is implemented in the numerical simulation software, as shown in Figures2, 3, and 4.Figure 2
Mach number with streamlines after 1 s.Figure 3
Mach number and streamlines after 5 s.Figure 4
Mach number and streamlines after 10 s.The static temperature of the air is 293.15 K and the duration of the simulation is 10 s. A fully compressible, turbulent flow was considered, and the governing equations are as follows.ConsiderContinuity Equation
(2)
∂
ρ
∂
t
+
∇
→
·
ρ
V
→
=
0
.Navier-Stokes Equation
(3)
ρ
·
∂
V
→
∂
t
+
V
→
∇
→
V
→
=
∇
→
·
-
p
-
2
3
μ
+
μ
T
∇
→
·
V
→
-
2
3
ρ
k
·
I
+
μ
+
μ
T
∇
→
V
→
+
∇
→
V
→
T
.Energy Equation
(4)
ρ
·
C
p
·
∂
T
∂
t
+
V
→
·
∇
→
T
=
∇
→
·
λ
∇
→
T
+
μ
·
∇
→
V
→
+
∇
→
V
→
T
-
2
3
·
∇
→
·
V
→
·
I
:
∇
→
V
→
+
V
→
·
∇
→
p
.Turbulent Kinetic Energy Conservation Equation
(5)
ρ
·
∂
k
∂
t
+
V
→
∇
→
k
=
∇
→
·
μ
+
μ
T
σ
k
∇
→
k
+
P
k
-
ρ
ϵ
.Turbulent Kinetic Energy Dissipation Rate Conservation Equation
(6)
ρ
·
∂
ϵ
∂
t
+
V
→
∇
→
ε
=
∇
→
·
μ
+
μ
T
σ
ϵ
∇
→
ϵ
+
C
ϵ
1
ϵ
k
P
k
-
C
ϵ
2
ρ
ϵ
2
k
.Turbulent Kinetic Energy Production Rate
(7)
P
k
=
μ
T
·
∇
→
V
→
:
∇
→
V
→
+
∇
→
V
→
T
-
2
3
·
∇
→
·
V
→
2
-
2
3
ρ
k
·
∇
→
·
V
→
.Turbulent Viscosity Equation
(8)
μ
T
=
ρ
C
μ
k
2
ϵ
.Perfect Gases State Equation
(9)
p
=
ρ
R
T
.In order to simulate this kind of flow, the High Mach Number Flow module of COMSOL was used. For further information, please see the CFD Module, High Mach Number Flow Interface, pages 282–304.Figures2, 3, and 4 show the Mach number fields with streamlines, respectively, after 1 s, 5 s, and 10 s from the air inlet in the chamber.As can be seen, a subsonic expansion of the gas occurs. The Mach number initially rises because of the mass flow rate increase as the valve opens. Then, while the vessel internal pressure rises along with the gas density, at a constant and fixed mass flow rate, the velocity decreases and therefore also the Mach number decreases. A great vortex, entrained by the jet of air, completely fills the vessel. This structure of the flow is manifestly very effective in mobilizing dust contained in the cavity.Figures5 and 6 present a comparison between numerical results obtained with the model presented and experimental results measured in STARDUST facility. The graph in Figure 5 shows the pressurization rate inside the chamber for the first 10 seconds. Please note that the experimental results are presented with error bars. The graph in Figure 6 shows the air velocity over time, calculated from the Ma number. The results are still not perfectly equal due to insufficiently fine grid and numerical viscosity effects but represent a first step towards the validation of the numerical model with the experimental data collected in STARDUST facility.Figure 5
Pressurization of the chamber in the first 10 s.Figure 6
Local air velocity inside the STARDUST facility in the first 10 s.The second case study is about tungsten powder removal by an air jet flowing inside a cavity containing the dust. Figure7 shows the physical two-dimensional domain (all lengths in mm).Figure 7
Physical domain.The air enters at 100 m/s from the left upper side and exits from the upper right side that is imposed to be at atmospheric pressure. Initially, the air is at rest at 1 atm and a layer of tungsten powder, 3 cm thick, lies at the bottom of the cavity. The diameter of the solid particles is 0.4µm and the tungsten density is 19500 kg/m3. The equations describing the behaviour of this two-phase flow are as follows.ConsiderAverage Density Equation
(10)
ρ
=
ρ
c
·
1
-
φ
d
+
ρ
d
·
φ
d
.Continuity Equation
(11)
∂
ρ
∂
t
+
∇
→
·
ρ
V
→
=
0
.Navier-Stokes Equation
(12)
ρ
·
∂
V
→
∂
t
+
V
→
∇
→
V
→
=
∇
→
·
-
p
-
2
3
μ
+
μ
T
∇
→
·
V
→
-
2
3
ρ
k
·
I
+
μ
+
μ
T
∇
→
V
→
+
∇
→
V
→
T
-
∇
→
·
ρ
·
c
d
·
1
-
c
d
·
u
→
s
l
i
p
-
D
m
d
1
-
c
d
·
∇
→
φ
d
φ
d
·
u
→
s
l
i
p
-
D
m
d
1
-
c
d
·
∇
→
φ
d
φ
d
T
+
ρ
·
g
→
.Mass Balance between the Two Phases Equation
(13)
ρ
c
-
ρ
d
·
∇
→
·
φ
d
·
1
-
c
d
·
u
→
s
l
i
p
-
D
m
d
·
∇
→
φ
d
-
∇
→
φ
d
·
u
→
d
+
ρ
c
·
∇
→
·
V
→
=
0
.Momentum Balance between the Two Phases Equation
(14)
u
→
d
=
V
→
+
1
-
c
d
·
u
→
s
l
i
p
-
D
m
d
·
∇
→
φ
d
φ
d
,
c
d
=
ρ
d
·
φ
d
ρ
,
D
m
d
=
μ
T
ρ
·
σ
T
,
3
·
C
d
·
ρ
c
4
d
d
·
u
→
s
l
i
p
·
u
→
s
l
i
p
=
-
ρ
-
ρ
d
·
-
∂
V
→
∂
t
-
V
→
·
∇
→
·
V
→
+
g
→
.Effective Viscosity Equation
(15)
μ
=
μ
c
1
-
φ
d
/
φ
m
a
x
2.5
·
φ
m
a
x
.Drag Coefficient Calculation (Schiller-Naumann Model)
(16)
C
d
=
24
R
e
p
·
1
+
0.15
·
R
e
p
0.687
,
for
R
e
p
<
1000
,
0.44
,
for
R
e
p
>
1000
.Particle-Based Reynolds Number
(17)
R
e
p
=
d
d
·
ρ
c
·
u
→
s
l
i
p
μ
.Energy Equation
(18)
ρ
·
C
p
·
∂
T
∂
t
+
V
→
·
∇
→
T
=
∇
→
·
k
∇
→
T
+
μ
·
∇
→
V
→
+
∇
→
V
→
T
-
2
3
·
∇
→
·
V
→
·
I
:
∇
→
V
→
+
V
→
·
∇
→
p
.Turbulent Kinetic Energy Conservation Equation
(19)
ρ
·
∂
k
∂
t
+
V
→
∇
→
k
=
∇
→
·
μ
+
μ
T
σ
k
∇
→
k
+
P
k
-
ρ
ϵ
.Turbulent Kinetic Energy Dissipation Rate Conservation Equation
(20)
ρ
·
∂
ϵ
∂
t
+
V
→
∇
→
ε
=
∇
→
·
μ
+
μ
T
σ
ϵ
∇
→
ϵ
+
C
ϵ
1
ϵ
k
P
k
-
C
ϵ
2
ρ
ϵ
2
k
.Turbulent Kinetic Energy Production Rate
(21)
P
k
=
μ
T
·
∇
→
V
→
:
∇
→
V
→
+
∇
→
V
→
T
-
2
3
·
∇
→
·
V
→
2
-
2
3
ρ
k
·
∇
→
·
V
→
.Turbulent Viscosity Calculation
(22)
μ
T
=
ρ
C
μ
k
2
ϵ
.Perfect Gases State Equation
(23)
p
=
ρ
R
T
.The full model takes into account the drift (or slip) speed between the air and the particle, and the Schiller-Naumann drag model was used. To simulate this flow, the multiphase module using the turbulent Euler/Euler transport model was used. For further details, see the CFD Module, Multiphase Flow, the Euler/Euler Model Interface, pages 360–371.Figures8, 9, and 10 show the dust volume fraction and the streamlines at 0 s, 0.5 s, and 1 s.Figure 8
Dust initial volume fraction.Figure 9
Dust volume fraction with streamlines after 0.5 s.Figure 10
Dust volume fraction with streamlines after 1 s.As can be seen, the vortex that forms in this driven cavity rapidly entrains the dust, so its fallout and the contamination of the surroundings of the release point of the powder have to be expected. The high mobilization speed is due to the velocity of the inlet air jet. The order of magnitude of the speed (100 m/s) is common when dust is blown by an air jet during an accident/terroristic attack. Figure11 shows the evolution in time of the mass content of dust in the cavity, normalised with respect to its initial value.Figure 11
Normalised dust content versus time.The mobilization speed has been confirmed by experiments carried out in the STARDUST-U facility of the Industrial Engineering Department of the University of Rome Tor Vergata. This kind of accident can occur in Tokamaks, when a LOVA happens or when hazardous dust is accidentally or voluntarily exposed to the wind. It is clear from Figure11 that almost the total amount of dust is resuspended in the first second after the air inlet in accordance with the peak of air velocities measured (see [19, 54, 57]) and the dust mobilized during the experimental campaigns (see [56, 58, 59]).Once the multiphase model has been developed and validated through a comparison with the experimental data, changing the boundary conditions (geometry, temperatures, flow rates, etc.), it can be used to simulate the release or mobilization of other dangerous compounds (like chemical agents or biological agents).


## 3. Toxic Gas Release

The diffusion of harmful gases in the atmosphere is the result of several causes: terroristic attack, failure of a vessel in a chemical plant, gaseous exhaust of a factory, and many others [1–5]. Here, a release of chlorine was considered as a case study. Referring to Figure 10, chlorine flows upwards in a vertical pipe, blown from the right by the wind. The wind speed is 5 m/s and the chlorine pressure at the base of the pipe is 1.5 atm. On the upper and left edges of the domain, atmospheric pressure was imposed. The simulation time is 10 s. Initially, the domain is filled with air; then chlorine starts to flow and to be released in the atmosphere.The domain is two-dimensional, as shown in Figure12, and all dimensions are expressed in meters.Figure 12
Physical domain.The gas is considered as a mixture of air and chlorine. The simulation is unsteady and was carried out with COMSOL Multiphysics. The flow is turbulent, compressible, and nonisothermal. The general equations describing the behaviour of the flow are as follows.ConsiderContinuity Equation
(24)
∂
ρ
∂
t
+
∇
→
·
ρ
V
→
=
0
.
Navier-Stokes Equation
(25)
ρ
·
∂
V
→
∂
t
+
V
→
∇
→
V
→
=
∇
→
·
-
p
-
2
3
μ
+
μ
T
∇
→
·
V
→
-
2
3
ρ
k
·
I
+
μ
+
μ
T
∇
→
V
→
+
∇
→
V
→
T
.
Energy Equation
(26)
ρ
·
C
p
·
∂
T
∂
t
+
V
→
·
∇
→
T
=
∇
→
·
k
∇
→
T
+
μ
·
∇
→
V
→
+
∇
→
V
→
T
-
2
3
·
∇
→
·
V
→
·
I
:
∇
→
V
→
+
V
→
·
∇
→
p
.
Turbulent Kinetic Energy Conservation Equation
(27)
ρ
·
∂
k
∂
t
+
V
→
∇
→
k
=
∇
→
·
μ
+
μ
T
σ
k
∇
→
k
+
P
k
-
ρ
ϵ
.
Turbulent Kinetic Energy Dissipation Rate Conservation Equation
(28)
ρ
·
∂
ϵ
∂
t
+
V
→
∇
→
ε
=
∇
→
·
μ
+
μ
T
σ
ϵ
∇
→
ϵ
+
C
ϵ
1
ϵ
k
P
k
-
C
ϵ
2
ρ
ϵ
2
k
.
Turbulent Kinetic Energy Production Rate
(29)
P
k
=
μ
T
·
∇
→
V
→
:
∇
→
V
→
+
∇
→
V
→
T
-
2
3
·
∇
→
·
V
→
2
-
2
3
ρ
k
·
∇
→
·
V
→
.
Turbulent Viscosity Equation
(30)
μ
T
=
ρ
C
μ
k
2
ϵ
.
Perfect Gases State Equation
(31)
p
=
ρ
R
T
.
Species Conservation Equation
(32)
∇
→
·
j
→
i
+
ρ
·
∂
ω
i
∂
t
+
V
→
·
∇
→
ω
i
=
R
i
.
i
th Species Diffusive Flux
(33)
j
→
i
=
-
ρ
D
i
m
+
ν
T
S
c
T
·
∇
→
ω
i
+
ρ
ω
i
D
i
m
·
∇
→
M
n
M
n
+
D
i
T
·
∇
→
T
T
.
i
th Species Mass Molecular Diffusion Coefficient
(34)
D
i
m
=
1
-
ω
i
∑
k
≠
i
x
k
/
D
i
k
.
Mean Molar Mass
(35)
M
n
=
∑
i
ω
i
M
i
-
1
,where ω
i is the mass fraction of the ith species, x
i is its molar fraction, S
c
T is the Schmidt turbulent number, and D
i
k is the binary diffusion coefficient between ith and kth species. This simulation was performed by using the Concentrated Species Interface and the Single-Phase Flow Interface. For further details, see the CFD Module, Single-Phase Flow Interface module, Theory for Turbulent Flow Interfaces, pages 136–172, and the Chemical Reactions Module, Chemical Species Transport Interfaces, and Transport of Concentrated Species, pages 196–211.Figures13, 14, and 15 show the chlorine mass fraction with streamlines at 1 s, 5 s, and 10 s.Figure 13
Chlorine mass fraction with streamlines after 1 s.Figure 14
Chlorine mass fraction with streamlines after 5 s.Figure 15
Chlorine mass fraction with streamlines after 10 s.As can be seen, an unsteady vortices system forms behind the pipe and the chlorine plume, and this gas is rapidly mixed in the air by the eddies, so the whole zone behind the release point of the gas is contaminated.


## 4. Conclusions

Expansion into vacuum, dust mobilization, and gaseous mixing were investigated by means of numerical simulations. It was found that this approach is well suited for studying this kind of safety problems and makes it possible to predict the transient evolution of the phenomena involved in CBRNe events. A further development of this technique will be the study of accidents/terroristic attacks occurring in three-dimensional domains under realistic conditions. It is in fact well known that turbulent flows are intrinsically unsteady and three-dimensional, rather than bidimensional. So, as turbulence plays a crucial role in the flow behaviour, especially with respect to pollutants transport, a more realistic 3D modelling of the flow is to be preferred. It is also known that a 3D domain allows movements in the third dimension which are not possible in 2D. In other words, a 3D simulation fits better the real flow than a 2D one. The drawback of this approach is the greater requirement of RAM and the higher number of processors.Building synthetic models on the basis of numerical results obtained for each class of problem under several conditions (speed of wind, mass flow rate of the pollutants, and so on) is intended in the future. The approach is the same as the one used by experimental researchers to determine empirical correlations. In such a way, the models can be used by the in-field operators just by inserting as inputs the above-mentioned conditions and will give as an output the pollutants concentration, which is the required data to evaluate the safety risks.



---

*Source: 6271853-2016-01-21.xml*