Guillaume Vignat
email: guillaume.vignat@centralesupelec.fr

Preethi Rajendram Soundararajan

Daniel Durox

Aymeric Vi

Antoine Renaud

S Ébastien Candel

A Joint Experimental and LES Characterization of the Liquid Fuel Spray in a Swirl Injector

The quality of liquid fuel spray injection determines to a large extent the steady state performance and dynamics of gas turbine and aero-engine combustors. The present investigation is focused on the detailed characterization of the liquid fuel spray in a single sector test rig targeted at aeroengine applications. The liquid fuel (heptane) is injected in a hollow cone spray pattern by a simplex atomizer and the injector comprises a radial swirler. Two features of the droplet distribution that are less commonly found in the technical literature are identified. First, the distributions of mean droplet diameters exhibit non-axisymmetric patterns, a lack of symmetry that is investigated for three types of swirlers differing by their swirl number and/or head loss. Second, it is found that the size-conditioned velocity distributions feature a single wide peak for small droplets and become bimodal for the largest droplets, with a first peak at low velocities, and a second one at higher velocities. The spray behavior analysis is complemented by making use of Large Eddy Simulations with Lagrangian Particle Tracking. Droplet injection is achieved with a model in which the initial size and velocity distributions are specified from experimental data in the atomizer near field. The initial spray interacts with the lateral injector surface and requires a droplet-wall interaction model accounting for the existence of a liquid film. Simulations do not retrieve the lack of rotational symmetry that is found experimentally indicating that this is not linked to the nature of the swirling flow. This is also consistent with further experiments with a different atomizer confirming that this is due to imperfections in the initial atomizer geometry. Another result is that certain swirler designs appear to be more robust to these atomizer imperfections. Simulations accounting for the liquid film yield a bimodal distribution for the droplets' axial velocity distribution which would not be obtained without this model indicating that it is important to represent the droplet-wall interaction, a feature that is not commonly found in the literature.

INTRODUCTION

Liquid-fueled aero-engine injectors have to handle a wide range of fuel flow rate from engine idle to full power conditions (typically 1:40) while ensuring good atomization in order to maintain combustor performance levels, in terms of pollutant emissions and operability (ignition, blow out, combustion instabilities...). To this end, many larger engines use several atomizers, or even several injectors, one of which, dubbed "pilot", is designed to handle lower flow rates and ensure reliable operation under low power settings [START_REF] Lefebvre | The role of fuel preparation in low-emission combustion[END_REF]. Pressure-swirl atomizers are often used for that purpose as they give access to the required performances by generating a spray of fine droplets in a hollow cone pattern with a large angle and a homogeneous circumferential liquid distribution even under low fuel pressure and flow rates [START_REF] Lefebvre | Gas turbine combustion : alternative fuels and emissions[END_REF].

The general objective of the present article is to investigate the complex interaction taking place inside the injector between the swirling flow, the droplet cloud produced by the pressure atomizer and the injector geometry. This question is examined by combining detailed experiments and simulations. On the experimental side attention is first focused on the impact of the injector design on the circumferential droplet size distribution. The velocity distribution is measured and its features are interpreted using Large Eddy Simulation (LES).

At this stage it is worth reviewing the literature dealing with these various issues. Early characterizations of pressure-swirl atomizers are reviewed by Lefebvre [START_REF] Lefebvre | Atomization and Sprays[END_REF]. Much of the work was focused on obtaining correlations between injector geometrical dimensions and performance such as discharge coefficients, cone angle and overall droplet diameter. As atomizers are integrated inside injectors, interactions between the air flow and the spray must be considered. This is accomplished for example by Rajamanickam & Basu [START_REF] Rajamanickam | Insights into the dynamics of spray-swirl interactions[END_REF] who study the interactions between the hollow cone liquid sheet formed by a simplex atomizer and a swirling air flow by systematically varying the air flow rate to investigate different flow regimes. Time resolved measurements show that the dominant mechanism for sheet break-up and primary atomization is driven by surface tension effects when the momentum of the liquid is large compared to that of the air. At higher air flows, the motion of the liquid is strongly coupled with that of the air leading to flapping and catastrophic break-up of the liquid sheet by vortices in the air flow. The angle of the spray also increases with increasing air flow rate, a result also reported in [5] where the spray angle of a simplex atomizer is shown to increase when the air flow rate and swirl number of the surrounding air coflow are augmented. The dynamics of the disperse phase is investigated for example in [START_REF] Hardalupas | Particle Dispersion in a Vertical round Sudden-Expansion Flow[END_REF] where particle dispersion takes place in a jet after a sudden expansion. Adding a swirling motion can lead to segregation depending on droplet size and associated Stokes number, and constitutes the working principle of cyclonic separation [START_REF] Dring | Particle trajectories in swirling flows[END_REF][START_REF] Kriebel | Particle Trajectories in a Gas Centrifuge[END_REF]. Detailed characterization of droplet size and velocity in laboratory scale aeroengine-type combustion chamber have been made possible by advances in laser diagnostics, in particular by the development of the Phase Doppler Anemometer (PDA). An early application of this instrument by McDonnell et al. [START_REF] Mcdonnell | Symmetry assessment of an air-blast at-omizer spray[END_REF] showed its potential for laboratory characterization of sprays formed by aero-engine injectors. More recently, PDA has been used by Lecourt et al. [10] to fully characterize an aeronautical injector in non-reactive and reactive conditions and provide data that was then used to initiate numerical simulations [START_REF] Jaegle | Eulerian and Lagrangian spray simulations of an aeronautical multipoint injector[END_REF][START_REF] Guedot | Modeling and Analysis of the Interactions of Coherent Structures with a Spray Flame in a Swirl Burner[END_REF][START_REF] Franzelli | Large Eddy Simulation of Swirling Kerosene/Air Spray Flame Using Tabulated Chemistry[END_REF]. These studies explored Eulerian and Lagrangian particle tracking approaches to model spray combustion. One lesson learned from these investigations is that a complex interaction takes place between the spray and the Precessing Vortex Core (PVC). This helicoidal coherent vortical structure that is commonly found in strongly swirling flows is also present in aeronautical injectors. This interaction subsequently perturbs the flame. Recent experiments by Renaud et al. [START_REF] Renaud | Experimental Study of the Precessing Vortex Core Impact on the Liquid Fuel Spray in a Gas Turbine Model Combustor[END_REF] on spray-PVC interaction indicate that the PVC attaches at the tip of the atomizer, and produces a periodic disturbance of the fuel supply. Although the PVC vanishes at the injector outlet, the perturbation on the spray still influences the flame. A more detailed, joint experimental and numerical study on the interaction between a spray formed by a prefilming atomizer and a PVC was carried out by Keller et al. [START_REF] Keller | Numerical and experimental investigation on droplet dynamics and dispersion of a jet engine injector[END_REF]. They compared several numerical simulation approaches and note that contrary to Unsteady Reynolds Averaged Navier Stokes (URANS), Large Eddy Simulation (LES) is better suited to predict interactions between droplets and vortical structures, leading to an improved modeling of fuel-air mixing. Their simulations also feature spray clustering: the large droplets gather and form a helical shaped volume of higher liquid loading and equivalence ratio, in phase with the PVC, somewhat at variance with Franzelli et al. [START_REF] Franzelli | Large Eddy Simulation of Swirling Kerosene/Air Spray Flame Using Tabulated Chemistry[END_REF] who find fuel loading perturbations that are delayed with respect to the PVC.

In the area of combustion dynamics, the response of droplets and of the fuel atomizer itself to acoustic perturbations is also found to play a key role in the flame response. Cold flow and spray response to acoustic perturbations of a swirl spray injector examined experimentally and with LES in [START_REF] Gajan | Investigation of Spray Behavior Downstream of an Aeroengine Injector with Acoustic Excitation[END_REF] indicate that a slow convective wave of fuel loading is formed. This is attributed to the modulation of the droplet transport by the air flow, and to the modulation of the atomization process in the air-blast atomizer. Kim et al. [START_REF] Kim | The behavior of liquid fuel sprays in acoustically-forced air swirler flows[END_REF] used LES in a similar cold flow study. Eckstein et al. [START_REF] Eckstein | Experimental study on the role of entropy waves in low-frequency oscillations in a RQL combustor[END_REF] studied a self-excited oscillation in a single sector, laboratory scale, kerosene fueled aeronautical burner. They observe a modulation in droplet diameter being released by their injector and note that subsequent modulation on evaporation and droplet acceleration delays may play a key role in the dynamics of their flame. The joint experimental and numerical (LES) investigation [START_REF] Tachibana | Experimental and numerical investigation of thermo-acoustic instability in a liquid-fuel aero-engine combustor at elevated pressure: Validity of large-eddy simulation of spray combustion[END_REF], carried out on a high pressure, kerosene fueled test rig, underline the importance of the delay introduced by droplet evaporation on the thermoacoustic feedback loop. Finally, perhaps the best demonstration of the impact of liquid fuel injection properties on combustion instabilities is provided by Lee et al. [START_REF] Lee | Slow" active control of combustion instabilities by modification of liquid fuel spray properties[END_REF], who use an actively controlled Nanomiser™to vary the fuel droplet size of their injector, allowing them to move from a stable to an unstable combustion regime of operation.

