MULTI-PHYSICS METHODOLOGY FOR PHASE CHANGE DUE TO RAPIDLY DEPRESSURISED TWO-PHASE FLOWS

We propose a multi-physics methodology to predict the phase change of a rapidly depressurised ﬂuid usually found in ﬁre suppression systems. The system consists in a pressurised ﬁre extinguishing agent (Novec-1230) that is injected into a compartment through a nozzle that produces the atomisation of the agent. A multi-scale approach is proposed where the detail of the atomisation process is captured through a high ﬁdelity multiphase phase ﬁeld method. Once the ﬁre extinguishing agent is atomised, a particle tracking method is used to model the phase change and the spreading of the agent into the compartment. An experimental campaign has been carried out to validate the numerical methodology. The experimental results provide the axial penetration and the spreading angle of the spray with time. Dif-fuse back-light illumination and Schlieren techniques are used to permit the visualisation of Novec-1230 in liquid and vapor state.


INTRODUCTION
Aeronautical on-board fire suppression systems are based on the interruption of the propagation of chain reactions typically found in aeronautical fuels. Fire suppression systems historically used hydrofluorocarbons (HFCs) as a fire protection fluid, the most well-known and used one being Halon (R13B1, CFBr3). However due to their high global warming potentials (GWPs), the global regulatory phase-down under the Montreal Protocol and the availability of proven, more sustainable alternatives, industry has been pushed towards more environmentally friendly alternatives.
Among these alternatives, we can find Novec-1230 Fire Protection Fluid [1]. The main difference between Halon and Novec-1230 is, in addition to the lower global warming potential of the latest, its higher molecular weight. As a result, the boiling point at ambient pressure of the latter is much higher (-57.8 • C compared to 49.2 • C at ambient pressure) which results in a more difficult vaporisation and dispersion. This difference is one of the reason to prevent the direct implementation of Novec-1230 as a new suppression agent, which possibly requires a redesign of the complete fire suppression system. In this work, a numerical multi-scale approach is proposed to predict the atomisation, the phase change and the spreading of the agent into the compartment. This methodology is applied to a generic fire suppression system and compared to experiments to show its validity.
The modelisation of multiphase flows is not a canonical problem, therefore different approaches can be found in the literature. Two different approximations will be considered here: Euler-Euler (EE) models and Euler-Lagrange (EL) models.
In EE models the two phases are treated mathematically as interpenetrating continua. The concept of volume fraction (or void fraction, normally denoted by φ) as the fraction of a particular infinitesimal control volume which is occupied by each phase (minimum two), is introduced. Volume fraction is supposed to vary continuously in space and time. Conservation equations for each phase (Navier-Stokes equations) together with physical constitutive relations (usually from experimental correlations) and a partial differential equation for the evolution of the volume fraction are posed. The solution of this set of equations describe the evolution of the flow. Different levels of complexity can be found to describe the EE models. Among the most simple is the Volume Of Fluid (VOF) model [2], which defines a single set of momentum equations shared by all phases, while the volume fraction is tracked throughout the domain following a convective equation. A similar approach is followed by phase-field methods [3], where a phase-field parameter which identifies each phase is defined. The evolution of this parameter is modeled through a partial different equation with more physical meaning than in VOF, see e.g. Cahn-Hilliard equation [4,5,6]. More complex models consider that each of the phases is described by its own density, velocity, temperature and pressure (see for example the Baer and Nunziato model [7,8]). These models are usually closed with mathematical relations for the interaction between the phases (e.g., drag or pressure relaxation).
On the contrary, EL methods consider that only one of the phases is continuous (Euler) while the other/s is/are dispersed (Lagrange). The continuous Navier-Stokes equations are solved for the Euler phase while the Lagrange particles are tracked (individually or as groups) [9]. Lagrange particles are used to represent objects that fall under the resolution level of the numerical grid. Liquid droplets are one of the most common examples. Two different levels of coupling between the Euler and Lagrange phases are usually considered. In the one-way coupling, the dispersed phase is transported by the continuous phase but the latter is not affected by the former. The two-way coupling also consideres the effect of the Lagrange into the Euler phase.
Due to the limitation imposed by current computers it is not possible to address the whole atomisation and vaporisation simulation with an EE high fidelity simulation. Therefore in this work a multi-scale approach is proposed where the atomisation process (near field) is captured by a high fidelity EE simulation while the vaporisation and dispersion of the agent (far field) is simulated with a EL model. For the near field, the proposed approach studies the primary atomisation of a liquid jet injected into air using a high fidelity Large Eddy Simulation (LES) EE model. It is based on an entropy stable conservative level set scheme stabilised with an entropy stable method [10] for the regions where the interface can be resolved, while it includes a subgrid contribution when the slip-velocity of the two-phase flow cannot be resolved with the given filter size. For the far field, an EL approximation is solved which includes a Eulerian description for the gas while the liquid is assumed to be dispersed and tracked as individual droplets. This model is only valid to simulate the flow outside the nozzle (after the fire extinction agent has been atomised). In this work, a two-way coupled model based on previous work in the simulation of fire suppression systems [11,12,13] is used to solve the droplets vaporisation.
Both models are coupled through a boundary condition for the EL model. Results of the accurate simulation of the atomisation process obtained with the LES EE model will be used as a boundary condition for the computation of the far field domain. The LES EE simulation is performed with ALYA [14] while the EL simulation is performed with ANSYS-Fluent [15].
Additionally, an experimental campaign is carried out to validate the numerical methodology. The facilities used for the experiments are the CMT-Motores Térmicos (UPV). Previous experiments with liquid and vapour phases axial penetrations as well as spreading angle were evaluated using different experimental techniques, separately or simultaneously [16,17]. In this work, two different techniques are used: Diffuse Backlight Illumination (DBI) technique [18] is employed to visualise the liquid phase, and the Schlieren technique [19] for the vapor phase. Spray penetration and cone angle are measured at up to four visualisation windows located at different axial positions.
The remaining of this text is organised as follows: first, in Section 2, the methodology is described. This section is divided into two parts: first, the experimental methodology and later the numerical methodology. Then, in Section 3, the numerical results are compared with experimental results. Finally, in Section 4 conclusions are presented.

