# Safety Analysis in Large Volume Vacuum Systems Like Tokamak: Experiments and Numerical Simulation to Analyze Vacuum Ruptures Consequences

Authors: A. Malizia; I. Lupelli; M. Richetta; M. Gelfusa; C. Bellecci; P. Gaudio
Journal: Advances in Materials Science and Engineering
 (2014)

---

## Abstract

The large volume vacuum systems are used in many industrial operations and research laboratories. Accidents in these systems should have a relevant economical and safety impact. A loss of vacuum accident (LOVA) due to a failure of the main vacuum vessel can result in a fast pressurization of the vessel and consequent mobilization dispersion of hazardous internal material through the braches. It is clear that the influence of flow fields, consequence of accidents like LOVA, on dust resuspension is a key safety issue. In order to develop this analysis an experimental facility is been developed: STARDUST. This last facility has been used to improve the knowledge about LOVA to replicate a condition more similar to appropriate operative condition like to kamaks. By the experimental data the boundary conditions have been extrapolated to give the proper input for the 2D thermofluid-dynamics numerical simulations, developed by the commercial CFD numerical code. The benchmark of numerical simulation results with the experimental ones has been used to validate and tune the 2D thermofluid-dynamics numerical model that has been developed by the authors to replicate the LOVA conditions inside STARDUST. In present work, the facility, materials, numerical model, and relevant results will be presented.


---

## Body

## 1. Introduction

Fusion power is a promising long term candidate to supply the energy needs of humanity [1]. From the safety view point nuclear fusion holds inherent and potential safety advantages over other energy sources. In the framework of the European Fusion Program, based mainly on magnetic confinement tokamak-type machines, the complete physical and technological basic demonstration of fusion was approached by the ITER1 (International Thermonuclear Experimental Reactor) engineering and conceptual design. In particular ITER will be the first challenge to demonstrate licensable fusion safety and environmental potential of fusion and thereby provide a good precedent for the safety of future fusion power reactors. ITER should be able to generate 500 MW, but the high cost of construction has led to work on a reduced scale option [2]. In magnetic confinement devices the plasma edge2 and surrounding material surfaces provide a buffer zone between the high temperature conditions in the plasma core and the normal “terrestrial’’ environment. The interaction between the plasma edge and the surrounding surfaces profoundly influences the conditions in the plasma core3 and is a key engineering issue. Robust solutions to issues of plasma-material interactions (PMIs) are required to realize a commercially attractive fusion reactor [2]. The plasma edge needs to provide good thermal insulation and prevent impurity influx from poisoning the burning plasma core. The wall has to withstand the intense heat load and particle flux from the plasma core, over months or years of operation, with little or no maintenance. The wall surface plays an important role in the recycling of hydrogen isotopes and in plasma fuelling. The approach to practical fusion reactors inevitably leads to an increase in plasma energy content, pulse duration, and cumulative run time [2]. Plasma physics effects and PMIs that are only partially observed or accessible in present day then experiments become important. Higher heat loads, more intense transient heating events (i.e., edge localized modes (ELMs)4 and disruptions5), and the predicted magnitude of plasma facing component (PFC) damage by melting and evaporation are critical issues [3].The orders of magnitude increase of the duty cycle6 in a next step device will lead to centimeter scale erosion of PFCs. This represents three to four orders of magnitude increase from present tokamaks7, a change that is much larger than the change in any of the core physics parameters needed for ignition [2]. The edge plasma and the wall are a strongly coupled system whose interactions range over an extraordinary width of scale, from eV scale atomic interactions to hundred megajoule disruptions. PMIs critically affect tokamak operation in many ways. Erosion by the plasma determines the lifetime of PFCs and creates a source of impurities8, which cool and dilute the plasma. Deposition of material onto PFCs alters their surface composition and can lead to long term accumulation of large in-vessel tritium inventories. Retention and recycling of hydrogen from PFCs affect fuelling efficiency, plasma density control, and the density of neutral hydrogen in the plasma boundary, which impacts particle and energy transport [2]. The primary driver for the interactions between core plasma, edge plasma, and wall is the power deposited or generated in the plasma core. The fraction of this power which is not radiated from the plasma core as bremsstrahlung9 or line radiation is transported across field lines to the plasma edge [2]. The plasma edge has a strong influence on the plasma core transport processes and thereby on the energy confinement time10 [4] and it plays an essential role in the transition from the low (L) to the enhanced (H) plasma confinement regime [5]. A schematic representation of the regions of the plasma and the boundary walls in a divertor tokamak is shown in Figure 1 that represents the poloidal11 cross section of a tokamak plasma (with a single magnetic null divertor configuration), illustrating the regions of the plasma and the boundary walls where important PMIs and atomic physics processes take place. [2].Figure 1
A schematic representation of the regions of the plasma and the boundary walls [2].The characteristic regions (Figure1) of the plasma inside the VV are(i)
the plasma core;
(ii)
the edge region just inside the separatrix12;
(iii)
the scrape-of-layer (SOL)13 plasma outside the separatrix;
(iv)
the divertor plasma region, which is an extension of the SOL plasma along field lines into the divertor chamber.The baffle14 structure is designed to prevent neutrals from leaving the divertor. In the private flux region below the X point15, the magnetic field surfaces are isolated from the rest of the plasma [2]. The outermost closed magnetic field surface is characterized by a zero in the poloidal magnetic field within the vessel known as an “X point.” This boundary is called the last closed flux surface (LCFS) or separatrix. Magnetic field surfaces inside the LCFS are closed, confining the plasma ions. The edge region, just inside the LCFS, contains significant levels of not fully ionized impurities and perhaps also neutral particles. Impurity line radiation and neutral particles transport some power from here to the wall. The remaining power, P
SOL, enters the region outside the LCFS either by conduction or by convection [2]. This region is known as the scrape of layer or “SOL” as here the power is rapidly “scraped of” by electron heat conduction along open field lines, which are diverted to intersect with material targets, in what is known as a “divertor”; see inset in Figure 1 [2].Cross-field transport of power continues in the SOL, either through conduction and/or through convection, and the relative rates of the cross-field versus parallel transport along magnetic field lines determine the power width of the SOL,λ
p. Since parallel transport is very rapid compared to cross-field transport rates, the SOL is very thin in the radial direction. The power e-folding width λ
p (the distance required to decrease by a factor 1/e) typically ranges between 3 and 10 mm at the outside midplane in present devices and is expected to be 5 to 30 mm wide in a device such as ITER [6]. Typically, density and temperature SOL e-folding widths are factors of two to five larger than power widths [7]. The interaction of the edge plasma with the PFCs is determined by plasma density, temperature, flows, power fluxes, and neutral fluxes and is most intense in the vicinity of the “strike point,” where the separatrix intersects with the divertor target plate (see inset in Figure 1). The plasma density and temperature determine the flux density and energy of plasma ions striking plasma wetted surfaces. These, in turn, determine the rate of physical sputtering, chemical sputtering, ion implantation, and impurity generation. In addition, the plasma conditions determine where eroded material is redeposited, and to what degree codeposition of tritium occurs. The power flow determines the level of active structural cooling required. Neutral fluxes also cause erosion and codeposition, even in areas that are not in direct contact with the plasma [2]. The wall is a major source of neutral atoms and molecules and hence of plasma ions. Outside the divertor, in the region between the SOL plasma and the vessel wall, the density of neutrals is relatively low, typically a factor of 30–300 smaller than in the divertor. However, when spatially integrated over the large chamber wall, the total flux (in particles per second) can be quite significant, in absolute number comparable to the number of recycling ions at the divertor plates. Energetic hydrogen from the plasma can be retained in the wall or return as neutrals to the plasma edge. Since the quantity of hydrogen in the wall is typically much greater than the quantity in a plasma, small relative changes in recycling of hydrogen from the wall strongly affect fuelling of individual discharges. Neutral particles coming from the wall readily cross magnetic field lines, increasing in energy as they repeatedly charge exchange into the boundary plasma until they are ionized. They then diffuse back out across field lines as ions carrying plasma power through convection. Since the divertor is likely to be opaque to recycling neutrals, this main chamber recycling goes on independently of the divertor [7, 8]. In addition to this cross-field flow, there may be additional parallel and poloidal flows in the SOL in the main chamber induced by a variety of forces; observations in JET [9, 10] and JT-60U16 [11] imply strong drift around the SOL from the outer to the inner divertor, for a B
x
∇
B direction that is downwards in the vessel [2]. The average ion leaves and reenters the plasma many times over the course of a discharge; this process is called recycling, and average recycling times are in the range of 1 to 10 ms. Typically, the flux of hydrogen isotopes from the wall is a larger global particle source into the plasma than external fuelling by gas puffing [12], pellet injection [13], or neutral beam injection [14].Wall released hydrogen (H, D, or T) can lead to an uncontrolled density rise, and high plasma performance in present machines is often only possible when the influx of hydrogenic ions and impurities from the wall is controlled by wall conditioning [2]. Hydrogen (protium) influx is undesirable as it dilutes the reactivity of deuterium-tritium plasma. Wall conditioning has been used since the early phases of fusion research to provide the necessary conditions for plasma production and high plasma performance [2, 15–23]. High temperature bake and specialized discharges are used to deplete hydrogen trapped in the wall and to either volatilize and pump out contaminants from the wall or to bind them chemically to the wall by gettering [2]. Many advances in plasma performance in the last decade have been achieved through improved wall conditioning, for example, boronization in Alcator C-Mod17 [24, 25] and DIII-D18 [26, 27], JT-60U [28], TEXTOR19 [29], START20 [30], and, most recently, the best performance has occurred on MAST21 [31] and NSTX22 [32]. Lithium conditioning in TFTR23 [33–35] produced a dramatic 64-time increase in the fusion triple product (
n
e
τ
E
T
i
) [2] (Figure 2). Liquid lithium is being explored as a potential plasma facing material in fusion reactors [36], and its use may open very attractive, stable, high beta “zero recycling’’ regimes [37]. Even when the plasma confinement does not improve, the use of wall conditioning permits for density and recycling control, greatly widening the operational space in JET [38], ASDEX-Upgrade24 [39], and during long duration (>1 min) plasma discharges on Tore Supra25 [40–42].Figure 2
A 64-time increase in the fusion triple productne
τ
ETi in TFTR with wall conditioning [2].The description above is related to the plasma interaction with walls; another important point for the safety issues concerns is the plasma interaction with the divertor26. The main function of the divertor system is to exhaust the major part of the alpha particle power as well as He and impurities from the plasma. It is the main interface component that, under normal operation, is located between the plasma and the material surfaces, it must tolerate high heat loads while at the same time providing neutron shielding for the vacuum vessel and magnet coils in the vicinity of the divertor. Although good progress has been made in the understanding of divertor plasma physics, there continue to be some uncertainties, and hence the divertor remains an experimental device, which it is anticipated will need to be replaced and upgraded several times during the life of ITER. To facilitate rapid replacement, remote maintainability of the divertor has been given a high priority [43].The main components of the divertor system are (see Figures3 and 4):(i)
a divertor cassette body, that is reusable to minimize activated waste and provides neutron shielding and a mechanical support for different possible arrangements of plasma interfaces;
(ii)
inner and outer vertical targets, which are the plasma-facing components (PFCs) which in their lower part interact directly with the scrape-of layer (SOL) plasma and in their upper part act as baffles for the neutrals;
(iii)
the private flux region (i.e., the space below the separatrix which has no flux line connections to the main plasma) PFC, which in turn consists of
(a)
a dome, located below the separatrix X-point, seeing mainly radiation and charge exchange (CX) neutrals; the dome additionally baffles neutral particles and protects the liner and the neutral particle reflector plates from the SOL plasma;
(b)
inner and outer neutral particle reflector plates that together with the lower ends of the vertical targets form a “V” shape that confines particles in the divertor channels to aid in reduction of peak heat flux by encouraging partial plasma detachment from the plate;
(c)
a semitransparent liner that protects the cassette body from direct line of sight of the plasma while allowing He and other impurities to be pumped away;
(iv)
support pads integrated into the cassette to provide locking and alignment of the divertor cassettes on the rails;
(v)
divertor to VV gas seals, to prevent backstreaming of gas from the divertor into the main plasma chamber;
(vi)
cooling pipe interfaces connecting the divertor cassettes to the radial cooling pipes at each divertor port;
(vii)
special diagnostic cassettes providing access for diagnostics;
(viii)
rails supporting the cassettes, part of the vacuum vessel.Three PFCs (inner and outer vertical targets and a private flux region assembly) are mounted on each cassette body in the hot cell by special semiautomatic tools [43].Figure 3
Cross section through the divertor. This shows the V-shaped strike zones and the pumping path as well as the attachment to the toroidal rails [43].Figure 4
Divertor cassette showing the large openings in the dome support. Also the transparent liner is shown allowing for pumping access from the private flux region [43].We define the divertor zone to be the volume below the X point, which also usually means below the divertor baffle structure; see Figure1. We assume that plasma ions strike the divertor plate, are recycled back into the plasma either as atoms or molecules, and are ionized/dissociated primarily in the divertor, which is not directly penetrating the confined plasma or the SOL outside the divertor. Neutrals are either prevented from leaving the divertor by the opacity of the divertor plasma itself (through ionization) or are simply blocked by the mechanical baffle structure. Such a divertor condition we call an “opaque divertor’’ [7]. In the “detached regime,’’27 the plasma temperature (density) in the divertor is significantly lower (higher) than in the SOL outside of the divertor, for example, at the outside midplane. Plasma detached regimes are present at moderate or high levels of collisionality in the SOL, which is sufficiently high in which parallel field temperature gradients can be sustained between the main SOL, where power enters by cross-field transport from the confined plasma and the sink of power, in or near the divertor. Detached conditions are usually found at moderate to high levels of the ratio of plasma density/input power in present machines [2, 7]. With typical plasma temperatures of ≈100 eV at the midplane and of ≈3 eV in a detached divertor, the corresponding ion energies will be about 500 and 15 eV, respectively. The difference in erosion rate between these two energies can be enormous since 15 eV is below the physical sputtering threshold for most materials [2]. We see from Figure 1 that the SOL width varies around the poloidal circumference, being the smallest at the outside midplane and expanding near the divertor, particularly near the X point. Typical expansion factors are ≈4 between the divertor plate and the midplane, reflecting the decrease in the poloidal field in the divertor compared with the midplane location [2]. This expansion has a very beneficial effect in that it effectively reduces the power flux on the divertor plates while increasing the plasma volume available for volumetric loss mechanisms such as impurity/hydrogenic radiation and volume recombination [7, 44]. Radiative processes (e.g., due to intrinsic or intentionally added impurities in the divertor region) can produce further reduction in power flux. Most present divertor machines have reduced the power flowing to the divertor plate by factors of >5, using radiative processes [2, 7]. Thus, through flux expansion and radiative processes, the divertor can deliver a factor >20 reduction in incident power density over what would be experienced by a (hypothetical) toroidal28 belt limiter5 normal to the SOL in the main chamber [45].The PFC materials are optimized for the plasma conditions in the specific regions. For example, near the baffle region of the divertor (Figure10), there lies the interface between the hot SOL plasma outside the divertor (>100 eV) and the high neutral density in the divertor. This region is expected to generate significant fluxes of energetic charge exchange5 atoms, which bombard the neighboring structures. These particles have energies >100 eV, and up to perhaps several keV, and cause physical sputtering of PFCs in the baffle region. For this reason, the present ITER design calls for a high Z material, tungsten, in this region, which has a low yield for physical sputtering [46]. In the tokamaks (the ones “ITER like”) the energy leaves the plasma in form of electromagnetic radiation and kinetic energy of particles. Plasma surrounding wall is irradiated by ions, charge-exchange neutrals, electrons, photons (
γ
,
X
,
UV
), and neutrons. All of them modify material properties, from the very surface to the bulk. Therefore, blanket materials must be compatible with ultrahigh vacuum, cryogenics (cry pumps), magnetohydro dynamics, neutron irradiation, and handling of high heat loads [47]. As a consequence, there are stringent requirements regarding the properties of plasma facing components (PFCs)as follows:(i)
high thermal conductivity;
(ii)
good thermomechanical properties and resilience to thermal shocks;
(iii)
nonmagnetic;
(iv)
low activation by neutrons and resistance to radiation damage;
(v)
low accumulation of hydrogen isotopes accompanied by low chemical affinity to hydrogen in order to avoid chemical erosion leading to the formation of volatile compounds;
(vi)
high affinity to oxygen towards formation of stable and nonvolatile oxides is also important for gettering oxygen impurity species in a reactor [47].Properties of no single element, compound, or alloy can satisfy all points of that list. Only few candidate materials for the plasma-facing wall are seriously considered, as(i)
stainless steel;
(ii)
beryllium;
(iii)
tungsten.Behavior of these elements under plasma conditions, that is, particle bombardment and high heat flux deposition, is very different [48]. Therefore, their distribution on the reactor wall is not accidental. Beryllium (safety problems are described in [49]) will be used for the main chamber wall, port limiters, and baffle, whereas the divertor dome and upper vertical target will be covered with tungsten tiles. A detailed distribution of tungsten and carbon in the divertor is shown in [47]. Tungsten is the preferred material for the divertor, except for the area near the strike points. Each of these three candidate materials has some inherent advantages and disadvantages, and their application depends on the specific operational requirements [2]. Because of its low physical sputtering yield and high threshold energy for physical sputtering,W is foreseen in regions where the significant flux of hydrogen neutrals, predominantly of low energy, can lead to large sputtering and unacceptably short lifetimes for materials such as carbon and beryllium. Choosing the proper material is important also to guarantee the integrity of fuel. The fuel failure modes determining fuel lifetime were reviewed and fuel integrity was analyzed and compared with the failure criteria. The correlations of the cladding characteristics were modeled starting by the well-known characteristics of austenitic modified 316 SS (PNC316), ferritic-martensitic steel (PNC-FMS) and oxide dispersion strengthened steel (PNC-ODS) that have been analyzed in [50]. Significant progress has been achieved in the international research effort on reduced activation ferritic/martensitic steels for fusion structural applications. Because this class of steels is the leading structural material for test blankets in ITER and future fusion power systems, the range of ongoing research activities is extremely broad [51–54].One of the main challenges for fusion reactors is the compatibility between a reactor-grade plasma and the materials facing the plasma (the first wall). Until 2011, JET used carbon composite (CFC) tiles for the first wall. However, from JET’s D-T experiments, it is obvious that carbon composites are not suitable for the tritium operation due to high carbon migration, leading to tritium deposition in walls. Therefore, the ITER design comprises a beryllium-clad first wall in the main chamber, while use of carbon tiles is limited to the region where the edge plasma is deflected onto the wall (divertor strike points) and tungsten tiles are to be used elsewhere on the divertor. Tungsten is very resistant to high temperatures (melting only at 3695 degrees Celsius) but it is a heavy element (proton number 74) that can pollute plasmas considerably; it is highly ionized in extreme plasma temperatures, which causes immense energy losses due to plasma radiation and dilutes the D-T fuel. Beryllium is a light element with a proton number just 4. However, it melts at just 1284 degree Celsius. The combination of beryllium and tungsten has never been tested in a tokamak, let alone in one with ITER-relevant geometry and plasma parameters like JET [55]. According to these considerations, the authors, in the new experiments of dust resuspension with STARDUST facility, will focus the attention mainly on tungsten and stainless steel.Dust is produced inside the vacuum vessel of a tokamak by interaction of the plasma with the components of the first wall and the divertor, as we have seen before. In the past, this has received little attention within the fusion community, mainly, because dust is both a safety and not an operational problem in existing tokamaks [2]. However, the ITER design has highlighted the fundamental need to deepen the understanding of the production of dust in tokamaks since this may directly affect the safe operation of a next step device [2]. The safety issues associated with tokamak dust [56–59] include radiological hazard (tritium and activation products), toxicity, and chemical reactivity with steam and air.A variety of “undesired’’ products result from PMIs (e.g., films, flakes, debris, and dust) (Figure5), and it is not clear how many of these should be considered a “dust’’ hazard. Small amounts of dust have been collected in various tokamaks during scheduled vacuum vessel events [60, 61].Figure 5
PMI undesired products.Flaking and breakup of films, resulting from redeposition of eroded and vaporized materials, is expected to be the primary source of dust in a device such as ITER. Arcing, which tends to release micron size particulates, may also play an important role, but investigations are needed to better quantify the effects [2]. If redeposited carbon layers are susceptible to arcing, released particulates will be saturated with deuterium-tritium, presenting safety and operational issues [2]. The radiological hazard and/or toxicity of the dust depend on the specific dust material, how much it is produced, and how well the dust is confined. Confinement depends on the dust particle size (e.g., particles larger than ≈100µm will not be transported to the environment) [2]. In vacuum vessel breach events, dust (beryllium, tungsten, or stainless steel) can also react chemically with steam and air oxidizing and producing hydrogen. In particular, accident scenarios that result in water or steam exposure of hot plasma facing materials (e.g., during loss of coolant accidents) are of concern for fusion reactor designs such as ITER which include beryllium as the first wall material.As we have seen before, a recognized safety issue for future fusion reactors fueled with deuterium and tritium is the generation of sizeable quantities of dust [62–64]. Several mechanisms resulting from material response to plasma bombardment in normal and off-normal conditions are responsible for generating dust of micron and submicron length scales [64]. Neutron activation and tritium uptake in the progenitor material and the particulate deposits themselves make the dust radioactive. In the past ITER the administrative limits for dust inside the VV have been established as 200 kg of carbon, 100 kg of beryllium, and 100 kg of tungsten, on the base of the estimated dust radioactivity. Now the safety limit for dust inside the ITER vacuum vessel (VV) has been proposed to the safety authorities to be 1000 kg, without any precision on the composition. However, there is also a “hot dust” limit (defined as dust on surfaces withT > 400°C) of 6 kg each of C, Be, and W [65]. These limits have been developed to avoid, also in the case of severe accident, the evacuation of the population from the area surrounding the plant [66]. The loss of coolant accidents (LOCA), loss of coolant flow accidents (LOFA), and loss of vacuum accidents (LOVA) [67] are types of accidents that may jeopardize components and plasma vessel integrity and cause dust mobilization risky for workers and public [68]. In the LOCA event the cooling tubes installed into plasma-facing components are broken and the cooling water enters into the vacuum vessel. Then the cooling water boils and evaporates because of the high temperature of the in-vessel components and the low pressure in the vacuum vessel. Consequently, the pressure in the vacuum vessel increases rapidly [69]. Then some safety devices such as VVPSS29 are supposed to operate. However, we could consider some malfunction of the VVPSS [70]. Dust does not strongly adhere to surfaces and is capable, in case of LOVA, of being mobilized, thereby potentially challenging containment barriers in postulated accident scenarios [71]. Additionally the particulate may be chemically reactive and/or toxic. Substantial uncertainties presently exist in estimates of dust production rate, amount, distribution, and mixed materials effects within fusion reactors, potentially imposing burdensome regulatory restrictions on reactor operation. Dust’s impact on fusion reactor safety is partly determined by its spatial distribution and how the particles move around during an accident. There are two primary safety concerns centered on dust mobilization. Significant dust quantities located at positions giving rise to efficient mobilization may contribute greatly to the calculated site source term for a given accident scenario. Secondly, mobilized dust greatly increases the total surface area available for chemical reactions that produce hydrogen during steam ingress accidents. Accidental introduction of air into the vacuum vessel may also lead to airborne dust concentrations capable of rapid oxidation or dust explosion. Accumulation of dust in certain locations that have potentially limited mobilization, such as grooves between tiles of divertor target plates, could display enhanced rates of chemical reactivity compared to bulk planar material. A map of the main risks due to dust mobilization is reported in Figure 6. Possible safety consequences from accidental dust mobilization prompted greater attention to dust in the safety analyses of high energy density machines. In this field by the priority has been given to the issue concerns the simulation of dust transport caused by a continuum phase like water or air ingress into the vacuum vessel (VV). In particular the specific approach includes(1)
scaled experiments to investigate mobilization of dust behavior during accidents; experiments should appropriately simulate geometry, flow conditions, temperature distributions, and structural components that effect mobilization of dust with different characteristics;
(2)
investigating chemical reactivity of dust with various characteristics and exposure configurations;
(3)
providing numerical model for analysis of dust mobilization; validation and benchmarking of the analysis codes are provided by experiment systems.Since 2006 we have made experimental campaigns to reproduce LOVA event using the facility STARDUST [72–81] (Small Tank for Aerosol Removal and Dust) developed by Quantum Electronics and Plasma Physics Group at University of Rome “Tor Vergata” in collaboration with ENEA Laboratories of Frascati. By means of this experimental facility the authors are able to analyse the influence on dust resuspension, of different types of LOVAs during operative conditions (wall temperature of 120–130°C). By the mean of STARDUST facility, the LOVAs have been reproduced with different pressurization rates of(1)
300 Pa/s which is the estimated consequence of a 0,02 m2 wide breach during the first seconds of a LOVA, as defined by General Safety and Security Report (GSSR) [82];
(2)
100 Pa/s and 500 Pa/s (lower and higher than the previous one).The flow rates that have been set to obtain these flow rates are listed in Table 1.Table 1
Flow rate setting.
Temperature
Pressurization rate
Flow rate
120–130°C
100 Pa/s
10 lt/min
300 Pa/s
26.5 lt/min
500 Pa/s
52 lt/minFigure 6
Map of risk due to dust mobilization.For the design and licensing of the nuclear fusion power plant computational fluid dynamics (CFD), codes are considered an essential tool [83]. Accuracy is also required for the CFD analysis of a LOVA in order to provide essential data to the engineering and designers to design the safety features. The present work deals with development and experimental validation of CFD model which simulates the complex thermofluid-dynamic field and gives some indication about internal hazardous dust mobilization during a vessel filling at near vacuum conditions for supporting first instant of LOVA event safety analysis. The development of the CFD model and experimental activities are carried out in strong correlation in order both to understand the capabilities of computational code and to predict correctly the characteristics of the flows during a LOVA event. STARDUST facility is selected both for the model development phase and for the validation case. The extraction of numerical data has been done in the same locations of experimental ones in order to make comparisons and validate the numerical 3D model that will be in able to reproduce different thermofluid-dynamic scenarios in case of LOVAs. A summary of the principal numerical and experimental results will be presented, analyzed, and discussed in the paper.


