Simulating schlieren and shadowgraph images from LES data

Geometrical optics ray-tracing is used to derive schlieren and shadowgraph images from large-eddy simulation (LES) data of a jet in supersonic crossflow and to compare with experimental data. Including the components of the optical system that forms the image in the simulation is found to be important. The technique produces images that replicate flow physics more faithfully than straight-line path integration and other techniques, and more efficiently than physical-optics techniques. Applications of these simulated images are demonstrated in supersonic flows. Time-correlated pairs of shadowgraph images taken from the LES using this technique are used in conjunction with an image-correlation velocimetry technique to compare the estimated convection velocity field in the LES to that of experiments of the same flow. Agreement between the two is good with a maximum variance of 5% by some metrics. This technique can aid in the validation of LES results, allowing quantitative comparison between experiment and simulation, and to extract information unattainable by experiment alone. Comparisons of simulated and experimental jet penetration into the supersonic freestream are also made.


Introduction
Schlieren and shadowgraph images are an important component of experimental techniques in fluid dynamics. As computational fluid dynamics evolves as an analysis tool with more powerful computers and more advanced models, the ability to both replicate and augment experimental results based on simulated flow fields by forward-modeling becomes even more relevant.
Experiments are often analyzed to infer the physics that yield their results. This approach can be described as inverse modeling: given the observations, infer what led to them. However, many observations are not invertible and cannot uniquely yield the information that led to them.
Forward modeling, on the other hand, aims to simulate observations by explicitly modeling both the phenomena studied and the measurement process. With sufficient knowledge of the physics and processes that led to the measurements, this forward-modeling approach allows more direct comparisons between theory and experiment in observation space.
The discussion below aims to illustrate this approach by describing a geometrical-optics technique for simulating schlieren and shadowgraph images from three-dimensional numerical data by tracing rays through the simulated flowfield, and processing these rays by including and modeling the optical components required for image formation using ray tracing, to a level of approximation deemed adequate.
Flow-field data from spatially and temporally resolved numerical simulations can augment information derived from experimental measurements, as well as contributing to simulation validation. Traditional validation methods have relied on pointwise velocity, temperature, and pressure measurements; 'eye-norm' comparisons of experimental schlieren and shadowgraph images to contours of simulated density-gradient fields; or inferred quantities that include velocity (e.g., Beresh et al. 2006;Santiago and Dutton 1997) and scalar fields (e.g., VanLerberghe et al. 2000;Lin et al. 2010), for example.
Illustrating an application of numerical images, this paper describes a technique for direct comparison of simulation results to experimental data in high-speed flows. This technique also extracts inferred convection velocityfield measurements from time-correlated, high-speed experimental and numerically simulated shadowgraph images, as well as other information.
In high-Reynolds number flows with a variable index of refraction (IoR), turbulent structures marked by refractive-index gradients at their interfaces are discernible in schlieren and shadowgraph systems. The interface geometry as well as the convection velocity of these structures can then be determined from sequences of two or more time-correlated images.
In the case studied here of a turbulent jet in supersonic crossflow, IoR-gradient interfaces are generated as the injected fluid entrains, disperses, and mixes with the crossflow fluid. Cymbalist (2016) describes a frame-averaged schlierenimage correlation velocimetry (SICV) technique suitable for high-speed complex turbulent flows. This paper compares experiment and computation by simulating the experimental measurement, as part of the experiment itself, and then processing both experimental and simulated images using the same SICV methodology (e.g., Burns et al. 2015).
IoR variations deflect/refract rays passing through them. Their curvature is given by Settles (2001) as where n( , t) is the local IoR field.
The deviation angle of these rays, , (the integral of the curvature in z) can be imaged using schlieren techniques, and their convergence and divergence (partial derivative of in x and y) can be imaged using shadowgraph techniques. Since the integral of n∕n is ln(n) , the deviation angle becomes x i = ∫ x i [ln(n)]dz ; a more compact notation.
Schlieren images are modeled in terms of the first derivative of the natural log of the IoR field, i.e., ln(n) , while shadowgraph images are typically modeled in terms of its second derivative (e.g., Liepmann and Roshko 1957;Vasiliev 1971;Settles 2001). The index of refraction is a function of density and composition, and for gases, n = 1 + , with ≪ 1 , such that gradients of ln(n) ≈ can frequently be approximated as gradients of density, , only.
A common method for approximating schlieren images in simulated flow-fields is to plot contours of |∇ | at a single plane in the flow (e.g., Srinivasan and Bowersox 2008; Kawai and Lele 2010;Chai and Mahesh 2011;Ferrante et al. 2011;Vuorinen et al. 2013). This approach highlights shocks and eddy structures. By plotting contours of ∕ y = y for a single planar slice, Wang et al. (2013) showed gradient directions for the jet center-plane shock system in their simulation.
To account for the 3D nature of turbulent flow, rays passing through near-unity IoR flows, like air, are often approximated as straight lines with intensity on an imaging focal plane calculated as the integral along these paths of y for schlieren (e.g., Svakhine et al. 2005;Guo and Zhang 2004) and ∇ 2 2D for shadowgraph (e.g., Svakhine et al. 2005). Yates (1993) documented this schlieren straight-line method for n∕n. Brownlee et al. (2011) accounted for ray refraction by the flow using a full ray equation to trace curved ray paths through an IoR field. They simulated schlieren and shadowgraph intensity from the concentration of rays or 'photons' reaching an image plane with a Gaussian filter for smoothing. Gori and Guardone (2018) built on this work to make use of multi-processor and graphical-processing unit capabilities at the same time, with an extension to non-ideal compressible fluid flows. Brownlee et al. (2011) acknowledged the need to incorporate the effect of the optical system for a faithful comparison between experimental and computational shadow images. We extend this ray-processing method by including the optical system that forms the image in the model, both ahead of and after the probed IoR field, in addition to tracing rays through the turbulent IoR field, as Brownlee et al. (2011) recommended. Anyoji and Sun (2007) and Sun (2007) modeled imageformation in simulated optical systems to assess the effect of optical elements. They applied Snell's Law at cell boundaries with a 2D computational flow-field acting as a transparency in the test section. In the case of optical elements, Snell's Law can provide a reasonable approximation. However, Snell's Law is best suited for interfaces between distinct media with different indices of refraction. Its application to continuously varying IoR fields can yield artifacts and low fidelity overall.
The following sections describe the experiment, numerical method used to model the flow field, and the numerical modeling of the schlieren and shadowgraph image formation. The resulting simulated schlieren and shadowgraph images are presented, compared with experimental results, as well as with images based on other techniques typically employed.

