Investigating 3-D Effects on Flashing Cryogenic Jets with Highly Resolved LES

For the development of upper stage rocket engines with laser ignition, the transition of oxidizer and fuel from the pure cryogenic liquid streams to an ignitable mixture needs to be better understood. Due to the near vacuum conditions that are present at high altitudes and in space, the injected fuel rapidly atomizes in a so-called flash boiling process. To investigate the behavior of flashing cryogenic jets under the relevant conditions, experiments of liquid nitrogen have been performed at the DLR Lampoldshausen. The experiments are accompanied by a series of computer simulations and here we use a highly resolved LES to identify 3D effects and to better interpret results from the experiments and existing 2D RANS. It is observed that the vapor generation inside the injector and the evolution of the spray in the combustion chamber differ significantly between the two simulation types due to missing 3D effects and the difference in resolution of turbulent structures. Still, the observed 3D spray dynamics suggest a suitable location for laser ignition that could be found in regions of relative low velocity and therefore expected low strain rates. Further, measured droplet velocities are compared to the velocities of notional Lagrangian particles with similar inertia as the measured droplets. Good agreement between experiments and simulations exists and strong correlation between droplet size and velocity can be demonstrated.


Introduction
The aspiration to develop environmentally friendly and lighter upper stage rocket engines leads to continuing attempts to replace hydrazine and its derivatives with more conventional fuels, such as methane, kerosene, or hydrogen.As these fuels are not hypergolic, fuel replacement strategies must be combined with a lightweight ignition technique, such as laser ignition (Manfletti and Kroupa 2013).Before ignition, the fuel is injected into near vacuum conditions present at high altitudes and in space, causing the fluid to transit into a superheated state and resulting in rapid and strong evaporation.This process is typically referred to as flash boiling or flash evaporation.A well-based understanding of this complex physical process is required to achieve a successful ignition under these extreme conditions.The modeling and simulation of flashing sprays poses a significant challenge due to the inherent multi-scale problem caused by bubble nucleation in the range of nanometers, bubble growth, breakup and coalescence with scales at a macroscale and the resulting droplets at a micrometer scale (Sher et al. 2008;Witlox et al. 2007;Rees et al. 2020).While DNS simulations of single aspects, such as bubble growth, are possible (Loureiro et al. 2020), resolving all scales within one simulation is not feasible with current computational resources.Hence, bubble nucleation, growth, breakup, and droplets are typically not resolved and have to be modeled.After the spray breakup, the droplets constitute a disperse phase which is typically simulated by Lagrangian particles, whereas the continuous phase, here the vapor, is modeled on an Eulerian grid (Yeoh and Tu 2010).However, prior to spray breakup, the liquid is the continuous phase and the vapor bubbles can be treated as Lagrangian particles.Further, a continuous and separate discrete phase cannot be defined during spray breakup.To avoid the challenge of switching the phases for the Lagrangian treatment and the difficult spray breakup region, an alternative approach where the full set of conservation equations are solved for each phase on an Eulerian grid is more feasible.Often, the additional assumption that vapor and liquid are closely coupled and move with the same velocity is applied, thus removing the need for the momentum exchange terms.Then, only one set of conservation equations for the vapor-liquid mixture has to be solved together with an additional transport equation for the volume fraction.The mixture properties are obtained through an appropriate mixing rule of the two phases.This so called one fluid model is commonly applied to solve two phase sprays and jets for flashing conditions (Schmidt et al. 2010;Guo et al. 2020;Lee et al. 2009;Lyras et al. 2017;Mohapatra et al. 2020;Saha et al. 2017;Neroorkar et al. 2011).Independent of the modeling framework, phase change and nucleation has to be modeled.A simple and early model to describe the phase change is the homogeneous equilibrium model (HEM), which assumes that the vapor and liquid are always in a thermodynamic equilibrium (Bilicki and Kestin 1990).However, this model requires the assumption of an infinite fast relaxation time towards equilibrium conditions, whereas the main aspect of flashing flows is that the liquid is in a superheated, meta-stable state, which is not in equilibrium (Sher et al. 2008).To account for the finite time required to reach equilibrium conditions, the HEM model has been expanded to include a relaxation time, resulting in the homogeneous relaxation model (HRM) (Bilicki and Kestin 1990).A description of the relaxation time, based on flashing water experiments, has been found by Downar-Zapolski et al. (1996).Despite the fit on flashing water experiments, the proposed model has shown wide applicability for various fluids and conditions with only minor modifications of the empirical coefficients (Saha et al. 2017;Lee et al. 2009;Guo et al. 2020).Further, the HRM model does not require detailed information about the interphase surface, which distinguishes it from other models such as the Hertz-Knudsen model (Persad and Ward 2016).The one fluid model, combined with the HRM model, has already been successfully applied to simulate flashing flows for various fluids and conditions, inside and outside the injector (Schmidt et al. 2010;Guo et al. 2020;Lee et al. 2009;Lyras et al. 2017;Mohapatra et al. 2020;Saha et al. 2017;Neroorkar et al. 2011).
In addition to the complex physics of the flashing process, compressibility effects may alter the flow dynamics as shocks can be observed in flashing sprays (Lamanna et al. 2014;Vieira and Simoes-Moriera 2007;Poursadegh et al. 2017).Similar to gaseous underexpanded flows, the evaporation inside the injector causes the vapor-liquid mixture to exit with a higher than ambient pressure, leading to an acceleration of the vapor into supersonic conditions, which is terminated by a shock system (Franquet et al. 2015).The presence of shocks was corroborated by experimental studies of a short converging nozzle with iso-octane as the working fluid (Vieira and Simoes-Moriera 2007).Further, shocks were observed for flashing acetone sprays by Lamanna et al. (2014) and also for propane using a gasoline direct injector (GDI) by Poursadegh et al. (2017).It is therefore essential to include the compressibility effects in the simulation of flashing sprays to capture the gas dynamics correctly, e.g., for GDI systems, the resulting shocks and their interaction is a crucial aspect to explain behaviors such as spray collapse (Gärtner et al. 2021;Guo et al. 2021;Lacey et al. 2017).
To understand the behavior of flashing cryogenic liquids, experiments with cryogenic nitrogen at conditions present at high altitudes have been performed at the DLR Lampoldshausen, providing shadowgraph images, droplet velocities, and droplet size distributions (Rees et al. 2019(Rees et al. , 2020)).A joint numerical-experimental study attempted to explain nearly motionless or slightly upstream floating dark structures that were observed in the experiment.2D RANS simulations identified a recirculation zone about 20D downstream of the injector, which seemingly explained the experimental observations.To investigate the effect of turbulence modeling as well as potentially missing 3D effects in existing 2D simulations, a 3D LES of the same setup is presented here.Previously, only shadowgraph images had been available for a qualitative comparison of the flow structures.With the more recently availability of experimental droplet measurements a quantitative comparison of the simulation and experiment is possible.This allows for a validation of our simulation models and an improved analysis of the flow behavior.This work discusses the effects and associated restrictions of 2D implementations within the injector, see Sect.4.2, as well as on the recirculation zone, see Sect.4.4.The flow behavior inside the combustion chamber and particularly the quantitative comparison of droplet velocities and simulation results are presented in Sect.4.3.In addition, one-way coupled Lagrangian particles, representing the measured droplet size distribution, are added to the simulation to investigate the correlation between droplet size and velocity in comparison to the LES results.A better understanding of the flow dynamics in flashing cryogenic jets is achieved, which may aid the identification of suitable ignition locations for combustible mixtures.These findings aid the development of novel upper stage rocket engines.

