Error analysis of supersonic air-to-air ejector schlieren pictures

The scope of this article is focused on general analysis of errors and uncertainties possibly arising from CFD-to-schlieren pictures matching. Analysis is based on classic analytical equations. These are firstly evaluated with the presumption of constant density gradient along the ray course. In other words, the deflection of light-ray caused by density gradient is negligible in compare to the cross size of constant gradient area. It is the aim of this work to determine, whether this presumption is applicable in case of supersonic air-to-air ejector. The colour and black and white schlieren pictures are carried out and compared to CFD results. Simulations had covered various eddy viscosities. Computed pressure gradients are transformed into deflection angles and further to ray displacement. Resulting computed lightray deflection is matched to experimental results.


Introduction
The basic idea of the Schlieren or knife-edge method for indicating small optical disturbances was described by Topler some eighty years ago, and thirty years later Rayleigh investigated the optical principles involved.Nowadays the Schlieren method has become a standard tool for engineers in high-speed aerodynamics.Unfortunately it involves some uncertainty which requires an analysis.This analysis is usually being overcome or neglected.Once the experiment is completed, schlieren pictures can be used for flow analysis itself or for validation of numerically computed flow field.Schlieren pictures provide relatively cheap information about flow field and so the CFD does.The matching of its results is than the logical step in applied research when moving from experimental field into CFD base.The validation of CFD is usually done in conjunction with other experiments which allows some level of quantification.

Analytical solution
The refractive index n of a gas was found by Gladstone and Dale to be related to its density  by the law defined as 1 . (1) So, given the refraction index, the density may be calculated.This property forms the basis of the schlieren and other optical methods of analysis of gaseous flow.Let's suppose a Cartesian co-ordinate system, the zdirection being the axis of the schlieren system in the disturbance region, and let the refractive index at any point be defined by , , =K.Take an incident ray parallel to the z-axis entering this region.Its direction cosines will be (0, 0, 1) and after passing through the system which contains disturbance, the incident ray is refracted with x-coordinate angle given as follows . (2) In y-coordinate, the situation is similar and refraction is defined by equation .
More detailed view about equations ( 2) and (3) may be found in ref. [1].In this article, we try to match the CFD simulations to experimental schlieren data.Additional equations are applied to post process computed data with the aim to simulate behaviour of light ray in numerically predicted flow field.To manage this issue, we need to perform some data handling.For the pressure based solver with presumption of ideal gas, the pressure corresponds with density via equation of state .(4) With the presumption of sufficiently fine step we can define density gradient as follows

,
Where the individual gas constant and is local temperature.From CFD simulations, we are able to enumerate a coordinate-oriented pressure derivation , , .With relation given as follows , (6) we can define coordinate oriented density variation as .(7) Refraction coefficient variation is linked with density by relation , (8) where is refraction index found for air with density .For this index we can use Edléns relation defined in [2].From equations (5-8) we can derive differential form of refraction index relation as follows .
In previous text, it has been found that light-ray passing through the disturbance region in z-coordinate direction is related to the refraction index by relations (2, 3) In these equations, has meaning of test section depth.In aerodynamic experiment the depth is usually required to be sufficiently large to eliminate the influence of third (z) dimension as much as possible.In contrast, the best way to capture the density variation with the equations (2, 3) is to minimize the L. It is because the large L allows the light-rays declined at very enter of test section to enter the region with different or / .Resulting declination than corresponds to curvilinear integral of (2, 3) along the light-ray path within test section instead to local value required.The scheme of this problem is shown in the figure 1.The light-ray is oncoming from left to test section.Once it enters the region with density gradient the ray is distorted towards higher density region.If the distortion angle and size of test section (L) is sufficiently large, then the ray enters to region with different density gradient and the ray is distorted in new angle .Outgoing ray angle than comes from the last diffraction within test section and from distortion on test section glass window.The glass influence may be neglected if the oncoming the glass is sufficiently small.To understand, whether this influence may be neglected and whether the presumption of constant density gradient along the ray path is feasible, the method of postprocessing CFD data was developed.This method is based on equations (2, 3) which is reformulated to the form .
Resulting ray angle in y-coordinate comes from last distortion within test section, for which , where If we accept the presumption of the constant density or pressure gradient within the region limited in zcoordinate, given as ∆ , then we can find from equation where Once the equation ( 12) defines the ray-position derivative, the actual position may be found by its integration what gives and All values of are needed to compute the increase of x and y coordinate of the ray regarding to deflection in interval by mean of equation and These equations define actual coordinate, for density gradient evaluation.