Experimental and numerical turbulent flow field
The schlieren and shadowgraph simulation technique proposed below can be used with any gas flow field with variable IoR and optical setup. In this case, the technique is applied to a transverse jet in supersonic flow with applications to combustion and propulsion in high-speed flows. The resulting flow is three dimensional, unsteady, and turbulent with important temporal and spatial flow field information that is not readily accessible via experimental diagnostics. The experimental setup and numerics used to produce the LES flow field are described in the following section. This sets the imaging technique's context and, in particular, the optical elements considered in its development.

Facility overview
The experiments are performed at the GALCIT supersonic shear-layer laboratory ( S 3 L ), a blowdown facility with test times on the order of seconds (e.g., Hall 1991;Slessor 1998;Bond 1999;Bergthorson et al. 2009;Bonanos et al. 2009).
In the experiments described here, a weak, underexpanded, sonic helium jet is injected into a supersonic nitrogen crossflow at M ∞ = 1.5 at ambient total temperature. The jet-to-crossflow momentum-flux ratio (e.g., Mahesh 2013), is near-unity, with a velocity ratio, u j ∕u ∞ ≈ 2 . The transverse jet is injected at a right angle through a circular orifice of diameter, d = 5.08 mm (0.2 �� ) . The edge of the injector in the lower plenum has a fillet of radius 1.59 mm ( 1∕16 �� ). Table 1 summarizes the flow parameters. Further details of the gas-delivery system are given in Cymbalist (2016) and the references above to the facility.
The optically accessible portion of the test section is detailed in Cymbalist (2016) and summarized here. The upper guidewall is slightly diverged ( div = 1.05 • ) to offset boundary-layer growth and avoid imposed streamwise pressure gradients. The distance between the upper and lower guidewalls is 32.2 mm ( 1.26 ′′ ) at the exit of the converging-diverging nozzle. A 30 • expansion ramp follows the quasi-uniform area section and accelerates the flow, isolating the upstream flow from downstream boundary conditions and other effects. The jet orifice is centered 95.25 mm ( 3.75 ′′ ) upstream of the ramp and supplied from the lower plenum. The width of the test section is 152.4 mm ( 6.0 ′′ ) throughout.
The sidewalls of the optically accessible portion of the test section are 50.8-mm ( 2.0 ′′ )-thick BK7 glass, lined with

Diagnostics and data processing
A high-speed camera captures time-correlated, high-resolution schlieren/shadowgraph image pairs using a folded-Z schlieren system, shown in Fig. 1, with a pulsed LED light source of small, but non-negligible extent that is accounted for, as described below. The schlieren system is based on a pair of 10-inch spherical mirrors (Hermanson 1985, SM1 and SM2 in Fig. 1), and uses high-quality 85 mm, f/1.4, lenses (L1 and L2) at full aperture to form the image on the camera sensor. At f/6.5, the spherical mirrors do not perfectly collimate the light passing through the test section, which does not emanate from a point source, however. A PCO.Dimax HD camera was modified to support double-framing with a Δt IF = 6 μs interframe time. A white CREE XM-L LED with a custom-designed electronics driver is pulsed on each side of the interframe (pulse duration, Δt LED = 0.6 μs ) to generate the images in each pair.
The system is capable of capturing up to 1500 imagepairs per second at HD-level resolution ( 1920 × 1080 pixels). Some 4000 image pairs are captured during a run, in addition to 100 prior calibration frames.
In general, shadowgraph images are created by focusing the camera some finite distance, g, away from the test section centerline so what is imaged is the deflection of rays at this location. The value of g controls the maximum sensitivity of the image. However, in flows with strong IoR gradients, such as shocks and supersonic mixing of gases with significantly different indices of refraction, contact shadowgraphs, where g → 0 , are commonly used (Chang and Chang 2000;Settles 2001;Discetti and Ianiro 2017), as sensitivity is not a limiting factor, while the precise location of shocks and other features is more important-something that becomes distorted as the distance g increases. For flows without strong gradients, a contact shadowgraph produces scant variation in illumination in the viewing plane. In the experiments described here, the camera was not refocused between schlieren and shadowgraph images. The shadowgraph images are, therefore, best described as focused contact shadowgraphs.
There is an effective continuum between schlieren and shadowgraph images in this case. A schlieren image with a knife edge sufficiently far from the standard 50% mark reverts to a focused contact shadowgraph, increasing contributions from second IoR derivatives.

