3D transient CFD simulation of an in-vessel loss-of-coolant accident in the EU DEMO fusion reactor

Accidental transients in fusion reactors, such as the in-vessel loss-of-coolant accident (LOCA) considered here, are generally analyzed using system-level codes. However, because of their lumped nature, such tools cannot predict local pressure and temperature values, which are instead of interest to assess the integrity of the vacuum vessel (VV) structures. It is then fundamental to prove that the system-level tools’ predictions are at least conservative. In this work, we analyze the helium flow inside the VV following a LOCA in the helium-cooled blanket of the EU DEMO, using a 3D transient computational fluid-dynamics (CFD) model implemented in the commercial STAR-CCM+ code. In view of the large pressure ratio, a hypersonic flow regime (Ma > 5) develops, with the formation of shock fronts requiring a high time and space resolution. The model is applied to compute the evolution of the pressure distribution inside the VV and then to compare it with the information provided by the system-level GETTHEM model. The predictions by the two models of the intervention time of the VV pressure suppression system are compared, showing that the 0D model is conservative.


Introduction
Safety is one of the main design drivers of the EU DEMO reactor [1] by the EUROfusion Consortium, since the preconceptual phase.
One of the design basis accidents is the in-vessel loss-ofcoolant accident (in-VV LOCA), occurring when a break in Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. a part the cooling system (e.g. due to a rupture of either an in-vessel component or an in-vessel pipe) causes a release of pressurized coolant inside the vacuum vessel (VV).
The subsequent pressurization of the VV must be carefully analyzed to confirm its confinement function of hazardous material is respected, as the VV is the first barrier for the containment of radioactive contaminants. Such analyses are generally carried out using system-level computational tools, such as RELAP [2], MELCOR [3][4][5] or GETTHEM [6][7][8], which adopt a lumped-parameter 0D/1D modelling of the relevant components. The information provided by such tools is however limited to volume-averaged values, and cannot predict local peaks (of pressure and/or temperature) on the VV surfaces. Due to the large computational effort required, 3D computational fluid-dynamics (CFD) modeling of transients of this kind is generally avoided, with few exceptions limited to smaller domains, such as experimental facilities, see e.g. [9,10], where a loss-of-vacuum accident is simulated numerically and experimentally in the STARDUST facility.
The failure criterion for the VV is instead based on a local value, i.e. the pressure on the gyrotron diamond windows, which should not overcome 2 bar. Moreover, as the pressure distribution may be strongly nonuniform, using the average value to predict the intervention delay of safety systems (e.g. burst disks or active valves) may introduce large under-or over-estimations, so that the volumeaveraged system-level approach needs at least to be proved conservative.
In the present work, a 3D transient CFD model of the EU DEMO VV is built and used to simulate for the first time the flow evolution following a LOCA, in order to explore the physics of the complex hypersonic flow resulting from the rupture. The CFD results are then compared with those obtained with GETTHEM, in order to quantitatively assess the limitations of a system-level modeling approach. The main physical phenomena are recalled in appendix A, while the solvers and algorithms developed to solve the problem and their verification and validation in a simple 2D steady-state case are described in appendixes B and C, respectively.

Scenario and geometry description
The accident considered for the present analysis initiates from a failure of a portion of the first wall (FW) of the EU DEMO helium-cooled pebble bed (HCPB) breeding blanket (BB) [11], either due to an unmitigated plasma disruption or a beam of runaway electrons causing the melting of the EUROFER97 FW structure.
The helium used as coolant (nominal parameters: 8 MPa, 300 ÷ 520 • C) is then suddenly released inside the plasma chamber, see figure 1, which pressurizes as the coolant expands. The main features of the resulting highly underexpanded jet are recalled in appendix A. To avoid overpressurization, the VV is equipped with a pressure suppression system, i.e. a low-pressure environment separated from the VV via one or more burst disks (BD), passively opening when the pressure reaches 1.5 bar, and one or more active 'bleed valves', useful in case of small leakages.
In the present analysis, the transient starts when the break occurs (instantaneously at t = 0 s) and ends when the volumeaveraged pressure reaches 1.5 bar (i.e. when a 0D model would predict the opening of the BD). The intervention of bleed valves is not considered here, since it was already shown to have a negligible effect in the case of a large break [6,13].
The fluid domain for a 3D analysis should include the whole free volume inside the VV (i.e. all the volume available for the coolant to expand). However, in view of the complexity of the problem, several simplifications are introduced in the geometry, as reported in the following, aiming at a reduction of the computational cost.
The different BB segments are separated by tight gaps (2 cm in non-operating conditions), but the flow in the gaps is expected to be negligible with respect to the flow in the main plasma chamber. This not only because of the gap size, but also because (a) the actual gap dimension may be very different in operating conditions, due to electromagnetic loads and thermal expansion (which may even completely close the gaps), and (b) part of the gap between the VV and BB will be filled with attachments and other equipment, which are not designed in detail yet.
In view of this assumption, a first simplification is to neglect the gaps and keep only the main plasma chamber as fluid domain. Then, as the up-down asymmetry is limited to the divertor region, to further reduce the size of the domain a simplified, up-down symmetric geometry is built, according to figure 2, keeping the same total volume (divided by four, due to up-down and left-right symmetry) of the VV (~3000 m 3 ).
As reference scenario, the worst case among those identified as design-basis accidents in [14] is chosen, involving a failed FW surface of 5 m 2 on the outboard (OB) BB. Considering the internal geometry of the HCPB, this corresponds to exposing 1031 cooling channels to the plasma chamber, so a total flow area of~0.32 m 2 is used [13,15]. In this model, for the sake of simplicity, this flow area is lumped in a circular inlet on the OB side, see figure 2, with the coolant entering the domain in the radial direction. This is indeed a further simplification of the problem, as in real conditions the coolant would flow in the chamber from several small channels spread across the 5 m 2 of failed FW. In addition, according to the FW failure mode (e.g. melting, sublimation, burst, …), the flow might not be totally in the radial direction but could have a nonnegligible initial toroidal component. However, we expect it to have essentially no impact on the main outcomes of this work, i.e. on capturing the local (wall) effects of an in-VV LOCA and on comparing the 3D model with a 0D one.
The suppression system is assumed to be connected to the plasma chamber via the neutral beam injector openings, with one BD having an area of 1 m 2 [7,15]. In this work, the worst case is assumed for the location of the break with respect to the BD, i.e. they are assumed at a toroidal distance of 180 • , see figure 2, because in this case the distance to be travelled by the pressure wave to reach the BD is maximized.