Experimental
A boarded fire extinguisher system may consist of only a bottle and a nozzle or on a complex arrangement of pipes and nozzles [20]. In order to simplify the test rig and the following numerical validation, a nozzle with swirling oriented in the axial direction has been used. The nozzle of Spraying Systems Co. ref. 1/4GG-SS3009, shown in Figure 4, is used. It has an orifice outlet diameter of 2mm and includes a swirler in order to add tangential component to the flow and so improve spray atomisation giving a spray angle of θ = 30 deg. The design of the test vessel is shown in Figure 2. It consists of a prism with an interior volume of 750x750x1500mm, required to avoid any interference of the vessel with the spray development. Girders and columns are made of 45x45mm Fasten basic profiles, welded to each other. Silicone and o-rings are used to seal the connections with windows and removable components (i.e. nozzle holder). The four profiles located in the middle of the vessel avoid bending (and eventually breaking) of the windows and walls. The four walls are made of stainless steel plates of 5mm thick, and in the same way as the profiles, are welded to the structure. Two large transparent windows, 750x1500mm in size and 19mm thick, are placed on two opposite sides of the prism. A metallic frame is employed to hold the window in place and ensure sealing of the vessel. One of the smaller plates (on the right hand side in Figure 2) is provided with a hole just in the center to assembly the nozzle holder, and other two connections for a pressure sensor and a thermocouple. The other smaller plate, located just in front, includes connections for another pressure sensor and the vacuum pump. This last connection is used also to re-circulate the ambient gas. Since the suppressant vapor phase is heavier than air, it may increase the discharge ambient density (up to 7 times if a saturated air and Novec-1230 gas mixture is achieved). A sketch of the complete test rig is depicted in Figure 3. A piston accumulator (Parker ref. A3ES0116L2ERC) is used to set the injection pressure (upstream pressure). A 20MPa pressurised nitrogen bottle is connected to the accumulator and, by means of a manual valve, allows increasing or decreasing the pressure of the piston. Before pressurising, the accumulator is filled up with the fire suppression agent (either water or Novec-1230). This accumulator is directly connected to the nozzle with a pipe of 10mm inner diameter, and two solenoid valves  Manometers are used in all deposits and containers to monitor the pressure, as well as thermocouples to monitor the temperature. A piezoelectric Kistler pressure sensor (SN 2144942/2013) is located at the pipe just upstream the nozzle. This pressure signal is utilised to measure the injected mass flow rate and therefore the total injected mass. The complete test rig is placed inside a climatic chamber with 30 kW of cooling/heating power capacity capable of reproducing altitude conditions (up to 3000 m) [21]. The cooling power of the chamber is used to cooldown the ambient air when it is required.
Optical instruments utilised in the measurements (described in following sections) do not allow to record the whole field of view of 750x1500mm. Therefore, up to four different positions are selected as fields of view, as depicted in Figure 4. Three of them are aligned with the axis of the nozzle, and the fourth one is moved some distance from it in order to verify the measured spray angle. As mentioned before, two different measurements techniques are used: Diffuse Backlight Illumination (DBI) and the Schlieren technique. DBI consists on the determination of the shape of the spray liquid phase based on the silhouette obtained by the obstruction of a beam of diffuse light with the jet [18]. For the setup utilised in this work, the light beam is produced by a high power continuous light source and goes through a diffuser which homogenises and smooths the background of the spray image. An example of the type of images obtained with this technique is shown in Figure 8 (top).