The present analysis is motivated by observations of combustion dynamical phenomena in an annular combustor "MICCA-Spray", an atmospheric test rig at EM2C laboratory used to investigate ignition and azimuthally coupled combustion dynamics. This system has been more recently converted to allow liquid fuel injection in order to be more representative of a small helicopter engine combustor [START_REF] Vignat | High Amplitude Combustion Instabilities in an Annular Combustor Inducing Pressure Field Deformation and Flame Blow Off[END_REF][START_REF] Prieur | Ignition dynamics in an annular combustor for liquid spray and premixed gaseous injection[END_REF]. To better understand flame dynamics in this rig, it is necessary to take a closer look at fuel injection in this configuration and obtain information that may apply more generally to swirl injectors comprising a pressure atomizer.

The aim of this article is to present detailed measurements of cold flow spray that are necessary inputs to high fidelity simulations. These data are obtained on a single sector featuring a central injector "SICCA-Spray" (a sister test rig to "MICCA-Spray"). This configuration will first be introduced, along with the design of the injector and the diagnostic tools that have been employed. Several swirlers are used to assess the impact of swirler design and in particular that of the swirl number on the circumferential distribution of fuel. The spray is then characterized using laser tomography and Phase Doppler Anemometry. These data are then combined with Large Eddy Simulations. This is accomplished by first presenting the numerical setup, validating calculations with experimental results and finally gathering insights from the LES and experimental measurements to discuss key aspects of the injected spray.

EXPERIMENTAL SETUP AND INJECTOR DESIGN

The "SICCA-Spray" atmospheric test rig (Fig. 1) comprises a plenum fed with compressed air. In the present unconfined cold flow experiments, the mixture of air and fuel created by the injector is dumped to the atmosphere. The air flow rate is controlled by two Bronkhorst EL-Flow mass flow controllers (relative accuracy of 1%). Liquid heptane is injected through the atomizer, and its flow rate is measured by a Bronkhorst CORI-Flow controller, with a relative accuracy of 0.5%. For clarity, we will designate in this paper the atomizer as the element through which the liquid fuel forms a spray (orange in Fig. 1), while "injector" is used to designate the whole system placed between the air plenum and the combustion chamber through which air is channeled, mixed with fuel and prepared to stabilize combustion. The injector consists in an axial channel that is fed with air from the plenum and conveys the flow to a tangential hole swirler. The flow is set in rotation and starts mixing with the fuel before being exhausted in the chamber or in the atmosphere through first a 5 mm long conical section followed by a 1 mm long cylinder. The outlet diameter of the injector is d in j = 8 mm. The atomizer (orange in Fig. 1(c), part number IJ2582-0-00) is manufactured by ADIndustries Hydraulics (formerly MicroMécanique Pyrénéenne, Oloron-Sainte-Marie, France) and is a pressure-swirl, simplex design with a manufacturer specified flow number of FN = 0.26 l.h -1 bar -0.5 and an atomizer final orifice radius of r at = 40 µm. Following the final orifice, the atomizer features a diverging cup with an included half angle of 59 deg and a final diameter of 1.5 mm. The tip of the atomizer is set 5.6 mm in recess from the dump plane.

The spray is characterized with a Dantec Dynamics FlowExplorer Fiber PDA 2 component Phase Doppler Anemometer that simultaneously measures fuel droplet diameters and velocities in a small probe volume (0.14 mm × 0.14 mm × 0.23 mm). Spray tomographies are obtained with a 600 mW Nd:YaG continuous laser operating at 532 nm wavelength. A Panasonic Lumix FZ38 camera is used for acquiring images.

SWIRLER DESIGN IMPACT ON FUEL SPRAY PAT-TERN

It is first interesting to examine the impact of the swirler design on the fuel spray circumferential pattern. A homogeneous circumferential distribution of the spray is usually desired to ensure symmetry of the flame and avoid formation of hot spots which may reduce the service life of the high pressure distributor and increase pollutant emissions [START_REF] Lefebvre | Gas turbine combustion : alternative fuels and emissions[END_REF]. This characteristic of the spray is however seldom studied. Mc-Donnell et al. [START_REF] Mcdonnell | Symmetry assessment of an air-blast at-omizer spray[END_REF] performed such a study on three methanolfueled air-blast atomizer used in a helicopter engine, and showed the advantage of PDA against the more traditional patternator for such a study. However, only manufacturing dispersion was considered.

In the present section, the spray characteristics are examined for six variants of the injection system. Two different atomizers are used: the first designated as AG was observed to form a circumferentially homogeneous spray. The second, AD, was damaged by repetitive electrical discharges during plasma assisted combustion experiments. Its circumferential spray pattern is skewed.

The sprays formed with the two atomizers are characterized for three different swirlers, dubbed 707, 712 and 716, designed by varying the diameter d th and the excentricity R 0 of the tangential holes (see Fig. 1(d)). Note that swirler 716 is the one being used in a recent investigation of azimuthal combustion instabilities [START_REF] Vignat | High Amplitude Combustion Instabilities in an Annular Combustor Inducing Pressure Field Deformation and Flame Blow Off[END_REF]. The operating point considered in what follows corresponds to a heptane mass flow rate of 500 g h -1 and an air flow rate of 2.3 g s -1 .

The aerodynamic characteristics of the different swirlers are summed up in Tab. 1. The head loss ∆P of the air flow between the upstream side of the injector and the atmosphere is measured using a Kimo MP111 differential pressure gauge with a 1% relative accuracy. Laser Doppler Velocimetry (LDV) yields flow velocities in the absence of fuel injection. In these experiments, the air is seeded with silicon oil droplets with a Mean Diameter of around 2.6 µm [START_REF] Durox | Flow seeding with an air nebulizer[END_REF]. Results are shown in Fig. 2. Swirler 707 and 712 are designed to have identical velocity profiles downstream of the injector,

S = 2R in j 0 U θ U x r 2 dr R in j 2R in j 0 U x 2 rdr (1)
are quoted in Tab. 1. The dynamics of the Precessing Vortex Core (PVC) is also investigated in the absence of spray using a constant temperature hot wire anemometer (Dantec-Dynamics miniCTA with 55P16 probe). The frequency f PVC of the m = 1 PVC (first order in azimuth, single helix PVC) is reported for all three swirlers in Tab. 1. Although the external, time-averaged aerodynamic profiles of swirlers 707 and 712 are quite similar, one observes that the frequency of the PVC is higher for swirler 707. The PVC frequency of swirler 716 is also close to that of swirler 712. However, this injection system features an intermittency between a first order (m = 1), single helix PVC, and a second order (m = 2), double helix PVC with a frequency twice that of the single helix PVC. Characterization of a double helix PVC in SICCA-Spray is presented using time-resolved particle tomography and LES in [START_REF] Vignat | Tomographie à très haute cadence de flammes turbulentes swirlées au moyen du dioxyde d'étain[END_REF].