## 2. Materials

The flow diagram for LOVAS reproduction inside STARDUST facility is showed in Figure7.Figure 7
STARDUST functional scheme for LOVAs.


### 2.1. Hardware Components Used

The hardware has been mounted and set in order to achieve both boundary conditions comparable with those expected in ITER and to reproduce different types of LOVA. The hardware used in the experiments are as follows.(i) External Wall Heaters. The walls can be heated until 300 degrees Celsius by the electric heaters wrapped around the tank, according to the position showed in the Figure 8.Figure 8
Position and powers of electrical resistances.(ii) Vacuum Pump. A vacuum pump (Alcatel Pascal 2021 SD) that allows obtaining values of vacuum until 35 Pa.(iii) Line for Compressed Air. The compressed air is fed inside the tank by an external compressor connected to STARDUST through an external cylindrical tank (stainless steel “lung”). The pressure output is almost constant and fixed at 2, 8 atm by a Fiac air compressor valve; model F0008, a valve to capture the water in the air is also present on the line.(iv) Pressure Gauges. There are three pressure gauges placed on the higher part of STARDUST (as showed in Figure 9): (a) Leybold Heraeus barometer; (b) BOC EDWARDS ASG-2000-NW16 pressure gauge; (c) Alcatel AP 1004 Pirani. The barometer is useful to visualize immediately the internal pressure values, but it is not electrically connected to PC and monitored via software. The pressure gauge BOC EDWARDS ASG-2000-NW16 is an instrument that allows measuring pressure in a range from 1 to 2000 mbar, with a precision of ±0.2%, and it is used to measure the pressure inside the chamber from 1 atm to 1000 Pa. The Alcatel AP 1004 Pirani is a gauge with a range from 1000 to 5*10−4 mbar (100000 to 5*10−2 Pa) and it is used to measure the pressure inside the chamber for values lower that 1000 Pa.Figure 9
From the left: barometer; Pirani; BOC Edwards ASG.Figure 10
Valve A and Valve B on STARDUST facility.(v) Electropneumatic Valve. The electropneumatic valve is a valve placed on the higher part of STARDUST and it allows for, when it is open, the air extraction by a vacuum pump.(vi) Flow Meter. The instrument used is a MKS mass-flow meter and controller 1559A type. This instrument is able to regulate a flow gas inlet and give in output the correct value of the flow gas through itself. This is a massive flow, so corrections due to temperature and pressure variation are not important. It allows fluxing compressed air in the range of 0–100 slm30. The flow meter allows for the pressurized air inlet through 2 valves. Valve A simulates a LOVA at equatorial port level of ITER; Valve B simulates a LOVA at divertor port level of ITER (Figure 10).(vii) Feed Valve. The feed valve is an external valve that allows compressed air to flow through the flow meter.(viii) Thermocouples J Type. In STARDUST J type thermocouples (iron/constantan) are used. They work from −210°C to +750°C. Inside STARDUST 4 thermocouples are placed in three different points of the wall, inside steel pockets (thermocouples 1, 2, and 4), and one (thermocouple 3) is free inside the tank to monitor the environmental temperature during the LOVA experiment (Figure 11).Figure 11
Thermocouples inside STARDUST.(ix) Pressure Transducers. To map the velocity flow values in the tank two XCE-093-2D [84] Kulite pressure transducers have been placed inside a ring support system. Some pressure detectors can work at temperature ranging from −55 to 273°C (the thermal operational range of these pressure transducers is higher than the one used in the previous experimental campaigns [75–79]). The array allows for the measurements simulating operative condition of the tokamak and then in hot walls conditions. We can evaluate velocity using an empirical equation (1) [85] shown below:
(1)
v
=
2
γ
R
T
¯
M
(
γ
-
1
)
[
(
P
Δ
+
P
s
P
s
)
(
γ
-
1
)
/
γ
-
1
]
,
where(i)
γ
:
(
c
p
/
c
v
) (~1.4 for air);
(ii)
R
: universal gas constant [8,314 J/(mol*K)];
(iii)
T
¯: mean temperature of thermocouples (K);
(iv)
M: air molecular mass (28,968 g/mol);
(v)
P
s
: static pressure measured by the pressure gauges (Pa);
(vi)
P
Δ
: differential pressure (
P
Δ
=
P
T
-
P
s
) measured with pressure transducers (Pa);
(vii)
P
T
: total pressure (Pa).(x) Acquisition and Control System Hardware. Currently, the acquisition and control of STARDUST are controlled by a homemade Labview software routine. The main components are COMPAQ DAQ system (National Instruments) that allows for the acquisition of 8 input outputs cards . In our case there are two input and two output cards. The input cards are NI 9219 card constituted by 4 universal analogical channels as input, at 24 bit; NI 9211 card, 4 thermocouples input channels. The output cards are NI 9485 card, 8 digital output channels, with SSR31; NI 9264 card, 16 voltage analogical output, channels at ±10 V. It allows managing all the hardware described above and of acquire all their data with a maximum frequency of 100 Hz.


### 2.2. Software Components Used

The software programs used in this activity are(i)
Labview, to manage the experiment control and acquire the data from the devices;
(ii)
Matlab, to elaborate the experimental data and compare them with the numerical ones;
(iii)
ANSYS-CFX, to implement 2D numerical simulations.Labview. The acquisition system has been developed in Labview; it is a dedicated software developed to execute experiments with STARDUST and acquire data with a controlled frequency with maximum values of 100 Hz. (Figure 12).Figure 12
Labview software for STARDUST control.Matlab. The homemade codes are used for experimental data elaboration and for the benchmark between numerical and experimental data obtained.ANSYS-CFX. The commercial CFD code ANSYS-CFX [86] has been chosen as the solver because it has the capability to treat the main physical phenomena occurring during a LOVA and it has been developing simulation solutions that are well suited to address engineering and regulatory concerns and standards in the global nuclear industry. The ANSYS CFX solver uses finite elements (cell vertex numeric) to discretize the domain. It focuses on the implicit pressure based coupled algebraic multigrid approach to solve the governing equations; this technique solves the fluid-dynamic equations (foru,v, w, p) as a single system.


## 2.1. Hardware Components Used

The hardware has been mounted and set in order to achieve both boundary conditions comparable with those expected in ITER and to reproduce different types of LOVA. The hardware used in the experiments are as follows.(i) External Wall Heaters. The walls can be heated until 300 degrees Celsius by the electric heaters wrapped around the tank, according to the position showed in the Figure 8.Figure 8
Position and powers of electrical resistances.(ii) Vacuum Pump. A vacuum pump (Alcatel Pascal 2021 SD) that allows obtaining values of vacuum until 35 Pa.(iii) Line for Compressed Air. The compressed air is fed inside the tank by an external compressor connected to STARDUST through an external cylindrical tank (stainless steel “lung”). The pressure output is almost constant and fixed at 2, 8 atm by a Fiac air compressor valve; model F0008, a valve to capture the water in the air is also present on the line.(iv) Pressure Gauges. There are three pressure gauges placed on the higher part of STARDUST (as showed in Figure 9): (a) Leybold Heraeus barometer; (b) BOC EDWARDS ASG-2000-NW16 pressure gauge; (c) Alcatel AP 1004 Pirani. The barometer is useful to visualize immediately the internal pressure values, but it is not electrically connected to PC and monitored via software. The pressure gauge BOC EDWARDS ASG-2000-NW16 is an instrument that allows measuring pressure in a range from 1 to 2000 mbar, with a precision of ±0.2%, and it is used to measure the pressure inside the chamber from 1 atm to 1000 Pa. The Alcatel AP 1004 Pirani is a gauge with a range from 1000 to 5*10−4 mbar (100000 to 5*10−2 Pa) and it is used to measure the pressure inside the chamber for values lower that 1000 Pa.Figure 9
From the left: barometer; Pirani; BOC Edwards ASG.Figure 10
Valve A and Valve B on STARDUST facility.(v) Electropneumatic Valve. The electropneumatic valve is a valve placed on the higher part of STARDUST and it allows for, when it is open, the air extraction by a vacuum pump.(vi) Flow Meter. The instrument used is a MKS mass-flow meter and controller 1559A type. This instrument is able to regulate a flow gas inlet and give in output the correct value of the flow gas through itself. This is a massive flow, so corrections due to temperature and pressure variation are not important. It allows fluxing compressed air in the range of 0–100 slm30. The flow meter allows for the pressurized air inlet through 2 valves. Valve A simulates a LOVA at equatorial port level of ITER; Valve B simulates a LOVA at divertor port level of ITER (Figure 10).(vii) Feed Valve. The feed valve is an external valve that allows compressed air to flow through the flow meter.(viii) Thermocouples J Type. In STARDUST J type thermocouples (iron/constantan) are used. They work from −210°C to +750°C. Inside STARDUST 4 thermocouples are placed in three different points of the wall, inside steel pockets (thermocouples 1, 2, and 4), and one (thermocouple 3) is free inside the tank to monitor the environmental temperature during the LOVA experiment (Figure 11).Figure 11
Thermocouples inside STARDUST.(ix) Pressure Transducers. To map the velocity flow values in the tank two XCE-093-2D [84] Kulite pressure transducers have been placed inside a ring support system. Some pressure detectors can work at temperature ranging from −55 to 273°C (the thermal operational range of these pressure transducers is higher than the one used in the previous experimental campaigns [75–79]). The array allows for the measurements simulating operative condition of the tokamak and then in hot walls conditions. We can evaluate velocity using an empirical equation (1) [85] shown below:
(1)
v
=
2
γ
R
T
¯
M
(
γ
-
1
)
[
(
P
Δ
+
P
s
P
s
)
(
γ
-
1
)
/
γ
-
1
]
,
where(i)
γ
:
(
c
p
/
c
v
) (~1.4 for air);
(ii)
R
: universal gas constant [8,314 J/(mol*K)];
(iii)
T
¯: mean temperature of thermocouples (K);
(iv)
M: air molecular mass (28,968 g/mol);
(v)
P
s
: static pressure measured by the pressure gauges (Pa);
(vi)
P
Δ
: differential pressure (
P
Δ
=
P
T
-
P
s
) measured with pressure transducers (Pa);
(vii)
P
T
: total pressure (Pa).(x) Acquisition and Control System Hardware. Currently, the acquisition and control of STARDUST are controlled by a homemade Labview software routine. The main components are COMPAQ DAQ system (National Instruments) that allows for the acquisition of 8 input outputs cards . In our case there are two input and two output cards. The input cards are NI 9219 card constituted by 4 universal analogical channels as input, at 24 bit; NI 9211 card, 4 thermocouples input channels. The output cards are NI 9485 card, 8 digital output channels, with SSR31; NI 9264 card, 16 voltage analogical output, channels at ±10 V. It allows managing all the hardware described above and of acquire all their data with a maximum frequency of 100 Hz.


## 2.2. Software Components Used

The software programs used in this activity are(i)
Labview, to manage the experiment control and acquire the data from the devices;
(ii)
Matlab, to elaborate the experimental data and compare them with the numerical ones;
(iii)
ANSYS-CFX, to implement 2D numerical simulations.Labview. The acquisition system has been developed in Labview; it is a dedicated software developed to execute experiments with STARDUST and acquire data with a controlled frequency with maximum values of 100 Hz. (Figure 12).Figure 12
Labview software for STARDUST control.Matlab. The homemade codes are used for experimental data elaboration and for the benchmark between numerical and experimental data obtained.ANSYS-CFX. The commercial CFD code ANSYS-CFX [86] has been chosen as the solver because it has the capability to treat the main physical phenomena occurring during a LOVA and it has been developing simulation solutions that are well suited to address engineering and regulatory concerns and standards in the global nuclear industry. The ANSYS CFX solver uses finite elements (cell vertex numeric) to discretize the domain. It focuses on the implicit pressure based coupled algebraic multigrid approach to solve the governing equations; this technique solves the fluid-dynamic equations (foru,v, w, p) as a single system.