Experiment
Experimental research is performed by supersonic air-to air ejector.More detailed view can be taken in reference [3][4][5][6].There is a complex aerodynamic including subsonic, transonic and supersonic flows within the ejector.See figure. 2 to clearly understand the operating of the two-dimensional ejector.The third dimension, not included in the figure 2, is the ejector's depth (z-coordinate) set to 80 mm.Primary air compressed to , which exceeds the critical isentropic pressure ratio in the throat defined as flows from the air supply to the primary nozzle in the middle left of the figure 2. Aerodynamic throat appears near the geometrical throat of the primary channel.Downstream, supersonic air velocity is induced in the diverging part of the nozzle.Isentropic Design Mach number given by section rate 2 was set to , 2.197 [7].After the supersonic air stream passes the trailing edge, viscosity and turbulence effects overwhelm secondary stream and the secondary flow is induced from the surrounding.As soon as the mass flow rate of the overwhelmed air gets big enough, the throat appears and secondary flow become supersonic downstream the throat.Design Mach number of secondary flow is , 1.609.Velocity of the primary and secondary flows are given only by ejector geometry and vary one to the other than.Shear layer divides the regions of various flow speeds in vertical direction.Because of non-zero trailing edge angle, both streams must adapt stream in respect to the new direction.This course-adaption happens across oblique shock waves [8], being inclined by the angle defined by downstream Mach number as Once the oblique shock wave is being induced it spreads towards the symmetry plane and ejector wall.Near the symmetry plane, the interaction of two oblique shocks takes place.In the near-wall region the interaction with boundary layer causes separation of subsonic boundary layer and deflected reflection of the shock might by observed.Afterwards, wall-reflected-shock crosses the region of variable flow velocity called shear layer.Within this interaction originally straight shock is s-like deformed as shown in the reference [5], [6].These interactions repeat many times thanks to the multiple reflections.Shock spreads downstream, than.As every shock causes the decrease of the stream Mach number, every subsequent shock weakens downstream until faded of.More detailed theoretic about supersonic flows may be found in the reference [7].All mentioned shock waves, their interactions, shear and boundary layers comes with large gradients of flow variable.These gradients spread only on very limited distance.Cross sectional size of shear layer derived from schlieren pictures and numerically established are shown in the figure 4.Two types of experimental setup were used.The classical schlieren gap and knife to take B/W pictures and to establish directional deflection of observed issue and the circular four quadrant color filter to colorize the segment, where the light ray were bent to.The ray deflected direction is given by the direction of density gradient.In our experiment we used the filter colorized to the blue (upper left), green (upper right), red (lower left), orange (lower right), seen in the direction of oncoming light rays.There is clear un-colorized spot in the middle of filter.This spot enables the unbent beams to remain uncolored.Electric spot flash-light is used to produce sufficient light intensity.With that light intensity, we were able to set exposure time of capturing DLSR camera to 1/1000 s.This enables take Schlieren pictures clear and bright enough.

CFD simulations
The and directions change shock-like across the wave.In conjunction with relative macro-scale of finite elements, the governing equations defined for FVM should not be valid in the near-shock regions.Further details may be found e.g. in [9].The strongest form of the oblique shock is the normal shock.Even though flow direction remains the same across, state variables and Mach number rapidly changes from supersonic to subsonic.Although theoretical bases of normal shock and its position in channel are well established in e.g.[7], experimental and numerical investigation remains uncompleted.From these reasons, we decided to investigate only the design regime of the ejector.The commercial code FLUENT was used to solve 2D RANS equations by finite volume technique.For all equations, convective terms are discretized using a second-order upwind scheme; inviscid fluxes are derived using a order flux splitting achieving the necessary upwinding and dissipation close to shocks.Diffusion terms are always cast into a central difference form.The discretized system is solved in a coupled way.We used k-omega turbulence model introduced by Wilcox in [10] an laminar model to complete the system of governing equations .The standard k-omega model has a very strong sensitivity to the free stream conditions.
That is why we used the shear stress transport (SST) modification of the k-omega model provided by Menter in [11].SST overcomes this deficiency because it is based on a blending approach.Indeed, to archive the different desired features, the standard k-epsilon model turned into k-omega formulation, is used from the wake region and outside; while the original model is used in the near wall region.To perform these features, a blending function F1 is designed to be one inside the boundary layer and to change gradually to zero in the wake region, if present.The resulting equation for  is as follows More detailed view and expressions of all terms in Equation ( 9) is provided by Menter in [13].The criterion for assessing the convergence was based on the root mean square of the density residues expressed by Where M is the number of grid points and  is the variable considered to check (mass, energy, momentum, etc.).Generally, computations are stopped when residuals fall below 1x10-6 and when the solution is was no longer changing.In addition, at convergence, the mass imbalance is checked on each inlet and outlet boundaries.All mentioned requirements were successfully meet in 50 000 iteration steps.Courant-Friedrichs-Lewy (CFL) number set to 1. Boundary conditions were set to pressure inlet/outlet with non-reflecting boundary condition enabled.Turbulent intensity at the primary nozzle inlet was set to 3% and gauge total pressure to p 01 =300 000 Pa. Secondary outlet is set to 98 000 Pa, what is the atmospheric pressure it experimetal lab.Outlet gauge pressure was set to p E =-70 000 Pa.These conditions applies for design regime.Symmetric computational domain consists of 650 000 quadrilateral elements with original mesh size 0.2 mm.Gradual boundary layer in ten levels from the size of 0.005 mm in wall adjacent cell to free stream value was modeled to keep the wall y + 1.