Simulation setup
In this section, the main aspects of the setup of the CFD model are described, while the details can be found in appendix B. The software used for the simulation is STAR-CCM+ v. 2019. 3 [16]. The boundary conditions adopted in this simulation are reported in figure 2: 'Stagnation inlet' is the boundary condition used to impose sonic or supersonic flow at an inlet in STAR-CCM+ [17].
Since a supersonic flow is expected, both the static pressure p s = 8 MPa (called 'supersonic static pressure' in STAR-CCM+) and the total pressure p t must be prescribed at the inlet, according to where γ is the ratio of the specific heats of the gas and Ma is the Mach number (ratio of fluid speed at the boundary to speed of sound in the medium). Since the inlet section is the minimum section, sonic flow conditions are expected, so Ma = 1 is imposed in equation (1). Also the total temperature T t is specified at the inlet, according to where T s is the static temperature (410 • C) and Ma = 1 for the reason above. The imposed supersonic static pressure of 80 bar leads to an initial pressure ratio equal to 800 (see section 3.1.2 below).
The assumption of constant inlet pressure is reasonable, as the time span of interest for this simulation is limited to few tens of milliseconds, whereas the timescale for the depressurization of the Primary Heat Transfer System is of some seconds [6,8,13].
No outlet boundaries have been considered, since the domain is closed (except for the inlet) until the burst disks break, and the phase of the transient following that moment is beyond the scope of this work.

Initial conditions.
Within the VV, for computational reasons, it is assumed the presence of He instead of vacuum. The He speed is set to zero (stagnation) at a static pressure of 0.1 bar and a temperature equal to 300 K.
Note that the pressure within the VV can change from 1 mPa during standby up to 1-10 Pa during the plasma operation. However, 0.1 bar was chosen in order to have a sufficiently low Knudsen number (Kn~10 −7 ), so that the continuum assumption holds. Recall that Kn is defined according to equation (3): where λ is the molecular mean free path, L is a characteristic length (i.e. the distance between the inlet and the wall in front of it), k B is the Boltzmann constant, T is the temperature of the system, σ is the particle diameter and p is the pressure of the system. The cross-check of this assumption asks for the development of different physical models able to cope with high Kn (such as the direct simulation Monte Carlo [18,19]), which is however beyond the scope of this work. This assumption may produce a shift in the timing (i.e. to rupture the BDs), but does not affect the comparison with the 0D model, which is run starting from the same initial condition.
Concerning the initial temperature condition, the actual value in such scenario is not known (as well as its distribution in the domain). However, the initial mass of He in the domain is small compared to the mass at the end of the transient, so that also the total (initial) energy of He in the domain is negligible. Therefore the initial temperature value is not expected to significantly affect the results. Indeed, different 0D analyses of in-vessel LOCAs found in literature assume different initial VV temperatures, in the range~290 K to~600 K (see e.g. [3,5,6].), but the VV temperature evolution is almost independent on the initial condition, see e.g. [6]. Furthermore, the comparison with the 0D model is consistent since the same initial conditions are chosen.