Frame-averaged schlieren image-correlation velocimetry
Sets of image pairs are post-processed using the frame-averaged cross-correlation algorithm detailed in Cymbalist (2016) to produce convection-velocity estimates. This technique shares elements with the method of Meinhart et al. (2000) for microfluids applied to schlieren and shadowgraph images, and extends the work of authors such as Jonassen et al. (2006) to intermittent turbulence. The process is outlined briefly here. Extracting convection velocity from multiple image pairs occurs in two steps. First, the images are processed to partially suppress quasi-stationary structures such as shocks and expansions, and background Moiré patterns including fringe patterns from the BK7 and Pyrex glass sheet contact surfaces. The mean intensity of a sequence of N images, Ī , is subtracted from the individual intensity fields from a pair, I A,n and I B,n , where n is the index of an image pairs in a sequence.
Next, sets of N image-pairs undergo a frame-averaged cross-correlation to determine the time-averaged displacement of coherent structures and the convection velocity field. The size of the A N tile, the B N template, and the number of images (N) used to extract the time-averaged signal depends on the flow, image resolution and quality, and image-pair spacing. For this particular flow, [32 × 16]-sized tiles and [64 × 32]-sized templates were extracted from N = 32 image pairs. The final cross-correlation signal is given as Greater detail about this method is documented in Cymbalist (2016, Appendix A).

Numerical method
The numerical simulation of the flow relies on the US3D computational framework: an unstructured, parallel, finite-volume compressible Navier-Stokes solver (e.g.,   et al. 1999). In the present work, second-order inviscid fluxes are used with second-order Crank-Nicholson time integration and the full-matrix point-relaxation method (Wright et al. 1996). This fully implicit method allows time steps to be dictated by considerations of accuracy, rather than stability. Turbulence is modeled using improved delayed detachededdy simulation (IDDES): the hybrid Reynolds-averaged Navier-Stokes (RANS)/LES method of Shur et al. (2008). This acts as wall-modeled LES with the Spalart-Allmaras one-equation turbulence model providing closure in regions of attached boundary layers, and the constant-coefficient Smagorinsky subgrid scale model for regions away from walls, using a Smagorinsky constant of 0.18.
A hexahedral simulation grid is built on a CAD model of the S 3 L test section using the LINK3D software from GoHypersonic. The grid extends into the lower plenum to capture the correct injected mass-flow rate and flow features when plenum total temperature and pressure are matched to experimental values (Peterson and Candler 2010), as in this case. The grid is clustered to resolve the boundary layer on all four walls with a base grid of ∼ 68 M cells, and ∼ 88 M after clustering.
Separate flush-wall simulations were compared to noninjecting experiments to provide the guidewall divergence required to match the near-zero streamwise pressure gradients in the experiments. In the simulation, this required 0.5 • compared to 1.05 • in the experiment since the RANS boundary layers on three of the walls grew more slowly than either the time-varying simulated boundary layer on the bottom wall, or the experimental boundary layers.
Both experiments and computations have shown the dependence of jet behavior on its interaction with the crossflow boundary layer (e.g., Fric and Roshko 1994;Subbareddy et al. 2006;Ferrante et al. 2011), even though the sensitivity for normal injection is expected to be lower than for inclined-jet injection (e.g., Peterson and Candler 2010). The outflow plane from a separate 3D RANS simulation of the converging-diverging nozzle that generates the Mach 1.5 crossflow for the S 3 L wind tunnel provides the basis of the inflow for the jet simulation described here.
Time-and spanwise-varying perturbations are added to the boundary layer on the bottom wall using the digital filter technique (Kline et al. 2003;Xie and Castro 2008;Touber and Sandham 2009) based on the work of Kartha et al. (2015).
The turbulence field is generated to match prescribed Reynolds stresses. Time-averaged statistics compare favorably with experimental results (Kartha et al. 2015). Figure 2 shows an example of the resulting flow field. Q-criterion isosurfaces (e.g., Chacin et al. 1996;Haller 2005) are colored by streamwise velocity in a half volume. The center mid-span plane plots contours of |∇ | , highlighting the jet and shear layer shock system in this plane. The figure illustrates the complexity of the LES flow responsible for the simulated schlieren and shadowgraph images.

Simulated shadowgraph and schlieren images
The method for simulating schlieren and shadowgraph images presented here can be applied to any optical system and variable-IoR gas-flow experiment. However, the optical elements and their arrangement, and the IoR field in question are specific to the experiment described in Sect. 2. In this section, first the ray-tracing method developed by Brownlee et al. (2011) and its implementation in the present paper is outlined. Then the extension to previous work by modeling the optical system is discussed, and finally the imageformation approach.