Results
To get knowledge about the absolute level of pressure, a pneumatic measurement was performed by the fast piezotransducers.Static pressure was measured at mixing chamber wall.Results are shown in the figure 5.         Experimental setup enables to establish the directional declination by relations tan and tan , where , are positions of schlier knife, where observed distubance is just beeing cuted out and f is the focal lenght of shlier objective.On our euipment Zeis 80 the f=1200 mm.Experimentally observed position of knife was established to 10.2 for shear layer ,which gives 0.5°.In both CFD simulations, laminar and k-ω SST, the declinatin can be found about 5°.Situation with x-orinted declination is quite similar, altought the is not any strictly vertical effect.The declination of shock waves moves within the interval ± 5°.Ecperimentally observed declination moves within ± 1°.During the run of experiments, the regime with strong closing shock within mixing chamber was captured.Unfortunately this regime is not captured numerically.This regime is shown in the figure.17

Conclusions
The numerical and experimental research of air flow in the supersonic air-to -air ejector was performed.
Experimental research was based on colorized and b/w schlieren pictures respectively.Pneumatic measurement of static pressure at mixing chamber wall gives general overview of pressure level and enables to set sufficient boundary conditions for CFD simulations.Position of shear layer between primary and secondary streams was observed and its position is matched to those numerically predicted by various eddy viscosity models.Shear layer position is derived from schlieren pictures scaled 20:1.Many unwanted phenomenon may influence the captured picture in such scale.For these reasons the risk analysis was performed by analytical equations.It was shown that both, k ω SST and laminar predicts ray declination in range about 5° for shear layer.Experimentally investigated declination angle moves about 1°.It is very likely that with incresed mesh density, the gradient would steepen even more.From these reasons we can assume that light ray distorted in regions with steepest density gradient experience such strong deflection that it enters the region of different density gradient.Hence, captured schlieren picture may get some inadequate phenomenon compared to infinitesimally thin situation.Thus, in consequent work, the analysis of ray positions will be established step-by-step with difference ∆l in zcoordinate direction.The variable density gradient should be captured, then.

Fig. 1 .
Fig. 1.Scheme of schlieren light-ray distortion within test section and its possible discretization by refraction.

Fig. 4 .
Fig. 4. Distribution of static pressure along mixing chamber wall.Simulation and experiment results.Experimental results are matched to numerical simulations performed with conformity to chapter 4.The position of shear layer measured from 20:1 scaled schlieren picture is shown in the figure5.Experiment is matched to numerical simulations with various eddy viscosity models.

Fig. 5 .
Fig. 5. Position of shear layer in supersonic ejector.Simulation and experiment results.Deflection of light-ray is computed at the base of equations (2, 3).There are general overview and scaled plots of deflection angle for k ω SST and laminar model in the figures 6, 7 and 8, 9 respectively.The x and y coordinate of distorted ray defines displacement relative to origal coordinate at the exit from test section (entering to second glass).These coordinates are computed at the base of equations (2), (3) and (14), (15).Results for k ω SST and laminar model respectively are shown in the figures 10-13.Full scale and croppe are included.

Fig. 14 .
Fig. 14.Black and white schlieren picture with knife in vertical position 10.2 -design regime.

Fig. 15 .
Fig. 15.Black and white schlieren picture with knife in horizontal position 5.5 -design regime.