Model and solvers
3.2.1. Model. Several features of the problem at hand call for a careful selection of the models.
Since the pressure ratio between the upstream conditions and those of the chamber is much larger than 2, a supersonic flow is expected. This means that the flow will for sure be compressible so that the pressure and velocity (or continuity and momentum) equations need to be coupled from the beginning to the energy equation and to the equation of state of the fluid.
A first choice must be made on the segregated vs. coupled approach to solve the set of Navier-Stokes equations (including the energy equation). In the case of compressible, supersonic flows, it is advisable to opt for a coupled approach (enabled with the 'coupled flow' and 'coupled energy' features in the selected software), which is for sure heavier from a computational point of view, but guarantees better stability of the computation, since the pressure-velocity-temperature fields are strongly coupled.
Concerning the equation of state, a simple ideal gas law has been chosen. This choice is due to the fact that the fluid entering the chamber can be considered ideal, being at high temperature. Once in the chamber, it will expand towards low pressure, remaining in the ideal gas range. However, in the high Mach regions, very low temperatures are expected. Therefore, in this region a cross-check of the validity of the model is envisaged.
The local Reynolds number computed according to equation (4): where ρ is the fluid density, v is the average speed and µ is the fluid dynamic viscosity, is >10 6 inside the jet region throughout the transient, so the flow regime is fully turbulent and, based on [20], the k-ω SST (Menter) [21] turbulence model has been chosen for the analysis. In addition, since an adaptive mesh is employed, a 'solution interpolation' model is enabled. This model is needed simply because the solution must be interpolated on the new mesh each time the mesh is adapted. Despite the easiness to understand the need for this model, the choice of the settings of this model is insidious. By default, the interpolation is performed on the computed quantities, such as pressure, temperature and the components of the velocity. This strategy guarantees the smoothness of the corresponding fields before and after the interpolation. However, it does not guarantee the conservation of the mass, momentum and energy. Since several mesh adaptations are expected, large errors would be introduced. Therefore, in order to keep the mass, momentum and energy conservation before and after remeshing the domain, the 'conservative interpolation' option is selected, which guarantees the conservation of mass, momentum and energy before and after the interpolation.

Solvers.
According to the complexity of the problem at hand, also the solver settings need careful tuning. The (modified) solvers and solver settings adopted in this paper are summarized below, while the details are collected in appendix B. Due to the nature of the phenomenon, i.e. propagation of shock fronts, a dedicated mesh adaptivity strategy is needed to optimize the mesh generation in order to follow such fronts. The alternative, sort of brute force approach, would be to mesh very finely the entire domain. The inlet region is meshed with a static mesh, see figure 3. The characteristic mesh size adopted in this region gives an average volume of each mesh polyhedron of~0.6 cm 3 . This means that, following a brute force approach, meshing the entire computational domain (~750 m 3 ) would end up with more than 1 billion cells. In addition, the problem at hand needs a very small time step to follow the propagation of several shock fronts, leading to a very large number of time steps to solve the entire transient. This leads to the impossibility of solving the problem with a static mesh and calls for the development of an algorithm to suitably adapt the mesh based on the solution throughout the transient.
Only the start-up of the simulation is performed with the static mesh with only the inlet region finely meshed.
After the simulation start-up, the meshing strategy is switched to adaptive, thus each 10 time steps the domain is re-meshed according to the refinement algorithm chosen. The local mesh refinement at the inlet is kept constant, since the inlet is a crucial location throughout the entire transient.
The mesh refinement algorithm is described in detail in appendix B.

Results
In this section, the main results obtained with the detailed 3D simulation are discussed.