Numerical Modeling
The accurate numerical solution of the two phase flow with phase change and the faithful capture of transonic effects is challenging.Typically, density-based solvers coupled with suitable Riemann solvers are chosen to compute transonic flows and to resolve the shock front.
However, these kinds of solvers are not well suited for subsonic regimes.To overcome this problem, a pressure based, hybrid approximate Godunov-type solver has been developed by Kraposhin et al. (2015Kraposhin et al. ( , 2018)).It can be applied to a wide range of Mach numbers as well as non-ideal gas conditions.This technique has also been applied to a two phase solver (UniCFD 2021).One disadvantage is that the solver requires the phases to be compressible and that phase change is currently not included.With the aim of simulating the complete injection process of pure liquid at the injector inlet, the vapor bubble-liquid mixture in the injector, and the spray in the chamber, not only subsonic and supersonic regions are crossed but also a transition of nearly incompressible fluid to the compressible mixture occurs.This is illustrated in Fig. 1.
Due to these reasons, a compressible, pressure-based, two-phase solver with phase change has been developed in OpenFOAM, which solves the compressibility of each phase separately.This solver has already been successfully applied to simulate flashing sprays in different conditions using a RANS based approach (Gärtner et al. 2020(Gärtner et al. , 2021)).Further, its capability to predict shock size and position correctly could be validated with flashing acetone experiments, for which unambiguous measurements of the shock front exist (Gärtner et al. 2020).In contrast to typical one fluid approaches, the solver uses separate energy transport equations for both phases.This separation of the two phases is required due to the significant temperature differences between the two phases at the same location (Gärtner et al. 2020).With this adaption the governing equation system is given by, (1) where the volume fraction is defined as = V L ∕V .The variables and p are density and pressure with the overbars denoting Reynolds averages.Further, u , K, t , h, k Eff , g repre- sent the mixture velocity vector, kinetic energy, turbulent viscous stress tensor, enthalpy, effective thermal conductivity and gravity, respectively.The tildes and double prime superscripts denote Favre averages and their fluctuations, whereas overbars denote Reynolds averages.The subscripts " L " and " G " refer to the fluid properties of the respective phase, subscripts "SG" and "SL" denote saturation conditions of the gas or liquid at the given pressure level, while properties without an index are volume-averaged quantities.The turbulent contributions of the sub-grid stresses in the momentum equation are included in the turbulent viscous stress t .For the energy equation, the effective conductivity k Eff = k t + k is the sum of laminar and turbulent conductivities.Here, the turbulent conductivity is calculated with the turbulent Prandtl number, Pr t , set for all simulations to 1.0.For the transport of the volume fraction the additional term ∇ ⋅ L �� L u �� is modeled with a gradient diffusion assumption (Anez et al. 2019), Here, t is the turbulent viscosity and Sc t the turbulent Schmidt number, which is set to Sc t = 0.5 (Tominaga and Stathopoulos 2007).A detailed derivation of the equation system is found in Gärtner et al. (2020).In the following the overbars and tildes will be omitted in the equations to improve the readability.
In the LES implementation, the turbulent viscosity, t , is modeled with the Wall-Adapting Local Eddy-viscosity model (WALE), which is based on the square of the velocity gradient tensor (Nicoud and Ducros 1999).This model accounts for the turbulent viscosity behavior at the wall and has shown to be superior to the wall damped Smagorinsky model when applied to variable density flows (Huang et al. 2021).
To solve the coupled equation system, the mass and momentum equation are combined to derive a pressure equation, (3) (5) The variable, a p represents the diagonal terms of the momentum equation and H(u) the off- diagonal component.The total derivative of the density requires the compressibility of the mixture, which is expressed by Brennen (2005), From Eq. ( 8), it is already apparent that the compressibility of the mixture is significantly higher than the individual components.This results in a lower speed of sound, such that velocities of 50 ms −1 and below can lead to supersonic conditions (Karplus 1958).This reduction of the speed of sound is the cause for the observed shock structures of flashing sprays for comparable low gas flow velocities.
The empirical homogeneous relaxation model (HRM) is chosen to model the mass transfer, ṁL .This model describes the change of mass fraction, , as being proportional to the difference to equilibrium conditions combined with a relaxation time, Θ , (Downar- Zapolski et al. 1996) Here, is the mass fraction, h L the enthalpy of the superheated liquid, h SG and h SL the satu- ration conditions at the current pressure p.The relaxation time is based on a constant factor Θ 0 , the void fraction, , and the normalized pressure difference, , with p s (T) denoting the saturation pressure and p c the critical pressure of the fluid.The high-pressure fit coefficient for Θ 0 and the exponents and of Downar-Zapolski et al.  (1996) have shown in previous works a good agreement with the flashing behavior of cryogenic liquid nitrogen Gärtner et al. (2020).Therefore, the high-pressure fit coefficients of Θ 0 = 3.84 × 10 −7 , = −0.54, and = −1.76 are selected for this work.