The circumferential spray droplet size distribution is determined by moving the laser probe volume of the PDA system on a cylindrical mesh. The Sauter mean droplet diameters d 32 are displayed in polar coordinates in Fig. 3 by placing a symbol at a distance from the center corresponding to the d 32 value. The plots correspond to three radii with respect to the injector axis. The patterns pertaining to the damaged atomizer AD and swirlers 707 and 712, are not axisymmetric. This is confirmed in the polar plots of the mean diameter d 10 corresponding to the same injection systems and operating conditions (not shown here). The circumferential Sauter mean droplet diameter distribution pertaining to swirler 716, with the damaged atomizer AD, is more regular. When atomizer AG is employed, the circumferential pattern of the spray is much more homogeneous for all three swirlers. A slight asymmetry is still present with swirler 716 and atomizer AG and found to be quite reproducible. From these observations, one concludes that the damaged atomizer AD is mainly responsible for the skewed circumferential distribution, while an undamaged atomizer AG creates a fairly homogeneous spray. One also finds that an injector with a higher level of swirl interacts with the spray such that the droplet pattern is more robust to defects in the atomizer. It is however not possible to conclude at this point about the physical processes that homogenize the spray distribution or create the slight asymmetry with swirler 716. The internal aerodynamics of swirler 716 are probably conducive to a higher level of turbulent mixing. It appears indeed from velocity data measured downstream of the injector (Fig. 2) that the turbulent intensity is higher with this swirler compared to 707 and 712. An- other possibility may be that the interaction between the flow and the conical liquid sheet formed at the tip of the atomizer plays a role in the azimuthal distribution of droplet sizes. The higher swirl number of swirler 716 is creating a stronger recirculation zone near the tip of the atomizer. These features may have a direct impact on the atomization process itself [START_REF] Rajamanickam | Insights into the dynamics of spray-swirl interactions[END_REF]. One reviewer also suggests that the PVC might be the cause of the asymmetry observed with swirler 716 and atomizer AG. Indeed, if the transition between the single and double helix PVC of this swirler happens preferentially at some azimuthal position, it may be the root of the deformed pattern that has been observed.

DETAILED EXPERIMENTAL CHARACTERIZA-TION OF THE FUEL SPRAY

The focus is now placed on the detailed analysis of PDA measurements obtained using atomizer AG and swirler 707. this region is quite low, as is well exemplified in the laser tomography image in Fig. 5. This is a common problem when computing Sauter mean diameter from PDA measurements as noted in [START_REF] Mcdonnell | Interlaboratory comparison of phase Doppler measurements in a research simplex atomizer spray[END_REF]. In general, the droplet diameter profiles have a bowl shape with larger droplets on the outside, which is expected due to segregation based on the Stokes number in swirling flows [START_REF] Dring | Particle trajectories in swirling flows[END_REF][START_REF] Kriebel | Particle Trajectories in a Gas Centrifuge[END_REF]. One may note however that the difference between the profiles of d 10 and those of d 32 decreases sharply with the distance from the backplane, indicating that an originally highly polydisperse distribution at z = 5 mm becomes nearly monodisperse at z = 20 mm and that this is probably due to a combination of swirl induced particle segregation and vaporisation. Finally, one may also note that in the tomographic image of Fig. 5, the signal is very intense near the backplane around r = r in j = 4 mm. It was observed in this system that a film of liquid droplets is formed on the conical convergent section of the injector, and that this film is being atomized in an air blast manner. Droplet size histograms are shown in Fig. 6 for different radial positions at z = 5 mm. These positions correspond to areas where large numbers of droplets are present (see Fig. 5). The distributions are highly polydisperse, and essentially feature three peaks: the first, at approximately d = 10 -15 µm appears small here in this volumic representation, but corresponds to the main peak of the droplet number distribution. The second is around d = 30 -40 µm while the third represents a small number of relatively large droplets (d 70 µm). One may suspect that these multiple peaks arise from the different atomization processes taking place in the experiment, and result: from the flow interacting with the atomizer spray; from the liquid film escaping from the conical convergent; and from air-blast by the swirling air flow in the injector endpiece.

To further analyze the spray characteristics, the distribution of axial velocity conditioned by the droplet diameter is shown in Fig. 7. For smaller droplets, the distribution is fairly wide and continuous, as could be expected as these droplets behave essentially like tracers in a highly turbulent flow. For larger droplets, the velocity distribution is much narrower with a sharp peak around 29 m s -1 . At r = 5 mm, a second peak seems to be present at a much lower velocity for larger droplets. This is confirmed by the velocity histograms in red in Fig. 8. These velocity histograms pertain to droplets with a diameter larger than 45 µm. At r = 4 mm the histograms exhibit a single peak around 29 m s -1 . At r = 5 and r = 6 mm the distribution becomes bimodal with a second peak at a much lower velocity. The velocity statistics for the air flow at the same positions are shown in black alongside those of larger droplets in Fig. 8. These histograms feature a single very wide peak at a velocity that notably differs from that of the large fuel droplets.

Such a bimodal velocity distribution for large droplets is also identified by Wang et al. [START_REF] Wang | Experimental Study of a Model Gas Turbine Combustor Swirl Cup, Part 11: Droplet Dynamics[END_REF] in the injector from the CFM56 series of engines. This injector is a two-stage counterswirling design. Fuel is injected on the centerline in a hollow cone pattern and subsequently hits the wall of a venturi, forming a film. The authors attribute the bimodality of the velocity distribution to two possible and distinct droplet histories: those formed by the atomizer and those formed by the film have different velocities. Given that filming is also observed in the "SICCA-Spray" experiments, a similar explanation for the bimodality of the velocity distribution can be invoked here: some droplets directly originate from the atomizer while others are formed in the liquid fuel film located on the convergent section of the injector. However, in the present case, the highest velocity values correspond to the droplets from the atomizer. Indeed, using a simple mass and momentum balance on the atomizer, one concludes that droplets are injected with an axial velocity in the order of 30 to 50 m s -1 and, given their large size, their velocity does not change much when they reach the measurement section. This will be confirmed and complemented with the large eddy simulations carried out in what follows.

LARGE EDDY SIMULATION OF THE GASEOUS AND DISPERSED PHASES Euler-Lagrange numerical setup

The present simulations are carried out on the injector equipped with swirler 707, using the fully compressible, LES-filtered, Navier-Stokes solver AVBP jointly developed by CERFACS and IFPEN [30]. Three gaseous species (O 2 , N 2 and C 7 H 16 ) are transported under a unity Lewis assumption. Integration is based on TTGC, a third order in space and time numerical scheme [START_REF] Colin | Development of High-Order Taylor-Galerkin Schemes for LES[END_REF], along with NSCBC bound-ary conditions [START_REF] Poinsot | Boundary conditions for direct simulations of compressible viscous flows[END_REF]. Boundaries within the injection system are modeled using a logarithmic law of the wall. The impact of subgrid scales on the resolved momentum is accounted for using the SIGMA model [START_REF] Baya Toda | A subgridscale model based on singular values for LES in complex geometries[END_REF]. The extent of the numerical domain is that of the schematic view in Fig. 1: the plenum is shortened, and the injector exhausts to an atmosphere with a slow coflow (1 m s -1 ) added for stability. The mesh comprises 28 million tetrahedral elements. Starting from an initial, uniform mesh size in the injector (∆x = 0.3 mm), the automatic mesh refinement method proposed in [START_REF] Daviller | A Mesh Adaptation Strategy to Predict Pressure Losses in LES of Swirled Flows[END_REF] is used to reduce the cell sizes down to ∆x 60 µm. Validation of gaseous phase velocity profiles is performed on the purely gaseous computation without droplet injection. Numerical results are compared to LDV measurements in Fig. 9 showing a very good agreement for the three velocity components and the velocity fluctuations.

The dynamics of the disperse liquid phase is represented with a Lagrangian formalism. The Schiller-Naumann correlation [START_REF] Schiller | A Drag Coefficient Correlation[END_REF] is used to compute the drag force on the particles. Evaporation is estimated using the model proposed by Abramzon and Sirignano [START_REF] Abramzon | Droplet vaporization model for spray combustion calculations[END_REF] with a constant value for Prandtl and Schmidt numbers: Pr = 0.976 and Sc = 1.343 [START_REF] Sanchez | Modeling the dispersion and evaporation of sprays in aeronautical combustion chambers[END_REF]. Dependence of evaporation on the droplet Reynolds number uses a Frossling correlation. Two-way coupling between the disperse and gaseous phase is taken into account using a second order linear interpolation. Gravity and shear effects are considered negligible compared to drag and the time integration for the trajectories uses a two-step Runge-Kutta scheme.

Particle-wall interaction is a critical aspect in this simulation, because some of the droplets hit the wall of the coni- cal end piece. Particles may splash on the wall, rebound, slip in an isolated manner or form a film. This simulation uses the film model proposed by Chaussonnet et al. [38]. This model imposes the velocity of particles close to the wall as a group based on the solution of a set of Saint Venant Eulerian equations for film thickness and momentum.

Liquid injection model