Flow field
A video showing the time evolution of the jet throughout the transient as well as the corresponding dynamic mesh adaptation is available in [22] (see supplementary data, available online at (stacks.iop.org/NF/60/126001/mmedia)); the noteworthy features of the transient are also reported and commented below.
The initial phase of the evolution of the jet is shown in figure 4. As discussed in appendix A, as soon as the high pressure fluid is free to expand, it develops an almost planar front, see figure 4(a), which is not yet affected by the finite shape of the aperture. This phase is short, as the jet starts spreading (or diffracting) on the sides and the leading front develops a spherical shape, see figure 4(b). Recall that the jet startup features a leading shock, followed by a second shock, i.e. the Mach disk: this is clearly visible in figure 4(c). Furthermore, the onset of the slipstream emanated from the triple point, i.e. the interception point of Mach disk, intercepting shock and reflected shock, is also visible. In this phase, the evolution of the jet is close to a free underexpanded jet startup, since the information about the confined geometry has not reached the jet yet.
Note that also the mesh has been plotted in figures 4(a)-(c), in order to appreciate the capability of the adopted algorithm to follow the jet propagation.
The evolution of the jet after the impact on the VV inboard surface is shown in figure 5. The leading shock is the first impacting on the walls, see figure 5(a), causing a sharp peak in the pressure on the wall surface, see next section. It reflects back, see figure 5(b), impacting on the Mach disk that was developing and propagating right behind it, see figure 5(c). Consequently, the Mach disk is strongly deformed, adapting to the shape of the VV inboard surface, without, however, ever 'touching' it. This is a positive effect, since the Mach shock is stronger than the leading one [23], i.e. it features stronger variations of pressure and density and it travels at a higher speed. In figures 5(d)-(f ), a first spreading of the jet in the lateral direction is visible as well as a subsequent shrinkage due to the presence of the VV outboard walls, i.e. the jet lateral boundaries shrink, moving away from the outboard wall (compare the jet lateral extension in figures 5(e) and (f )). The clear formation of the triple point can also be appreciated, from which a strong slipstream starts flowing away from the jet core. Note also that, outside the jet boundary, a strong recirculation takes place. This justifies the adoption of a viscous turbulent flow model rather than an inviscid (Euler) approximation, which would have worked well if only the jet evolution was of interest.
The further evolution of the jet is shown in figure 6. It can be seen that the jet shrinks, i.e. the Mach disk moves far away from the inboard surface, still keeping all the macroscopic characteristic features discussed in appendix A. This is due to the increase of the average pressure in the volume, which reduces the pressure ratio (which is in turn the parameter that defines the Mach disk location, see appendix C). The average pressure reaches 1.5 bar at t = 53 ms, see figure 9, which reduces the pressure ratio from 800 to about 50. From [24], a pressure ratio of 50 would lead at steady-state to a Mach disk positioned at about 4.5 m from an exit with diameter 0.6 m. This is actually not far from the situation shown in figure 7, where the Mach number distribution at t = 53 ms is shown. Indeed, the Mach disk is located at 4.3 m from the inlet (see figure 7(b)), i.e. less than 5% different from the expected value. However, the case at hand is transient, as the average pressure in the VV is increasing and the jet has undergone an impingement on the inboard VV wall, thus it is not easily comparable to a steady state result. Additionally, the curvature of the top part of the VV induces recirculation of the He back towards the inlet, see figure 7(a). Nevertheless, the increase of the pressure in the volume leads to a reduction of the jet dimensions, as expected.

Pressure field
The pressure distribution at t = 53 ms is reported in figure  8(a). The inlet portion has been cut away since it has very high values (up to 80 bar, which is the inlet pressure boundary      condition) and they are not of interest for the present discussion. The strong expansion of the jet is evident: it goes from 80 bar at the inlet down to 0.1 bar (assumed initial pressure of the chamber for this exercise) in~2.5 m. A large portion of the jet is at very low pressure, which is consistent with the corresponding very high speed, since the energy content switches from being mainly internal at the inlet to mainly kinetic close to the shock front. After the shock, close to the inboard wall of the chamber, the pressure rises again, see also figure 8(b), to~3 bar. However, downstream of the shock, and in general outside the jet, the pressure distribution varies strongly with time, see also the next section.

Pressure evolution.
The pressure can significantly vary as the pressure waves pass by at a given location. The interest is to monitor the pressure on the walls of the outboard portion of the chamber, since the weakest point, e.g. the apertures for the plasma heating devices (including the gyrotron diamond windows), are located there, as well as on the regions close to the rupture, where the highest value of pressure is reached, e.g. on the wall in front of the rupture itself.
Note that the pressure evolution in such locations is independent of the relative position of the rupture and BD (as the latter is still closed), and consequently these results depend on the position only relatively to the break location.
The pressure evolution on the outboard side of the chamber (on the equatorial plane) has been monitored at four different points on the equatorial plane, see figure 9(a): assuming a cylindrical reference system centered in the center of the chamber and the azimuthal coordinate starting from the position of the rupture, the four monitored points are at 20 • , 80 • , 140 • , 180 • from the inlet, respectively. It can be noticed that, as the pressure wave travels in the chamber, the pressure on the selected points raises subsequently, see the 'A' peaks in figure  9(a). Once the pressure reaches the opposite side with respect to the inlet, it reflects back, leading to a second pressure rise, see the 'B' peaks in figure 9(a). This is the reflection of the information that the domain is closed, which reaches the inlet region (which is the last portion of the domain that gets the information on how big itself domain is) after 30-40 ms. This is consistent with a rough estimation, taking into account the length which the wave needs to travel from the inlet to the BD and back to the inlet, i.e. L~2πR, and the speed of sound, i.e. c~2200 m s −1 , which gives a transit time of L/c~34 ms.
The points on the chamber walls where the largest pressure values are reached are close to the inlet, see figure 9(b). The maximum pressure is reached in front of the rupture and it spikes up to~3.5 bar after only 2 ms; this is consistent considering the distance from the inlet L IB~6 m and the (average) speed of sound of~3000 m s −1 . The pressure wave travels along the inboard wall but with a progressively lower intensity: after few milliseconds, the peak on the inboard surface 20 • from the point in front of the rupture is~2 bar. The pressure at the top of the chamber reaches a very small peak in the beginning and stays at about 0.5 bar until the global pressure rise begins, i.e. at t~35 ms.
After the first 30-40 ms, the pressure evolution within the entire volume and in particular on the wall (except for the region close to the inlet) starts following the average pressure increase, which is no more linked to the travelling pressure wave, i.e. to local effects, but to the overall pressurization of the domain.