LES Case Setup
The computational domain consists of the injector with a conus as well as the combustion chamber, see Fig. 2.Here we define, the origin of the coordinate system to be at the injector outlet and the z-axis corresponding with the axial direction of the injector.The injector has to be included for a faithful approximation of the well defined boundary conditions of the experiment, which provide mass flow rate, injection pressure as well as temperature at the injector inlet (Rees et al. 2020).The injection and boundary conditions are listed in Tables 1 and 2. Further, it can be assumed that at the injector inlet the nitrogen is still pure liquid.The initial conditions assume that the inlet conditions prevail in the entire injector.At the start of the simulation, the domain size can be reduced to only include the injector and a small part of the combustion chamber.With simulation progress the domain size is adapted until the final dimensions of 35 mm diameter and a length of L = 40 mm are reached.To check the mesh dependency of the solution, three meshes with 40, 67, and 80 million cells are tested, see Table 3.While no significant changes of velocity, vapor distribution, and LES quality (defined by the ratio of resolved to total turbulent kinetic energy) could be found in the chamber, a mesh refinement of the boundary layer at the injector wall (leading to the setup with 80 m cells) is needed to prevent spurious flows at block mesh boundaries, see Fig. 3.The largest mesh is therefore used for the averaging of the LES field.In analogy to experimental data, time is reported as time after start of injection (SOI).
To resolve the shocks present for this flashing flow (Gärtner et al. 2020) with an LES suitable accuracy, the WENO scheme is chosen to discretize the momentum, energy and volume fraction transport equation (Gärtner et al. 2020).To ensure stability of the equation system a maximum Courant number of Co = 0.85 is set, which results in a time step of about 2E-8 s.Further, the OpenFOAM PIMPLE algorithm is used with five outer iterations to reach convergence of the governing equation system.
Fig. 3 Visualization of the mesh and velocity magnitude at the transition from injector cone to pipe for the medium and fine mesh.For the medium mesh the spurious flow, which appear at the block mesh boundary between the circumferential blocks around the core, is circled.The white lines represent the block mesh boundaries