A critical aspect of spray simulations in a Lagrangian framework is the specification of the initial position, size and velocity distributions for the particles. There are various models in the literature that describe pressure-swirl atomizers such as the FIMUR model [START_REF] Sanjosé | Fuel injection model for Euler-Euler and Euler-Lagrange large-eddy simulations of an evaporating spray inside an aeronautical combustor[END_REF]. However these models do not yield velocity distributions that are consistent with what is known from the present measurements. Indeed, it can be experimentally observed, as shown in Fig. 10, that the atomizer and the air flow interact, leading to a drastic change in injection properties when air flow is present. Such drastic changes in the near atomizer flow field caused by addition of a swirling air flow addition are also reported in the literature: the term "explosive breakup" is coined for this phenomenon by Hopfinger & Lasheras [START_REF] Hopfinger | Explosive breakup of a liquid jet by a swirling coaxial gas jet[END_REF], who studied a round water jet in a swirling air flow. Above a certain critical value of the swirl number, the plain jet underwent a violent radial expansion, forming a hollow cone and dramatically improving atomization. In [START_REF] Muthuselvan | Experimental study of spray breakup phenomena in small-scale simplex atomizers with andwithout air swirl[END_REF], the authors analyze the spray formed by a simplex atomizer in a swirling air flow, and observe a transition similar to that of Fig. 10 from a closed spray to a hollow cone as air flow is increased. They observe that this transition is affected by a strong hysteresis. Thus, a dynamical injection model accounting for coupling with the air flow would be required. This is not attempted in the present simulations and a simpler uncoupled model is adopted in which the initial spray characteristics are adjusted by making use of experimental data obtained in the injector near field. Four main parameters must be determined: the position of fuel injection, the initial droplet size distribution, their velocity and injection angle.

The initial droplets are injected on a ring of inner ra- dius 27 µm and outer radius 54 µm. The droplet size distribution is the one measured at z = 5 mm. This choice is driven by two reasons. First, there is at present no accurate model for the determination of droplet size distribution at the atomizer's exit, and there is no direct optical access to perform in situ experimental measurements. Second, several physical phenomena (secondary atomization, splashing, dropletdroplet interactions) are not accounted for in the simulation because of the lack of fully predictive models. By imposing the experimental distribution, one circumvents these issues, making the assumption that the corresponding physical processes only take place in the immediate vicinity of the atomizer and are characterized by short time scales. Simulations with diameter distributions from measurements at several different radial positions as input were performed, but using a droplet diameter PDF averaged over all experimentally available radial positions gave the best results. The injected distribution is corrected for evaporation between the atomizer and the PDA measurement section assuming that it follows a d 2 law. This correction is important for the smaller droplets, for which this simplistic evaporation law is acceptable as they behave as tracers in an isothermal flow [START_REF] Lefebvre | Atomization and Sprays[END_REF]. For this correction, a constant value of the d 2 law coefficient λ d2 = 1.65 × 10 -8 m 2 s -1 given by [START_REF] Chauveau | An analysis of the droplet support fiber effect on the evaporation process[END_REF] is used. The time of flight of particles from the atomizer to the measurement point is estimated using LES. The histogram of the corrected distribution is shown in Fig. 11. One may observe that droplets initially smaller than approximately 5 µm are fully evaporated and cannot be recovered. In practice, this has a negligible impact as these droplets represent a small fraction of the total injected mass, and this will not bias the comparison between experiment and LES at the PDA measurement section. One side advantage is that the tracking of these droplets would have induced a high CPU cost.

To determine the axial injection velocity by taking advantage of experimental data, it is instructive to examine to travel from the atomizer to the PDA measurement volume, the increase in axial velocity for such a particle is less than 3 m s -1 . Thus, the velocity of the larger droplets is roughly the same as their injection velocity. It is then reasonable to assign an initial value 29 m s -1 to all the initial droplets.

The injection angle of each droplet is randomly chosen by normally distributing its values with a mean θ 0 and a standard deviation δθ. Initial values for these parameters are deduced from the tomographic images of the spray with air flow (Fig. 10). A parametric study of several LES is then conducted by varying θ 0 and δθ around these initial values to find the optimal set of parameters, θ 0 = 37 deg and δθ = 8 deg.

Validation of the simulation

Volumetric cumulative distribution function of droplet diameter are compared between LES and PDA measurements at two locations in Fig. 12. At r = 5 mm, the agreement is excellent. At r = 4 mm, the proportion of large droplets is overpredicted by about 20% by volume. Mean and Sauter mean diameters are compared with PDA measurements in Fig. 13. The agreement is fairly good, especially in the regions characterized by higher liquid volume flux that are of higher practical interest. Some discrepancies are observed in the inner recirculation region at z = 5 mm. This could be expected as the number of droplets in this region is very low, also leading to high statistical uncertainties. Overprediction of d 32 around r = 4 mm are mainly due to the film model and this topic is discussed in Appendix A.

The probabilities of axial velocity conditioned by the droplet diameter are gathered in Fig. 7 for the points located at r= 4, 5 and 6 mm and at 5 mm from the outlet. The bimodality of the velocity distribution for large droplets observed experimentally at r = 5 mm is retrieved albeit overestimated. As is seen experimentally, this bimodality is also not observed at r = 4 mm. The agreement between LES and experimental measurements is also excellent at r = 6 mm. The locations of the velocity peaks are also well retrieved for all positions and for both peaks.

Discussion

Backward-forward Lagrangian particle tracking is employed in Fig. 14 to examine droplet histories from the LES data. Only droplets larger than d > 40 µm are considered. The following process is used. Droplet trajectories are projected to the (rz) plane. All droplets that pass in the vicinity of the control volume (0.5 mm diameter torus around the point of interest at r = 5 mm, z = 5 mm) are then extracted. In this manner, the history of droplets whose behavior is examined in the center image of Fig. 7 can be analyzed.

Figure 14(a) presents the trajectories of all 516 large particles observed to come in the vicinity of (r, z) = (5, 5) mm during 5 ms of simulation, along with the geometry of the injector (in black). Two populations can be distinguished. Some droplets hit the conical convergent wall of the injector, others do not, and those two groups are well separated. The dynamics of 20 particles that do not collide with the wall is examined in more details in Fig. 14(b). Their axial velocity at every point in the trajectory is shown by their color. These particles have a high initial axial momentum communicated at injection. They are slightly accelerated after the outlet of the injector and reach velocities of the order of 30 m s -1 in the vicinity of the point of interest (r, z) = (5, 5) mm. Those particles stay closely together over a long distance.

Figure 14(c) is focused on particles that have formed a film on the convergent section. As they hit the wall at a close to normal angle, their initial momentum is mostly dissipated Additional materials to this paper present experimental measurements that have been performed on the spray formed when using injector 712 at the same operating point. Results and conclusion are similar to those obtained with swirler 707. The same characterization is also presented for swirler 716 at a slightly differing operating point corresponding to Flame Describing Function measurements in [START_REF] Prieur | Experimental determinations of Flame Describing Functions of swirling spray flames[END_REF]. The fuel flow rate is nearly identical to that used for the present characterization of swirlers 707 and 712, but the air flow rate is 15% higher. Similar features, and in particular the bimodality of the velocity distribution are also found. However, the probability of slower velocities is much higher with swirler 716. This may be caused by a wider spray angle which is observed on laser illumination in the near field of the atomizer.

CONCLUSION

This investigation is concerned with issues of technical importance in the engineering of swirling injectors for jet engine applications. The study is carried out on a laboratory scale injector designed to emulate a small helicopter engine injector. This unit features a pressure swirl simplex atomizer creating a hollow cone shaped spray of heptane droplets. The geometry is such that some of the droplets in the spray hit the injector end piece inner surface forming a liquid film. Using an experimental characterization of the spray in the form of polar plots, it is shown that asymmetries in the spray can be identified, and that these features are linked to manufacturing variability or in the present case to a damaged atomizer. Tests with swirlers having different levels of rate of rotation indicate that injector flows characterized by a higher swirl number tend to smooth out irregularities in the spray pattern. A second item of considerable importance in the numerical simulation of liquid spray injectors is that of size and velocity distribution of the droplets. It is shown that suitable choices of the size and velocity distributions imposed in the atomizer near field are needed to obtain a good match between experimental data and results of simulations. The injected distributions are derived by optimizing the input set of parameters determined from laser tomographic visualizations of the spray and Phase Doppler Anemometry. This allows high fidelity LES simulations of the spray. It is shown in particular that an experimentally observed bimodality in the velocity distribution for large droplets is retrieved and that this feature can be linked to the droplet histories. Droplets forming the liquid film are slowed down and form a second, low velocity peak in the distribution. It is worth noting that the bimodal character of the velocity distribution in the spray of aeronautical injection systems was observed some time ago by McDonnell et al. [START_REF] Mcdonnell | Symmetry assessment of an air-blast at-omizer spray[END_REF], and that their interpretation is validated by the present work.