Temperature field
The He temperature distribution at t = 53 ms is shown in figure 10(a). Note that both very high and very low temperatures are reached. Where the jet is expanding more, e.g. inside the jet itself, the He cools down to few tens of K. In these same regions, due to the very low temperatures, the Mach number reaches the maximum values, see again figure 7(a), as also hinted by equation (2).
The temperature profile in figure 10(b) shows that the largest temperature values are reached in front of the inboard wall. This is expected, since there the He has a low speed, but it is being pressurized, see the pressure profile in figure 8(b), by the expanding jet. The pressure and temperature behavior are qualitatively similar, since the jet expansion leads to an energy conversion from internal (high pressure, high temperature) close to the inlet region to kinetic (high speed) close to the jet front. This is also confirmed combining the isentropic relations reported in equations (1) and (2), leading to  Temperatures in the range of 1200-1500 K as those shown in figure 10(a) can be reached, but they should not be a problem for the integrity of the FW. This is because it features a 2mm thick tungsten armor whose melting point is above 3000 K and could show recristallization only above 1500 K [25].

3D CFD vs. 0D system-level model
The results of the 3D model are then compared with those obtained with a 0D model developed with the GETTHEM code, which was already applied to the analysis of this kind of transients [6,8]. The GETTHEM model adopted in this work is a modified version of that adopted in [8], fixing the mass flow rate at the inlet of the VV to the value computed by the CFD model, in order to guarantee that the total He mass in the system at the end of the transient is the same in both models; a sketch of the model is reported in figure 11. The GETTHEM simulation is stopped at the same condition of the CFD one, i.e. when the average pressure in the VV reaches 1.5 bar.
The time evolution of the pressure in the VV as computed by the GETTHEM and CFD models is reported in figure 12, together with the time evolution of the average pressure on the BD surface as computed by the CFD model.
The main quantity of interest here is the time needed to reach the BD rupture threshold: as the GETTHEM model computes a single pressure value for the entire VV, the rupture would happen when the average pressure value reaches the 1.5 bar threshold, i.e. after 49 ms. Conversely, the 3D model can predict the local pressure evolution (i.e. averaged only on the BD surface). In this case, the threshold is first reached after 19 ms, i.e. when the first pressure wave reaches the BD (see section 4.2.1 above), and the pressure stays above the threshold for 4.4 ms, see figure 12 (green curve). A time span of 4 ms is a lower bound of the time needed to cause the rupture of most commercially available BDs [26]. In the case of rupture, the remaining part of the computed transient clearly has little physical meaning, as the fluid would start flowing outside of the VV towards the suppression tank, causing the (average) pressure in the VV (magenta curve in figure 12) at least to rise with a smaller slope (and eventually to reduce). To analyze this phase of the transient in detail, a model of the suppression tank should be connected to the 3D model, which is however beyond the scope of this work. In the current absence of details about the BD which will be employed in the EU DEMO, it can also be postulated that the BD does not rupture before the pressure goes back below 1.5 bar (i.e. before 23.4 ms), so it is still interesting to observe the CFD prediction of the pressure on the BD surface (green curve in figure  12) after this time. In this case, the pressure on the BD surface overcomes the BD rupture threshold again at about the same time predicted by the GETTHEM model (blue curve in figure 12).
Based on the above, it can then be concluded that: • in case the BD does not rupture after the first burst in local pressure, i.e. before 23 ms in figure 12, the GETTHEM prediction of the average pressurization in the VV is supported by/consistent with the CFD one; indeed, in terms of the average pressure in the VV, the difference between the two models is marginal since, after the first few milliseconds of the transient, the pressurization rate is roughly the same for both; • in case the BD ruptures as a consequence of the first burst in local pressure, the GETTHEM prediction of the pressurization is conservative, as the pressure inside the VV will increase after the rupture at a lower rate than predicted by GETTHEM.
As a final remark, it is important to highlight that the local pressure level may be much higher than the average one (as highlighted in section 4.2.1 above), so the loading of other critical components (such as the gyrotron  diamond windows) should be analyzed with detailed (3D) models.

