Comparison of cryogenic and non-cryogenic droplet impact dynamics at low velocities using OpenFOAM

The goal of this study is to analyse, through a newly developed framework in OpenFOAM, the impact dynamics of cryogenic droplets. Interaction of cryogenic liquid droplets and jets with solid surfaces is of interest to a wide range of applications. To emphasise the differences between non-cryogenic and cryogenic cases, the same impact condition (based on nondimensional numbers) is simulated for both fluids. First, the impact of a water droplet on a flat, solid surface is investigated. Then, choosing the same Weber, Ohnesorge and Reynolds numbers, a Liquid Oxygen (LO2) droplet into gaseous Nitrogen (N2) hitting a solid surface is simulated. The ambient temperature and pressure are below the oxygen critical point, limiting the investigation at a sub-critical regime. The algebraic Volume of Fluid (VoF) method is employed with adaptive mesh refinement on the interface region. Numerical treatments to improve the interface description are implemented as well. The results obtained for both cases are investigated, comparing the droplet morphology evolution for the two fluids. Differences are observed mainly in the receding stage, once the droplet has reached the maximum spreading, with the receding stage of the cryogenic case characterised by a faster dynamic.


Introduction
The interest in the study of cryogenic applications has recently increased among the academic and industrial community, due to the potential use of cryogenic liquids to a range of fields including energy storage, cooling, quantum computing and even biomedical applications. Interaction of cryogenic liquid droplets and jets with solid surfaces is one of the more intriguing areas within the cryogenic field that however only limited research has been performed, mostly from the experimental side. For example an interesting application involving impact dynamics is in cryosurgery, where the extreme cold produced by Liquid Nitrogen (LN 2 ) vaporisation is used to destroy abnormal tissue. LN 2 is applied directly to the cancer cells with a spraying device. Various physical phenomena occur during the impact of a non-cryogenic liquid droplet on a solid substrate, such as the spreading, fingering, air entrapment, coalescence, shedding, solidification, bouncing, and splashing. Several reviews are present in literature that classify these different phenomena (see for example [1]). Physically these phenomena are mostly controlled by the fluid properties, the surface properties and the impact velocity. How these parameters control the impact of cryogenic fluids has not been fully investigated yet.
Numerical simulations are fundamental to highlight the role of the different physics behind droplet impact since they can provide insight into processes that experimentation is hard. Several numerical approaches have been proposed in the literature over the years for the modelling of multiphase interfacial phenomena. Great effort has been dedicated to correctly account for the transport of the interface and the consideration of the capillary forces. The role of the capillary forces is important both in cases with low velocity impact where capillary forces dominate throughout the process as well as to high impact case where although in the initial stages of the impact inertia is dominant as the droplets spread and slow down, capillarity becomes significant [2,3,4]. As regards the cryogenic droplets impact, when they are the subject of investigation the focus is mainly on phenomena relevant to heat transfer and on phase change. Rebelo et al. [5] studied the evaporation process of single liquid nitrogen droplets when submerged into a non-cryogenic liquid, showing the influence of the initial droplet size and of the surface tension on the evaporation rate. Van Limbeek et al. [6] studied the impact of a single liquid nitrogen drop on a smooth sapphire prism through high-speed frustrated total internal reflection imaging. Varying the prism temperature and impact velocity of the drops a phase diagram of the impact characteristics was observed. They showed how also in these conditions the cooling power of a drop is strongly related to the wetting behaviour of the impacting drops.
The aim of our study is to focus on the wetting behaviour of cryogenic droplets, providing novel insight into whether the current literature established for non-cryogenic conditions can be applied to the case of cryogenic drops. In this regards, the same impacting conditions, characterised by the Weber, Ohnesorge and Reynolds numbers, are applied to a cryogenic and non-cryogenic case. First, the impact of a water droplet on a flat, solid surface is investigated. Then, the impact behaviour onto a solid surface of a Liquid Oxygen (LO 2 ) droplet moving into gaseous N 2 is simulated at the same conditions. The algebraic Volume of Fluid (VoF) method is employed with adaptive mesh refinement on the interface region. Numerical treatments to improve the interface description are implemented as well. The simulations have been performed in OpenFOAM with a newly developed code. The results obtained for both cases are analysed, focusing on the comparison the droplet morphology evolution.