APPENDIX A. SPRAY PREDICTION AT THE INJEC-TOR OUTLET

It is interesting to examine the nearfield of the injector outlet. In this region and more specifically around r = 4 mm, Fig. 13 indicates that the calculated d 32 is over-predicted. This can be understood by examining the backward-forward lagrangian tracking displayed in Fig. 15 for particles pass-ing around r = 4 mm, z = 5 mm corresponding to the injector outlet. Large particles with high Stokes numbers observed at this position originate from the film boundary condition and have joined the film by hitting the small vertical wall found at the end of the conical convergent. These large particles conserve a high level of axial momentum that leads to clustering around r = 4 mm further downstream, while smaller particles will disperse more readily. In the experiment, the 90 • corner is not sharp, but presents a small, non uniform chamfer due to the manufacturing process. If a more realistic, ragged edge geometry were adopted for this corner, a more realistic dispersion of the large particles would be achieved by the film model. A secondary edge atomization model would also be useful in this respect, but would require an improved knowledge of the spray at injection. This clustering effect probably also plays a role in the underestimation of the diameters in the outer flow region.

Fig. 1 :

 1 Fig. 1: From left to right: ( ( (a a a) ) ) Schematic view of the "SICCA-Spray" test rig. The air plenum is not shown (see however [21] for more details). Approximate position of the laser sheet (green), the field of view of the camera (blue dotted square) and the laser beams of the PDA system (red) are shown. ( ( (b b b) ) ) True color photographs of "SICCA-Spray" under reactive conditions. The "M-shaped" flame is visible in blue. A Nd:YaG laser sheet is used to obtain a droplet tomography in two directions, illustrating the hollow cone shape of the spray. ( ( (c c c) ) ) Exploded view of the injector and its components. In purple, the main body, in orange the liquid fuel atomizer, in translucide teal the tangential air swirler and in gold the injector outlet, which is flush mounted in the chamber backplane. Adapted from [22]. ( ( (d d d) ) ) Schematic view of the swirler seen from above indicating the main dimensions of this component. ( ( (e e e) ) ) Shadowgraph visualization of a typical spray near the tip of the atomizer.

  (a) Mean axial velocity. (b) Mean tangential velocity. (c) Mean radial velocity. (d) Mean turbulent velocity.

Fig. 2 :

 2 Fig. 2: Mean air flow velocity profiles measured 5 mm downstream of the dump plane for three injection systems (see Tab. 1).

  (a) Atomizer AD, Swirler 707. (b) Atomizer AD, Swirler 712. (c) Atomizer AD, Swirler 716. (d) Atomizer AG, Swirler 707. (e) Atomizer AG, Swirler 712. (f) Atomizer AG, Swirler 716.

Fig. 3 :

 3 Fig. 3: Circumferential distribution of the Sauter mean droplet diameter d 32 .For each data point the azimuthal angle on the graph corresponds to the azimuthal angle in the experiment, the color of the point to the radial coordinate of the measurement point, and the radial position in the graph corresponds to the value of d 32 (the gray circles correspond to 5, 10, 15... µm as indicated in the graphs when this is possible). Measurements acquired at z = 2.5 mm.

Figure 4 Fig. 4 :

 44 Fig. 4: Radial profiles of mean droplet diameter d 10 and Sauter mean diameter d 32 obtained at z = 5, 10 and 20 mm above the backplane (red dotted lines in Fig. 5). Swirler 707. Error bars indicate the 95% statistical confidence interval obtained from bootstrapping.

Fig. 5 :Fig. 6 :

 56 Fig. 5: Laser tomography of the spray above the dump plane, in the presence of flow. Swirler 707. Red dotted lines: PDA measurements heights.

  (a) r = 4 mm. (b) r = 5 mm. (c) r = 6 mm.

Fig. 7 :

 7 Fig.7: Probability of axial velocity conditioned by the droplet diameter at z = 5 mm and at radial positions r = 4, 5 and 6 mm. On the top row, experimental results measured by PDA. On the bottom row, simulation (LES) results that will be discussed in the LES section. In both cases, circumferential averaging is performed.

Fig. 8 :

 8 Fig. 8: In red, axial velocity histogram for large droplets (diameter greater than 45 µm) at 3 radial positions (from top to bottom, r = 4, 5 and 6 mm, all taken at z = 5 mm). In black, air axial velocity histogram from LDV measurements.

Fig. 9 :

 9 Fig. 9: Comparison between LES velocity profiles and LDV measurements at z = 5 mm. The plots show mean axial and tangential velocities as well as overall RMS velocity.

1 Fig. 10 :

 110 Fig. 10: Tomographic visualization of the spray at the tip of the atomizer. The spray is illuminated by a laser sheet just above the outlet of the atomizer. In this experiment, the end plate of the atomizer (golden disc in Fig. 1(c)) is removed to gain optical access. The laser sheet crosses the atomizer on its centerline, and the spray image is obtained with a long exposure time. Left image : fuel injection without air flow. Right image : with air flow.

Fig. 11 :

 11 Fig. 11: Measured (blue) and injected (brown) droplet distribution. The injected distribution is corrected for evaporation.

Fig. 12 :

 12 Fig. 12: LES and PDA volumetric cumulative distribution function of droplet sizes at r = 4 (top) and r = 5 mm (bottom). z = 5 mm.

Fig. 13 :

 13 Fig. 13: Comparison of mean (black) and Sauter mean (red) diameters from the LES with PDA measurements. Three measurement sections are considered (z = 5, 10 and 20 mm). Azimuthal averages are performed in both cases. On the right axis, the radial distribution of the axial liquid volume flux is shown as a blue solid line.

Fig. 14 :

 14 Fig.14: Backward -forward Lagrangian particle tracking of all droplets that pass inside the control volume (0.5 mm diameter torus around the point of interest at = 5 mm, z = 5 mm). In black, geometry of the upper part of the injector. On the right, trajectories are colored by axial velocity.

Table 1 :

 1 Dimensions and main aerodynamic characteristics of the three swirlers with an air mass flow rate of 2.3 g s -1 . See Fig.1for the definition of dimensions d th and R 0 . Swirler 716 yields a significantly different flow pattern: due to the higher rotation rate, the angle of the air jet is wider, the level of turbulence and the recirculating velocities are higher, and the diameter of the recirculation zone is augmented. Experimental swirl numbers calculated at z = 5 mm using the conventional formula[START_REF] Beer | Combustion aerodynamics[END_REF]

	Swirler	d th	R 0	∆P	S	f PVC (kHz)
		(mm)	(mm)	(kPa)
	707	4	4.6	2.79 0.62	3.95
	712	3	2.28	3.45 0.6	3.44
	716	3.5	4.7	4.4	0.78	3.47

Acknowledgements

This work benefited from the support of SafranTech (contract NF5Z-5100), and of project FASMIC ANR-16-CE22-0013 of the French National Research Agency (ANR). This work was granted access to the HPC resources of CINES under the allocation 2019-A0072B10315 attributed by GENCI. The authors wish to thank J.M. Dupays, E. Lo Schiavo, D. Laera and L.Y.M. Gicquel for their help in the use of AVBP and the film model boundary condition.


Nathalie Revol
email: nathalie.revol@inria.fr

Luis Benet
email: benet@icf.unam.mx

Luca Ferranti
email: luca.ferranti@uwasa.fi

Sergei Zhilin
email: szhilin@gmail.com

Testing interval arithmetic libraries, including their IEEE-1788 compliance

Keywords: Unit tests for interval arithmetic libraries, Test cases for interval arithmetic, Testing IEEE 1788-2015 compliance

As developers of libraries implementing interval arithmetic, we faced the same difficulties when it came to testing our libraries. What must be tested? How can we devise relevant test cases for unit testing? How can we ensure a high (and possibly 100%) test coverage? In this paper we list the different aspects that, in our opinion, must be tested, giving indications on the choice of test cases. Then we examine how several interval arithmetic libraries actually perform tests. Next, we introduce two frameworks developed specifically to gather test cases and to incorporate easily new libraries in order to test them, namely JInterval and ITF1788. Not every important aspects of our libraries fit in these frameworks and we list extra tests that we deem important, but not easy, to perform.

Introduction