Conclusions and perspective
A CFD model to simulate an in-vessel LOCA from a heliumcooled BB for the EU DEMO has been developed.
The model is 3D and transient, and is able to capture the hypersonic flow expected to develop inside the plasma chamber in view of the extremely high pressure ratio.
The model, as well as the chosen mesh adaptation strategy, have been first validated against experimental data from the literature in a simpler 2D steady-state hypersonic flow scenario, and then applied to the real EU DEMO in-VV LOCA scenario. The model is able to resolve the strong velocity, pressure and temperature gradients developing in the shock fronts, and allows evaluating the local pressure evolution at different locations on the VV walls.
The model results have been compared with those obtained with a 0D code, which is the approach normally adopted for accidental transient analyses. The 0D code is able to capture with acceptable accuracy the average pressure evolution; however, it necessarily loses important information on local peak values, which are of interest for the integrity of the containment barriers. In addition, the 0D code prediction of the intervention time of the Burst Disk is shown to be equal to or larger than that obtained with the 3D code.
In perspective, some simplifications introduced in the model shall be removed, namely concerning the domain and in particular adding a detailed model of the inlet region. Moreover, the analysis of the first instants of the transient shall be carried out with a different physical model, able to resolve the He distribution before the pressure level of 0.1 bar is reached; the two models could then be combined to perform analyses of the entire transient and have accurate estimates of the accident timing. In addition, the CFD model shall be connected with a system-level model of the Primary Heat Transfer System where the break occurs, to provide the CFD with boundary conditions more representative of the actual problem, in particular, to continue the transient beyond the intervention of the Burst Disk. Finally, a similar analysis for the release from a water-cooled BB is planned, where, however, additional complications will be present, due to flashing and two-phase flow.

Appendix A. Description of the phenomenon
The phenomenon to be studied in the scenario presented in section 2 is the evolution of an underexpanded jet. This situation is typically found at the exit of a convergent-divergent nozzle when the exit pressure is larger than the ambient pressure [27]. Depending on the pressure ratio, i.e. the ratio between the high-pressure environment and the or low-pressure one (η = p high /p amb ), the flow pattern changes significantly. It is useful to start describing briefly the phenomena occurring in the so-called moderately underexpanded jets (η <~3) [24].
The fluid at the exit plane is supersonic and the only way a supersonic flow has to expand is through the so-called Prandtl-Meyer expansion fan [28], which is a series of isentropic expansion waves. After going through the expansion fan, the flow is at the correct pressure, i.e. p amb ; however, the other effect associated with the expansion is also a change in direction of the jet velocity, which is deflected away from the exit axis, see figure 13(a). At the same time, the expansion waves emanated from the edge of the exit are reflected at the 'nozzle' exit centerline and the flow (at ambient pressure, but with the 'wrong' direction) goes again through an expansion fan, recovering the initial direction i.e. parallel to the nozzle centerline but expanding further, i.e. below p amb . The reflected expansion fan is now reflected at the jet boundary, i.e. the boundary at which p amb is present, thus converting it into compression waves, see figure 13(a), which may coalesce into an oblique shock [29] (also called intercepting shock in this context [24]), i.e. a shock that is not normal to the flow direction. Regardless of the shock type, a shock always compresses the flow. Indeed, the oblique shock compresses again the fluid at the correct pressure, though deflecting it again away from the centerline. The oblique shock then reflects at the centerline and the flow, going through the reflected oblique shock is again deflected in the correct direction, however being compressed again (above p amb ). At this point, the flow conditions are equivalent to those at the nozzle exit: a flow parallel to the nozzle exit and at a pressure higher than p amb , which then must expand through expansion fan and the pattern repeats in the same way. In a real fluid, where viscous effects are involved, there are dissipations, thus this pattern (also called 'shock cell') cannot replicate indefinitely.
For highly underexpanded jets (2 < η < 4), which is also the case simulated in appendix C, the reflection of the oblique (or intercepting) shock does not happen at the centerline, but it intersects a normal shock (also called Mach disk) and the oblique shock is then reflected [24]. This phenomenon is also referred to as Mach reflection [29,30], i.e. when the deflection angle of the oblique shock is large, a regular reflection is not possible and a normal shock is formed, which intersects the oblique shock. The point at which the intercepting shock and the Mach disk intersect, and from which the reflected oblique shock is emanated, is called triple point. From the triple point, a slipstream is emanated, which is a shear layer that separates the (subsonic) flow downstream of the Mach disk from the (supersonic) flow upstream of the reflected shock.
In case of very highly underexpanded jet (η > 4), which is the most interesting case for this work, the number of shock cells diminishes as the pressure ratio increases, leading to the limit condition in which only a single large shock cell is present and in which the Mach disk becomes curved, and thus cannot be strictly considered a normal shock anymore, see figure 13(b). The unique cell is also called 'barrel' and, therefore, the Mach disk is also named 'barrel shock'. This is the case in this work, i.e. η ≫ 4, thus a single shock cell is expected.
In addition, in the case examined here, an obstacle is present in front of the jet, i.e. the inboard portion of the VV. In the case of free jet impingement, three different regions can be identified [31]: a first region where the effects of the presence of the obstacle are negligible, so called free jet regime; a second region (called impingement regime) where the interaction of the jet with the obstacle causes a change in the flow direction and a third region (wall jet regime) where the flow becomes essentially parallel to the obstacle surface. The case of a supersonic jet impingement is characterized by the same features [32].
So far, the discussion has been focused on stationary phenomena (underexpanded jets and jet impingement). However, the scenario to be analyzed is focused on the initial portion of a LOCA, therefore it will include the initial development of the jet as soon as the rupture takes place, i.e. the ambient at high pressure comes in contact with that at low pressure. For this reason, it can be useful to briefly describe the expected relevant phenomena during the initial stages of the jet evolution, which have been investigated both experimentally and numerically in the literature [23,33].
Different stages can be distinguished in the jet evolution, see also section 4.1: first, a so-called leading shock starts propagating with a planar front in front of the rupture. The presence of the walls close to the rupture induces a diffraction of the leading shock (for this reason it is also called 'diffracting shock') [23]. The solid boundary of the rupture leads also to the formation of an oblique shock emanated from the solid boundary, which will soon become the intercepting shock. It intersects a forming shock which propagates behind the leading shock and, at the early stages, does not reach the centerline. This shock will propagate and grow, forming the Mach disk. From the intersection points, slip lines starts being emanated, as already discussed in the steady-state features of the underexpanded jet.
At the intersection between the intercepting shock and Mach disk, so-called vortex rings are formed. These are regions of strong flow instabilities, turbulence and vorticity. These vortices are originated by the shear between the subsonic surrounding environment and the supersonic jet, which drives the so-called Kelvin-Helmholtz instabilities and, entraining the surrounding fluid, ends up in the formation of vortices.
In order to solve for these local phenomena, a numerical model able to capture these small scale effects is needed, as for example large eddy simulations [34] or even direct numerical simulation. However, the larger scale, main phenomena have been shown to be captured well even solving the Navier-Stokes equations in the inviscid limit, i.e. Euler equations. This is because the jet region main features are well approximated  with an isentropic expansion (still being multi-dimensional and transient, thus asking for a numerical solution of the problem). Neglecting the small-scale perturbations allows also to treat the problem as axisymmetric (in case of axisymmetric geometry) [33].
The jet evolution continues and the leading shock starts dissipating, due to its continuous expansion in the quiescent medium and it asymptotically travels at sonic speed [33]. On the other hand, the Mach disk becomes stronger even if it is also expanding, because the compression waves travelling from the aperture accumulate on the Mach disk itself.
At this point, the further evolution of the jet is not of interest here because of the presence of the obstacles that induces deformation in the jet and it limits the jet expansion. In addition, the expansion is attenuated by the growing average pressure in the VV volume, which reduces the pressure ratio, thus decreasing the strength (and the extension) of the jet itself.