Single-pass Schlieren
Light rays are deflected when trespassing through a medium with density changes. Schlieren imaging technique makes use of this phenomenon and allows detecting the spray contour by filtering or discarding some of the deflected light [19]. A single-pass setup was used in this work. A white-light source was placed at the focal length of a collimating lens. After this lens, the now parallel light beams pass through the testing region. The deviated beam is collected by another lens, being a diaphragm located at its focal distance, just before the high-speed camera [18]. As before, an example of the type of images obtained with this technique is shown in Figure 8 (bottom). Vapor phase is barely visible to the eye in the first field of view, but for further distances it covers the whole field of visualisation.

Numerical
Regarding the numerical methodology, two different approaches are used depending on the zone of interest. As mentioned before, two zones are considered: near field and far field. For the near field, the proposed approach studies the primary atomisation of a liquid jet injected into air using a high fidelity Large Eddy Simulation (LES) EE model. For the far field, an EL approximation is solved which includes a Eulerian description for the gas while the liquid is assumed to be dispersed and tracked as individual droplets. Both models are coupled through boundary condition for the EL model. Results of the accurate simulation of the atomisation process obtained with the LES EE model are used as a boundary condition for the computation of the far field domain.

Near field: Euler-Euler model
The near field corresponds to the interior of the nozzle and the primary atomisation zone. No phase change is considered here. In this work, the proposed Euler-Euler (EE) approach to study primary atomisation for air and Novec-1230 using LES is based on an entropy stable conservative level set method. This model is implemented in Alya [14] solver. The modeling approach includes an extension of the conservative level set equation stabilised with an entropy stable method [10] for the regions where the interface can be resolved, while it includes a subgrid contribution when the slip-velocity of the two-phase flow can not be resolved with the given filter size [22]. The system of equations is given by the conservation of momentum, liquid volume fraction φ, and liquid/gas interface Σ in the incompressible limit, with ∇ · u = 0. The system is closed by the equation of state for a liquid/gas mixture given by ρ = φρ l + (1 − φ) ρ g , where the density of the gas ρ g and liquid ρ l phases can be obtained by equations of state for gas and liquid respectively, but are assumed to be constant on this study. The liquid/gas interface density Σ represents the quantity of liquid/gas interface per unit of volume and it is used to provide the characteristics and sizes of the liquid droplets in the dense part of the spray. The conservation equations will be presented below using spatially filtered quantities for LES. A decomposition of Σ following Chesnel at al. (2011) [22] is used and the subgrid surface density Σ that evolves similar to Σ is solved instead. This separation is based on the existence of a minimum surface density Σ min due to the presence of liquid and a subgrid surface density Σ that evolves similar to Σ = Σ min + Σ , where Σ min is obtained from empirical correlations based on DNS [22] and given by with ∆ being a characteristic length scale (filter size in standard LES with implicit filtering) and α = 2.4. Therefore, the transport equation for the liquid-gas interface can be expressed in terms of Σ rather than Σ.
The complete set of governing equations is presented here using filtered quantities and reads: where ν is the minimum between the upwind viscosity and the stable entropy. The turbulent flux τ φ and the slip velocity Σ ( u − u Γ ) are both modelled using a gradient assumption closure [22,23]. The source terṁ Σ int represents the production/destruction of surface density by the mean shear, turbulence and liquid structure interactions [22] and is modeled bẏ The characteristic time scale τ b can be related to a turbulent time scale as in the dense zone of the spray most of the liquid structures break-up and coalescence by turbulence [22] and the same approach will be used here.
The equilibrium value of surface density Σ crit is obtained from a critical Weber number expressed in terms of balance of turbulent kinetic energy and takes the form of The spatial discretisation is based on a low dissipation finite element method with an explicit algorithm of the fractional step to solve the velocity/pressure coupling [24] and an entropy stable conservative level set [10] for the volume fraction equation. The main difference with Guermond work [10] is the use of a cell-based viscosity instead of an edge-based version. The temporal discretisation is based on a third order energy preserving Runge-Kutta scheme.