Many libraries implement interval arithmetic, from XSC [START_REF] Hofschuster | C-XSC 2.0 -A C++ Library for Extended Scientific Computing[END_REF] and FiLib [START_REF] Lerch | filib++, a Fast Interval Library Supporting Containment Computations[END_REF] among the pioneers, to Octave/Interval [START_REF] Heimlich | Interval arithmetic in GNU Octave[END_REF], ValidatedNumerics.jl [START_REF]ValidatedNumerics[END_REF] and Moore [START_REF] Mascarenhas | Moore: Interval arithmetic in modern C++[END_REF] for the more recent ones. In 2010, a comparison of several libraries [START_REF] Revol | Standardized Interval Arithmetic and Interval Arithmetic Used in Libraries[END_REF] indicated that the underlying approach of the definition of interval arithmetic was, more or less, different and specific to each library. It was thus impossible to compare the results obtained by different libraries, as so many differences impacted the computations. To enable comparisons, it was decided to standardize interval arithmetic. A collective effort, launched in 2008, led to the standardization of interval arithmetic, specified in the IEEE 1788-2015 standard [START_REF]IEEE Standard for Interval Arithmetic[END_REF].

The development of IEEE 1788-2015 compliant libraries, and thus the comparison of their results, was therefore made possible. However, before even considering comparisons between libraries, the developers of an interval arithmetic library need to check their developments, to ensure that their library behaves correctly, in particular with respect to IEEE 1788-2015 compliance. How do they usually proceed? Through intensive tests. Formal proof that each operation is correct is the next step, but it will not be considered in what follows.

Another use of tests is, for users, at installation time, to check that the newly installed library behaves properly on one's platform. This is not a 100% guarantee of correctness, but it increases the user's confidence.

Developers test for many aspects of their libraries, not only for IEEE 1788-2015 compliance, but also for instance the handling of exceptional or invalid inputs and so on. Unit tests are very common: for each operation or function, one prepares a set of inputs and checks whether the library computes the expected output, or at least an enclosure of it. Rapidly, the developers realize that, in order to get sufficient confidence in the implementation of a given operation, in order to test a high enough coverage of the code (or 100% of it), a huge number of cases must be tested. This implies that, for each implemented operation, a large part of the development time is devoted to devising these test cases. We are developers of different libraries, we faced the same situation and we elaborated similar tests. We decided to put our expertise in common and to share it. The preliminary result of this collaboration is the content of this article.

In what follows, we propose to list, in Section 2, the different features that must be tested. We then compare, in Section 3, this list of recommended tests with the tests that are actually included in some well-known libraries (for which we have access to their test suites). Our goal is to devise test cases for each situation we have identified and to share them: in Section 4 we introduce two existing frameworks that, on the one hand, ease the integration of a new library to be tested, and that, on the other hand, offer a significant set of test cases, as well as a convenient mechanism to add new ones. As no setup is perfect, in Section 5, we then give a list of other types of tests deemed useful in our development of interval arithmetic libraries, but still missing. Our short-to midterm-goal is to enrich and share test suites, and enhancing the capabilities of the frameworks introduced in Section 4.

What must be tested?

Testing floating-point arithmetic has been studied thoroughly, as testified by Nelson Beebe's web page for software [START_REF] Beebe | IEEE 754 floating-point test software[END_REF] and more generally by the huge amount of bibliographic references and of programs available. Similarly, in this article, we gather and organize the important features that must be tested in order to assess the correctness, compliance and possibly quality of a library implementing interval arithmetic.

General remarks about unit tests

Let us denote by L the tested library, and intervals using boldface: x, y. We focus here on unit tests, for a function, denoted by f . A unit test case is a pair composed of the input argument x and the expected output y.

First, the output y must be the tightest representable interval enclosing f (x); otherwise a very accurate library could compute z such that z y and still z ⊃ f (x). To ensure the tightness property for y, typically, one computes the endpoints y and ȳ of y = [y, ȳ] using a precision higher than the computing precision of L. Let us illustrate this with L using Binary64 floating-point numbers, on the (simple) example of the exponential function. We assume that the given floating-point implementation of exp does not provide correct rounding, but that, for any precision q, it satisfies RD q (exp(x)) ≤ exp(x) ≤ RU q (exp(x)) for RD q rounding downwards and RU q rounding upwards, both in precision q. Given x = [x, x], to compute the infimum y of y, one gets the approximations RD q (exp(x)) and RU q (exp(x)) of exp(x) in high precision q, and finally round them downwards in the target precision = 53. If RD 53 (RD q (exp(x))) and RD 53 (RU q (exp(x))) are equal, then they are the sought value y for the infimum of exp(x).

Once the test cases are devised, what conclusions can be drawn from the comparison between y and z the result computed by L? Requiring equality may be too demanding. Inclusion is required, but we want to dismiss, for instance, an implementation of sin that returns [-1, 1] for any argument. If the accuracy of the function is given in the specification of the library, how can it be used and checked?

A last general recommendation is to incorporate the following procedure: if (x, y) is given as a pair of input and output, pick at random (a reasonably large number of) values x ∈ x and check whether f (x) ∈ y with f (x) computed by the underlying arithmetic, and whether f ([x, x]) ⊂ y. Those tests are disconnected from the knowledge of the implementation of f and may hit a zone not considered (forgotten) in the development of f , such as an overlooked quadrant for a trigonometric function.

Tests common to all interval arithmetic libraries

-Easy test cases. First, easy cases are tested: these test cases are both easy to devise and easy to compute. They constitute preliminary tests, at an early stage of the implementation of the function, to identify and correct the most obvious bugs. These test cases are chosen to cover, roughly speaking, the various possible magnitudes of the arguments, but without any exhaustivity. They also span positive and negative values, if the domain and range of the function so permit. -Special and exceptional values. A time-consuming task, when implementing a function, is the handling of special and exceptional values. In floating-point arithmetic, 0 is a special value, because it has two representations, +0 and -0, that is, signed zeros. It is thus valuable to test whether one gets the same result, independently of the sign of the representation of 0, when 0 is an endpoint of either an input argument or of the output. 1 is not a special value for floating-point arithmetic. However, it is often considered as a special value, for instance for the logarithmic functions and the arc-cosine: its image is the special value 0.

Infinities are also special values in floating-point arithmetic. If they are supported by the library, one must have test cases containing infinities as their endpoints, both for inputs and for outputs. If infinities are not supported by the library, it is even the more so important to test input arguments that yield overflow, such as the addition or the multiplication of

[0, MAXREAL] by [0, MAXREAL] or [0.75 * MAXREAL, MAXREAL] by [0.75 * MAXREAL, MAXREAL]. Another example is exp[x, x]
) with x such that exp(x) rather small and x such that exp(x) overflows, or both x and x such that exp(x) overflow. Again, when floating-point arithmetic is used, NaN is an exceptional value, that is the result of an invalid computation. How is a NaN handled when it occurs as the endpoint of an input interval? The answer should be that, even if it has no mathematical meaning, in practice NaNs propagate. Some libraries use NaN to denote the emptyset, as is the case in Intlab [START_REF] Rump | INTLAB -INTerval LABoratory[END_REF]. Care must then be taken for the union or intersection of two intervals.

-Cornercases. Another family of test cases is designed in order to test the difficulties encountered during the implementation, such as difficult-to-round values, or values at the boundaries of the domain of the function, or close to points whose image is an extremum. A typical example is a floating-point input close to an integer multiple of π 2 , for trigonometric functions. The various tools for testing floating-point arithmetic can be a source of inspiration. Indeed, with floating-point arithmetic, as there is a sudden change in absolute error between two consecutive floating-point numbers when the binade changes, it is worth exploring several binades, both for the inputs and for the output values. This is particularly true for mid-rad representation of the intervals. With floating-point arithmetic, another delicate zone concerns subnormal numbers. Test cases should contain intervals [0, x] and [x, x] such that x and x are subnormals -for the logarithmic functions, or intervals [x, x] such that their exponential has one or two subnormal endpoints.

-Functions specific to intervals. Functions such as the union (convex hull of the union) or the intersection must be tested when one, or both endpoints are special or exceptional values, or when one argument is the empty set, if it is supported by the library. Functions such as the midpoint or radius have been thoroughly studied in [START_REF] Goualard | How do you compute the midpoint of an interval?[END_REF], test cases can be taken from this reference. -Input and Output. Finally, I/O is the place where the most unexpected things can happen: much creativity is needed to cover a large variety of input values given as strings.

Tests about IEEE-1788 compliance