Test section ray tracing
The interaction between a collimated light beam and the testsection flow field is modeled using geometrical optics, which assume an infinitesimal wavelength, . This assumption is reasonable in the majority of the flow-field regions since physical length scales are significantly greater than for visible light (shocks with the thickness of a few mean free paths may be an exception).
Geometrical optics treats light as rays, for which the equation is approximated by where is the position vector of the light ray, s is its arc length, and n is the local index of refraction along the ray path (e.g., Born and Wolf 1970).
In addition to the negligible-wavelength assumption, paraxial rays (rays with small deviation angles and close to the optical axis) and a short transit time of light through the optical system relative to changes in the flow (i.e., frozen flow) are assumed.
To trace the ray through a known IoR field, given in Cartesian coordinates, Eq. (4)  where ∕ s is a unit vector in the direction of the ray. For most gases at moderate pressure, n ≈ 1 and so treating as a direction vector is a reasonable approximation. Two partial differential equations can then be written (Eq. 6) that are updated at each point in space to trace the ray through the test section. This is done here with a central differencing scheme for second-order accuracy: The IoR is calculated based on partial densities in the simulated flow field using the following equation: where i is the partial density of the ith gas, and K i is its Dale-Gladstone constant, allowing for full modeling of the composition and density-dependence of the local IoR.

Initial conditions and light-source extent
For each ray traced through the test section, an initial position, , and direction, , are provided. The unstructured mesh used in the LES poses a challenge for ray tracing. This is mitigated by projecting the LES solution onto a uniform 3D Cartesian mesh for the purposes of ray tracing. The Cartesian coordinates used in the analysis are then x = streamwise, y = transverse, and z = spanwise.
This regular mesh has a resolution equivalent to that of the well-resolved LES grid center plane (1356 × 308 × 400 cells for the full test section). The CFD grid has greater resolution in the boundary layers and at the center of the jet. Higher resolution is not adopted in the ray-tracing grid. Only scant mixing and weak IoR gradients are encountered in the jet-core region, which is not optically interesting.
Boundary-layer gradients on vertical walls have a negligible effect on rays passing through them, while vertical gradients in the floor and ceiling boundary layers are so strong that, even with the reduced resolution, rays are found to be deflected outside the finite lens aperture, as detailed below. In this case, a greater boundary-layer resolution in the raytracing grid would have had limited effect on the formation of the boundary-layer image.
Ray deflection in turbulent regions is dominated by the geometry of interfaces between freestream and turbulent regions, and influenced considerably less by interactions with IoR fluctuations internal to the turbulence, as shown in Dimotakis et al. (2001). As a consequence, the smallest eddies captured by LES suffice for schlieren and shadowgraph ray-tracing modeling for image formation.
In terms of resolution, we are limited by that of the underlying CFD grid: increasing the resolution of the optical grid or the number of rays of light used to trace the image reduces the width of the Gaussian kernel used to reconstruct the image at the end and does not materially change the definition of the shadow image since the additional rays are traveling through the same subset of the flow field.
Doubling the number of incident light rays used to trace the image took approximately twice as long with no discernible difference in the resulting image.
Consider an example of a plane wave incident on the initial (x, y)-plane of the test-section mesh (equivalent to a point source on the optical axis). For each cell in the initial (x, y)-plane, 16 evenly spaced rays are distributed giving the initial conditions = [x, y, 0] and = [0, 0, 1] . Trials with a higher areal density of incident rays did not produce a more detailed numerical image. The experimental optical system has a light source with a small but finite extent: a 2.1 × 2.1 mm square LED, whose extent cannot be ignored, as discussed below. To simulate schlieren images faithful to the experiments this extended source is discretized into 25 points, as shown in Fig. 3. Initial conditions are set by tracing the principal ray passing through the center of SM1 in Fig. 1 originating from each source location, as shown in the left panel in Fig. 4. In the context of geometrical optics, polarization effects are neglected, and the optical system is simplified by unfolding it into a linear system, neglecting the 45 • mirrors and treating the spherical mirrors as thin lenses of equivalent optical power.
The right panel of Fig. 4 shows the parameters describing a typical ray. Here, and are angles of elevation and azimuth, respectively, and Δx and Δy are offsets from the The 25 sets of rays are combined to form the final image intensity. The red points in Fig. 3 contribute full intensity to the final image, the ones in orange (at the edge of the LED) contribute 50% intensity, while corner points (green) contribute only 25% intensity to represent contributions of the associated tile areas on the LED. The full test section images shown in Fig. 11 required tracing ∼ 166.3 M rays.