Far field: Euler-Lagrange model
The far field focuses into the vaporisation process solved through a Euler-Lagrange model (EL). In the EL approximation the ambient gas is solved following an Eulerian description while the liquid, assumed to be dispersed, is tracked as individual droplets. This model is only valid to simulate the flow outside the atomiser (after the fire extinction agent has been atomised), therefore a valid inflow boundary condition should be provided in ANSYS-Fluent. The explanation will be decomposed into two parts. First we will explain the modifications required in the Euler solver to take into account the Lagrange (drops) part. Second we will review the equations that model the evolution of each of the droplets.

Euler solver
The Euler equations are described with the extra terms coupled with the Lagrange equations.
Momentum equation: where now ρ refers to the gas density and u its velocity, p is the pressure (obtained from the incompressibility condition ∇ · u = 0), f b represents the momentum transferred from particles to the gas (see Eq. 13), µ is the gas viscosity and g is the gravitational acceleration.
Energy equation: where energy E is based on the specific heat and the temperature T , k is the thermal conductivity and q b represents the energy transferred from the particles to the gas (see Eq. 18).

Species transport equation:
An additional transport equation for mass fractions is included to track the evolution of the fire extinction agent (vapor) in the gas phase, where ρ is the density of the air, Y 1 is the mass fraction of fire extinction agent vapour in the gas, D lg is the diffusion coefficient of fire extinction agent vapour in the gas andṁ evap is the source term due to evaporation (see Eq. 17).

Lagrange solver
The Lagrange solver (tracking of the particles) is described here. It is important to notice that, to address computational efficiency, not every droplet is tracked as the computational effort would be too high. Instead, a smaller number of "superdroplets" or parcels are tracked, where each parcel represents many individual real droplets with the same diameter and thermophysical properties [11].
Momentum equation: where u p is the particle velocity, g the gravity, C d its drag coefficient (function of its Reynolds number), A p,c the particle cross-sectional area and m p the particle mass. Gas velocity, u, and gas density, ρ, are obtained from the Euler part of the solver. The momentum equation solved by ANSYS-Fluent takes into account buoyancy effects, which are negligible for liquid droplets moving through the air. Being the velocity of the particles known, the particle position, x p , is determined from the equation, Based on a sphere, the drag coefficient is computed as: where µ(T ) is the dynamic viscosity of air at temperature T and r p the particle radius.
By summing the forces transferred from each particle in a grid cell and dividing by the cell volume, V , the momentum transferred from particles to the gas (see Eq. 7) is obtained: Energy equation: Over the course of a time step of the gas phase solver, the droplets in a given grid cell evaporate to form part of the gas phase. The evaporation rate is a function of the liquid equilibrium vapour mass fraction, Y 2 , the local gas phase vapour mass fraction, Y 1 , the (assumed uniform) droplet temperature T p and the local gas temperature, T . The mass and energy transfer can be described by the following set of equations: where B m is the Spalding mass number given by: Here, m p is the mass of the droplet, A p,s the surface area of the liquid droplet, k m the mass transfer coefficient (to be discussed below), ρ the gas density, c p the liquid specific heat, k h the heat transfer coefficient between the droplet and the gas and h v the latent heat of vaporisation of the liquid. Again, by summing for all the particles in a cell we get the source term for the advection equation, see Eq.9, which is:ṁ where V is the cell volume. As far as the source term for the energy equation, see Eq.8, it reads: where h l is the liquid specific enthalpy, and h v is the latent heat of vaporisation of the liquid.
In order to evaluate Eq. 14, the vapour mass fraction of the extinction agent in the gas Y 1 is obtained from Eq. 9 while the liquid equilibrium vapour mass fraction is obtained by assuming chemical equilibrium between phases defined by the Clausius-Clapeyron equation: where X 2 is the equilibrium vapour volume fraction, W is the molecular weight of the gaseous species (fire extinction agent), W a is the molecular weight of air, R is the universal gas constant and T b is the boiling temperature of the fire extinction agent at standard atmospheric pressure.
Mass and heat transfer coefficients between liquid and gas are described with analogous empirical correlations. The mass transfer coefficient, k m is described by: being Sh the Sherwood number, D lg the binary diffusion coefficient between the fire extinction agent vapour and the surrounding gas (air), L is a length scale (equal to the droplet diameter), Re D is the Reynolds number of the droplet (based on the droplet diameter, D, and the relative air-droplet velocity) and Sc is the Schmidt number ν/D lg (ν = µ/ρ is the kinematic viscosity). An analogous relationship exists for the heat transfer coefficient: where Nu is the Nusselt number, k is the thermal conductivity of the gas and Pr is the Prandtl number. These correlations are used by ANSYS-Fluent. It should be noticed that the exchange of mass and energy between liquid droplets and the surrounding gases is computed droplet by droplet (or parcel by parcel). After the temperature of each droplet is computed, the appropriate amount of vaporised liquid is added to the given mesh cell, and the cell gas temperature is reduced slightly based on the energy lost to the droplet. The integration methods of these equations may create numerical instabilities or inaccuracies [13].