Results and Discussion
In the following section the results of the LES are presented and compared to 2D RANS simulations and the experimental data.The section is split into four parts, firstly validating the phase change model by comparison of the mass flow rate with experiments, secondly investigating the effect of turbulent eddies on the vapor distribution in the injector, thirdly the dynamic behavior of the spray in the combustion chamber and lastly investigating the dark structures in the shadowgraph images at locations at about 20 to 30D downstream of the injector.

Validation of HRM Model Parameters
For a quantitative comparison of the experiment to the simulation, the mass flow rate through the injector is chosen, as it develops freely and is subject to the known pressure boundaries.In addition, the mass flow rate is directly affected by the evaporation rate inside the injector.Therefore it is a suitable quantity to validate the HRM model parameters.To include the turbulent inflow condition, a time-variant input, sampled from a highly resolved LES of the pipe and conus leading up to the injector, is used.To fix the velocity as well as the pressure at the inlet boundary is numerically unstable.Hence, while solving the momentum predictor step, the pressure boundary at the inlet is treated as zero gradient.
The pressure equation then uses the total pressure boundary condition and does not force the velocity flux, calculated from the predictor step, to the velocity boundary.This mimics the fixed flux extrapolated pressure boundary condition of OpenFOAM without relying on a Neumann boundary condition.This process is required as the pressure information of the outlet boundary cannot travel upstream to the inlet due to the supersonic shock developing in the chamber.
The mass flow rate is measured at the inlet patch of the simulation using the liquid density at the given pressure and temperature value.After the end of the start up phase of the simulation, the mass flow rate is constant over time at 10.4 g s −1 , which is slightly below the target value of 10.6 g s −1 .However, considering the comparable large experimental uncertainty of 1.77 g s −1 it is well within the expected bounds.Further, the transition radius of cone to injector pipe affects the separation zone and, with it, the mass flow rate.Compared to the simulation, minor manufacturing uncertainties could also lead to a slight difference between the simulation and the target value.