Propagator
On exiting the test section, each ray has a coordinate in relation to the optical axis, offsets Δx and Δy , and a direction vector described by the angles of elevation, , and azimuth, , (Sect. 3.2). The angles can be found by rearranging Eq. (8).
Optical-transfer matrices are used to propagate these rays through the remainder of the optical system (Fig. 1). Gauss (1840) first proposed two linear simultaneous equations relating the height and angle of an output ray to its input that he solved using the method of Gaussian brackets. Sampson (1913) translated this into a matrix method for optics, later widely adopted (e.g., O'Neill 1963;Halbach 1964;Gerrard and Burch 1975).
The extension for skew rays, in addition to meridional rays, is used (Attard 1984). Equations 9-11 show Attard's matrix relations for three elements of interest here.
Propagation through constant IoR, along a distance, d: Transmission through an interface (Snell's law): Transmission through a lens: Figure 5 illustrates the effect of Eq. (11). The transfer matrix relates the incoming and outgoing ray position and direction.
In the case of a thin lens, for example, the ray position is not shifted between the two sides of the lens, but the ray is redirected depending on the lens focal length and the distance between the ray and the optical axis. The finite aperture of each lens is included in the model by cutting off rays that are incident outside it. The cutoff provided by the (here horizontal) knife edge is also treated as an optical element. In the same manner, as in the experiment, a process of trial and error sets the height of the cutoff, with the goal of a 10% overall reduction in illumination at the image plane to match image properties in the experiments. While a 50% knife edge setting may be a standard recommendation for schlieren, experimental images in this case were empirically optimized for contrast since the experimental data were mainly limited by the energy available in each light pulse. The experiment generated thousands of images per second. A less than 50% knife Fig. 4 Left: method for determining direction and offset at the test section initial plane caused by the extended light source. Right: geometry used to describe the offset from the optical axis, Δx and Δy , and the angle of elevation, , and azimuth, , of the ray Fig. 5 Example of the application of an optical-transfer matrix for a lens. The lens changes the angle of elevation and azimuth of the ray, but not its offset from the optical axis edge position is generally considered to reduce sensitivity. However, schlieren sensitivity was not a significant issue in this experiment because of the large refractive index gradients in helium and nitrogen interfaces. In the experimental optical setup, the first spherical mirror represents the limiting aperture for the field of view, while the last two confocal lenses provide the limiting apertures for ray deflection. The maximum deflection of rays from the optical axis that may be supported is 1.5 • : well within the paraxial approximation. Importantly, rays in certain high-gradient regions, such as the lower boundary layer, are found to be deflected entirely outside the lens aperture and, as a consequence, do not contribute to image formation.

Image formation
The end result of the ray tracing and propagation is a set of coordinates of rays incident on the focal-plane imaging array. While this may grant an overall impression, it does not render a proper image.
To produce an intensity pattern, as required for an image, each ray is convolved with a Gaussian kernel. A Gaussian beam, as in a laser, comes close to representing a ray (Gerrard and Burch 1975). Gaussian convolution smooths the image, as shown in Fig. 6 by comparing the left panel to the right, and also approximates diffraction effects.

Test case and code verification
Replicating schlieren and shadowgraph images for computational data requires a numerical model for the propagation of light through the optical system. The complex flow field created by the jet in supersonic crossflow (cf. Fig. 2) makes a direct assessment of the model difficult. In this section, a simple shadowgraph test case and simple flow-field are used instead to verify and validate the image-formation method.
Four methods of modeling intensity at the image (focal) plane are compared: the exact analytical solution to Fresnel diffraction; Fresnel diffraction calculated using Fast Fourier Transform (FFTs); geometric ray tracing; and the straight-line path-integral method [ ∇ 2 2D ln(n) ] described in Sect. 1.
The two light-propagation models-physical and geometric optics-derive from Maxwell's equations through simplifying approximations. Physical or Fourier optics considers the wave nature of light, providing a more complete description of electromagnetic propagation than geometrical optics that treats light as thin beams or rays.
Fourier optics accurately treats the highly oscillatory functions in Maxwell's equations, posing an ongoing challenge (e.g., Bruno and Kunyansky 2001). Convergence of the 2D FFTs in this analysis requires hundreds of thousands of points in each dimension, rendering it intractable with current computing power (Chaubell et al. 2009). The simplest optical system that can reproduce an image of an IoR field is a shadowgraph consisting of collimated light passing through the field, and a screen to form the image. A Fourier-optics analysis represents this as a convolution between the test-section phase disturbance and a spherical wave (e.g., Papoulis 1968;Goodman 2005).
A simple IoR field, n(x, y, z), for which an analytical Fourier optics solution can be derived, provided a test. The IoR field is piecewise linear, i.e., as illustrated in Fig. 7. Here, n 0 is the refractive index of the undisturbed field, is the maximum deviation of n from the undisturbed value, and the parameters Δz and b are the width and height of the disturbed IoR field, respectively. was set at 0.0001, which is approximately 40% of the maximum difference in refractive index between N 2 and He at visible wavelengths and standard temperature and pressure.
An incident wave can be expressed as a complex function  where is the position vector, A( ) the wave amplitude, and ( ) its phase, which is the solution of the Helmholtz equation for a monochromatic wave. When transmitted through a transparent medium, amplitude variations are negligible, but the phase changes are not. The phase delay is proportional to the optical path length difference ( S ) caused by passage through the region of varying IoR where n 0 is the unperturbed refractive index, and k 0 = 2 ∕ 0 is the wavenumber of light in a medium with n = n 0 .
An image is formed at a screen or focal-plane array by variations of the incident light intensity. In physical optics, intensity is expressed as I( ) = |U( )| 2 . Geometric optics relies on the distribution of rays on the image plane, as described in Sect. 3. Figure 8 shows the image-plane intensity resulting from the test-case IoR field, calculated by each of the four methods. Figure 8, top left, is the analytical solution to the Fresnel diffraction theory and provides the basis for comparison. Figure 8, top right, is the result of taking the requisite convolution via FFTs. The close match to the analytical solution required 100,000 points to converge to an adequate approximation, even for the simple case. This method is, therefore, not practical. Figure 8, bottom left, is produced by the ray-tracing method described in Sect. 3.1 and required ∼ 5000 rays to converge with a criterion for convergence based on accepting a 0.5% variance in the outcome as a function of the number of rays. To assess this, images were computed with 5120 rays, and 10,240 rays, with convergence achieved with 5120. The fringing characteristic of diffraction, especially owing to the sharp test-case IoR gradient transitions, is not captured by this method, as expected, but the principal phenomena are well represented on a mesh 20 times coarser than that required for the FFT solution. Figure 8, bottom right, is produced by the straight-line integral method. Locations where IoR field gradients change sharply are discontinuities in all derivatives and as such should map to an infinitesimal point of infinite magnitude. Given the discretized nature of a computation, the spikes have a finite width and depth, and only approximate the analytical result. We conclude that, in balancing between computing power and accuracy, geometric ray-tracing offers an advantage for problems of the type of interest here.
From the test cases, we would expect the ray tracing technique to capture the salient flow features. The maximum illumination and its absence are reduced in comparison with the analytical solution, although in accord with the FFT. A notable limitation of ray tracing comes in regions with infinitely sharp, non-differentiable interfaces as in this test case example, which despite being extreme admits a simple analytical solution. Actual phenomena investigated are smoother and do not entail such singularities.

Results
In this section, the importance of modeling the optical system is illustrated. Simulated images created using geometric ray tracing and modeling of the optical system are presented and compared to experimental images. An example of an application of these images is shown using the SICV algorithm for quantitative comparisons between the simulation and experiment.

Modeling the optical system
Including the optical system in the image-simulation process can be credited for improved image fidelity. A total of ∼ 10% of the rays used to trace the full test section images are eliminated by finite lens apertures, without additional filtering. The spherical mirrors also limit the test-section field of view. The confocal lenses cut off a portion of the rays in both the lower boundary layer (noted in Sect. 5.2) and the corner of the expansion fan.
If the finite aperture of the lenses is ignored, most of the lower boundary-layer rays would be approximately returned to their original location. Those from nearest the wall where density gradients are strongest are diverted into the crossflow, to form a new line above the boundary layer. Both the lack of near-wall darkness and the presence of a bright line would be unphysical.
In addition to finite-aperture effects, the finite light-source extent plays an important role, particularly in modeling 3. An extended light source may be thought of as a distribution of point sources across the vertical and horizontal extent of the source (Settles 2001). As described above, we discretize the light source into 25 separate sources. This has the effect of illuminating every point of the test section with bundles of multiple rays.
Since the source and the knife edge are conjugate pairs, the source image is reconstructed at the knife edge. Some portion of undisturbed rays from each discretely modeled source are cut off by the knife edge while rays that passed through the varying IoR field will have been deflected above or below the knife edge. The final effect at the image plane is a uniform dimming of the undisturbed portions of the test section image, with additional dimming or brightening of areas in accord with ray deflections. Figure 9 shows schlieren images created with rays originating from three different light-source points: (top) the top right corner, (middle) the optical axis, (bottom) the bottom right corner. In each case, the rays have passed through different regions of the test section and were differently deflected. The three images are surprisingly different considering the small distance that separates the ray origins (maximum ∼ 2 mm) and respective angles of incidence through the flow.
The 1 mm separating the optical axis from the top and bottom of the light source represents a maximum additional angle at the initial plane of the test section, s,max = 0.04 • . While this is small, it is significant when compared with the predominant exit direction of rays from the test section, most of which travel approximately straight. Figure 10 plots the probability density function (PDF) of exit angles of elevation and azimuth for the optical axis source as a 2D contour plot (left) and as a surface (right). Both plots are color coded in terms of log(PD) , where PD is the probability density.
The left panel in Fig. 10 includes a black square representing the extent in each direction of s,max . This region is small, but the surface plot in the right panel of Fig. 10 shows the sharp spike contained within it. Indeed, 56% of all rays exiting the test section are bounded by this box, and s,max ∕ ≥ 10% for 92% of all rays. Figure 10 demonstrates the importance of including the small spatial extent of the light source in modeling schlieren images. Individual slices in Fig. 9 appear artificial with large, stark unilluminated regions. The combined effect yields more realistic images that capture the flow features observed in the experimental schlieren and shadowgraph images.
There was no attempt to capture the absolute intensity of every ray and, as a consequence, a uniform illumination level was left adjustable to compare with experiment. The average illumination level in the undisturbed freestream portion of the computed images was matched to that of the same region in the experimental images. Figure 11 shows examples of experimental shadowgraph and schlieren images, and their counterparts resulting from the process described in Sect. 3 for the full test section. For comparison, schlieren and shadowgraph images created using the straight-line path-integral method described in Sect. 1 are shown in Sect. 5.3.

Numerically simulated images
The computational schlieren (Fig. 11, upper middle) and shadowgraph (Fig. 11, bottom) provide a good rendition of the experimental schlieren and shadowgraph images. The improved correlation with experiment noted in the previous section can be discerned in the region of the interaction between the shock and the jet, particularly in the schlieren image. Additionally, near-wall rays in the lower boundary layer are deflected entirely outside the lens aperture resulting in dark regions with no light, even in the shadowgraph with no knife edge filter. Omitting the finite aperture size results  There are differences between the simulated and experimental images resulting from both the modeling of the optics, as well as the modeling of the flow-field itself. Notably, the expansion-ramp shear layer in the experiment is at ∼ 7 • below the horizontal, while in the simulation the angle increases to ∼ 14 • , and the shear layer reattaches on the bottom wall before exiting. However, in both the experiment and simulation, the expansion region isolates the main test section region from conditions downstream, with the preramp region the section of interest. The simulations register a faster initial growth, but penetration of the jet at later stages matches the experiment well, and conditions in the region of interest show a generally good match between simulation and experiment.
A final difference to note is the presence of significantly more weak shocks in the experiment than the simulation. These are caused by small stationary geometrical discontinuities on wall boundaries, and to a lesser extent, eddies within the boundary layer. Large scales in the latter, however, convect with velocities close to the freestream and subsonic with respect to the freestream and do not generate Mach waves. The simulation included a ridge at the location of leading edge of the injector insert a few cells high and the same width to mimic this weak shock, but is otherwise mathematically smooth, leading to vastly reduced shock activity. It should be noted that these shocks cause dark patches on the upper surface of the jet body in the shadowgraph experimental images where the shocks pass through this interface. Additionally, smaller scale eddies within the boundary layer that convect sufficiently slowly to be supersonic with respect to the freestream are assumed to be part of the wall model in the simulation rather than present as effective roughness like the outer eddies, and generate no Mach waves. Figure 12 compares images produced by the method proposed in this paper and by the straight-line path-integral method. The intensity in a schlieren image produced by a horizontal knife edge is calculated as I ∝ ∫ y (ln n)dz , as described in Sect. 1. Shadowgraph images created by the straight-line integral and ray-tracing techniques are more similar than the two schlieren images. However, the straightline integral shadowgraph images are less smooth, with artificially sharp contrasts between regions of light and dark, consistent with the test case shown in Sect. 4. The same can be said for the straight-line integral schlieren image. Second, straight-line integral images do not capture some important flow physics evident in both experimental and raytraced images. Notably, the near-wall boundary layer in the straight-line integral images does not register the complete absence of light, characteristic of the other images. In the shadowgraph image, the near-wall region is in fact lighter than that above it. The straight-line integral images also do not capture the interaction between the reflected shocks and the jet body.

Straight-line path integral images vs. ray tracing
Finally, accounting for only transverse IoR gradients means that the bow shock in the schlieren image is incorrectly rendered. There is a very faint darkening ahead of the bright white intensity of the bow shock, but this is dominated by the high-intensity region behind it.

Extension to the straight-line path integral schlieren technique
Line-of-sight methods best illustrate shocks and expansions with gradients principally in line with the cutoff. Yates (1993) computed schlieren images, for example, use a horizontal knife edge with all shocks at less than 45 • . Pagendarm and Post (1995), on the other hand, show a conical flow with vertical knife edge where leading-edge attached shocks are at or below 45 • , i.e., principally out-of-line with the knife-edge. Here, the simulated image shows the bright center of the shock from the experimental image, but does not clearly capture the dark regions on either side of this dark core. In the case shown here, the bow shock and density gradients are principally out-of-line with the knife edge and not adequately captured by a straight-line computation that considers gradients in y. An extension to the straight-line integral method for schlieren images is proposed: the presence of a horizontal knife edge does not occlude the contribution of the second derivative of the IoR field in the direction parallel to the knife edge. As discussed in Sect. 2.2, there is an effective continuum behavior between shadowgraphs and schlieren images. A knife edge far below a 50% level reverts to a contact shadowgraph, while one at 50% is a traditional schlieren. The schlieren image intensity can then be estimated as where has units of 1/length. For a focused shadowgraph optical system, as here, is large since it is a function of the focal length of the lens ahead of the knife edge, the unobscured height of the source image at the knife edge, and the distance between the phase object in the test section and the imaging screen. Combining the contributions of the two derivatives required some trial and error to set the value of for best match with experimental images. The resulting image is shown in Fig. 13. The proposed extension to the straight-line integral schlieren intensity calculation improves correspondence with experiment. Intensity rendering in the shock region is improved and regions such as the upper edge of the near-injector jet are darkened. Despite this improvement, this technique still does not capture the shock-jet body interaction nor the near-wall boundary-layer shadow that is rendered well by the raytracing method described here.
If lacking the computational resources to engage in the ray tracing-based method detailed in this paper, the proposed extension to the straight-line path method can render improved images over considering only the first derivative of density perpendicular to the knife edge, for example.
Straight-line integral images took approximately 2 min to produce on a standard desktop with an i7 processor and 32 GB of RAM. Ray-tracing through the test section took 3 h on the same machine, propagating through the optical system took 15 s and applying the Gaussian filter took 0.5 h for 26.4 M rays. Ray tracing takes substantially longer than the straight-line integral method, but is more physically faithful. Adopting the octree method of Brownlee et al.  Gori and Guardone (2018) also leverage the GPU for image generation and rendering. The principal efficiency is over computing images with Fourier optics as a further increase in physical fidelity. This latter would require greater than an order-of-magnitude mesh refinement, following the test problem in Sect. 4.

Simulated image application
The SICV algorithm outlined in Sect. 2.3 (Cymbalist 2016) is applied to experimental images and those created from the numerical simulation. A schlieren or shadowgraph image is the projection of a 3D field onto a 2D surface. The output of the algorithm is both temporally and spatially averaged across the span of the flow. In the following analysis, shadowgraph images were used. Figure 14 plots the time-and spanwise-averaged streamwise velocity of imaged convected turbulent structures from the experimental (top) and computational (bottom) jets. The crossflow outside the jet body and the boundary layers offers no signal to track. Velocity estimates very close to the jetorifice ( − 0.5 < x∕d < 1.5 centered around the jet orifice axis) were deemed unreliable and ignored. The simulations have a lower resolution than the experimental images, and scant detail is captured in the initial regions of the jet body in the numerical images. For regions downstream of x∕d ∼ 3 , however, the convection velocity field in the experiments and simulations is adequately discerned and shows similar trends. The jet body is initially faster than the crossflow, as expected for a helium jet, but as entrainment and mixing occur, the jet body decelerates to below the crossflow velocity, re-accelerating to approximately match it towards the end of the domain.
To help quantify the comparison between the experimental and computational convection velocities, an average convection velocity field, u c , at a given streamwise location is determined by binning all measurements in a region −d∕2 < x < d∕2 for each integer x/d. Figure 15 compares the binned SICV convection velocity estimated from the experimental and computational images using the same SICV algorithm. Uncertainty bars denote the standard deviation of data within each bin. Closer to the injector, there is a wider range of velocity estimates within each bin because of high velocity gradients. For a weak jet ( J ∼ 1 ), the deflection into the crossflow has a contribution from pressure differences between the up-and downstream faces of the jet entering the crossflow (e.g., Hasselbrink and Mungal 2001). If the jet fluid exits the injector at a higher speed than the crossflow, pressure effects in weak jets can lead the jet wake to be significantly faster streamwise than the crossflow, despite a normal jet contributing no streamwise component of momentum (e.g., Schetz and Billig 1966).
While u c,com < u c,exp at all locations, the maximum difference is 5% and the correspondence between the experiment and simulation is good, both quantitatively and qualitatively. The rapid deceleration of the estimated convection velocity to below the freestream velocity is captured in both experiment and simulation, as well as the gradual re-acceleration. Figure 15 also includes freestream readings from the experiment, extracted from Mach-wave angle estimation, and the simulation based on the average of fluid within the volume and outside the boundary layer containing less than 0.1% helium.

Conclusions
Geometric ray tracing and optical-system modeling was used to derive improved-fidelity simulated schlieren and shadowgraph images based on LES data of a normal jet in supersonic crossflow. The method provides a good approximation of experimental images, correctly capturing several flow features observed experimentally, such as the jet-shock interaction and inferred velocity fields.
Elements of the experimental schlieren system modeled include the finite extent of the light source, the light collimation, the varying-IoR flow-field, relay lenses, a knife edge cutoff, and image formation at the focal-plane imaging array. Previous work by Brownlee et al. (2011) focused on modeling light propagation through the varying IoR gas field. Here, that work is extended by incorporating other elements of the experimental schlieren system in the system modeling, on both sides of the disturbed IoR field, to improve quantitative image formation. Brownlee et al. (2011) note that the finite light-source extent and lens relays may be important in computing faithful reproductions of experimental images. This is confirmed by modeling the optical system alongside the optical-propagation through the IoR field. It is shown here to have a significant effect on the fidelity of images produced, relative to images produced by ray tracing of the gas flow alone from a point source. The forward-modeling approach described in Average convection velocity estimated using SICV from experimental ( U c,exp , light blue circular symbols) and computational ( U c,comp , red square symbols) images. Note large origin offsets. Jet convection velocity is compared to freestream readings from the experiment (Mach-wave angle estimation, U shk , scattered circles) and simulation (average over the central half of the volume, U ∞,comp , orange squares) this paper improves comparisons between theory and experiment in observation space.
Large refractions in the near-wall boundary-layer region orient rays found to be occluded by finite-aperture optics downstream in the imaging system. Failure to model this optical element results in a lack of some experimentally registered image features, and spurious unphysical highintensity streaks.
Despite only ∼ 2 mm, the light-source extent has a significant effect on rays cut off by the schlieren knife edge. Almost 60% of rays exiting the test section are deflected by an angle less than or equal to that created by an offset subtended by the light-source extent, which, as a consequence, cannot be ignored. Contributions from each of 25 discrete origins used to model the extended source are found to be discernibly different; their superposition improves image simulation than from a single point source.
Both this ray-tracing technique, and, to a lesser degree, the extension to the straight-line path integral technique improve analysis of the computational images in two key quantitative metrics: the representation of the intensity contribution of the boundary layer, and the improvement of the representation of shocks in computational schlieren images. In both cases, using a simpler image-generation method leads to quantitatively inaccurate intensity levels in these regions.
The application of such simulated images for quantitative analysis was also demonstrated. Using 32 pairs of time-correlated images with the schlieren image correlation velocimetry algorithm allowed estimates of convection-velocity maps for both simulation and experiment. In the jet body region of interest, a good match is found between the average convection velocity in the experiment and simulation with a maximum variance of 5% between the two, also providing a useful validation of the numerical simulation of the flow.
This imaging technique, in combination with SICV, can aid in quantitative comparisons between experiment and simulation with comparisons drawn directly between convection velocity fields computed with the same method, extending information that can be extracted from experiments.