Coupling between Euler-Euler solver and Euler-Lagrange solver
As previously mentioned, none of the previous approaches is well suited to solve the complete physical phenomena (discharge, atomisation and vaporisation of the fire suppression agent). Therefore in this work a multi-scale approach is proposed where the atomisation process (near field) is captured by the high fidelity EE simulation while the vaporisation and dispersion of the agent (far field) is simulated with the EL model. Figure 5 represents the modelisation of the problem with the two approaches. On the one hand, the blue box represents the two-phase LES solver domain, where the atomisation process takes place. On the other hand, the red dashed box represents the Euler-Lagrange solver domain where the liquid can be considered dispersed and the phase change occurs. Between the two domains, the dashed black line represents the solid cone injection for the Euler-Lagrange solver. Following [25] both models are coupled through boundary condition for the EL model. Results of the accurate simulation of the atomisation process obtained with the LES EE model will be used as a boundary condition for the computation of the far field domain.

Figure 5: Sketch of the coupling between EE and EL models
This methodology proposes solving a two-phase Euler problem where the mesh resolution can capture the details of the flow and transition to a Euler-Lagrange model where the liquid volume fraction (that represents the droplets) falls below the mesh resolution limit. Then processing the results and using that as a boundary condition for the Euler-Lagrange simulations allows to adress the vaporisation of the droplets. The main parameter controlling the transition from the two-phase Euler model to the Euler-Lagrange model is the liquid volume fraction,φ, that represents the volume occupied by the liquid with respect to the cell volume. The transition is complete when the liquid volume fraction becomes lower than 0.1. The transition zone is composed of the computational cells that form the border with the dense zone (i.e. the zone where the liquid volume fraction is greater than 0.1). The diameter of the droplets is calculated as the Sauter mean diameter, given by:d The number of droplets per generated parcel is obtained from mass conservation: where V cell is the volume of each transitional cell. Finally, the Euler-Lagrange simulation uses an extended domain and injects the particles with the droplet diameter and velocity values extracted from the transition zone with a solid cone model.

RESULTS
The reference application case corresponds to the atomisation of a rounded jet of Novec-1230 issuing from a single injector nozzle. The experimental facility is located at CMT and the operating conditions of the experiment are detailed in Table 1.
Time  The remaining of this section includes experimental and numerical results.