Effect of Turbulence Modeling in the Injector
To compare the effect of dimensionality and the associated effect on the representation of turbulent structures in the injector, the 3D LES results are compared to a 2D RANS simulation using the same numerical settings and a comparable mesh resolution.Figure 4a shows the instantaneous liquid volume fraction for the LES (top) and the RANS result (bottom).The separation zone at the inlet is larger for the LES and more vapor is generated and then distributed further into the center of the injector by turbulent eddies.A comparison of the mean volume fraction (solid line) and velocity (dashed line) as a function of the normalized radial position is shown in Fig. 4b.The radial distribution of the mean properties confirm the finding from the instantaneous view that more liquid is evaporated and transported further into the injector center.Even though the vapor distributions vary between the two simulation types, the overall mass flow varies only by 3.6% due to an increased velocity of the LES in the near wall region.This shows that resolving the large eddies affects the vapor and velocity distribution significantly inside the injector and that a simpler RANS model is not sufficient to study the injector behavior in detail.

Spray Dynamics in the Combustion Chamber
Due to the optical dense area close to the injector exit and the challenging experimental conditions, it is not possible to get quantitative data of the spray break-up process after the injector exit and only qualitative shadowgraph images are available for comparison, see Fig. 5b.Further downstream of the injector, however, where a more dilute regime can be observed, quantitative measurements of droplet characteristics are possible.Droplet velocities and size distributions are provided downstream of the injector in Rees et al. (2020) and allow for a qualitative and quantitative comparison of simulations with measurements in the combustion chamber.In Fig. 5a the mean velocity distribution at the center plane of the spray is compared to the measured droplet velocities and Fig. 5b shows an experimental shadowgraph image.In the slip stream predictions are quite good but the velocities in the spray center do not match in regions close to the injector.After 20 mm the simulation results are in better agreement with the measured droplet velocities.Close to the injector outlet less than 100 Doppler bursts were detected, hence favoring larger and faster droplets.Due to the low statistics in this region the points in the range of 10 ≤ y ≤ 15 and 0 ≤ z ≤ 10 are excluded from the following analysis, marked by a red square in Fig. 5a.We note that the one fluid method imposes the assumption that the two phases have no relative velocities.The droplet size distributions provided by the PDA measurements suggest however, that the droplets generated by the spray may be too large to fully follow the gas flow.To investigate the behavior of the droplets, passive particles are injected into the simulation at positions where the liquid volume fraction is between 1 to 5%.The particles are initialized with the local velocity value and initial droplet sizes are sampled from the measured size distribution at y = 0 and z = 10 mm that is depicted in Fig. 6.After injection the particle size is fixed and only the drag force is acting on the particles given by the standard sphere drag model (Putnam 1961) The particle velocity field is time and space averaged, with an averaging interval of 0.7 ms and spatial averaging in circumferential direction.For the averaging procedure all particles are tracked and given a unique ID to avoid duplicate counts of the same particle in one averaging bin. Figure 7 shows that the particle velocities match well with the experimental data for the values in the slip stream and at 10 mm downstream of the injector.In addition, the velocity direction of the particle data is displayed as white arrows with the same scaling for the velocity magnitude as used for the experimental data set, showing good agreement of the predicted and measured particle flow directions.Positions relatively ( 13) ) are excluded.The data quality of the experiment at these positions is, to reiterate, poor as only a low two digit number of Doppler bursts were detected (Rees et al. 2020), and it is likely that a bias for larger particles with higher inertia and larger radial velocities exist for this part of the domain.
A quantitative comparison of the experimentally measured velocity magnitude with the LES and particle data at the axial positions z∕D = [5, 10, 20, 30] is given in Fig. 8. Due to the circumferential averaging of the particle data, only results for the positive y-axis are presented for the particles.Further, as experimental points at y∕D = [10, 15] and z∕D = [5, 10] are excluded, only one measurement point remains closest to the injector at z∕D = 5 .The line plots confirm the finding of Fig. 5a, that the LES results (blue line) match well in the slipstream and further downstream, whereas close to the injector and the center line, they significantly deviate.In contrast, the particle data (red line) matches well with the experimental measurements, as observed in Fig. 7.
The correlation of droplet size and velocity is also visible in Fig. 9, which shows the axial velocity of the one fluid solution (black line), the mean particle velocity (simulation -orange line) and the mean droplet velocity (experiments -orange symbols) along the centerline taking averages from the entire particle and droplet populations.Note that the simulated values are an integrated average of all particles found in a radius of 2 mm around the centerline.This is due to the circumferential averaging, which is needed as the the zero volume of the centerline itself prevents sampling and a finite cross-section needs be defined to have a significant number of particles in each averaging bin.The figure also shows a conditioned velocity for the particle population with diameters larger than 10 μ m (blue line).Measured data also include droplet velocities averaged from droplets larger than 10 μ m only and they are indicated here by the blue square symbols.The conditioned (large droplet) velocities of the experiments agree very well with the corresponding (large particle) simulated data.Due to the larger mass and momentum of these droplets they decelerate more slowly than the smaller droplets.The computed velocity of the total particle population, however, decreases much faster and reaches the value of the gas-phase velocity after 20 mm.The difference between the mean particle velocity and the measured droplet values at 15 mm and 20 mm is likely to be due to a bias for larger droplets in the measurements.The experiments recorded fewer than 100 droplets at measurement locations close to the injector with z ≤ 20 mm and a bias for large droplets was observed (Rees et al. 2020).This bias can also be observed in the droplet size distributions at the axial locations of 10, 15, and 20 mm downstream of the injector as shown in Fig. 10.For the locations at 10 mm to 20 mm the experimentally measured mean droplet size increases which is due to the mentioned bias in the measurement, whereas the particle size distribution is shifted slightly towards smaller particles.However, this shift in particle size distribution is not due to evaporation, as the particle size is fixed at the injection, but due to larger particles having a higher inertia and thus continuing to travel radially outwards out of the averaging volume used for circumferential averaging.
These results show that droplets with larger diameters do not follow the gas flow perfectly and that they reach the gas flow velocity only with a delay at about 20 mm to 40 mm downstream of the injector.The exact position depends on the droplet size.Further, larger droplets will travel beyond the extent of the shock size which also explains why the one fluid simulation does not predict a transport of liquid into the regions close to the injector but far away from the centerline (marked as area B in Fig. 5b).The good agreement between the size conditioned particles and the experimental data suggests that the underlying velocity field is correct.As they are only affected by the drag force, that the gas flow exerts on the particles.Further, one-way coupling as effectuated here, seems justified as the volume fraction is about 1% after the spray break-up and the small droplets will not significantly affect the gas flow velocity.Therefore, these results do not invalidate the onefluid solution used to calculate the underlying velocity field but rather demonstrate, that only the motion of the smaller droplets and the gas phase is captured accurately.Capturing the dynamics of larger droplets will, however, require a more detailed approach, such as a combined Euler-Lagrange approach.