## 3. Methods

### 3.1. Experimental Protocol Implemented

The LOVA experimental protocol developed by the mean of STARDUST facility is showed in Figure13.Figure 13
LOVA flowchart experimental protocol.These experiments are useful for the evaluation of the capabilities of the code for modeling local gas velocity field at low pressure condition. During the present experimental campaign, focused on both the development and validation of the CFD model of LOVA, valve A is exclusively used in order to provide preliminary results of velocity fields for symmetric configuration. During the experiments the thermofluid-dynamic condition of a LOVA inside STARDUST is monitored by the use of the hardware and software described above in order to have experimental data of the internal pressure, temperature, flow rate, and velocity behavior with a frequency of 100 Hz. These data are important to give the proper boundary conditions for the numerical simulations and to implement benchmark with the numerical data in order to validate a numerical model from a thermofluid-dynamic point of view. It is important to underline that the LOVAs reproduced are those described in Table1. The pressure transducers placed inside STARDUST to measure the velocity flow fields are called PT 455 and PT 461 (Figure 14), and they have been placed inside a ring support system (Figure 15).Figure 14
Pressure transducer.Figure 15
Ring pressure transducers support system.The support has been positioned at several distances from the valves air inlet positions:(i)
4 cm;
(ii)
23 cm;
(iii)
45,5 cm;
(iv)
68, 5 cm.The pressure transducers are placed inside the tank, on this support, at the following quotes considering (Table 2, Figure 16).Table 2
Pressure transducers positions.
Probes coordinates
Name
Plane
X [m]
Y [m]
Z [m]
A
Equatorial
0.01
0.25
0
B
0.23
0.25
0
C
0.45
0.25
0
D
0.69
0.25
0
E
0.69
0.25
−0.25
F
0.455
0.25
−0.25
G
0.23
0.25
−0.25
H
Symmetry
0.01
0.49
0
I
0.23
0.49
0
L
0.455
0.49
0
M
0.69
0.49
0
N
0.69
0.01
0
O
0.45
0.01
0
P
0.23
0.01
0
Q
0.01
0.01
0
Directions
Name
Plane
Point
Line 1
Equatorial
B–G
Line 2
C–F
Line 3
D–E
Line 4
Symmetry
H–Q
Line 5
I–P
Line 6
L–O
Line 7
M-NFigure 16
Pressure transducers position inside STARDUST.The pressure transducer 455 has been placed on points A, B, C, and D with the sensitive element on Valve A (facing the flow field air flux) as in Figure17 that represents the direction of transducers on line I-P positions B (PT 455) and P (PT 461).Figure 17
Pressure transducers orientations.The pressure transducer 461 has been placed on the other points (by rotating the support) with the sensitive element on the lid (in front of the flow field air reflux). For each position, the experiments have been repeated for each configuration reported in Table2 at operative conditions (T
walls
: 120–130°C) and with both the initial pressure conditions considered (100 Pa and 1000 Pa). All the parameters are acquired and they stored when the internal pressure equals the external ones. At these time instants the acquisition is stopped and the results are automatically stored.


### 3.2. Numerical Simulation

The numerical simulations have been carried to reproduce the flow field of the same experimental campaign scenarios used during the laboratory tests; the comparison between numerical data and experimental one is used for CFD code validation. The aim of these simulations was to guarantee the accuracy required for the analysis of a LOVA in order to provide a sufficient margin in the design of safety features. They have used those features to add useful information undetectable with the experiments (like the radial velocity into the tank or the streamlines of the flow). In the nearly future they might be used to address the phenomenological questions of the LOVA that will result as part of the licensing of future nuclear fusion reactor without the need to make conservative assumptions in analysis or to apply potentially nonconservative assumptions that lead to accident termination.


#### 3.2.1. Analytic Unsteady Solution of Pressurization Process

The first order of an analytic study is carried out. It is aimed at understanding the fundamental phenomena occurring during a LOVA event, in particular the pressurization of the vessel. The flow field and, consequently, the dust mobilization are considered effects of this phenomenon. An analytic solution of pressure history of adiabatic and isothermal pressurization of a confined volume has been adopted [87] in order to allow for a preliminary comparison between the developed CFD model and a simplified analytical model known in the scientific literature [87].The boundary condition has been identified as STARDUST itself and the following assumptions have been adopted to simplify the system.(i)
During transient pressurization, it is assumed that thermodynamic equilibrium exists inside the vessel and with respect to the surroundings.
(ii)
The gas (air) is ideal with constant specific heat. The equation of state for an ideal gas (2) is used to calculate the pressure in the chamber:
(2)
P
t
+
Δ
t
=
m
t
+
Δ
t
R
T
t
+
Δ
t
V
.
(iii)
The air in the vessel is at rest.
(iv)
Spatially uniform thermodynamic properties.
(v)
Flow through the inlet is 1D and isentropic.
(vi)
The cross-sectional diameter of the inlet is constant.An implicit successive substitution method is used to solve the governing equations.(1)
Initially(
t
=
0
) internal P
,
T, and m (the mass of air into the chamber) are set according to the experimental data provided by STARDUST experiment.
(2)
Begin a new time step(
t
=
t
+
Δ
t
) and update variables in the current time step; temperature at the current time step is calculated from the internal energy as follows:
(3)
T
t
+
Δ
t
=
u
t
+
Δ
t
c
v
.
(3)
Begin the iteration loop for the current time step and make a guess of the variable values. Usually, in unsteady computation, the values prevailing at the previous time step are the best guess.
(4)
Calculate resident mass,m
t
+
Δ
t, from the mass conservation equation
(4)
m
t
+
Δ
t
=
m
t
+
m
˙
t
+
Δ
t
Δ
t
,
where the mass flow rate can be expressed by the relationship between the stagnation pressure of the supply line and the static pressure in the vessel:
(5)
m
˙
=
2
γ
(
γ
-
1
)
P
S
ρ
S
(
p
P
S
)
2
/
γ
(
1
-
p
P
S
)
(
γ
-
1
)
/
γ
.
(5)
Calculate internal energy from the energy conservation equation(6)
(
m
·
u
)
t
+
Δ
t
=
(
m
·
u
)
t
+
m
˙
t
+
Δ
t
h
s
Δ
t
.
(6)
Calculate the pressure in the chamber from the equation of state.
(7)
Calculate the fractional difference in pressure between successive iterations.
(8)
If the fractional difference is larger than a prespecified small number known as the convergence criterion (convergence criterion for the present problem is set to 10–6), the iteration loop is repeated (steps 3 to 7) until the solution is converged.The adiabatic case, as expected, presents the highest slope, while the slope is lower for the isothermal condition. The variation of initial internal pressure does not seem to influence significantly the pressurization process. The internal temperature seems to increase the slope of the pressurization curve.


#### 3.2.2. Governing Equations of Mathematical Model

The governing equations are applied to an infinitesimally small control volume located in a moving fluid [87]. The set of governing equations can be cast as a single equation called the transport equation (7), Table 3 (or, for a Newtonian fluid, Navier-Stokes equations (NSEs)), for the property (or generic field) Φ
:
(7)
∂
(
ρ
m
Φ
)
∂
t
︸
Rate
of
change
in
the
amount
of
extensive
property
in
the
control
volume
+
∂
(
ρ
m
u
j
Φ
)
∂
x
j
︸
Convection
term
=
∂
∂
x
j
[
Γ
eff
∂
Φ
∂
x
j
]
︸
Diffusion
term
+
S
Φ
︸
Net
source
term
.
In these equations ω
k is the species mass fraction for a given species k, c
p
m is constant-pressure specific heat, h is the enthalpy, p is the pressure, T is the temperature, u
i is the velocity, and Q
′
′
′ is internal heat generation rates. The suffix m refers to the fluid mixture [88]; for a single component fluid, the suffix may be dropped and the equation of mass transfer becomes irrelevant. Similarly, the suffix eff indicates effective values of mass diffusivity D, viscosity μ, and thermal conductivity λ; in turbulent flows, however, the transport properties assume values which are higher compared to the fluid’s ones, as the effective transport properties turn out to be properties of the flow rather than of the fluid [88].Table 3
Generalized representation of transport equations.
Equation
Φ
Γ
eff
S
Φ
Conservation of mass for the mixture
1
0
0
Equation of mass transfer for speciesk
ω
k
ρ
m
D
eff
R
k
Momentum equations
u
i, i = (1, 2, 3)
u
i
μ
eff
-
∂
p
/
∂
x
i
+
ρ
m
B
i
+
S
u
i
Energy equation(enthalpy form)
h
λ
eff
/
c
p
m
Q
′
′
′
Energy equation (temperature form)
T
λ
eff
/
c
p
m
Q
′
′
′
/
c
p
mThe rate of change (or time derivative) term is to be invoked only when a transient phenomenon is under consideration.The termρ
m
Φ denotes the amount of extensive property available in a unit volume.The convection second term takes into account the transport ofΦ due to bulk motion; this first-order derivative term is relatively uncomplicated but assumes considerable significance when stable and convergent numerical solutions are to be economically obtained [88]. The net source term implies an algebraic sum of sources and sinks of Φ; therefore, in a chemically reacting flow (combustion), a given specie k may be generated via some chemical reactions and destroyed (or consumed) via some others and R
k will comprise both positive and negative contributions. In addition, some chemical reactions may be exothermic and endothermic, making positive and negative contributions to Q
′′′ [88]. Similarly, the term B
i in the momentum equations may represent a buoyancy force, a centrifugal and/or Coriolis force, an electromagnetic force, and so forth. Sometimes, B
i may also represent resistance forces. The terms S
u
i represent viscous terms arising from Stokes’s stress laws. The transport equations described above must be extended with constitutive equations of state for density and for enthalpy; in this case, the standard Redlich-Kwong real gas model is used because it is considered as one of the most accurate models. Equation (8) of state is written as
(8)
p
=
R
T
ν
-
(
0.08664
R
T
c
/
p
c
)
-
(
0.42747
R
2
T
c
2
/
p
c
)
(
T
/
T
c
)
-
0.5
ν
(
ν
+
(
0.08664
R
T
c
/
p
c
)
)
,
where ν is the specific volume and p
c and T
c are critical pressure and temperature. Therefore, enthalpy (9), entropy (10), and the specific heat capacities at constant volume (11) and at constant pressure (12) can be calculated as(9)
h
(
T
,
v
)
=
∫
ν
ref
∞
(
T
(
d
p
d
T
)
V
-
p
)
d
ν
T
ref
+
∫
T
ref
T
c
V
0
d
T
-
∫
v
∞
(
T
(
d
p
d
T
)
V
-
p
)
d
ν
T
+
e
(
T
ref
,
v
ref
)
︸
+
p
v
e
(
p
,
T
)
(10)
s
(
T
,
v
)
=
∫
ν
ref
∞
(
d
p
d
T
)
V
d
ν
T
ref
+
∫
T
ref
T
c
p
0
T
d
T
-
R
ln
⁡
(
p
p
ref
)
-
∫
v
∞
(
d
p
d
T
)
V
d
ν
T
+
s
(
T
ref
,
v
ref
)
,
(11)
c
v
=
(
∂
u
∂
T
)
v
=
∂
∫
T
ref
T
(
c
p
0
(
T
)
-
R
)
d
T
∂
T
-
2
(
0.42747
R
2
T
c
2
/
p
c
)
(
T
/
T
c
)
-
0.5
(
0.08664
R
T
c
/
p
c
)
T
×
log
⁡
(
1
+
(
0.08664
R
T
c
/
p
c
)
v
)
,
(12)
c
p
=
c
v
+
v
T
(
(
-
(
∂
p
/
∂
T
)
v
v
(
∂
p
/
∂
v
)
T
)
2
×
(
-
1
v
(
∂
p
/
∂
v
)
T
)
-
1
)
,where c
V
0 and c
p
0 are the zero pressure ideal gas specific heat capacities. The greatest impediment to obtain physically accurate solutions is offered by the diffusion and the net source (S
Φ) terms of the transport equation because both of them require empirical information. In laminar flows, the diffusion term represented by the second-order derivative offers no difficulty because Γ is a fluid property and can be accurately determined (via experiments) isolating the flow under consideration. In turbulent (or transitional) flows, determination of Γ
eff requires considerable empirical support. This is named as turbulence modeling (see next paragraph).


#### 3.2.3. Turbulence Modelling