Governing Equations
Our suggested framework is designed for simulations of compressible and immiscible fluids. The inteface tracking of the two-phase flow is performed with the VoF method as basis using the Open Source code OpenFOAM [7]. The native solver compressibleInterFOAM was used as starting point and modifications were made to account for the interface treatment. Being an one-fluid method, the continuity, energy and momentum equation for the mixture are solved: where p d = p − ρg · x is the peziometric pressure, f b the body force and the mixture stress tensor τ is defined as To track the single phase, the transport equation for the volume fraction is solved as only two phases are present, the volume fractions obey the algebraic relationship α 1 + α 2 = 1. The volume fraction is transported by the mixture velocity u. The transport properties of the mixture are calculated from the single phase properties through the volume fraction, for example the viscosity is defined as µ = µ 1 α + µ 2 (1 − α) while the density is defined as When dealing with bounded scalars such as α, it is fundamental to respect the boundedness of the variable. For the α transport equation, Eq. (3), the convective term can be source of unboundedness. Among the different strategies to mitigate the issue of the boundedness, in this study the Multidimensional Universal Limiter with Explicit Solution (MULES), proposed by Weller [8], is adopted. The MULES limiter consists of an explicit solver which is based on the Flux Corrected Transport in order to maintain the boundedness of the scalar. A detailed description of the algorithm is beyond the scope of this study, and more details about its implementation can be found in literature [9]. Despite being derived in an incompressible framework, the MULES scheme can be directly applied to the compressible case. To be suitable to this scheme, the volume fraction equation is rearranged and the contribution from the compressibility is accounted numerically as a source term. This allows to hold a similar numerical methodology for the transport of the volume fraction between the compressible and incompressible solvers. The methodology presented is pressure-based, which means that the pressure is solved, with the phase densities obtained algebraically through the equation of state.
As the gravity force is included in the pressure gradient through p d , the only body force present is the surface tension force. The surface tension force is treated as a pressure gradient across the liquid-gas interface and is calculated per unit volume based on the Continuum Surface Force (CSF) model proposed by Brackbill et al. [10].
where the mean interface curvature is defined as κ = ∇ · n I and the normal to the interface is defined as n I = ∇α/|∇α|.
As well known in literature [11], the evaluation of the interface properties in algebraic VOF can results in an imbalance of the pressure and surface tension force at the interface, thus in the presence of spurious currents. To mitigate this problem, the interface region is subject to numerical treatment to calculate the interface properties and the surface tension force. These treatments consists in a smoothing operation of the volume fraction for the calculation of the curvature, and on the sharpening for the evaluation of the surface tension force location. More details about their implementation into the method and their influence can be found in Tretola and Vogiatzaki [12].

Set-Up
The numerical domain and initial configuration is shown in Fig. 1. The domain has an extension of L X = L Y = 4 × d and L Z = 2.5 × d. On the x and y direction 160 cells are used, while 100 on the z direction. A level of refinement of 3 is employed, which means that the cell is halved up to 3 times. The refinement is operated in the region where 0.005 < α < 0.995. At the initial time step, the spherical droplet is set 0.1 × d 0 away from the solid wall. The velocity field within the droplet is initialised with the initial velocity of the moving droplet u 0 . The top and lateral boundaries are treated as open boundaries, while on the bottom boundary the no-slip boundary condition is applied for the velocity fields, while a Neumann boundary condition is satisfied for the other fields. Applying this condition for the volume fraction means to neglect the influence of the contact angle on the boundary, as the focus of this study is only on the influence of the cryogenic fluids on the impact behaviour. Moreover exact values or guidelines for cryogenic contact angles do not currently exist. In a future study a sensitivity analysis of the contact angle influence will be performed.
The conditions investigated are summarised in Table 1. The cryogenic case and the noncryogenic case have the same non dimensional parameters, characterised by Weber, Reynolds and Ohnesorge numbers. Although these numbers are maintained the same in terms of absolute numbers, at cryogenic conditions, the gaseous phase is characterised by lower viscosity and higher density compared to the non-cryogenic case and the surface tension is lower. This also results in a different impact velocity.