Recirculation Zone
In the shadowgraph images of the experiment nearly motionless or slightly upstream floating structures can be found at about 20D to 30D downstream of the injector (Gärtner et al. 2020), see marked area A in Fig. 5b.A 2D RANS study of this case has shown that a recirculation zone develops in this region which could explain the observed behavior.The 3D LES, however, does not develop such a pronounced recirculation zone as can be seen in Fig. 11a, which shows the instantaneous pressure field in the center plane for the 3D and 2D simulations together with the velocity displayed as vector glyphs and colored by the axial velocity component.To improve clarity of presentation the pressure and velocity values are clipped.Instead of a recirculation zone, smaller regions with low axial velocity are present.They have been marked as area 'C' in Fig. 11a.The reason for this mismatch between the two simulation types can be found in the 3D effects and circumferential inhomogeneity that cannot be captured in the 2D simulation.In the 2D simulations, the recirculation originates from a velocity eddy that is created by the alternating low and high pressure pattern in the slip stream which causes the flow to rotate in a counter clockwise direction.This eddy then increases in magnitude by the large shear velocities between the slip stream and core velocity pushing it towards the center line of the spray, where it stays stationary for a short period.While the low/high pressure pattern also exists in the 3D LES, it is not as pronounced and distinctly separated as in the 2D RANS simulation.Further analysis of the pressure pattern in a cross-section at 8 mm downstream of the injector Fig. 11 a Pressure and axial velocity displayed as vector glyphs for 3D and 2D simulations.b cross sectional view of the instantaneous velocity and pressure for the 3D LES 8 mm downstream of the injector helps to explain the differences.We note a clear circumferential variation of pressure and velocity and fluid elements transporting mass into the low pressure regions can "escape" in lateral direction such that large eddies, as observed in the 2D simulations, cannot build up.This 2D effect seems very similar to the well known inverse energy cascade that is characteristic for 2D shear layer simulations, where a roll-up of the initial eddies can be observed.Thus, the existence of the rather strong recirculation zone in 2D appears to be related to the simplified representation of turbulence and the omission of the vortex stretching effect that leads to a decay of large vortex structures, prevents roll-up and the propagation of negative velocities towards the centerline.As a result we observe in the 3D LES small separated regions of low axial velocity in the range of 15 to 30D, but no contiguous recirculation area is present.In terms of combustor development, we can still identify the spray center at around 20D as a suitable location for mixture ignition as the low axial velocities will reduce strain on the ignition kernel, prevent immediate extinction and allow a fully burning flame to stabilize itself.