The IEEE 1788-2015 standard [START_REF]IEEE Standard for Interval Arithmetic[END_REF] was developed to enable comparisons of interval methods and their implementation. It mandates operations in Section 9; it provides "hooks" to integrate several flavors of interval arithmetic and it defines the so-called "set-based" flavor along with its set of recommended operations in Section 10; it mandates that intervals are decorated with flags that sum up the history of their computation (such as "every operation and function involved in this computation is defined and continuous on its arguments") in Section 8. It also specifies the handling of exceptions in Section 12.1. -testing the different accuracies: tightness/accurate/valid: we already discussed the difficulty to validate a computed result, to compare it with the expected output, when it is not the tightest possible one. The standard defines three possible levels of accuracy. Test cases must include all three possibilities for each pair of input and output, and check if the corresponding levels are available in the library. They must include the description of this level, the testing framework must be able to check whether the result is tight or valid (both are easy), or accurate (by computing the enclosing interval that is only slightly larger, as defined by the standard). -exceptional behavior: exceptions have already been mentioned above. The standard defines the different exceptions that must be signalled, the test cases must include these exceptions and check whether the correct exception is signalled. -standard specific I/O system: the standard defines precisely the various forms of inputs and outputs. At least one test case per form must be present in the test case. If one wants to test every possibility, with at most 3 values for each variable field lengths, one obtains a thousand test cases.

Tests specific to some libraries

Some libraries must need specific test cases, corresponding to their specificities.

-Representations of intervals. For instance, a library that uses the midrad representation of intervals must test that operations and I/O operate correctly on this representation. There are libraries that keep track of the openness or closedness of each endpoint: test cases must be provided to check that this information is correctly computed. If each possibility is tested, this can multiply by 4 the number of test cases. Some libraries do not employ floating-point types to represent intervals, either in inf-sup or mid-rad representation. Exact rational numbers can be used instead, as in MPRIA [START_REF]MPRIA: GNU Multi-Precision Rational Interval Arithmetic Library[END_REF] and JInterval [START_REF] Nadezhin | JInterval Library: Principles, Development, and Perspectives[END_REF]. Other libraries also allow complex intervals, that is, they use complex arithmetic as their underlying arithmetic. -Precision. A library that accomodates several floating-point formats, such as Binary32 and Binary64, must test that it operates soundly with both formats and with combination of both. Moreover, a library that uses arbitrary precision, such as MPFI [START_REF] Revol | Motivations for an Arbitrary Precision Interval Arithmetic and the MPFI Library[END_REF] or ARB [START_REF] Johansson | Arb: a C library for ball arithmetic[END_REF] must include tests for largely varying precisions, without claiming exhaustivity. -Directed roundings. In order to guarantee the inclusion property, libraries rely on directed roundings. They must have tests which are specific to how directed rounding is achieved. If it is achieved by changing the rounding mode, ideally it should be tested that this is thread safe and that using multiple threads doesn't lead to undefined behavior. If it is achieved via software, by emulating directed rounding in round to nearest (as it is done e.g. by default in IntervalArithmetic.jl), then they very likely rely heavily on the use of Error-Free Transforms, or EFT in short. Thus, hard to round cases and corner cases for EFT (when it overflows, underflows, returns NaN) should also be tested.

On tests timing

An important issue related to tests is time. To illustrate this question, let us mention that MPFI currently implements about 30 test cases for each function, and running the tests (make check) takes 20 seconds of user's time on a reasonably fast PC. If all these recommended test cases were implemented, there would be at least a few hundred test cases, at most a few thousands. The testing time will be multiplied according by a factor between 10 and 100. Shall we stick to our limited number of test cases and risk to have not enough tests performed, or shall we risk to have a more complete test cases coverage and no test performed at all, because of the time it takes?

Need for a unified framework for testing interval arithmetic libraries

While there are several shared features between all interval arithmetic libraries, they are all also somehow unique in design choices (how to handle direct rounding, tightness vs speed, etc.). This call for a unified framework to test interval arithmetic libraries. This has several benefits. First, a standard framework makes it easier to compare libraries against each other. It would also offer a tool for developers, to easily verify whether their library is compliant with the standard or not. Finally, an important use case of interval arithmetic is to perform rigorous computations, that can be used as mathematical proofs. While several theorems for rigorous computations have been proved, they all rely on the assumption that the underlying implementation of interval arithmetic is correct. For this reason, a unified testing framework, developed through inputs from several interval libraries developers, would increase the reliability of all libraries. Based on the above discussion, we identify the following requirements in a unified interval testing framework: an enclosure of the expected result. A unified framework should allow the user to switch between testing tightness and validity, ideally also allowing to set a threshold on how much wider the returned interval is allowed to be. 6. Self-validated: it should be possible for the user to verify that the results expected by the tests are indeed correct.

Interval arithmetic libraries and their test sets

During decades, several interval arithmetic libraries have been developed, each with its own design choices. Since we are interesting in testing, here we only list some of the ones for which tests are freely available.

1. MPFI [START_REF] Revol | Motivations for an Arbitrary Precision Interval Arithmetic and the MPFI Library[END_REF]: C library for arbitrary precision interval arithmetic. Based on MPFR. Unit tests for each function: some "easy" values, exceptions, exact cases. Currently not IEEE-1788 compliant.

2. JInterval [START_REF] Nadezhin | JInterval Library: Principles, Development, and Perspectives[END_REF]18]: Java library providing preliminary standard IEEE P1788 compliant exact implementation of interval arithmetic operations and controllable arbitrary precision for elementary functions. Use the set of rational numbers extended by -∞ and +∞ as a basic number type for intervals endpoints. On top of this implementation, tightest approximations of interval arithmetic and elementary functions are supported for float types binary16, binary32, binary64, binary128, etc.

3. libieeep1788 [START_REF] Nehmeier | libieeep1788: A C++ Implementation of the IEEE interval standard P1788[END_REF]21]: C++ implementation of the preliminary IEEE P1788 standard. Relies on MPFR. Unit tests for IEEE-1788 compliance.

4.

Octave interval package [START_REF] Heimlich | Interval arithmetic in GNU Octave[END_REF][START_REF] Heimlich | GNU Octave interval package[END_REF]: GNU Octave interval arithmetic package compliant with the standard. Test sets derived from libieeep1788, MPFI, C-XSC, FILIB plus several new tests for most of the standard functions. These tests constitute the current ITF1788 testsuite. Relies on MPFR for both arithmetic operations and elementary functions.

5.

ValidatedNumerics.jl [START_REF] Sanders | [END_REF][START_REF]ValidatedNumerics[END_REF]: Julia set of packages for interval arithmetic and applications. Currently not conformant with the standard. By default, direct rounding is handled via the software emulator RoundingEmulator.jl. However, it also supports changing rounding modes. Correctly rounded elementary functions are computed via CRlibm when possible, and using MPFR as fallback. For Binary64 bounds, uses the ITF1788 testsuite. Also has hardcoded tests for other non-standard functionalities (complex intervals, interval boxes, other number formats) 6. GAOL [START_REF] Goualard | Gaol: Not just another interval library[END_REF]: C++ library not compliant with the standard, e.g. lacks decorated intervals. Tight bounds on arithmetic operations using round upwards mode only. Supports CRlibm and IBM Accurate Portable Mathematical Library for elementary functions computations. Uses its own unit tests.

Existing frameworks

4.1 Why they are needed?

As a motivating example, consider Kuliamin's paper [START_REF] Kuliamin | Standardization and Testing of Implementations of Mathematical Functions in Floating Point Numbers[END_REF] about testing floatingpoint mathematical functions. The paper defines several criteria for devising tests, such as considering exceptional cases (overflow, underflow, NaN), inputs out of domain, special values (signed zero, subnormal numbers) etc. However, this leads to more than 15,000 values to test for the exponential function alone. When interval arithmetic is tested, this means that each endpoint must take 15,000 values: the combinatorics (even if one restricts the tests to proper intervals, where the left endpoint is not greater than the right endpoint) are in defavor of developing individually, manually, such intensive tests for each new library.

A natural way to increase effectiveness of individual libraries developers is to unite their efforts in order to accumulate tests in a properly structured common shared database of tests equipped by appropriate tools for tests usage. That is why the frameworks briefly described in the next two subsections are worth our attention.

JInterval P1788 Test Launcher

P1788 Test Launcher [START_REF] Nadezhin | JInterval Library: Principles, Development, and Perspectives[END_REF]19] is an interval arithmetic libraries testing application based on the JInterval library. The features of JInterval entail a freedom for P1788 Test Launcher in the selection of computation modes corresponding to a tested library.