Appendix B. CFD solvers and mesh refinement algorithm
In this appendix, the details behind the choice or modification of the numerical solvers as well as how the mesh refinement algorithm has been built are described.

Solvers
The most stable transient solver available is chosen ('implicit unsteady'). Furthermore, it has been equipped with a time step adaptation scheme. This is necessary because the most critical time range of the simulation is the start-up, where the cell faces of the inlet are at the upstream pressure, velocity and temperature, while the centroid of the adjacent cells in the domain are at the VV initial conditions. Therefore, a very small time step, such as 10 −8 s, is used at the beginning of the transient. However, such time step is not necessary throughout the transient. In particular, a parameter which qualifies the current time step is the Courant number Co = v ∆t/ ∆x, where v is the local speed, ∆t is the time step and ∆x is the local grid size (also known as Courant-Friedrichs-Lax, CFL, number). The strategy for the time step adaptation is therefore based on Co. In particular, if both the mean and maximum Co in the domain are below 0.5 and 5, respectively, the time step is increased by a factor 1.1, otherwise it is halved. This allows an automatic selection of the most suitable time step throughout the transient. In average, after the start-up phase, the resulting time step is around 2 µs.
The space discretization is carried out using finite volumes with an hybrid MUSCL 3rd order/Central Difference scheme [35], as recommended in [36]. Based on the Mach number, in the shock region, the scheme is limited based on Weighted Essentially Non Oscillatory (WENO) [37] principles, which is commonly employed as shock-capturing method and lowers the overall scheme order from 3rd to 2nd in the shock region.
Concerning the solver settings, the default options are not suitable for hypersonic flows [36]. In particular, the cycle employed in the algebraic multi-grid (AMG) solver is by default a V-type. However, an F-type cycle is better even if it is more expensive since it requires more internal cycles among the different grid levels adopted in the solver. Furthermore, it is suggested not to employ acceleration method, such as the default Bi-Conjugate Gradient, leading to a slower convergence but ensuring better robustness.   The pseudo-time step employed in the inner iterations of each time step is governed by the solver CFL number specified by the user in the coupled implicit solver settings. Note that this CFL number is not related to the Courant number adopted for the time step adaptation. The solver CFL number is tuned at the beginning of the simulation, starting from a low value (0.25) and ramping it linearly to 50 in 500 iterations. This slows down the convergence of the first time steps, but makes the internal iterations more stable.