A main key for the flow fields’ simulations is turbulence; it consists of fluctuations in the flow field in time and space [89]. It is a complex process, mainly because it is three-dimensional, unsteady, and multiscale. It can have a significant effect on the characteristics of the flow, and it occurs when the inertia forces in the fluid become significant compared to viscous forces (high Reynolds number). Two main approaches have been used to deal with this problem.(i)
Time-averaged simulations, where no eddy structures (vortices) are resolved; they employ a time averaged simulation with respect to the turbulence, and only mean statistics are predicted. The most common example is the Reynolds-averaged technique, where all flow variables are individually time averaged and the resulting equations are generally called the Reynolds-averaged Navier-Stokes (RANS) equations.
(ii)
Eddy-resolved simulations, where some or all of the eddy structures are resolved. This category includes large-Eddy simulations (LES) where the turbulence is only resolved up to some threshold wave numbers, beyond which a subgrid scale is employed; they are more accurate than the time-averaged techniques in terms of turbulent diffusion and other statistical flow features [90]. However, eddy-resolved formulations require a harsh price in terms of higher computational resources, both CPU memory and time.Reynolds-Averaged Navier-Stokes (RANS) Equations. Because turbulence consists of random fluctuation of various flow properties we use a statistical approach to this problem. The nonlinearity of Navier-Stokes equations leads to the appearance of momentum fluxes that act as apparent stresses throughout the flow; these momentum fluxes are unknown. We then derive equations for these stresses and the resulting equation includes additional unknown quantities. This problem is named “closure problem,” that is establishing a sufficient number of equations for all of the unknowns variables. Turbulence models try to solve a modified set of transport equations by introducing averaged and fluctuating components of velocity:
(13)
u
i
=
u
¯
i
+
u
i
′
=
1
Δ
t
∫
t
t
+
Δ
t
u
i
+
u
i
′
,
where u
,
u
¯, and u
′ are used to represent the instantaneous, average, and fluctuating terms, respectively. Δ
t is a time scale that is relatively large compared to the turbulent fluctuations but relatively small compared to the time scale to which the equations are solved. For simplicity, density fluctuations are considered negligible. Different from the Navier-Stokes equations, the momentum and scalar transport equations contain turbulent flux terms in addition to the molecular diffusive fluxes; these are the Reynolds stresses, ρ
u
ı
′
u
ȷ
′
-, that arise from the nonlinear convective term in the unaveraged equations and highlight the fact that convective transport due to turbulent velocity fluctuations enhances mixing over and above caused by thermal fluctuations at the molecular level. At high Reynolds numbers, turbulent fluxes are much larger than the molecular fluxes, as turbulent velocity fluctuations used to be larger than the mean free path of thermal fluctuations. The Reynolds averaged energy equation is
(14)
∂
ρ
(
h
+
(
1
/
2
)
u
i
u
i
+
(
1
/
2
)
u
ı
′
¯
2
︷
Turbulent
kinetic
energy
)
︷
Total
entalpy
,
h
tot
∂
t
-
∂
p
∂
t
+
∂
∂
x
j
(
ρ
u
j
h
tot
)
=
∂
∂
x
j
(
λ
∂
T
∂
x
j
-
ρ
h
u
ȷ
′
-
)
+
∂
∂
x
j
[
u
i
(
τ
i
j
-
ρ
h
u
ı
′
u
ȷ
′
-
)
]
︷
Viscous
work
term
+
Q
′′′
.
This equation contains an additional turbulence flux term ρ
h
u
ȷ
′
- compared with the instantaneous equation. Similarly, the Φ may be divided into an average component, Φ
¯, and a time varying component, φ
′, as
(15)
∂
(
ρ
Φ
)
∂
t
+
∂
(
ρ
u
j
Φ
)
∂
x
j
=
∂
∂
x
j
[
Γ
∂
Φ
∂
x
j
-
ρ
h
φ
′
-
]
+
S
Φ
,
where ρ
h
φ
′
- is the Reynolds flux. Turbulence models close the Reynolds averaged equations by providing models for the computation of the Reynolds stresses and Reynolds fluxes, based on the eddy viscosity hypothesis or the eddy diffusivity hypothesis, depending on the factor in which the stresses and fluxes are considered linearly dependent of (see [88]). Subject to these hypotheses, the Reynolds averaged energy equation becomes
(16)
∂
ρ
(
h
+
(
1
/
2
)
u
i
u
i
+
(
1
/
2
)
u
ı
′
¯
2
︷
Turbulent
kinetic
energy
)
︷
Total
entalpy
,
h
tot
∂
t
-
∂
p
∂
t
+
∂
∂
x
j
(
ρ
u
j
h
tot
)
=
∂
∂
x
j
(
λ
∂
T
∂
x
j
-
μ
t
P
r
t
∂
h
∂
x
j
)
+
∂
∂
x
j
[
u
i
(
τ
i
j
-
ρ
u
ı
′
u
ȷ
′
-
)
]
︷
Viscous
work
term
+
Q
′′′
.
Similarly, the Reynolds averaged transport equation for additional variables (nonreacting scalars) becomes
(17)
∂
(
ρ
Φ
)
∂
t
+
∂
(
ρ
u
j
Φ
)
∂
x
j
=
∂
∂
x
j
[
(
Γ
Φ
+
μ
t
σ
Φ
)
∂
Φ
∂
x
j
]
+
S
Φ
.Two-Equation Turbulence Models. As a good compromise between numerical effort and computational accuracy, two-equation turbulence models, where both the velocity and length scale are solved using separate transport equations, are widely used, like the RNG k
-
ε model, a derivation of the k
-
ε model; in the k
-
ε turbulence model, in fact, k
(
L
2
,
T
-
2
) is the turbulence kinetic energy and is defined as the variance of the velocity’s fluctuations and ε
(
L
2
,
T
-
2
) is the turbulence eddy dissipation (the rate at which the velocity’s fluctuations dissipate) [71, 91–93]. The k
-
ε model assumes that the turbulence viscosity is linked to the turbulence kinetic energy and dissipation via the relation:
(18)
μ
t
=
C
μ
ρ
k
2
ε
,
where C
μ is an empirical constant. The values of k and ε come directly from the differential transport equations for the turbulence kinetic energy (19) and turbulence dissipation rate (for full directional buoyancy model):
(19)
∂
(
ρ
k
)
∂
t
+
∂
(
ρ
u
j
k
)
∂
x
j
=
∂
∂
x
j
[
(
μ
+
μ
t
σ
k
)
∂
k
∂
x
j
]
+
μ
t
(
∂
u
i
∂
x
j
+
∂
u
j
∂
x
i
)
∂
u
i
∂
x
j
-
2
3
∂
u
k
∂
x
k
(
ρ
k
+
3
μ
t
∂
u
k
∂
x
k
)
︷
Turbulence
production
due
to
viscous
forces
,
P
k
-
μ
t
ρ
σ
ρ
g
i
∂
p
∂
x
i
︷
Buoyancy
producion
term
,
P
k
b
-
ρ
ε
,
(20)
∂
(
ρ
ε
)
∂
t
+
∂
(
ρ
u
j
ε
)
∂
x
j
=
∂
∂
x
j
[
(
μ
+
μ
t
σ
ε
)
∂
ε
∂
x
j
]
+
ε
k
(
C
ε
1
P
k
-
C
ε
2
ρ
ε
+
C
ε
1
·
max
⁡
(
0
;
P
k
b
)
sin
⁡
(
ϕ
)
)
.
The RNG-based k
-
ε turbulence model is derived from the instantaneous Navier-Stokes equations, using a mathematical technique called “renormalization group’’ (RNG) method; the transport equations for turbulence generation and dissipation are the same as those for the standard k
-
ε model, but the model constants differ, and the constant C
ε
1 is replaced by the function C
ε
1
RNG. The transport equation for turbulence dissipation becomes(21)
∂
(
ρ
ε
)
∂
t
+
∂
(
ρ
u
j
ε
)
∂
x
j
=
∂
∂
x
j
[
(
μ
+
μ
t
σ
ε
)
∂
ε
∂
x
j
]
+
ε
k
[
(
1.42
-
P
k
/
ρ
C
μ
RNG
ε
(
1
-
(
1
/
4.38
)
P
k
/
ρ
C
μ
RNG
ε
)
1
+
β
RNG
(
P
k
/
ρ
C
μ
RNG
ε
)
3
/
2
)
︷
C
ε
1
RNG
×
P
k
-
C
ε
2
RNG
ρ
ε
+
C
ε
1
RNG
·
max
⁡
⁡
(
0
,
P
k
b
)
·
sin
⁡
(
ϕ
)
[
(
1.42
-
P
k
/
ρ
C
μ
R
N
G
ε
(
1
-
(
1
/
4.38
)
P
k
/
ρ
C
μ
R
N
G
ε
)
1
+
β
R
N
G
(
P
k
/
ρ
C
μ
R
N
G
ε
)
3
/
2
)
︷
C
ε
1
R
N
G
]
.For all those turbulence models, compressibility corrections must be taken into account.


#### 3.2.4. The LES Model

The LES approach is more general than the RANS approach and avoids the RANS dependence on boundary conditions for the large-scale eddies. The concept behind large-Eddy simulations is that only the large, energy-carrying scales of motion are computed exactly and the effect of the subgrid scales, which tend to be more isotropic and easier to parameterize than the large ones, is modeled. This is possible because energy transfer between motions at different scales occurs through vortex-stretching processes, so that the energy transfer is predominantly from larger to smaller scales since the motions at the largest scales derive their energy from the forces driving the flow field, whereas the relatively strong viscous forces at the smallest scales disperse the small scale motions in the energy of molecular motions [71]. The governing equations for LES are obtained by filtering the time-dependent Navier-Stokes equations in the physical space. The filtering process effectively filters out the eddies whose scales are smaller than the filter width or grid spacing used in the computations, so that the resulting equations govern the dynamics of the large eddies [88]. The filtering process can be resumed in the following equation:
(22)
Φ
-
(
x
)
=
1
Δ
∫
V
G
(
x
-
x
′
Δ
)
Φ
(
x
′
)
d
x
′
,
where the filter function G has to satisfy the following properties:(1)
G
(
-
x
)
=
G
(
x
);
(2)
∫
-
∞
∞
G
(
x
)
d
x
=
1;
(3)
G
(
x
)
→
0 as |
x
|
→
∞ sufficiently fast so that all moments ∫
-
∞
∞
G
(
x
)
x
n
d
x
=
1
(
n
≥
0
) exist;
(4)
G
(
x
) is localized in (
-
1
/
2
,
1
/
2
).Since the flow field of interest is compressible, a density-weighted filtering (Favre filtering) is adopted to avoid subgrid scale terms in conservation of mass.With a little modification to the definition, a Favre filtered variable can defined as(23)
Φ
~
(
x
)
=
1
ρ
¯
Δ
∫
V
G
(
x
-
x
′
Δ
)
ρ
Φ
(
x
′
)
d
x
′
,
(24)
u
k
=
U
k
+
u
k
=
ρ
u
k
-
ρ
¯
+
u
k
(
Velocity
)
,
s
i
j
=
S
i
j
+
u
i
j
=
ρ
s
ı
ȷ
-
ρ
¯
+
u
k
(
Shear
Strain
)
,
σ
i
j
=
Σ
i
j
+
u
i
j
=
ρ
σ
ı
ȷ
-
ρ
¯
+
u
k
(
Shear
Stress
)
,
ρ
=
ρ
+
ϱ
=
ρ
¯
+
ϱ
(
Density
)
,
p
=
P
+
p
=
p
¯
+
p
(
Pressure
)
,
T
=
T
+
τ
=
ρ
T
-
ρ
¯
+
τ
(
Temperature
)
,
e
=
E
+
e
=
ρ
e
-
ρ
¯
+
e
(
Internal
Energy
)
,
h
=
H
+
h
=
ρ
h
-
ρ
¯
+
h
(
Enthalpy
)
,
where the instantaneous values of variables are denoted inbold, the filtered values in capital letters, and the subgrid components in small letters [71].The Favre-filtered NSEs can be written as(25)
∂
(
ρ
)
∂
t
+
∂
(
ρ
U
k
)
∂
x
k
=
0
,
∂
(
ρ
U
i
)
∂
t
+
∂
(
ρ
U
i
U
j
)
∂
x
j
=
-
∂
P
∂
x
i
+
∂
(
μ
[
2
·
(
1
/
2
)
(
(
∂
U
i
/
∂
x
j
)
+
(
∂
U
j
/
∂
x
i
)
)
︷
Favre-filteres
shear
strain
,
S
i
j
-
(
2
/
3
)
δ
i
j
S
k
k
]
︷
Favre-filteres
shear
stress
,
Σ
i
j
)
∂
x
j
+
∂
(
-
ρ
(
(
u
ı
u
ȷ
~
)
-
U
i
U
j
)
)
︷
Subgrid-scale
stress
tensor
τ
i
j
∂
x
j while the Favre-filtered energy equation (26) can be expressed as
(26)
∂
(
ρ
E
)
∂
t
+
∂
(
ρ
U
k
E
)
∂
x
k
=
∂
∂
x
k
[
λ
¯
C
v
∂
E
∂
x
k
]
-
P
∂
U
k
∂
x
k
+
Σ
i
j
S
i
j
-
∂
(
(
ρ
U
k
E
-
)
-
ρ
U
k
E
)
︷
SGS
Heat
Flux
∂
x
k
-
[
(
p
∂
u
k
∂
x
k
-
P
∂
U
k
∂
x
k
)
]
︷
SGS
Pressure
Dilatation
+
(
(
σ
ı
ȷ
s
ı
ȷ
-
)
-
Σ
i
j
S
i
j
)
.
︷
SGS
Viscous
Dissipation
The subgrid-scale stress tensor needs to be modeled. It can be decomposed further, with respect to the filter dimensions, into contributions from various levels of interactions between the scales of motion (see [71, 88]), in the table that follows (Table 4) the main differences between RANS and LES.Table 4
Comparison between RANS and LES.
RANS
LES
The velocity field equations’ quantities are averaged with respect to time.
The equation of velocity field quantities is averaged with respect to space.
The eddy structures are not solved, and therefore the model has to be calibrated with respect to the type of turbulent current under consideration, wherefrom their behaviour depends.
The eddy structures within the inertial range of the turbulent spectrum are solved by using a sufficiently fine grid, while the subgrid scales of motion are parameterized and modeled, as they can be considered isotropic.
The motion can also be solved in two dimensions, if it is averagely static.
The motion can only be solved in three dimensions, because of the constant presence of the time derivatives due to the unstatic behaviour of the solved eddy structures.
The equations’ number is cut within a predetermined order, and the unknown variables belonging to a higher order are modeled with approximate relations.
The equations’ number is cut with respect to the position within the inertial range, so the difference between what is computed and what is modeled is energetic rather than geometric.
The closure problem is introduced by the Reynolds’s stresses’ tensor and is generally solved by using the turbulent cinematic viscosity, which often causes blunders.
The closure problem is introduced by the subgrid stresses’ tensor and is solved by subgrid models, which have no significant influence on the solution’s precision.
High efficiency, especially for the k-e based closure.
Extra efficiency and computational precision.
Low computational effort and CPU memory required.
High computational effort and CPU memory required.


#### 3.2.5. Simulation Main Parameters

Three-Dimensional Domain Discretization. Some flow fields’ simulation issues, often those very remarkable in order to analyze and understand the real flow behavior and therefore to guarantee the pursued safety standards, need a much finer resolution of the simulations tools to be taken into account and identified. This is a consequence of the fact that these are too advanced to be detected at a rough level of the simulations; these issues are often related to situations where the 3D aspects of the flow and the geometrical effects have a significant influence, like the turbulence; that is why, a 3D axial-symmetric domain was considered in performing the CFD simulations. The spatial discretization consists of hexahedral elements for both the vessel domain and the horizontal inlet pipe. The 3D axial-symmetric domain was discretized (Table 5) using three structured grids: coarse (22.644 cells), fine (249.084 cells), and extra fine (2.264.402 cells). Grids are generated with an automatic procedure programmed in Gambit preprocessor. Considering the LES, the computational grid must be chosen such that the separation between the resolved and the subgrid scales occurs in the inertial subrange of the energy spectrum (where the statistics of motion are unaffected by the gross details of the flow field, unlike the large-scale motions and by the coefficient of viscosity, unlike the small-scale motions, but depend only on the average viscous dissipation of energy that determines the energy transfer across these scales). Accordingly, the fine grid size has been chosen to be one order of magnitude larger than that of the smallest scales.Table 5
Grid dimension.
Grid
Number of cells
Coarse
22.644
Fine
249.084
Extra fine
2.264.402Boundary Conditions. The first condition that has been considered is that a symmetry plane into two specular halves, only one of which has been analyzed, while the flow has been mirrored on the other side, as it has been done during the experimental campaign, has divided the tank. The normal velocity component at the symmetry plane boundary is set to zero u
n
=
0 and the scalar variable gradients normal to the boundary are also set to zero d
ϕ
/
d
n
=
0. At the inlet surface pressure and temperature are imposed. The total pressure (
P
=
101326
Pa
) is specified as an inlet boundary condition and the ANSYS-CFX solver computes the static pressure needed to properly close the boundary condition. The temperature is set to the environmental value (
T
=
296
K
). No-slip boundary conditions are assumed at all the walls. The velocity of the fluid at the wall is set to zero for all simulations, so the boundary condition for the velocity becomes u
wall
=
0. All the walls in the model, except for the heated wall in MC and AC experiments where a temperature condition is set, include the effect of heat transfer through the walls with the environment. A convection heat transfer coefficient (27) is considered for all the walls of the vessel [93]:
(27)
Nu
-
=
[
0.825
+
0.387
Ra
1
/
6
[
1
+
(
0.492
/
Pr
)
9
/
16
]
8
/
27
]
.
Heat flux at the wall boundary is calculated using q
w
=
h
c
(
T
b
-
T
n
w
), where h
c is a specified heat transfer coefficient (from Nusselt number), T
b is the specified boundary temperature (i.e, outside the fluid domain), and T
n
w is the temperature at the internal near-wall boundary element center node. The conditions of the environment are equal to T
=
296
K and P
=
101325 Pa (ambient conditions). The heated wall temperatures are maintained for both cases constant at fixed temperature, T
wall = 25°C for MC and T
wall = 110°C for AC, during the whole transient. As internal initial conditions, two different initial values, 100 and 1000 Pa, are considered in order to verify the effect of pressure on the pressurization of the vessel during the transient. The exact value of initial internal condition (due to the experimental difficulty of setting the internal pressure and temperature at the set point value) is provided by STARDUST data. The operating fluid used in the simulations is dry air assumed initially at rest. The density of the fluid is evaluated using the standard Redlich-Kwong real gas and the other thermodynamic properties are considered variable according to internal database.Time stepping. The selection of an appropriate time step size is essential in order to obtain good convergence. A too large time step will lead to poor convergence and a too small time step will lead to larger computational time. Within a given time step the transport equations are solved during an iteratively process until the convergence criteria are met for all equations. In these simulations a constant time step of 1 × 105 s is used. Each time step is considered to be converged when the maximum residual value is no higher than 1 × 10−6.


## 3.1. Experimental Protocol Implemented

The LOVA experimental protocol developed by the mean of STARDUST facility is showed in Figure13.Figure 13
LOVA flowchart experimental protocol.These experiments are useful for the evaluation of the capabilities of the code for modeling local gas velocity field at low pressure condition. During the present experimental campaign, focused on both the development and validation of the CFD model of LOVA, valve A is exclusively used in order to provide preliminary results of velocity fields for symmetric configuration. During the experiments the thermofluid-dynamic condition of a LOVA inside STARDUST is monitored by the use of the hardware and software described above in order to have experimental data of the internal pressure, temperature, flow rate, and velocity behavior with a frequency of 100 Hz. These data are important to give the proper boundary conditions for the numerical simulations and to implement benchmark with the numerical data in order to validate a numerical model from a thermofluid-dynamic point of view. It is important to underline that the LOVAs reproduced are those described in Table1. The pressure transducers placed inside STARDUST to measure the velocity flow fields are called PT 455 and PT 461 (Figure 14), and they have been placed inside a ring support system (Figure 15).Figure 14
Pressure transducer.Figure 15
Ring pressure transducers support system.The support has been positioned at several distances from the valves air inlet positions:(i)
4 cm;
(ii)
23 cm;
(iii)
45,5 cm;
(iv)
68, 5 cm.The pressure transducers are placed inside the tank, on this support, at the following quotes considering (Table 2, Figure 16).Table 2
Pressure transducers positions.
Probes coordinates
Name
Plane
X [m]
Y [m]
Z [m]
A
Equatorial
0.01
0.25
0
B
0.23
0.25
0
C
0.45
0.25
0
D
0.69
0.25
0
E
0.69
0.25
−0.25
F
0.455
0.25
−0.25
G
0.23
0.25
−0.25
H
Symmetry
0.01
0.49
0
I
0.23
0.49
0
L
0.455
0.49
0
M
0.69
0.49
0
N
0.69
0.01
0
O
0.45
0.01
0
P
0.23
0.01
0
Q
0.01
0.01
0
Directions
Name
Plane
Point
Line 1
Equatorial
B–G
Line 2
C–F
Line 3
D–E
Line 4
Symmetry
H–Q
Line 5
I–P
Line 6
L–O
Line 7
M-NFigure 16
Pressure transducers position inside STARDUST.The pressure transducer 455 has been placed on points A, B, C, and D with the sensitive element on Valve A (facing the flow field air flux) as in Figure17 that represents the direction of transducers on line I-P positions B (PT 455) and P (PT 461).Figure 17
Pressure transducers orientations.The pressure transducer 461 has been placed on the other points (by rotating the support) with the sensitive element on the lid (in front of the flow field air reflux). For each position, the experiments have been repeated for each configuration reported in Table2 at operative conditions (T
walls
: 120–130°C) and with both the initial pressure conditions considered (100 Pa and 1000 Pa). All the parameters are acquired and they stored when the internal pressure equals the external ones. At these time instants the acquisition is stopped and the results are automatically stored.


## 3.2. Numerical Simulation

The numerical simulations have been carried to reproduce the flow field of the same experimental campaign scenarios used during the laboratory tests; the comparison between numerical data and experimental one is used for CFD code validation. The aim of these simulations was to guarantee the accuracy required for the analysis of a LOVA in order to provide a sufficient margin in the design of safety features. They have used those features to add useful information undetectable with the experiments (like the radial velocity into the tank or the streamlines of the flow). In the nearly future they might be used to address the phenomenological questions of the LOVA that will result as part of the licensing of future nuclear fusion reactor without the need to make conservative assumptions in analysis or to apply potentially nonconservative assumptions that lead to accident termination.


### 3.2.1. Analytic Unsteady Solution of Pressurization Process

The first order of an analytic study is carried out. It is aimed at understanding the fundamental phenomena occurring during a LOVA event, in particular the pressurization of the vessel. The flow field and, consequently, the dust mobilization are considered effects of this phenomenon. An analytic solution of pressure history of adiabatic and isothermal pressurization of a confined volume has been adopted [87] in order to allow for a preliminary comparison between the developed CFD model and a simplified analytical model known in the scientific literature [87].The boundary condition has been identified as STARDUST itself and the following assumptions have been adopted to simplify the system.(i)
During transient pressurization, it is assumed that thermodynamic equilibrium exists inside the vessel and with respect to the surroundings.
(ii)
The gas (air) is ideal with constant specific heat. The equation of state for an ideal gas (2) is used to calculate the pressure in the chamber:
(2)
P
t
+
Δ
t
=
m
t
+
Δ
t
R
T
t
+
Δ
t
V
.
(iii)
The air in the vessel is at rest.
(iv)
Spatially uniform thermodynamic properties.
(v)
Flow through the inlet is 1D and isentropic.
(vi)
The cross-sectional diameter of the inlet is constant.An implicit successive substitution method is used to solve the governing equations.(1)
Initially(
t
=
0
) internal P
,
T, and m (the mass of air into the chamber) are set according to the experimental data provided by STARDUST experiment.
(2)
Begin a new time step(
t
=
t
+
Δ
t
) and update variables in the current time step; temperature at the current time step is calculated from the internal energy as follows:
(3)
T
t
+
Δ
t
=
u
t
+
Δ
t
c
v
.
(3)
Begin the iteration loop for the current time step and make a guess of the variable values. Usually, in unsteady computation, the values prevailing at the previous time step are the best guess.
(4)
Calculate resident mass,m
t
+
Δ
t, from the mass conservation equation
(4)
m
t
+
Δ
t
=
m
t
+
m
˙
t
+
Δ
t
Δ
t
,
where the mass flow rate can be expressed by the relationship between the stagnation pressure of the supply line and the static pressure in the vessel:
(5)
m
˙
=
2
γ
(
γ
-
1
)
P
S
ρ
S
(
p
P
S
)
2
/
γ
(
1
-
p
P
S
)
(
γ
-
1
)
/
γ
.
(5)
Calculate internal energy from the energy conservation equation(6)
(
m
·
u
)
t
+
Δ
t
=
(
m
·
u
)
t
+
m
˙
t
+
Δ
t
h
s
Δ
t
.
(6)
Calculate the pressure in the chamber from the equation of state.
(7)
Calculate the fractional difference in pressure between successive iterations.
(8)
If the fractional difference is larger than a prespecified small number known as the convergence criterion (convergence criterion for the present problem is set to 10–6), the iteration loop is repeated (steps 3 to 7) until the solution is converged.The adiabatic case, as expected, presents the highest slope, while the slope is lower for the isothermal condition. The variation of initial internal pressure does not seem to influence significantly the pressurization process. The internal temperature seems to increase the slope of the pressurization curve.


### 3.2.2. Governing Equations of Mathematical Model

The governing equations are applied to an infinitesimally small control volume located in a moving fluid [87]. The set of governing equations can be cast as a single equation called the transport equation (7), Table 3 (or, for a Newtonian fluid, Navier-Stokes equations (NSEs)), for the property (or generic field) Φ
:
(7)
∂
(
ρ
m
Φ
)
∂
t
︸
Rate
of
change
in
the
amount
of
extensive
property
in
the
control
volume
+
∂
(
ρ
m
u
j
Φ
)
∂
x
j
︸
Convection
term
=
∂
∂
x
j
[
Γ
eff
∂
Φ
∂
x
j
]
︸
Diffusion
term
+
S
Φ
︸
Net
source
term
.
In these equations ω
k is the species mass fraction for a given species k, c
p
m is constant-pressure specific heat, h is the enthalpy, p is the pressure, T is the temperature, u
i is the velocity, and Q
′
′
′ is internal heat generation rates. The suffix m refers to the fluid mixture [88]; for a single component fluid, the suffix may be dropped and the equation of mass transfer becomes irrelevant. Similarly, the suffix eff indicates effective values of mass diffusivity D, viscosity μ, and thermal conductivity λ; in turbulent flows, however, the transport properties assume values which are higher compared to the fluid’s ones, as the effective transport properties turn out to be properties of the flow rather than of the fluid [88].Table 3
Generalized representation of transport equations.
Equation
Φ
Γ
eff
S
Φ
Conservation of mass for the mixture
1
0
0
Equation of mass transfer for speciesk
ω
k
ρ
m
D
eff
R
k
Momentum equations
u
i, i = (1, 2, 3)
u
i
μ
eff
-
∂
p
/
∂
x
i
+
ρ
m
B
i
+
S
u
i
Energy equation(enthalpy form)
h
λ
eff
/
c
p
m
Q
′
′
′
Energy equation (temperature form)
T
λ
eff
/
c
p
m
Q
′
′
′
/
c
p
mThe rate of change (or time derivative) term is to be invoked only when a transient phenomenon is under consideration.The termρ
m
Φ denotes the amount of extensive property available in a unit volume.The convection second term takes into account the transport ofΦ due to bulk motion; this first-order derivative term is relatively uncomplicated but assumes considerable significance when stable and convergent numerical solutions are to be economically obtained [88]. The net source term implies an algebraic sum of sources and sinks of Φ; therefore, in a chemically reacting flow (combustion), a given specie k may be generated via some chemical reactions and destroyed (or consumed) via some others and R
k will comprise both positive and negative contributions. In addition, some chemical reactions may be exothermic and endothermic, making positive and negative contributions to Q
′′′ [88]. Similarly, the term B
i in the momentum equations may represent a buoyancy force, a centrifugal and/or Coriolis force, an electromagnetic force, and so forth. Sometimes, B
i may also represent resistance forces. The terms S
u
i represent viscous terms arising from Stokes’s stress laws. The transport equations described above must be extended with constitutive equations of state for density and for enthalpy; in this case, the standard Redlich-Kwong real gas model is used because it is considered as one of the most accurate models. Equation (8) of state is written as
(8)
p
=
R
T
ν
-
(
0.08664
R
T
c
/
p
c
)
-
(
0.42747
R
2
T
c
2
/
p
c
)
(
T
/
T
c
)
-
0.5
ν
(
ν
+
(
0.08664
R
T
c
/
p
c
)
)
,
where ν is the specific volume and p
c and T
c are critical pressure and temperature. Therefore, enthalpy (9), entropy (10), and the specific heat capacities at constant volume (11) and at constant pressure (12) can be calculated as(9)
h
(
T
,
v
)
=
∫
ν
ref
∞
(
T
(
d
p
d
T
)
V
-
p
)
d
ν
T
ref
+
∫
T
ref
T
c
V
0
d
T
-
∫
v
∞
(
T
(
d
p
d
T
)
V
-
p
)
d
ν
T
+
e
(
T
ref
,
v
ref
)
︸
+
p
v
e
(
p
,
T
)
(10)
s
(
T
,
v
)
=
∫
ν
ref
∞
(
d
p
d
T
)
V
d
ν
T
ref
+
∫
T
ref
T
c
p
0
T
d
T
-
R
ln
⁡
(
p
p
ref
)
-
∫
v
∞
(
d
p
d
T
)
V
d
ν
T
+
s
(
T
ref
,
v
ref
)
,
(11)
c
v
=
(
∂
u
∂
T
)
v
=
∂
∫
T
ref
T
(
c
p
0
(
T
)
-
R
)
d
T
∂
T
-
2
(
0.42747
R
2
T
c
2
/
p
c
)
(
T
/
T
c
)
-
0.5
(
0.08664
R
T
c
/
p
c
)
T
×
log
⁡
(
1
+
(
0.08664
R
T
c
/
p
c
)
v
)
,
(12)
c
p
=
c
v
+
v
T
(
(
-
(
∂
p
/
∂
T
)
v
v
(
∂
p
/
∂
v
)
T
)
2
×
(
-
1
v
(
∂
p
/
∂
v
)
T
)
-
1
)
,where c
V
0 and c
p
0 are the zero pressure ideal gas specific heat capacities. The greatest impediment to obtain physically accurate solutions is offered by the diffusion and the net source (S
Φ) terms of the transport equation because both of them require empirical information. In laminar flows, the diffusion term represented by the second-order derivative offers no difficulty because Γ is a fluid property and can be accurately determined (via experiments) isolating the flow under consideration. In turbulent (or transitional) flows, determination of Γ
eff requires considerable empirical support. This is named as turbulence modeling (see next paragraph).


### 3.2.3. Turbulence Modelling

A main key for the flow fields’ simulations is turbulence; it consists of fluctuations in the flow field in time and space [89]. It is a complex process, mainly because it is three-dimensional, unsteady, and multiscale. It can have a significant effect on the characteristics of the flow, and it occurs when the inertia forces in the fluid become significant compared to viscous forces (high Reynolds number). Two main approaches have been used to deal with this problem.(i)
Time-averaged simulations, where no eddy structures (vortices) are resolved; they employ a time averaged simulation with respect to the turbulence, and only mean statistics are predicted. The most common example is the Reynolds-averaged technique, where all flow variables are individually time averaged and the resulting equations are generally called the Reynolds-averaged Navier-Stokes (RANS) equations.
(ii)
Eddy-resolved simulations, where some or all of the eddy structures are resolved. This category includes large-Eddy simulations (LES) where the turbulence is only resolved up to some threshold wave numbers, beyond which a subgrid scale is employed; they are more accurate than the time-averaged techniques in terms of turbulent diffusion and other statistical flow features [90]. However, eddy-resolved formulations require a harsh price in terms of higher computational resources, both CPU memory and time.Reynolds-Averaged Navier-Stokes (RANS) Equations. Because turbulence consists of random fluctuation of various flow properties we use a statistical approach to this problem. The nonlinearity of Navier-Stokes equations leads to the appearance of momentum fluxes that act as apparent stresses throughout the flow; these momentum fluxes are unknown. We then derive equations for these stresses and the resulting equation includes additional unknown quantities. This problem is named “closure problem,” that is establishing a sufficient number of equations for all of the unknowns variables. Turbulence models try to solve a modified set of transport equations by introducing averaged and fluctuating components of velocity:
(13)
u
i
=
u
¯
i
+
u
i
′
=
1
Δ
t
∫
t
t
+
Δ
t
u
i
+
u
i
′
,
where u
,
u
¯, and u
′ are used to represent the instantaneous, average, and fluctuating terms, respectively. Δ
t is a time scale that is relatively large compared to the turbulent fluctuations but relatively small compared to the time scale to which the equations are solved. For simplicity, density fluctuations are considered negligible. Different from the Navier-Stokes equations, the momentum and scalar transport equations contain turbulent flux terms in addition to the molecular diffusive fluxes; these are the Reynolds stresses, ρ
u
ı
′
u
ȷ
′
-, that arise from the nonlinear convective term in the unaveraged equations and highlight the fact that convective transport due to turbulent velocity fluctuations enhances mixing over and above caused by thermal fluctuations at the molecular level. At high Reynolds numbers, turbulent fluxes are much larger than the molecular fluxes, as turbulent velocity fluctuations used to be larger than the mean free path of thermal fluctuations. The Reynolds averaged energy equation is
(14)
∂
ρ
(
h
+
(
1
/
2
)
u
i
u
i
+
(
1
/
2
)
u
ı
′
¯
2
︷
Turbulent
kinetic
energy
)
︷
Total
entalpy
,
h
tot
∂
t
-
∂
p
∂
t
+
∂
∂
x
j
(
ρ
u
j
h
tot
)
=
∂
∂
x
j
(
λ
∂
T
∂
x
j
-
ρ
h
u
ȷ
′
-
)
+
∂
∂
x
j
[
u
i
(
τ
i
j
-
ρ
h
u
ı
′
u
ȷ
′
-
)
]
︷
Viscous
work
term
+
Q
′′′
.
This equation contains an additional turbulence flux term ρ
h
u
ȷ
′
- compared with the instantaneous equation. Similarly, the Φ may be divided into an average component, Φ
¯, and a time varying component, φ
′, as
(15)
∂
(
ρ
Φ
)
∂
t
+
∂
(
ρ
u
j
Φ
)
∂
x
j
=
∂
∂
x
j
[
Γ
∂
Φ
∂
x
j
-
ρ
h
φ
′
-
]
+
S
Φ
,
where ρ
h
φ
′
- is the Reynolds flux. Turbulence models close the Reynolds averaged equations by providing models for the computation of the Reynolds stresses and Reynolds fluxes, based on the eddy viscosity hypothesis or the eddy diffusivity hypothesis, depending on the factor in which the stresses and fluxes are considered linearly dependent of (see [88]). Subject to these hypotheses, the Reynolds averaged energy equation becomes
(16)
∂
ρ
(
h
+
(
1
/
2
)
u
i
u
i
+
(
1
/
2
)
u
ı
′
¯
2
︷
Turbulent
kinetic
energy
)
︷
Total
entalpy
,
h
tot
∂
t
-
∂
p
∂
t
+
∂
∂
x
j
(
ρ
u
j
h
tot
)
=
∂
∂
x
j
(
λ
∂
T
∂
x
j
-
μ
t
P
r
t
∂
h
∂
x
j
)
+
∂
∂
x
j
[
u
i
(
τ
i
j
-
ρ
u
ı
′
u
ȷ
′
-
)
]
︷
Viscous
work
term
+
Q
′′′
.
Similarly, the Reynolds averaged transport equation for additional variables (nonreacting scalars) becomes
(17)
∂
(
ρ
Φ
)
∂
t
+
∂
(
ρ
u
j
Φ
)
∂
x
j
=
∂
∂
x
j
[
(
Γ
Φ
+
μ
t
σ
Φ
)
∂
Φ
∂
x
j
]
+
S
Φ
.Two-Equation Turbulence Models. As a good compromise between numerical effort and computational accuracy, two-equation turbulence models, where both the velocity and length scale are solved using separate transport equations, are widely used, like the RNG k
-
ε model, a derivation of the k
-
ε model; in the k
-
ε turbulence model, in fact, k
(
L
2
,
T
-
2
) is the turbulence kinetic energy and is defined as the variance of the velocity’s fluctuations and ε
(
L
2
,
T
-
2
) is the turbulence eddy dissipation (the rate at which the velocity’s fluctuations dissipate) [71, 91–93]. The k
-
ε model assumes that the turbulence viscosity is linked to the turbulence kinetic energy and dissipation via the relation:
(18)
μ
t
=
C
μ
ρ
k
2
ε
,
where C
μ is an empirical constant. The values of k and ε come directly from the differential transport equations for the turbulence kinetic energy (19) and turbulence dissipation rate (for full directional buoyancy model):
(19)
∂
(
ρ
k
)
∂
t
+
∂
(
ρ
u
j
k
)
∂
x
j
=
∂
∂
x
j
[
(
μ
+
μ
t
σ
k
)
∂
k
∂
x
j
]
+
μ
t
(
∂
u
i
∂
x
j
+
∂
u
j
∂
x
i
)
∂
u
i
∂
x
j
-
2
3
∂
u
k
∂
x
k
(
ρ
k
+
3
μ
t
∂
u
k
∂
x
k
)
︷
Turbulence
production
due
to
viscous
forces
,
P
k
-
μ
t
ρ
σ
ρ
g
i
∂
p
∂
x
i
︷
Buoyancy
producion
term
,
P
k
b
-
ρ
ε
,
(20)
∂
(
ρ
ε
)
∂
t
+
∂
(
ρ
u
j
ε
)
∂
x
j
=
∂
∂
x
j
[
(
μ
+
μ
t
σ
ε
)
∂
ε
∂
x
j
]
+
ε
k
(
C
ε
1
P
k
-
C
ε
2
ρ
ε
+
C
ε
1
·
max
⁡
(
0
;
P
k
b
)
sin
⁡
(
ϕ
)
)
.
The RNG-based k
-
ε turbulence model is derived from the instantaneous Navier-Stokes equations, using a mathematical technique called “renormalization group’’ (RNG) method; the transport equations for turbulence generation and dissipation are the same as those for the standard k
-
ε model, but the model constants differ, and the constant C
ε
1 is replaced by the function C
ε
1
RNG. The transport equation for turbulence dissipation becomes(21)
∂
(
ρ
ε
)
∂
t
+
∂
(
ρ
u
j
ε
)
∂
x
j
=
∂
∂
x
j
[
(
μ
+
μ
t
σ
ε
)
∂
ε
∂
x
j
]
+
ε
k
[
(
1.42
-
P
k
/
ρ
C
μ
RNG
ε
(
1
-
(
1
/
4.38
)
P
k
/
ρ
C
μ
RNG
ε
)
1
+
β
RNG
(
P
k
/
ρ
C
μ
RNG
ε
)
3
/
2
)
︷
C
ε
1
RNG
×
P
k
-
C
ε
2
RNG
ρ
ε
+
C
ε
1
RNG
·
max
⁡
⁡
(
0
,
P
k
b
)
·
sin
⁡
(
ϕ
)
[
(
1.42
-
P
k
/
ρ
C
μ
R
N
G
ε
(
1
-
(
1
/
4.38
)
P
k
/
ρ
C
μ
R
N
G
ε
)
1
+
β
R
N
G
(
P
k
/
ρ
C
μ
R
N
G
ε
)
3
/
2
)
︷
C
ε
1
R
N
G
]
.For all those turbulence models, compressibility corrections must be taken into account.


### 3.2.4. The LES Model

The LES approach is more general than the RANS approach and avoids the RANS dependence on boundary conditions for the large-scale eddies. The concept behind large-Eddy simulations is that only the large, energy-carrying scales of motion are computed exactly and the effect of the subgrid scales, which tend to be more isotropic and easier to parameterize than the large ones, is modeled. This is possible because energy transfer between motions at different scales occurs through vortex-stretching processes, so that the energy transfer is predominantly from larger to smaller scales since the motions at the largest scales derive their energy from the forces driving the flow field, whereas the relatively strong viscous forces at the smallest scales disperse the small scale motions in the energy of molecular motions [71]. The governing equations for LES are obtained by filtering the time-dependent Navier-Stokes equations in the physical space. The filtering process effectively filters out the eddies whose scales are smaller than the filter width or grid spacing used in the computations, so that the resulting equations govern the dynamics of the large eddies [88]. The filtering process can be resumed in the following equation:
(22)
Φ
-
(
x
)
=
1
Δ
∫
V
G
(
x
-
x
′
Δ
)
Φ
(
x
′
)
d
x
′
,
where the filter function G has to satisfy the following properties:(1)
G
(
-
x
)
=
G
(
x
);
(2)
∫
-
∞
∞
G
(
x
)
d
x
=
1;
(3)
G
(
x
)
→
0 as |
x
|
→
∞ sufficiently fast so that all moments ∫
-
∞
∞
G
(
x
)
x
n
d
x
=
1
(
n
≥
0
) exist;
(4)
G
(
x
) is localized in (
-
1
/
2
,
1
/
2
).Since the flow field of interest is compressible, a density-weighted filtering (Favre filtering) is adopted to avoid subgrid scale terms in conservation of mass.With a little modification to the definition, a Favre filtered variable can defined as(23)
Φ
~
(
x
)
=
1
ρ
¯
Δ
∫
V
G
(
x
-
x
′
Δ
)
ρ
Φ
(
x
′
)
d
x
′
,
(24)
u
k
=
U
k
+
u
k
=
ρ
u
k
-
ρ
¯
+
u
k
(
Velocity
)
,
s
i
j
=
S
i
j
+
u
i
j
=
ρ
s
ı
ȷ
-
ρ
¯
+
u
k
(
Shear
Strain
)
,
σ
i
j
=
Σ
i
j
+
u
i
j
=
ρ
σ
ı
ȷ
-
ρ
¯
+
u
k
(
Shear
Stress
)
,
ρ
=
ρ
+
ϱ
=
ρ
¯
+
ϱ
(
Density
)
,
p
=
P
+
p
=
p
¯
+
p
(
Pressure
)
,
T
=
T
+
τ
=
ρ
T
-
ρ
¯
+
τ
(
Temperature
)
,
e
=
E
+
e
=
ρ
e
-
ρ
¯
+
e
(
Internal
Energy
)
,
h
=
H
+
h
=
ρ
h
-
ρ
¯
+
h
(
Enthalpy
)
,
where the instantaneous values of variables are denoted inbold, the filtered values in capital letters, and the subgrid components in small letters [71].The Favre-filtered NSEs can be written as(25)
∂
(
ρ
)
∂
t
+
∂
(
ρ
U
k
)
∂
x
k
=
0
,
∂
(
ρ
U
i
)
∂
t
+
∂
(
ρ
U
i
U
j
)
∂
x
j
=
-
∂
P
∂
x
i
+
∂
(
μ
[
2
·
(
1
/
2
)
(
(
∂
U
i
/
∂
x
j
)
+
(
∂
U
j
/
∂
x
i
)
)
︷
Favre-filteres
shear
strain
,
S
i
j
-
(
2
/
3
)
δ
i
j
S
k
k
]
︷
Favre-filteres
shear
stress
,
Σ
i
j
)
∂
x
j
+
∂
(
-
ρ
(
(
u
ı
u
ȷ
~
)
-
U
i
U
j
)
)
︷
Subgrid-scale
stress
tensor
τ
i
j
∂
x
j while the Favre-filtered energy equation (26) can be expressed as
(26)
∂
(
ρ
E
)
∂
t
+
∂
(
ρ
U
k
E
)
∂
x
k
=
∂
∂
x
k
[
λ
¯
C
v
∂
E
∂
x
k
]
-
P
∂
U
k
∂
x
k
+
Σ
i
j
S
i
j
-
∂
(
(
ρ
U
k
E
-
)
-
ρ
U
k
E
)
︷
SGS
Heat
Flux
∂
x
k
-
[
(
p
∂
u
k
∂
x
k
-
P
∂
U
k
∂
x
k
)
]
︷
SGS
Pressure
Dilatation
+
(
(
σ
ı
ȷ
s
ı
ȷ
-
)
-
Σ
i
j
S
i
j
)
.
︷
SGS
Viscous
Dissipation
The subgrid-scale stress tensor needs to be modeled. It can be decomposed further, with respect to the filter dimensions, into contributions from various levels of interactions between the scales of motion (see [71, 88]), in the table that follows (Table 4) the main differences between RANS and LES.Table 4
Comparison between RANS and LES.
RANS
LES
The velocity field equations’ quantities are averaged with respect to time.
The equation of velocity field quantities is averaged with respect to space.
The eddy structures are not solved, and therefore the model has to be calibrated with respect to the type of turbulent current under consideration, wherefrom their behaviour depends.
The eddy structures within the inertial range of the turbulent spectrum are solved by using a sufficiently fine grid, while the subgrid scales of motion are parameterized and modeled, as they can be considered isotropic.
The motion can also be solved in two dimensions, if it is averagely static.
The motion can only be solved in three dimensions, because of the constant presence of the time derivatives due to the unstatic behaviour of the solved eddy structures.
The equations’ number is cut within a predetermined order, and the unknown variables belonging to a higher order are modeled with approximate relations.
The equations’ number is cut with respect to the position within the inertial range, so the difference between what is computed and what is modeled is energetic rather than geometric.
The closure problem is introduced by the Reynolds’s stresses’ tensor and is generally solved by using the turbulent cinematic viscosity, which often causes blunders.
The closure problem is introduced by the subgrid stresses’ tensor and is solved by subgrid models, which have no significant influence on the solution’s precision.
High efficiency, especially for the k-e based closure.
Extra efficiency and computational precision.
Low computational effort and CPU memory required.
High computational effort and CPU memory required.


### 3.2.5. Simulation Main Parameters

Three-Dimensional Domain Discretization. Some flow fields’ simulation issues, often those very remarkable in order to analyze and understand the real flow behavior and therefore to guarantee the pursued safety standards, need a much finer resolution of the simulations tools to be taken into account and identified. This is a consequence of the fact that these are too advanced to be detected at a rough level of the simulations; these issues are often related to situations where the 3D aspects of the flow and the geometrical effects have a significant influence, like the turbulence; that is why, a 3D axial-symmetric domain was considered in performing the CFD simulations. The spatial discretization consists of hexahedral elements for both the vessel domain and the horizontal inlet pipe. The 3D axial-symmetric domain was discretized (Table 5) using three structured grids: coarse (22.644 cells), fine (249.084 cells), and extra fine (2.264.402 cells). Grids are generated with an automatic procedure programmed in Gambit preprocessor. Considering the LES, the computational grid must be chosen such that the separation between the resolved and the subgrid scales occurs in the inertial subrange of the energy spectrum (where the statistics of motion are unaffected by the gross details of the flow field, unlike the large-scale motions and by the coefficient of viscosity, unlike the small-scale motions, but depend only on the average viscous dissipation of energy that determines the energy transfer across these scales). Accordingly, the fine grid size has been chosen to be one order of magnitude larger than that of the smallest scales.Table 5
Grid dimension.
Grid
Number of cells
Coarse
22.644
Fine
249.084
Extra fine
2.264.402Boundary Conditions. The first condition that has been considered is that a symmetry plane into two specular halves, only one of which has been analyzed, while the flow has been mirrored on the other side, as it has been done during the experimental campaign, has divided the tank. The normal velocity component at the symmetry plane boundary is set to zero u
n
=
0 and the scalar variable gradients normal to the boundary are also set to zero d
ϕ
/
d
n
=
0. At the inlet surface pressure and temperature are imposed. The total pressure (
P
=
101326
Pa
) is specified as an inlet boundary condition and the ANSYS-CFX solver computes the static pressure needed to properly close the boundary condition. The temperature is set to the environmental value (
T
=
296
K
). No-slip boundary conditions are assumed at all the walls. The velocity of the fluid at the wall is set to zero for all simulations, so the boundary condition for the velocity becomes u
wall
=
0. All the walls in the model, except for the heated wall in MC and AC experiments where a temperature condition is set, include the effect of heat transfer through the walls with the environment. A convection heat transfer coefficient (27) is considered for all the walls of the vessel [93]:
(27)
Nu
-
=
[
0.825
+
0.387
Ra
1
/
6
[
1
+
(
0.492
/
Pr
)
9
/
16
]
8
/
27
]
.
Heat flux at the wall boundary is calculated using q
w
=
h
c
(
T
b
-
T
n
w
), where h
c is a specified heat transfer coefficient (from Nusselt number), T
b is the specified boundary temperature (i.e, outside the fluid domain), and T
n
w is the temperature at the internal near-wall boundary element center node. The conditions of the environment are equal to T
=
296
K and P
=
101325 Pa (ambient conditions). The heated wall temperatures are maintained for both cases constant at fixed temperature, T
wall = 25°C for MC and T
wall = 110°C for AC, during the whole transient. As internal initial conditions, two different initial values, 100 and 1000 Pa, are considered in order to verify the effect of pressure on the pressurization of the vessel during the transient. The exact value of initial internal condition (due to the experimental difficulty of setting the internal pressure and temperature at the set point value) is provided by STARDUST data. The operating fluid used in the simulations is dry air assumed initially at rest. The density of the fluid is evaluated using the standard Redlich-Kwong real gas and the other thermodynamic properties are considered variable according to internal database.Time stepping. The selection of an appropriate time step size is essential in order to obtain good convergence. A too large time step will lead to poor convergence and a too small time step will lead to larger computational time. Within a given time step the transport equations are solved during an iteratively process until the convergence criteria are met for all equations. In these simulations a constant time step of 1 × 105 s is used. Each time step is considered to be converged when the maximum residual value is no higher than 1 × 10−6.


## 3.2.1. Analytic Unsteady Solution of Pressurization Process

The first order of an analytic study is carried out. It is aimed at understanding the fundamental phenomena occurring during a LOVA event, in particular the pressurization of the vessel. The flow field and, consequently, the dust mobilization are considered effects of this phenomenon. An analytic solution of pressure history of adiabatic and isothermal pressurization of a confined volume has been adopted [87] in order to allow for a preliminary comparison between the developed CFD model and a simplified analytical model known in the scientific literature [87].The boundary condition has been identified as STARDUST itself and the following assumptions have been adopted to simplify the system.(i)
During transient pressurization, it is assumed that thermodynamic equilibrium exists inside the vessel and with respect to the surroundings.
(ii)
The gas (air) is ideal with constant specific heat. The equation of state for an ideal gas (2) is used to calculate the pressure in the chamber:
(2)
P
t
+
Δ
t
=
m
t
+
Δ
t
R
T
t
+
Δ
t
V
.
(iii)
The air in the vessel is at rest.
(iv)
Spatially uniform thermodynamic properties.
(v)
Flow through the inlet is 1D and isentropic.
(vi)
The cross-sectional diameter of the inlet is constant.An implicit successive substitution method is used to solve the governing equations.(1)
Initially(
t
=
0
) internal P
,
T, and m (the mass of air into the chamber) are set according to the experimental data provided by STARDUST experiment.
(2)
Begin a new time step(
t
=
t
+
Δ
t
) and update variables in the current time step; temperature at the current time step is calculated from the internal energy as follows:
(3)
T
t
+
Δ
t
=
u
t
+
Δ
t
c
v
.
(3)
Begin the iteration loop for the current time step and make a guess of the variable values. Usually, in unsteady computation, the values prevailing at the previous time step are the best guess.
(4)
Calculate resident mass,m
t
+
Δ
t, from the mass conservation equation
(4)
m
t
+
Δ
t
=
m
t
+
m
˙
t
+
Δ
t
Δ
t
,
where the mass flow rate can be expressed by the relationship between the stagnation pressure of the supply line and the static pressure in the vessel:
(5)
m
˙
=
2
γ
(
γ
-
1
)
P
S
ρ
S
(
p
P
S
)
2
/
γ
(
1
-
p
P
S
)
(
γ
-
1
)
/
γ
.
(5)
Calculate internal energy from the energy conservation equation(6)
(
m
·
u
)
t
+
Δ
t
=
(
m
·
u
)
t
+
m
˙
t
+
Δ
t
h
s
Δ
t
.
(6)
Calculate the pressure in the chamber from the equation of state.
(7)
Calculate the fractional difference in pressure between successive iterations.
(8)
If the fractional difference is larger than a prespecified small number known as the convergence criterion (convergence criterion for the present problem is set to 10–6), the iteration loop is repeated (steps 3 to 7) until the solution is converged.The adiabatic case, as expected, presents the highest slope, while the slope is lower for the isothermal condition. The variation of initial internal pressure does not seem to influence significantly the pressurization process. The internal temperature seems to increase the slope of the pressurization curve.


## 3.2.2. Governing Equations of Mathematical Model

The governing equations are applied to an infinitesimally small control volume located in a moving fluid [87]. The set of governing equations can be cast as a single equation called the transport equation (7), Table 3 (or, for a Newtonian fluid, Navier-Stokes equations (NSEs)), for the property (or generic field) Φ
:
(7)
∂
(
ρ
m
Φ
)
∂
t
︸
Rate
of
change
in
the
amount
of
extensive
property
in
the
control
volume
+
∂
(
ρ
m
u
j
Φ
)
∂
x
j
︸
Convection
term
=
∂
∂
x
j
[
Γ
eff
∂
Φ
∂
x
j
]
︸
Diffusion
term
+
S
Φ
︸
Net
source
term
.
In these equations ω
k is the species mass fraction for a given species k, c
p
m is constant-pressure specific heat, h is the enthalpy, p is the pressure, T is the temperature, u
i is the velocity, and Q
′
′
′ is internal heat generation rates. The suffix m refers to the fluid mixture [88]; for a single component fluid, the suffix may be dropped and the equation of mass transfer becomes irrelevant. Similarly, the suffix eff indicates effective values of mass diffusivity D, viscosity μ, and thermal conductivity λ; in turbulent flows, however, the transport properties assume values which are higher compared to the fluid’s ones, as the effective transport properties turn out to be properties of the flow rather than of the fluid [88].Table 3
Generalized representation of transport equations.
Equation
Φ
Γ
eff
S
Φ
Conservation of mass for the mixture
1
0
0
Equation of mass transfer for speciesk
ω
k
ρ
m
D
eff
R
k
Momentum equations
u
i, i = (1, 2, 3)
u
i
μ
eff
-
∂
p
/
∂
x
i
+
ρ
m
B
i
+
S
u
i
Energy equation(enthalpy form)
h
λ
eff
/
c
p
m
Q
′
′
′
Energy equation (temperature form)
T
λ
eff
/
c
p
m
Q
′
′
′
/
c
p
mThe rate of change (or time derivative) term is to be invoked only when a transient phenomenon is under consideration.The termρ
m
Φ denotes the amount of extensive property available in a unit volume.The convection second term takes into account the transport ofΦ due to bulk motion; this first-order derivative term is relatively uncomplicated but assumes considerable significance when stable and convergent numerical solutions are to be economically obtained [88]. The net source term implies an algebraic sum of sources and sinks of Φ; therefore, in a chemically reacting flow (combustion), a given specie k may be generated via some chemical reactions and destroyed (or consumed) via some others and R
k will comprise both positive and negative contributions. In addition, some chemical reactions may be exothermic and endothermic, making positive and negative contributions to Q
′′′ [88]. Similarly, the term B
i in the momentum equations may represent a buoyancy force, a centrifugal and/or Coriolis force, an electromagnetic force, and so forth. Sometimes, B
i may also represent resistance forces. The terms S
u
i represent viscous terms arising from Stokes’s stress laws. The transport equations described above must be extended with constitutive equations of state for density and for enthalpy; in this case, the standard Redlich-Kwong real gas model is used because it is considered as one of the most accurate models. Equation (8) of state is written as
(8)
p
=
R
T
ν
-
(
0.08664
R
T
c
/
p
c
)
-
(
0.42747
R
2
T
c
2
/
p
c
)
(
T
/
T
c
)
-
0.5
ν
(
ν
+
(
0.08664
R
T
c
/
p
c
)
)
,
where ν is the specific volume and p
c and T
c are critical pressure and temperature. Therefore, enthalpy (9), entropy (10), and the specific heat capacities at constant volume (11) and at constant pressure (12) can be calculated as(9)
h
(
T
,
v
)
=
∫
ν
ref
∞
(
T
(
d
p
d
T
)
V
-
p
)
d
ν
T
ref
+
∫
T
ref
T
c
V
0
d
T
-
∫
v
∞
(
T
(
d
p
d
T
)
V
-
p
)
d
ν
T
+
e
(
T
ref
,
v
ref
)
︸
+
p
v
e
(
p
,
T
)
(10)
s
(
T
,
v
)
=
∫
ν
ref
∞
(
d
p
d
T
)
V
d
ν
T
ref
+
∫
T
ref
T
c
p
0
T
d
T
-
R
ln
⁡
(
p
p
ref
)
-
∫
v
∞
(
d
p
d
T
)
V
d
ν
T
+
s
(
T
ref
,
v
ref
)
,
(11)
c
v
=
(
∂
u
∂
T
)
v
=
∂
∫
T
ref
T
(
c
p
0
(
T
)
-
R
)
d
T
∂
T
-
2
(
0.42747
R
2
T
c
2
/
p
c
)
(
T
/
T
c
)
-
0.5
(
0.08664
R
T
c
/
p
c
)
T
×
log
⁡
(
1
+
(
0.08664
R
T
c
/
p
c
)
v
)
,
(12)
c
p
=
c
v
+
v
T
(
(
-
(
∂
p
/
∂
T
)
v
v
(
∂
p
/
∂
v
)
T
)
2
×
(
-
1
v
(
∂
p
/
∂
v
)
T
)
-
1
)
,where c
V
0 and c
p
0 are the zero pressure ideal gas specific heat capacities. The greatest impediment to obtain physically accurate solutions is offered by the diffusion and the net source (S
Φ) terms of the transport equation because both of them require empirical information. In laminar flows, the diffusion term represented by the second-order derivative offers no difficulty because Γ is a fluid property and can be accurately determined (via experiments) isolating the flow under consideration. In turbulent (or transitional) flows, determination of Γ
eff requires considerable empirical support. This is named as turbulence modeling (see next paragraph).


## 3.2.3. Turbulence Modelling

A main key for the flow fields’ simulations is turbulence; it consists of fluctuations in the flow field in time and space [89]. It is a complex process, mainly because it is three-dimensional, unsteady, and multiscale. It can have a significant effect on the characteristics of the flow, and it occurs when the inertia forces in the fluid become significant compared to viscous forces (high Reynolds number). Two main approaches have been used to deal with this problem.(i)
Time-averaged simulations, where no eddy structures (vortices) are resolved; they employ a time averaged simulation with respect to the turbulence, and only mean statistics are predicted. The most common example is the Reynolds-averaged technique, where all flow variables are individually time averaged and the resulting equations are generally called the Reynolds-averaged Navier-Stokes (RANS) equations.
(ii)
Eddy-resolved simulations, where some or all of the eddy structures are resolved. This category includes large-Eddy simulations (LES) where the turbulence is only resolved up to some threshold wave numbers, beyond which a subgrid scale is employed; they are more accurate than the time-averaged techniques in terms of turbulent diffusion and other statistical flow features [90]. However, eddy-resolved formulations require a harsh price in terms of higher computational resources, both CPU memory and time.Reynolds-Averaged Navier-Stokes (RANS) Equations. Because turbulence consists of random fluctuation of various flow properties we use a statistical approach to this problem. The nonlinearity of Navier-Stokes equations leads to the appearance of momentum fluxes that act as apparent stresses throughout the flow; these momentum fluxes are unknown. We then derive equations for these stresses and the resulting equation includes additional unknown quantities. This problem is named “closure problem,” that is establishing a sufficient number of equations for all of the unknowns variables. Turbulence models try to solve a modified set of transport equations by introducing averaged and fluctuating components of velocity:
(13)
u
i
=
u
¯
i
+
u
i
′
=
1
Δ
t
∫
t
t
+
Δ
t
u
i
+
u
i
′
,
where u
,
u
¯, and u
′ are used to represent the instantaneous, average, and fluctuating terms, respectively. Δ
t is a time scale that is relatively large compared to the turbulent fluctuations but relatively small compared to the time scale to which the equations are solved. For simplicity, density fluctuations are considered negligible. Different from the Navier-Stokes equations, the momentum and scalar transport equations contain turbulent flux terms in addition to the molecular diffusive fluxes; these are the Reynolds stresses, ρ
u
ı
′
u
ȷ
′
-, that arise from the nonlinear convective term in the unaveraged equations and highlight the fact that convective transport due to turbulent velocity fluctuations enhances mixing over and above caused by thermal fluctuations at the molecular level. At high Reynolds numbers, turbulent fluxes are much larger than the molecular fluxes, as turbulent velocity fluctuations used to be larger than the mean free path of thermal fluctuations. The Reynolds averaged energy equation is
(14)
∂
ρ
(
h
+
(
1
/
2
)
u
i
u
i
+
(
1
/
2
)
u
ı
′
¯
2
︷
Turbulent
kinetic
energy
)
︷
Total
entalpy
,
h
tot
∂
t
-
∂
p
∂
t
+
∂
∂
x
j
(
ρ
u
j
h
tot
)
=
∂
∂
x
j
(
λ
∂
T
∂
x
j
-
ρ
h
u
ȷ
′
-
)
+
∂
∂
x
j
[
u
i
(
τ
i
j
-
ρ
h
u
ı
′
u
ȷ
′
-
)
]
︷
Viscous
work
term
+
Q
′′′
.
This equation contains an additional turbulence flux term ρ
h
u
ȷ
′
- compared with the instantaneous equation. Similarly, the Φ may be divided into an average component, Φ
¯, and a time varying component, φ
′, as
(15)
∂
(
ρ
Φ
)
∂
t
+
∂
(
ρ
u
j
Φ
)
∂
x
j
=
∂
∂
x
j
[
Γ
∂
Φ
∂
x
j
-
ρ
h
φ
′
-
]
+
S
Φ
,
where ρ
h
φ
′
- is the Reynolds flux. Turbulence models close the Reynolds averaged equations by providing models for the computation of the Reynolds stresses and Reynolds fluxes, based on the eddy viscosity hypothesis or the eddy diffusivity hypothesis, depending on the factor in which the stresses and fluxes are considered linearly dependent of (see [88]). Subject to these hypotheses, the Reynolds averaged energy equation becomes
(16)
∂
ρ
(
h
+
(
1
/
2
)
u
i
u
i
+
(
1
/
2
)
u
ı
′
¯
2
︷
Turbulent
kinetic
energy
)
︷
Total
entalpy
,
h
tot
∂
t
-
∂
p
∂
t
+
∂
∂
x
j
(
ρ
u
j
h
tot
)
=
∂
∂
x
j
(
λ
∂
T
∂
x
j
-
μ
t
P
r
t
∂
h
∂
x
j
)
+
∂
∂
x
j
[
u
i
(
τ
i
j
-
ρ
u
ı
′
u
ȷ
′
-
)
]
︷
Viscous
work
term
+
Q
′′′
.
Similarly, the Reynolds averaged transport equation for additional variables (nonreacting scalars) becomes
(17)
∂
(
ρ
Φ
)
∂
t
+
∂
(
ρ
u
j
Φ
)
∂
x
j
=
∂
∂
x
j
[
(
Γ
Φ
+
μ
t
σ
Φ
)
∂
Φ
∂
x
j
]
+
S
Φ
.Two-Equation Turbulence Models. As a good compromise between numerical effort and computational accuracy, two-equation turbulence models, where both the velocity and length scale are solved using separate transport equations, are widely used, like the RNG k
-
ε model, a derivation of the k
-
ε model; in the k
-
ε turbulence model, in fact, k
(
L
2
,
T
-
2
) is the turbulence kinetic energy and is defined as the variance of the velocity’s fluctuations and ε
(
L
2
,
T
-
2
) is the turbulence eddy dissipation (the rate at which the velocity’s fluctuations dissipate) [71, 91–93]. The k
-
ε model assumes that the turbulence viscosity is linked to the turbulence kinetic energy and dissipation via the relation:
(18)
μ
t
=
C
μ
ρ
k
2
ε
,
where C
μ is an empirical constant. The values of k and ε come directly from the differential transport equations for the turbulence kinetic energy (19) and turbulence dissipation rate (for full directional buoyancy model):
(19)
∂
(
ρ
k
)
∂
t
+
∂
(
ρ
u
j
k
)
∂
x
j
=
∂
∂
x
j
[
(
μ
+
μ
t
σ
k
)
∂
k
∂
x
j
]
+
μ
t
(
∂
u
i
∂
x
j
+
∂
u
j
∂
x
i
)
∂
u
i
∂
x
j
-
2
3
∂
u
k
∂
x
k
(
ρ
k
+
3
μ
t
∂
u
k
∂
x
k
)
︷
Turbulence
production
due
to
viscous
forces
,
P
k
-
μ
t
ρ
σ
ρ
g
i
∂
p
∂
x
i
︷
Buoyancy
producion
term
,
P
k
b
-
ρ
ε
,
(20)
∂
(
ρ
ε
)
∂
t
+
∂
(
ρ
u
j
ε
)
∂
x
j
=
∂
∂
x
j
[
(
μ
+
μ
t
σ
ε
)
∂
ε
∂
x
j
]
+
ε
k
(
C
ε
1
P
k
-
C
ε
2
ρ
ε
+
C
ε
1
·
max
⁡
(
0
;
P
k
b
)
sin
⁡
(
ϕ
)
)
.
The RNG-based k
-
ε turbulence model is derived from the instantaneous Navier-Stokes equations, using a mathematical technique called “renormalization group’’ (RNG) method; the transport equations for turbulence generation and dissipation are the same as those for the standard k
-
ε model, but the model constants differ, and the constant C
ε
1 is replaced by the function C
ε
1
RNG. The transport equation for turbulence dissipation becomes(21)
∂
(
ρ
ε
)
∂
t
+
∂
(
ρ
u
j
ε
)
∂
x
j
=
∂
∂
x
j
[
(
μ
+
μ
t
σ
ε
)
∂
ε
∂
x
j
]
+
ε
k
[
(
1.42
-
P
k
/
ρ
C
μ
RNG
ε
(
1
-
(
1
/
4.38
)
P
k
/
ρ
C
μ
RNG
ε
)
1
+
β
RNG
(
P
k
/
ρ
C
μ
RNG
ε
)
3
/
2
)
︷
C
ε
1
RNG
×
P
k
-
C
ε
2
RNG
ρ
ε
+
C
ε
1
RNG
·
max
⁡
⁡
(
0
,
P
k
b
)
·
sin
⁡
(
ϕ
)
[
(
1.42
-
P
k
/
ρ
C
μ
R
N
G
ε
(
1
-
(
1
/
4.38
)
P
k
/
ρ
C
μ
R
N
G
ε
)
1
+
β
R
N
G
(
P
k
/
ρ
C
μ
R
N
G
ε
)
3
/
2
)
︷
C
ε
1
R
N
G
]
.For all those turbulence models, compressibility corrections must be taken into account.


## 3.2.4. The LES Model

The LES approach is more general than the RANS approach and avoids the RANS dependence on boundary conditions for the large-scale eddies. The concept behind large-Eddy simulations is that only the large, energy-carrying scales of motion are computed exactly and the effect of the subgrid scales, which tend to be more isotropic and easier to parameterize than the large ones, is modeled. This is possible because energy transfer between motions at different scales occurs through vortex-stretching processes, so that the energy transfer is predominantly from larger to smaller scales since the motions at the largest scales derive their energy from the forces driving the flow field, whereas the relatively strong viscous forces at the smallest scales disperse the small scale motions in the energy of molecular motions [71]. The governing equations for LES are obtained by filtering the time-dependent Navier-Stokes equations in the physical space. The filtering process effectively filters out the eddies whose scales are smaller than the filter width or grid spacing used in the computations, so that the resulting equations govern the dynamics of the large eddies [88]. The filtering process can be resumed in the following equation:
(22)
Φ
-
(
x
)
=
1
Δ
∫
V
G
(
x
-
x
′
Δ
)
Φ
(
x
′
)
d
x
′
,
where the filter function G has to satisfy the following properties:(1)
G
(
-
x
)
=
G
(
x
);
(2)
∫
-
∞
∞
G
(
x
)
d
x
=
1;
(3)
G
(
x
)
→
0 as |
x
|
→
∞ sufficiently fast so that all moments ∫
-
∞
∞
G
(
x
)
x
n
d
x
=
1
(
n
≥
0
) exist;
(4)
G
(
x
) is localized in (
-
1
/
2
,
1
/
2
).Since the flow field of interest is compressible, a density-weighted filtering (Favre filtering) is adopted to avoid subgrid scale terms in conservation of mass.With a little modification to the definition, a Favre filtered variable can defined as(23)
Φ
~
(
x
)
=
1
ρ
¯
Δ
∫
V
G
(
x
-
x
′
Δ
)
ρ
Φ
(
x
′
)
d
x
′
,
(24)
u
k
=
U
k
+
u
k
=
ρ
u
k
-
ρ
¯
+
u
k
(
Velocity
)
,
s
i
j
=
S
i
j
+
u
i
j
=
ρ
s
ı
ȷ
-
ρ
¯
+
u
k
(
Shear
Strain
)
,
σ
i
j
=
Σ
i
j
+
u
i
j
=
ρ
σ
ı
ȷ
-
ρ
¯
+
u
k
(
Shear
Stress
)
,
ρ
=
ρ
+
ϱ
=
ρ
¯
+
ϱ
(
Density
)
,
p
=
P
+
p
=
p
¯
+
p
(
Pressure
)
,
T
=
T
+
τ
=
ρ
T
-
ρ
¯
+
τ
(
Temperature
)
,
e
=
E
+
e
=
ρ
e
-
ρ
¯
+
e
(
Internal
Energy
)
,
h
=
H
+
h
=
ρ
h
-
ρ
¯
+
h
(
Enthalpy
)
,
where the instantaneous values of variables are denoted inbold, the filtered values in capital letters, and the subgrid components in small letters [71].The Favre-filtered NSEs can be written as(25)
∂
(
ρ
)
∂
t
+
∂
(
ρ
U
k
)
∂
x
k
=
0
,
∂
(
ρ
U
i
)
∂
t
+
∂
(
ρ
U
i
U
j
)
∂
x
j
=
-
∂
P
∂
x
i
+
∂
(
μ
[
2
·
(
1
/
2
)
(
(
∂
U
i
/
∂
x
j
)
+
(
∂
U
j
/
∂
x
i
)
)
︷
Favre-filteres
shear
strain
,
S
i
j
-
(
2
/
3
)
δ
i
j
S
k
k
]
︷
Favre-filteres
shear
stress
,
Σ
i
j
)
∂
x
j
+
∂
(
-
ρ
(
(
u
ı
u
ȷ
~
)
-
U
i
U
j
)
)
︷
Subgrid-scale
stress
tensor
τ
i
j
∂
x
j while the Favre-filtered energy equation (26) can be expressed as
(26)
∂
(
ρ
E
)
∂
t
+
∂
(
ρ
U
k
E
)
∂
x
k
=
∂
∂
x
k
[
λ
¯
C
v
∂
E
∂
x
k
]
-
P
∂
U
k
∂
x
k
+
Σ
i
j
S
i
j
-
∂
(
(
ρ
U
k
E
-
)
-
ρ
U
k
E
)
︷
SGS
Heat
Flux
∂
x
k
-
[
(
p
∂
u
k
∂
x
k
-
P
∂
U
k
∂
x
k
)
]
︷
SGS
Pressure
Dilatation
+
(
(
σ
ı
ȷ
s
ı
ȷ
-
)
-
Σ
i
j
S
i
j
)
.
︷
SGS
Viscous
Dissipation
The subgrid-scale stress tensor needs to be modeled. It can be decomposed further, with respect to the filter dimensions, into contributions from various levels of interactions between the scales of motion (see [71, 88]), in the table that follows (Table 4) the main differences between RANS and LES.Table 4
Comparison between RANS and LES.
RANS
LES
The velocity field equations’ quantities are averaged with respect to time.
The equation of velocity field quantities is averaged with respect to space.
The eddy structures are not solved, and therefore the model has to be calibrated with respect to the type of turbulent current under consideration, wherefrom their behaviour depends.
The eddy structures within the inertial range of the turbulent spectrum are solved by using a sufficiently fine grid, while the subgrid scales of motion are parameterized and modeled, as they can be considered isotropic.
The motion can also be solved in two dimensions, if it is averagely static.
The motion can only be solved in three dimensions, because of the constant presence of the time derivatives due to the unstatic behaviour of the solved eddy structures.
The equations’ number is cut within a predetermined order, and the unknown variables belonging to a higher order are modeled with approximate relations.
The equations’ number is cut with respect to the position within the inertial range, so the difference between what is computed and what is modeled is energetic rather than geometric.
The closure problem is introduced by the Reynolds’s stresses’ tensor and is generally solved by using the turbulent cinematic viscosity, which often causes blunders.
The closure problem is introduced by the subgrid stresses’ tensor and is solved by subgrid models, which have no significant influence on the solution’s precision.
High efficiency, especially for the k-e based closure.
Extra efficiency and computational precision.
Low computational effort and CPU memory required.
High computational effort and CPU memory required.


## 3.2.5. Simulation Main Parameters

Three-Dimensional Domain Discretization. Some flow fields’ simulation issues, often those very remarkable in order to analyze and understand the real flow behavior and therefore to guarantee the pursued safety standards, need a much finer resolution of the simulations tools to be taken into account and identified. This is a consequence of the fact that these are too advanced to be detected at a rough level of the simulations; these issues are often related to situations where the 3D aspects of the flow and the geometrical effects have a significant influence, like the turbulence; that is why, a 3D axial-symmetric domain was considered in performing the CFD simulations. The spatial discretization consists of hexahedral elements for both the vessel domain and the horizontal inlet pipe. The 3D axial-symmetric domain was discretized (Table 5) using three structured grids: coarse (22.644 cells), fine (249.084 cells), and extra fine (2.264.402 cells). Grids are generated with an automatic procedure programmed in Gambit preprocessor. Considering the LES, the computational grid must be chosen such that the separation between the resolved and the subgrid scales occurs in the inertial subrange of the energy spectrum (where the statistics of motion are unaffected by the gross details of the flow field, unlike the large-scale motions and by the coefficient of viscosity, unlike the small-scale motions, but depend only on the average viscous dissipation of energy that determines the energy transfer across these scales). Accordingly, the fine grid size has been chosen to be one order of magnitude larger than that of the smallest scales.Table 5
Grid dimension.
Grid
Number of cells
Coarse
22.644
Fine
249.084
Extra fine
2.264.402Boundary Conditions. The first condition that has been considered is that a symmetry plane into two specular halves, only one of which has been analyzed, while the flow has been mirrored on the other side, as it has been done during the experimental campaign, has divided the tank. The normal velocity component at the symmetry plane boundary is set to zero u
n
=
0 and the scalar variable gradients normal to the boundary are also set to zero d
ϕ
/
d
n
=
0. At the inlet surface pressure and temperature are imposed. The total pressure (
P
=
101326
Pa
) is specified as an inlet boundary condition and the ANSYS-CFX solver computes the static pressure needed to properly close the boundary condition. The temperature is set to the environmental value (
T
=
296
K
). No-slip boundary conditions are assumed at all the walls. The velocity of the fluid at the wall is set to zero for all simulations, so the boundary condition for the velocity becomes u
wall
=
0. All the walls in the model, except for the heated wall in MC and AC experiments where a temperature condition is set, include the effect of heat transfer through the walls with the environment. A convection heat transfer coefficient (27) is considered for all the walls of the vessel [93]:
(27)
Nu
-
=
[
0.825
+
0.387
Ra
1
/
6
[
1
+
(
0.492
/
Pr
)
9
/
16
]
8
/
27
]
.
Heat flux at the wall boundary is calculated using q
w
=
h
c
(
T
b
-
T
n
w
), where h
c is a specified heat transfer coefficient (from Nusselt number), T
b is the specified boundary temperature (i.e, outside the fluid domain), and T
n
w is the temperature at the internal near-wall boundary element center node. The conditions of the environment are equal to T
=
296
K and P
=
101325 Pa (ambient conditions). The heated wall temperatures are maintained for both cases constant at fixed temperature, T
wall = 25°C for MC and T
wall = 110°C for AC, during the whole transient. As internal initial conditions, two different initial values, 100 and 1000 Pa, are considered in order to verify the effect of pressure on the pressurization of the vessel during the transient. The exact value of initial internal condition (due to the experimental difficulty of setting the internal pressure and temperature at the set point value) is provided by STARDUST data. The operating fluid used in the simulations is dry air assumed initially at rest. The density of the fluid is evaluated using the standard Redlich-Kwong real gas and the other thermodynamic properties are considered variable according to internal database.Time stepping. The selection of an appropriate time step size is essential in order to obtain good convergence. A too large time step will lead to poor convergence and a too small time step will lead to larger computational time. Within a given time step the transport equations are solved during an iteratively process until the convergence criteria are met for all equations. In these simulations a constant time step of 1 × 105 s is used. Each time step is considered to be converged when the maximum residual value is no higher than 1 × 10−6.


## 4. Results

### 4.1. Experimental Results

The experiments have been conducted in operative conditions (wall temperature of 120–130°C) in order to understand(1)
the behavior of thermofluid-dynamic parameters at different pressurization rates;
(2)
the behavior of thermofluid-dynamic parameters at different initial pressure rates inside the tank.In this section the several experimental setups will be singularly discussed and then compared in order to understand the main difference between them and their influence on dust resuspension (strictly connected to flow field variation).These analyses are implemented for the first 5 seconds from the beginning of LOVA because it has been already demonstrated [74] that the dust resuspension phenomena are relevant in the first 2–4 seconds after the LOVA.The elaborated data have been analysed by grouping the points in which the pressure transducers have been placed into four longitudinal lines or, otherwise, into four sections, following symmetrical criteria illustrated in the following paragraphs.The experimental trends of(i)
internal pressure and
(ii)
flow ratehave been collected in the mean time because these are essential data to give boundary condition necessary to develop the simulations.The main result of experimental data analysis is that the temperature is almost constant in the first 5 seconds, while internal pressure increases linearly (Figure18) and the flow rate increases asymptotically in this time range (Figure 19). The velocity values analysed show that velocity trend is characterized by a rapidly increase of velocity for the first 0,5 s (Figure 20), while it decreases more slightly after the third second; the trend of the velocity measured by the pressure transducer 461 (placed in the lateral points) is generally more gently (Figure 21).Figure 18
Internal pressure trend, 100Pa_100 Pa/s_cold experiment_points C-F.Figure 19
Flow rate trend, 100 Pa_100 Pa/s_cold experiment_points C-F.Figure 20
Velocity trend measured by 455, 100 Pa_500 Pa/s_hot experiment_point B.Figure 21
Velocity trend measured by 461, 100 Pa_500 Pa/s_hot experiment_point P.Finally, the most relevant velocity values in all the points have been compared and will be analysed in the following paragraphs.The Maximum Velocity Values along the Four Flow Lines of STARDUST. In this paragraph, the experimental data will be analysed considering the division of the points in which the pressure transducers have been placed into four longitudinal lines (Figure 22):(i)
Central line: points A-B-C-D;
(ii)
Lateral line: points G-F-E;
(iii)
Superior line: points H-I-L-M;
(iv)
Inferior line: points Q-P-O-N.This division allows observing the behavior of the air flow after its inlet, towards the lid of the tank and back.Figure 22
Division of the half tank into four lines.The following evidences have been observed.(i)
In the central line, the maximum velocity values tend to decrease from the inlet point to the lid. In the other lines the velocity values tend to decrease in the same direction. It means that the inlet airflow, which is initially curbed, is accelerated by the impact on the lid and inverts its direction, giving rise to a recirculation (see Figure23).
(ii)
The maximum velocity values are 20% to 50% higher in hot experiments rather than in cold experiments (the spread is larger when the pressurization rate is 500 Pa/s, for reasons explained in the following point), as clarified in Figure24, where the data collected from pressure transducer 461 (which faces the lid of the tank) are represented.
(iii)
In hot experiments, the pressurization rate does not have a relevant influence on the velocity value; anyway, the velocity tends to increase slightly with the pressurization rate.
(iv)
When the vacuum is broken at 1000 Pa the maximum velocity values are 30% lower than those at 100 Pa. It is a consequence of the lowest differences of pressure between the compressed air and the vacuum vessel of STARDUST. The trends of internal pressure inside STARDUST, considering an initial internal pressure both of 100 Pa and 1000 Pa, are comparable so the experiments with an initial internal pressure of 1000 Pa will not be reproduced for thermofluid-dynamic analysis.
(v)
When the velocity values are relatively high (and, generally, when the pressurization rate is 300 Pa/s or 500 Pa/s) the values measured in the section next to the lid of the tank are averagely too high in the central line and too low in the others, going against the general trends of the respective lines (see Figure24). It might be related to the fact that in this very section the data registered by the pressure transducers are affected by the turbulence and the velocity increase generated by the impact of the air on the lid, or even by the presence of the measuring instrument’s cables, that might break the flow.Figure 23
Comparison between hot and cold experiments at 100Pa_100 Pa/s for the lateral line.Figure 24
Maximum velocity values central line 100Pa_300 Pa/s_hot experiment.Analysis of the Maximum Velocity Values into the Four Flow Cross Sections of STARDUST. In this section, the experimental data will be analyzed considering the division of the points in which the pressure transducers have been placed into four sections (Figure 25):(i)
Section one: points A-H;
(ii)
Section two: points B-G-I;
(iii)
Section three: points C-F-L.Quarters of section have been considered taking into account the symmetry of the tank; the superior quarters have been analyzed.Figure 25
Division of the half tank into four sections.The analysis of these areas is necessary to observe the timing of the maximum velocity values in the various sections.The following evidences have been observed.(i)
In section one the velocity maximum is first verified in A and then in H (with a delay of about 0,4 s), because in H it is due to the flow coming back from the lid, while in A it is due to the inlet flow.
(ii)
In section two the velocity maximum in I and G is delayed of≈0,3 s compared to the maximum in B; this might not be verified if the velocity is relatively high.
(iii)
In sections three and four, the maximum velocity values are generally contemporary in the whole section. It can be a consequence of the closeness to the lid and so to the recirculation point.After 1 second from the beginning of the flow inlet, the following evidences have been observed.(i)
In section one, in experiments with pressurization rate of 100 Pa/s, the velocity in A is about 20% higher than that in H, while in experiments with pressurization rate of 300 Pa/s or 500 Pa/s it is about 85% higher. This is probably because at 100 Pa/s the pressure differences and the velocity are lower, so that the airflow is spread in all directions, losing its compactness.
(ii)
The velocity values measured in the central and lateral points of the section become comparable in section three, where the inlet flow and the flow proceed from the lid get melted.


### 4.2. Numerical Results

The CFD applications show stable numerical results for modeling local gas velocity field at low pressure conditions. Velocities remain very high for all the transients; Figure26 shows the streamlines. The axial velocity values, measured on the symmetry plane of STARDUST (see Figure 26), decrease rapidly with the distance. Considering a range of distances from 25 to 75 cm from Valve A, reduction of almost the 25% of the axial air velocities has been measured.Streamlines.
(a)
(b)
(c)
(d)The calculated data show a satisfying global agreement with experimental data for all performed experiments and simulation (from 1/3 to 2/3 of the vessel about 33%). The facility area with higher axial velocity values is the one close to the air inlet that corresponds to the expansion jet flow area. A recirculating flow can be observed close to the right side of the cylindrical wall (negative axial velocity).The code is able to assess the radial velocity that cannot be detected during the experimental tests. Observing the velocity profile inside the vessel, we can see a reduction of axial velocity at point for 500 Pa/s pressurization rate case. This is caused by the typical compressible flow effect named choked effect. In the case of the upstream air pressure and vacuum conditions downstream of an inlet section (orifice), both the air velocity and the mass flow rate become limited when sonic velocity is reached through the inlet section and the expansion inside the vessel becomes limited. This effect can be appreciated in Figure27 where the maximum numerical velocity is plotted versus the pressurization rate imposed. This means that the highest velocities, that is, the driver parameter for dust mobilization inside the vessel, are reached for low pressurization rate around 300 Pa/s. The calculations show that the air velocities always increase when the vessel is heated (AC). The air velocity reaches its peak value, 770 m/s for 100 Pa/s AC case, at about 0.7 s.(i)
The code predicts the oscillating phenomenon that is not observed in the experimental curve probably because of the low response speed of pressure probes. The presence of pressure oscillations is due to the inertia phenomena involved in this kind of problem and could cause the release of dust from the vessel, contributing to the source term in an ITER loss of vacuum scenario.
(ii)
The numerical results with heat transfer at walls show that the pressurization process has an adiabatic unsteady behavior, while the isothermal unsteady model overestimation is clearly visible. From these results it can be assumed that for pressurization transient the effects related to heat transfer through the walls to the environment can be neglected. The filling time is also underestimated for adiabatic unsteady case. The obtained results of the pressurization rate and filling time show good agreement with experimental data.
(iii)
The difference between the finest grid and the coarser grid shows an unstable LES solution for the coarse grid that seems prone to instability in low flow regions. We can claim that no great improvements seem to be made by using the extra fine grid. The fine grid is selected as a suitable compromise between accuracy and computational efficiency (1.6 day of CPU time for 1 s of physical time).
(iv)
From the streamlines analysis, after vacuum breakdown, the jet generates vortex structures that roll up several times and the inlet gas particles follow larger ones which brush and interact with all boundaries. It has relevant safety implication because the mobilized dust, which follows the streamline during them motion, can escape or be trapped in VV asperity and then accumulated. This high dust concentration zones increase the probability of chemical reaction because they increase reactive surface area that may contribute to hydrogen production and potential explosion.Figure 27
Simulation velocity trend at different pressurization rates.


## 4.1. Experimental Results

The experiments have been conducted in operative conditions (wall temperature of 120–130°C) in order to understand(1)
the behavior of thermofluid-dynamic parameters at different pressurization rates;
(2)
the behavior of thermofluid-dynamic parameters at different initial pressure rates inside the tank.In this section the several experimental setups will be singularly discussed and then compared in order to understand the main difference between them and their influence on dust resuspension (strictly connected to flow field variation).These analyses are implemented for the first 5 seconds from the beginning of LOVA because it has been already demonstrated [74] that the dust resuspension phenomena are relevant in the first 2–4 seconds after the LOVA.The elaborated data have been analysed by grouping the points in which the pressure transducers have been placed into four longitudinal lines or, otherwise, into four sections, following symmetrical criteria illustrated in the following paragraphs.The experimental trends of(i)
internal pressure and
(ii)
flow ratehave been collected in the mean time because these are essential data to give boundary condition necessary to develop the simulations.The main result of experimental data analysis is that the temperature is almost constant in the first 5 seconds, while internal pressure increases linearly (Figure18) and the flow rate increases asymptotically in this time range (Figure 19). The velocity values analysed show that velocity trend is characterized by a rapidly increase of velocity for the first 0,5 s (Figure 20), while it decreases more slightly after the third second; the trend of the velocity measured by the pressure transducer 461 (placed in the lateral points) is generally more gently (Figure 21).Figure 18
Internal pressure trend, 100Pa_100 Pa/s_cold experiment_points C-F.Figure 19
Flow rate trend, 100 Pa_100 Pa/s_cold experiment_points C-F.Figure 20
Velocity trend measured by 455, 100 Pa_500 Pa/s_hot experiment_point B.Figure 21
Velocity trend measured by 461, 100 Pa_500 Pa/s_hot experiment_point P.Finally, the most relevant velocity values in all the points have been compared and will be analysed in the following paragraphs.The Maximum Velocity Values along the Four Flow Lines of STARDUST. In this paragraph, the experimental data will be analysed considering the division of the points in which the pressure transducers have been placed into four longitudinal lines (Figure 22):(i)
Central line: points A-B-C-D;
(ii)
Lateral line: points G-F-E;
(iii)
Superior line: points H-I-L-M;
(iv)
Inferior line: points Q-P-O-N.This division allows observing the behavior of the air flow after its inlet, towards the lid of the tank and back.Figure 22
Division of the half tank into four lines.The following evidences have been observed.(i)
In the central line, the maximum velocity values tend to decrease from the inlet point to the lid. In the other lines the velocity values tend to decrease in the same direction. It means that the inlet airflow, which is initially curbed, is accelerated by the impact on the lid and inverts its direction, giving rise to a recirculation (see Figure23).
(ii)
The maximum velocity values are 20% to 50% higher in hot experiments rather than in cold experiments (the spread is larger when the pressurization rate is 500 Pa/s, for reasons explained in the following point), as clarified in Figure24, where the data collected from pressure transducer 461 (which faces the lid of the tank) are represented.
(iii)
In hot experiments, the pressurization rate does not have a relevant influence on the velocity value; anyway, the velocity tends to increase slightly with the pressurization rate.
(iv)
When the vacuum is broken at 1000 Pa the maximum velocity values are 30% lower than those at 100 Pa. It is a consequence of the lowest differences of pressure between the compressed air and the vacuum vessel of STARDUST. The trends of internal pressure inside STARDUST, considering an initial internal pressure both of 100 Pa and 1000 Pa, are comparable so the experiments with an initial internal pressure of 1000 Pa will not be reproduced for thermofluid-dynamic analysis.
(v)
When the velocity values are relatively high (and, generally, when the pressurization rate is 300 Pa/s or 500 Pa/s) the values measured in the section next to the lid of the tank are averagely too high in the central line and too low in the others, going against the general trends of the respective lines (see Figure24). It might be related to the fact that in this very section the data registered by the pressure transducers are affected by the turbulence and the velocity increase generated by the impact of the air on the lid, or even by the presence of the measuring instrument’s cables, that might break the flow.Figure 23
Comparison between hot and cold experiments at 100Pa_100 Pa/s for the lateral line.Figure 24
Maximum velocity values central line 100Pa_300 Pa/s_hot experiment.Analysis of the Maximum Velocity Values into the Four Flow Cross Sections of STARDUST. In this section, the experimental data will be analyzed considering the division of the points in which the pressure transducers have been placed into four sections (Figure 25):(i)
Section one: points A-H;
(ii)
Section two: points B-G-I;
(iii)
Section three: points C-F-L.Quarters of section have been considered taking into account the symmetry of the tank; the superior quarters have been analyzed.Figure 25
Division of the half tank into four sections.The analysis of these areas is necessary to observe the timing of the maximum velocity values in the various sections.The following evidences have been observed.(i)
In section one the velocity maximum is first verified in A and then in H (with a delay of about 0,4 s), because in H it is due to the flow coming back from the lid, while in A it is due to the inlet flow.
(ii)
In section two the velocity maximum in I and G is delayed of≈0,3 s compared to the maximum in B; this might not be verified if the velocity is relatively high.
(iii)
In sections three and four, the maximum velocity values are generally contemporary in the whole section. It can be a consequence of the closeness to the lid and so to the recirculation point.After 1 second from the beginning of the flow inlet, the following evidences have been observed.(i)
In section one, in experiments with pressurization rate of 100 Pa/s, the velocity in A is about 20% higher than that in H, while in experiments with pressurization rate of 300 Pa/s or 500 Pa/s it is about 85% higher. This is probably because at 100 Pa/s the pressure differences and the velocity are lower, so that the airflow is spread in all directions, losing its compactness.
(ii)
The velocity values measured in the central and lateral points of the section become comparable in section three, where the inlet flow and the flow proceed from the lid get melted.


## 4.2. Numerical Results

The CFD applications show stable numerical results for modeling local gas velocity field at low pressure conditions. Velocities remain very high for all the transients; Figure26 shows the streamlines. The axial velocity values, measured on the symmetry plane of STARDUST (see Figure 26), decrease rapidly with the distance. Considering a range of distances from 25 to 75 cm from Valve A, reduction of almost the 25% of the axial air velocities has been measured.Streamlines.
(a)
(b)
(c)
(d)The calculated data show a satisfying global agreement with experimental data for all performed experiments and simulation (from 1/3 to 2/3 of the vessel about 33%). The facility area with higher axial velocity values is the one close to the air inlet that corresponds to the expansion jet flow area. A recirculating flow can be observed close to the right side of the cylindrical wall (negative axial velocity).The code is able to assess the radial velocity that cannot be detected during the experimental tests. Observing the velocity profile inside the vessel, we can see a reduction of axial velocity at point for 500 Pa/s pressurization rate case. This is caused by the typical compressible flow effect named choked effect. In the case of the upstream air pressure and vacuum conditions downstream of an inlet section (orifice), both the air velocity and the mass flow rate become limited when sonic velocity is reached through the inlet section and the expansion inside the vessel becomes limited. This effect can be appreciated in Figure27 where the maximum numerical velocity is plotted versus the pressurization rate imposed. This means that the highest velocities, that is, the driver parameter for dust mobilization inside the vessel, are reached for low pressurization rate around 300 Pa/s. The calculations show that the air velocities always increase when the vessel is heated (AC). The air velocity reaches its peak value, 770 m/s for 100 Pa/s AC case, at about 0.7 s.(i)
The code predicts the oscillating phenomenon that is not observed in the experimental curve probably because of the low response speed of pressure probes. The presence of pressure oscillations is due to the inertia phenomena involved in this kind of problem and could cause the release of dust from the vessel, contributing to the source term in an ITER loss of vacuum scenario.
(ii)
The numerical results with heat transfer at walls show that the pressurization process has an adiabatic unsteady behavior, while the isothermal unsteady model overestimation is clearly visible. From these results it can be assumed that for pressurization transient the effects related to heat transfer through the walls to the environment can be neglected. The filling time is also underestimated for adiabatic unsteady case. The obtained results of the pressurization rate and filling time show good agreement with experimental data.
(iii)
The difference between the finest grid and the coarser grid shows an unstable LES solution for the coarse grid that seems prone to instability in low flow regions. We can claim that no great improvements seem to be made by using the extra fine grid. The fine grid is selected as a suitable compromise between accuracy and computational efficiency (1.6 day of CPU time for 1 s of physical time).
(iv)
From the streamlines analysis, after vacuum breakdown, the jet generates vortex structures that roll up several times and the inlet gas particles follow larger ones which brush and interact with all boundaries. It has relevant safety implication because the mobilized dust, which follows the streamline during them motion, can escape or be trapped in VV asperity and then accumulated. This high dust concentration zones increase the probability of chemical reaction because they increase reactive surface area that may contribute to hydrogen production and potential explosion.Figure 27
Simulation velocity trend at different pressurization rates.


## 5. Discussion

Numerical and experimental results using initial pressure equal to 100 Pa in hot conditions, at different pressurization rates, will be discussed because they are considered most relevant for the validation case. Unless specified, the qualitative results must be considered valid for the nondiscussed tests.The following evidences have been observed.(i)
The average pressure inside the vessel during the simulations shows a linear trend. In the figures, which follow Figures28 and 29, the trend with different turbulence models and grids is, respectively, showed.
(ii)
Figures28 and 29 show that the influence of the turbulence model on pressure trend is negligible.
(iii)
The curves obtained in the performed simulations demonstrate that the LES computations are consistent and show generally better agreement with the experiments than RANS.
(iv)
The fine grid and the extra fine grid do not present substantial differences, so the first one has been chosen in order to reduce the time machine.By using the results obtained with this numerical model we have analyzed all the velocity trends with an initial pressure of 100 Pa and 1000 Pa, with different pressurization rates (100 Pa/s; 300 Pa/s and 500 Pa/s). We observed a good match between numerical and experimental data for all the points in which the PTs have been placed.Figure 28
Average pressure time trend for different turbulence models (100 Pa, AC).Figure 29
Average pressure time trend for different grids (100 Pa, AC).Figures30 and 31 show a comparison graph of velocity for point A at 100 Pa/s and 300 Pa/s with an initial pressure of 100 Pa/s.(i)
The axial velocity evaluated close to the air injection increases from zero to the maximum in the first second and then decreases. The CFD applications show stable numerical results for modeling local gas velocity field at low pressure conditions. Velocities remain very high for the entire transient; Figure30 shows the axial velocity magnitude history for the initial and final parts of the vessel. Assuming that the axial velocity profiles are symmetric with respect to the symmetry plane the velocity’s direction falls relatively rapidly with distance from the injection section along the axis (from 1/4 to 3/4 of the vessel about 25%). The calculated data show a satisfying global agreement with experimental data for all performed experiments and simulation.
(ii)
The temperature inside the vessel seems to be almost constant with a slight negative slope for the operative cases.
(iii)
For the experiments in operative conditions the average temperature decreases from initial internal temperature (384 K) to 305 K (12%) in about 18 sec. Results appear to be concordant, as simulations satisfactorily reproduce experimental results only for MC case. Agreement seems to be poor for AC case. This poor agreement may be caused by thermal high thermal inertia of TJ thermocouples.The calculated data show a satisfying global agreement with experimental data for all performed experiments and simulations.Figure 30
Velocity comparison at 100 Pa/s.Figure 31
Velocity comparison at 300 Pa/s.


## 6. Conclusion and Future Developments

The past experimental campaign has showed the validation of the 2D model to reproduce flow fields condition comparable to those obtained in STARDUST with a LOVA of 300 Pa/s. The authors demonstrate that, by means of STARDUST facility, it is possible to reproduce simulations of LOVA with different pressurization rates. It is possible to conduct experiments at 100 Pa/s, 300 Pa/s, and 500 Pa/s in order to consider also a little smaller and a little larger break in the ITER’s vessel walls and it is possible to reproduce different boundary conditions (like initial pressure and internal temperatures). The most relevant evidences that have been deducted from the data analysis are that the air flows through the central part of the vessel towards the lid and then back along the tank’s walls, after being reflected by the lid showing an air velocity higher (≈30%) in hot experiments rather than in cold experiments, which means that the problems caused by dust mobilization are higher in case of an accident occurring during the standard functioning of the reactor. By the experiments and simulations it has been showed also that the air velocity is lower (≈30%) when the initial pressure is 1000 Pa rather than 100 Pa and that the influence on the air velocity of the pressurization rates is much more remarkable in maintenance conditions, when the maximum is reached at 300 Pa/s, while in hot experiments the maximum is generally reached at 500 Pa/s. Within these experimental campaigns numerical 3D simulations have also been carried on by the mean of the ANSYS-CFX code in order to add missing information to the experimental scenario and to develop and validate a CFD code that can be a useful resource to handle LOVA problems without the experimental support. The data extracted by the simulations at the same points where the pressure transducers had been placed during the experimental campaign showed a satisfying agreement with the experimental ones and confirmed the evidences highlighted by the data analysis, adding information about the radial velocity, otherwise undetectable, allowing for a complete mapping of the flow field and velocity inside the tank. A numerical 3D model has been created in order to reproduce any fluid-dynamic conditions subsequent to LOVAs accidents. The most important future developments will develop a new geometry configuration of STARDUST, which has been already projected, making the choice to change the injection section on STARDUST in order to assimilate it close to a sector of a toroidal geometry and to avoid the wall effects of the lids.



---

*Source: 201831-2014-05-05.xml*