P1788 Test Launcher loads dynamic libraries (.so/.dll) with third party implementation of P1788 and checks the results obtained from a library with the tightest results computed internally using JInterval. Operations and functions from a tested library are called through the unified adaptor wrapper interface which must be preliminary implemented for the library being checked. The structure of the wrapper interface is adopted from C++ templates of libieeep1788 library. So, technically, the implementation of the interface in most of cases boils down to overloading methods of libieeep1788 classes describing setbased flavor of interval arithmetic for Binary64 value set as a rule. Wrappers for the following interval arithmetic libraries are already included in P1788 Test Launcher: Profil/BIAS, boost/interval, C-XSC, Filib, libieeep1788, libMoore, MPFI.

The launcher reads tests from plain text files of simple human-readable format and writes the results computed using tested library and JInterval and their relation into a plain text report. Test set included to the Launcher consists of over 14,000 tests which partly originated from libieeep1788, Filib, libMoore while the others are original.

ITF1788 -Interval Test Framework for IEEE 1788

The idea behind ITF1788 is to simplify the development of unit tests for an interval arithmetic library independently of the programming language it is written in, the testing framework used and library peculiarities, such as custom function names. ITF1788 is a Python engine for converting pre-calculated tests from domain specific language, called the Interval Test Language (ITL), to the code of unit tests according to the configuration describing the syntax of the language, the testing framework and the interval library.

The notation of ITL is easy to read and write and covers all notions of IEEE-1788 necessary to test an implementation. The following small example of test description from [START_REF] Kiesner | ITF1788: An Interval Testframework for IEEE 1788[END_REF] gives a general understanding of ITL composition. One can easily customize ITF1788 to support new programming languages, test frameworks and library specific features. The developer simply has to prepare several YAML configuration files. Additional flexibility can be obtained using implementation of Python callback functions for proper modifications of ITF1788 output.

Original ITF1788 engine and accompanying test sets for IEEE-1788 compliance testing have been developed by M. Kiesner, M. Nehmeier and J. Wolff von Gudenberg [START_REF] Kiesner | ITF1788: An Interval Testframework for IEEE 1788[END_REF]

What is missing? A roadmap for testing

A first important requirement, which is perhaps too obvious but still worth mentioning, is that a standard framework for tests must be open. It should also be portable across architectures and languages, and be properly documented. The ITF1788 test suite has developed an interesting approach to solve these particular issues. This is achieved by introducing a domain specific language for storing the unit tests, and allowing the user to pass the syntactic information necessary of the language through the configuration of some files [START_REF] Kiesner | ITF1788: An Interval Testframework for IEEE 1788[END_REF].

The unit tests of the ITF1788 suite cover numerous cases, including some corner cases, allow modularity and allow to test the tightness (equality assertion) or correctness (asserting containment). Both of these aspects are desirable and should be included with certain redundancy. For instance, to avoid that an implementation of sin trivially returns [-1, 1] and passes all containment assertions, tests that involve the radius of the returned interval (to be less than two) can be introduced that break the trivial implementation.

Yet, the unit tests of the ITF1788 suite are restricted to the binary64 floating point format. In our view, it is desirable to have the unit test coverage enlarged to include unit tests for quadruple and extended precision, including different precision values (e.g., 256, 512, 1024). Similarly, concrete tests should be included for other IEEE754 numeric formats, such as binary32.

Other unit tests should choose randomly (floating-point) values x ∈ x and evaluate that the inclusion f ([x, x]) ⊆ f (x) holds true. While these tests are certainly redundant, in particular if the implementation is correct, they serve to check values not explicitly covered in previous tests, in particular worst-case accuracy. Such tests may be particularly interesting during the early stages of development, considering specially the power function and the transcendental functions, in particular for extended precision intervals.

In other communities, also a priori interested in reliable computations, "friendly competitions" are organized, to compare what each tool offers: methods, limitations, new aspects, etc. The code of the tools must be publicly available, and the benchmarks must be reproducible. See for instance https://cps-vo.org/group/ARCH/FriendlyCompetition. Organizing such a friendly competition would be a sane motivation for developers, to update and increase the depth of unit tests.

Conclusions

This paper focused on testing strategies for interval arithmetic libraries. During the years, different libraries have adopted different design choices and as such comparison and reproducibility has been a challenge. The IEEE 1788-2015 standardization is a big step to solve this issues, but broad adoption of the standard is a slow process. Moreover, testing is critical in software development and, due to its rigorous nature, even more central in interval arithmetic. It is thus fundamental to have an extensive unified testing framework. In this work, we sketched a possible structure for such framework. Hopefully, developers of interval libraries will join forces and share experience to develop an extensive open-source cross platform testing framework.

/

  * Testing the addition function */ testcase addition.test { add [ -1.0, 1.0 ] [ empty ] = [ empty ]; add [ 1.0, 2.0 ] [ 3.0, infinity ] = [ 4.0, infinity ]; add [ 1.0, infinity ] [ -infinity, 4.0 ] = [ entire ]; // using hexadecimal notation add [ 0X1.FFFFFFFFFFFFP+0 ] [ 0X1.999999999999AP-4 ] = [ 0X1.0CCCCCCCCCCC4P+1, 0X1.0CCCCCCCCCCC5P+1 ]; } /* Testing the division function */ testcase division.test { div [ empty ] [ empty ] = [ empty ]; div [ -30.0, 15.0 ] [ entire ] = [ entire ]; } ITF1788 inputs ITL-files and converts tests to the code of unit tests for specified language, test framework and library.

  3 and the grammar for I/O in Section 13. Last but not least, it offers the possibility to return either the tightest possible result, or simply an accurate one, or a valid one in Section 12.10.1. "Compressed arithmetic" will not be mentioned, as it is not yet widely used. Testing IEEE 1788-2015 compliance means testing all of these aspects:-testing flavor compliance: each test case must be accompanied with the indication of the corresponding flavor: how can it be specified? If several flavors are tested (once they are defined by a new revision of the standard), either the testing framework or the library must offer a mechanism to change the current flavor in use.This is needed to test for instance ∅ +[START_REF] Hofschuster | C-XSC 2.0 -A C++ Library for Extended Scientific Computing[END_REF][START_REF] Lerch | filib++, a Fast Interval Library Supporting Containment Computations[END_REF] or

exp([-∞, 0]) in the set-based flavor and then

[START_REF] Heimlich | Interval arithmetic in GNU Octave[END_REF][START_REF] Lerch | filib++, a Fast Interval Library Supporting Containment Computations[END_REF]

+

[START_REF]ValidatedNumerics[END_REF][START_REF] Hofschuster | C-XSC 2.0 -A C++ Library for Extended Scientific Computing[END_REF]

in the (future?) Kaucher flavor.

-required functions / recommended functions: the standard mandates that some functions are implemented by the library, such as cosh, but only recommends other ones, such as expm1 for the set-based flavor. The tests must allow for executing functions if they exist in the library but not crashing when a recommended function is not available. -testing decorations: the set-based flavor defines 5 decorations. To test every possible meaningful combination of decorations for inputs and output, one must devise close to 25 test cases for unary operations and up to 125 test cases for binary operations. Care must be taken to exhaust the complete list of possibilities.

1 .

 1 Modularity: tests should be structured in a modular way, to allow users and developers choosing what parts to test. For example, tests for decorated intervals should be separated by tests for bare intervals. Test for recommended functions should be separated by tests for required functions. 2. Completeness: there should be tests for all functions required and recommended by the standards. The tests should contain both some simple normal cases and cornercases, such as exceptional behavior (handle NaN, overflow, underflow) and hard to round cases for finite precision data types. 3. Support for different number systems: allow to test when the bounds of the interval are Binary64, Binary32, arbitrary precision floats etc. 4. Test for tightness: verify that the computed interval is the smallest valid interval. 5. Test for validity: achieving tight bounds can be challenging or at least computationally expensive. Libraries at early stages or libraries with a focus on efficiency may prefer to test only that the computed interval is valid, i.e.

  . Later O. Heimlich contributed a lot to the original project and to its own fork of ITF1788[27]. The latter now contains:-configurations for programming languages C++, Octave, Julia, Python; -configurations for test frameworks Boost Test Library, test frameworks for Octave, Julia and Python);

	-almost 10,000 tests in ITL.
	Most of the tests accumulated in ITF1788 are derived from unit tests of C-XSC,
	FILIB, MPFI, libieeep1788.

-plugins for libraries:

• libieeep1788, C-XSC, GAOL, Ibex (in C/C++); • interval package, INTLAB toolbox (in GNU Octave); • pyIbex (in Python); • ValidatedNumerics (in Julia).