Conclusion
High order LES of flashing cryogenic nitrogen are performed.The simulations are carried out by a compressible, one fluid-two phase solver developed in OpenFOAM using the homogeneous relaxation model to describe phase change.The phase change model and its parameters have been validated by comparison with the measured mass flow rate through the injector.The effect of resolving turbulent eddies on the evaporation rate in the injector is investigated by comparing the results to 2D RANS simulations.The results show that more vapor is generated for the 3D LES due to turbulent eddies at the wall, a larger separation zone at the injector inlet and that more vapor is transported toward the center of the injector.The spray velocities in the combustion chamber of the one fluid solver match the measured droplet velocities well in the slip stream of the shock but are much lower than the measured values in the spray center just downstream of the shock front.The measured droplet size distributions indicate that not all droplets are small enough to follow the gas flow perfectly and additional one-way coupled particles are injected and tracked for a more meaningful comparison between experiments and simulations.These particles are injected with the droplet size distribution measured closest to the injector and the local velocity of the one fluid solution.The tracked particles provide good approximations of the experimental values in the spray center and in the slip stream, showing that the droplets follow the gas flow with a certain delay of about 10 to 20 D depending on the droplet size.This also matches with the velocity drop between 20 to 30D in the experimental data.Additional analysis of shadowgraph images revealed motionless structures which were explained by a recirculation zone found in previous 2D RANS simulations.This recirculation zone could not be reproduced by 3D LES and it is likely that the strong recirculation in 2D is wrongly caused by an inverse energy cascade since vortex stretching effects are omitted here.However, smaller eddies and zones of low axial velocity continue to be present in the 3D simulation in the same region, and the region can be identified as a suitable location for spray ignition in upper stage rocket engines burning conventional fuels.

Fig. 1
Fig.1Simulation regime of sub-to supersonic as well as incompressible to compressible flow for the case of flashing cryogenic liquid nitrogen.The background color is the liquid volume fraction and the black lines are velocity magnitude iso contours

Fig. 2
Fig. 2 Sketch of the injector with the computational domain marked by a red boundary.The spray with its shock system is outlined by the density gradient

Fig. 4 a
Fig. 4 a Instantaneous liquid volume fraction in the injector for the 3D (top) and 2D (bottom) simulation.The sample line for Fig. 4 is marked by a black line.b Mean liquid volume fraction (solid line) and velocity (dashed line) plotted at x∕D = 2.6 after the injector inlet for the 2D RANS and 3D LES

Fig. 5 aFig. 6
Fig. 5 a Mean velocity of the LES compared to the experimentally measured droplet velocities marked as circles with a scaled velocity vector in black.The white arrows represent the flow direction of the LES solution.b Experimental shadowgraph image of the spray, adapted from Gärtner et al. (2020) with permission of Elsevier

Fig. 7 Fig. 8 Fig. 9
Fig. 7 Velocity of the one-way coupled particles of the simulation projected onto a 2D plane compared to the measured droplet velocities.The white arrows denote the velociy magnitude and direction of the particles in the simulation and th black arrows represent magnitude and direction of the experimental measured droplets

Fig. 10
Fig. 10 Droplet size distribution measured in the experiment at z = 10, 15, 20mm compared to the particle size distribution in the simulation

Table 1
Injection conditions of

Table 3
Domain size and mesh resolution -Mesh size Cell size (μm)