Results and discussion
The operating point investigated falls into the rebound regime (based on the categorisation of non-cryogenic fluids). Figure 2 shows the evolution over time of the droplet as the iso-surface α = 0.5. The initial dynamics of the spreading are similar for the two cases and follow a non-splashing regime, as expected at this operating point [13]. At the initial stages, for both cases the droplet at the spreading phase is forming a circular crown at the side. Once the maximum spreading diameter has been reached, differences start arising between the cases: for the non -cryogenic case the droplet starts receding keeping a disk shape, until it reaches the centre of the disk where a liquid jets formation is observed. For the cryogenic case, the receding stage is different. The shape is not a disk and the jet formation at the centre happens earlier than the noncryogenic case.

Water-Air
To analyse the source of these differences, the velocity distribution plot is included. To compare the two cases, the velocity is normalised with the impact velocity u 0 for each case, see Table  1. Figure 3 shows the velocity distribution at the mid section, normalised by the corresponding impact velocity for each case, at different time steps. Looking at the mid-section distribution, the main differences in the velocity field are observed. The spreading is very similar for both the cases, with the velocity peak of |u| ≈ 2u 0 present at 0 |u|/u 0 Figure 3. Velocity distribution at the mid section. Velocity normalised by the corresponding impact velocity u0 for each case, see Table 1.
the lamella extending at the wall. Once the maximum spreading is reached, the cryogenic case recedes faster. The edges are larger for the cryogenic case and when the droplets recedes to the centre, resulting into a large blob with higher velocity at the tip. The non-cryogenic case instead presents a higher velocity in the gaseous phase, right above the droplet, which emerges from the receding in an elongated shape. Figure 4 shows the spreading at the wall r w over the time for the cryogenic and non-cryogenic case. The spreading is normalised by the initial droplet d 0 . A schematic representation of the spreading is showed as well.  The spreading stage is identical for both the cases, which reaches the maximum spreading at the same instant. Once the peak is reached, the cryogenic case shows a faster decay than the non-cryogenic. In both the cases the decay is linear, with a slope of ≈ 0.37 for the cryogenic case and of ≈ 0.29 for the non-cryogenic case.
To have a further analysis of the morphology, the centre of mass and mass velocity of the droplet are investigated, defined respectively as: • Centre of mass of the droplet C d : • Droplet Velocity u d : Figure 5 shows the trend over time of the droplet centre of mass and velocity magnitude and components for both the cases. The centre of mass shows a faster and higher increase after (d) reaching the minimum for the cryogenic case, with the final centre of mass slightly larger than the initial one. As regards the droplet velocity, after reaching a minimum, the velocity increases up to a maximum value similar for both the cases (|u d | ≈ 0.3u 0 ) but not at the same time, with the cryogenic case reaching this maximum earlier than the non-cryogenic one.

Conclusions
In this study, the wettability behaviour for the cryogenic droplet is observed. The heat transfer is not considered, focusing only on the morphology of the impact. First the shape of liquid droplet during the impact is analysed. The differences between the cases arise during the receding stage, where the cryogenic droplet recedes faster towards the centre than the non-cryogenic one. As a results a larger blob is rising after the liquid edge reaches the centre of the droplet.
To quantify the comparison, the comparison then focused on the spreading length r w and on the droplet shape parameters as the centre of mass and velocity, C d and |u d |. The normalised maximum spreading at the wall is similar for both the cases and it is reached at the same instant. After the peak, the spreading for the cryogenic droplet decays faster than the noncryogenic case. The decaying stage can be represented by a linear relation, with the slope of −0.37 for the cryogenic case and of −0.29 for the non-cryogenic case. This is reflected on the centre of mass of the droplet. After the impingement, the centre of mass rises faster over the initial value, while for the non-cryogenic case C d rises slower and at a lower value. The different behaviour observed can be due to the role of the gaseous phase. At cryogenic conditions, the gaseous phases are characterised by lower viscosity and higher density compared to the noncryogenic case. This difference results in the receding stage, once the initial kinetic energy of the droplet is dissipated. In the future investigations a wide range of operating conditions will be studied, introducing also various contact angles.