Experimental
This section is divided in two parts. First, the transient state the spray penetration and the mass flow during the time evolution is analysed. Later, in the steady state the spray angle and atomisation is studied. The results are shown with the both techniques mentioned before: DBI and Schlieren.
3.1.1. Transient state: Spray penetration Figure 6 depicts the spray penetration (a) and the mass flow rate (b), corresponding to the field of view #1 (see Figure 3). The shape of the curve in Figure 6a changes from a parabolic trend to a straight line because in the first moments of the injection process the spray accelerates as it is injected. Afterwards, the spray reaches the steady state in which the pressure stabilises and the injection mass flow rate is constant as shown Figure 6b . This effect can be seen in the sequence of Figure 7. From the start of injection until approximately 12ms there is a narrow liquid core corresponding o the fluid being accelerated in the injection system pipes, pushed out through the nozzle without the momentum provided by the swirler. After that, the swirl becomes effective, providing the spray with a significant radial velocity and hence promoting the atomisation process. Then, the steady state (Figure 7e) is reached. Unfortunately, the transient narrow liquid core always reaches the end of the field of view before the spray widens, for all injection conditions tested. This is an issue because the stabilisation of the injection occurs almost at the end of the visualisation window, explaining why in the penetration curves shown next there is a fast change in the slope at the end.

Steady state: Spray atomization
At the steady state, the spray atomisation can be analysed to obtain the spray angle. It is clearly observed in Figure 8 how the liquid droplets reach the end of the third field of view, a distance of about 1 m from the orifice outlet. There is a significant difference (15%) between the spray angle values measured by DBI (θ = 25deg) and Schlieren techniques (θ = 30deg) for Novec-1230, which means that liquid droplets are evaporating as they are injected.

Comparison between experimental and numerical results
In this section the numerical results from both models EE and EL are compared with experimental results.

EE model
The unstructured mesh employed to describe the injector system is composed by 73 million tetrahedrons that include the internal flow in the nozzle and the discharged atmosphere. Three levels of refinement are considered to characterise the internal flow and the near field after the expansion, the jet penetration up to 15D and the surrounding air with 0.1, 0.3 and 1.5 mm respectively (see Figure 9).  During stable operation, the liquid jet core penetrates with 30deg spray angle for about ten diameters and the jet breaks up into filaments and large droplets for an extension larger than the refinement region of 15D as shown in Figure 11. We obtain similar results in terms of spray penetration and angle as in experimental results (see views #1 in Figure 8). The dynamics of the spray can be visualised by an isocontour of liquid volume fraction φ=0.5 as shown in Figure 12 where droplet formation also occurs near the nozzle exit and in the presence of the liquid core. Figure 12: Iso-contour of liquid volume fraction φ=0.5 From the atomisation point of view, it is around 12D to 18D where the primary break up takes place and the liquid core is broken and where the extraction of liquid to compute droplets is obtained. Figure 13 shows the volumetric region from the continuous flow field where the droplets are extracted and used for the Lagrange model based on the Sauter mean diameter in the following section. Figure 13: Iso-surface of liquid of volume fraction coloured by the Sauter mean diameter

EL model
The far field simulation focuses into the vaporisation process solved through a Euler-Lagrange model (EL). For this simulation, the computational domain is a cylinder with a diameter of 2m and a height of 4m. An unstructured grid with 30,000 elements is generated with ANSYS-Meshing module. A refined area is considered to capture correctly the Eulerian solution (see Figure 14). The air is solved through the compressible Navier-Stokes equations. While the liquid is tracked with the Discrete Particle Model (DPM). This model is only valid to simulate the flow outside the atomiser (after the fire extinction agent has been atomised), therefore a valid inflow boundary condition should be provided in ANSYS-Fluent. The liquid particles are injected with  Figure 15 with DBI experimental results at steady state, where the evaporated particles are not visible. The spray angle of 30deg is well-predicted (see view #1), this angle is loosed at certain distance where the particle are evaporated (see view #2 and #4). Additionally, the numerical computations gives information about the distance where the Novec-1230 is completely evaporated. The vaporization distance is 2m from the outlet nozzle. In contrast with the DBI comparison, Figure 16 shows the Novec-1230 concentration in Fluent compared with Schlieren experimental results, where the density changes is captured. This comparison shows how the Novec-1230 evaporated maintain the 30deg along the domain, which it is captured with both experimental and numerical.

CONCLUSIONS
In this work, we have presented a multiphase methodology able to compute, at an affordable cost, the in-