Mesh refinement algorithm
In order to describe the mesh refinement algorithm, various quantities need to be introduced.
First, the representative dimension of a given polyhedral cell δ c in the current mesh is defined according to equation (5), see [17].
where V c is the volume of a given cell. The maximum and minimum target cell size, δ M and δ m , respectively, are defined according to equations (6) and (7), respectively. Note that δ M is the minimum between 0.5 m and three times the current size of the cell, i.e. δ M tends to coarsen the cell dimension, up to the limit of 0.5 m. The key quantity to capture the front propagation is the Mach number variation ∆Ma, see equation (8). The choice of the Mach number stems from the need of including the gradient of both velocity and temperature in identifying the locations in the domain that need mesh refinement.
A possible strategy to smoothly go from δ M to δ m and vice versa is to link the target dimension of the cell to the parameter on which the refinement depends, i.e. depending on ∆Ma, the cells are refined or coarsened within the range limited by δ m and δ M , according to equation (9), where δ c,ref is the cell size after remeshing.
Finally, the target cell size δ c,target , see equation (10), is chosen in the following way: in case in a given cell the quantity ∆Ma is above 0.1, the cell needs to be refined; the target cell size will then be the maximum between δ c,ref and δ m , in order to limit the minimum cell size in the mesh. In case ∆Ma is lower than 0.1, the target cell size is the minimum between δ c,ref and δ M , which limit the maximum size of the target cell achievable in the mesh.
The resulting mesh, at t = 53 ms, is reported in figure 3 where it is showed the capability of the present algorithm of following the propagating shock front and suitably refine the mesh there as well as coarsen the mesh where the front had already passed and no relevant gradients are present anymore. This allows keeping the total number of cells always below 10 million. In particular, the first~5 ms are crucial and heavier from the computational point of view because the jet is expanding, therefore, as the front grows, more and more cells are needed to cover the entire front. However, once the shock front has hit the wall in front of the inlet, it stops expanding, reaching a steady dimension, therefore the cell count does not grow anymore. In addition, as the average volume pressure increases, the gradients at the front are less and less severe, requiring fewer cells in the refined region. Thus, the cell count falls below 2.5 million after~30 ms, allowing a faster computation.

Appendix C. Model verification and validation in a 2D steady-state case
The model presented in section 3.2 is verified here against another published model [38] and validated against experimental data also published [39,40].
The selected case is a 2D steady state under-expanded jet with a pressure ratio~3. Even though the case of interest for us in this paper is clearly 3D, transient and with a larger pressure ratio, the validation on a simpler case should allow testing all the relevant features (model, solvers, mesh strategy) to be then employed in the more complex case.

Setup
The boundary conditions are the following, see figure 14: 1. Inlet: fixed supersonic and total pressure, total temperature 2. Side outlet: reference pressure and temperature gradient set to zero 3. Top outlet: fixed ambient pressure and temperature 4. Bottom boundary: symmetry axis 5. Wall: no-slip condition Initial condition: zero velocity, temperature equal to 300 K and pressure equal to 1 bar.

Model
The list of features of the model adopted in the simple case is the following: • 2D axisymmetric • Steady • Turbulent (k-ω SST [21]) • Ideal gas (air) • Coupled flow • Coupled energy The discussion on the models and on the mesh adaptivity strategy can be found in section 3.2.1 and 3.3, respectively, as well as in appendix B.
Since the simulation considered here is steady state, the update of the mesh is performed after 10 iterations rather than 10 time steps.

Results
The results computed for an inlet pressure equal to 3.25 bar are compared in figure 15. The comparison between the velocity and Mach number field computed in [38] and with STAR-CCM+ with the setup described above shows qualitatively the same behavior. In particular, the feature expected from an under-expanded jet is retrieved, see appendix A. For example, the location of the Mach disk is the same for both the reference and the setup and software adopted here.
Concerning the quantitative comparison of the results, a difference of about 10% in both the velocity and Mach maximum values is found. Note, however, that, for example, the meshing strategy is quite different. Therefore, a further step in the validation is undertaken and described below.
Experimental results for this setup are available [39,40], concerning, for example, the pressure profile across the first Mach disk. The profile of the static pressure computed with STAR-CCM+ shows a good agreement with the experimental data and with the simulation results of [38], see figure 16(a). The adaptive mesh adopted, see figure 17, in the current setup allows having a better description of the sharp gradient in correspondence of the Mach disk with respect to the simulation results of [38] where a static mesh is employed, since the discontinuity corresponding to the shock is better approximated, i.e. with less numerical diffusive effects.
A parametric scan varying the inlet pressure is performed. The quantity monitored in the parametric study is the position of the Mach disk with respect to the inlet, see figure 16(b), as a function of the pressure ratio, defined as the ratio of the exit (high) pressure and the ambient (low) pressure. The agreement with the experimental data available is very good, i.e. with less than 5% error, and this confirms that the models adopted as well as the adaptive mesh strategy are sound and their application to the transient analysis of an under-expanded jet should be reliable.