Vol.:(0123456789)
1 3
Experiments in Fluids (2019) 60:134
https://doi.org/10.1007/s00348-019-2774-6
RESEARCH ARTICLE
Simulating schlieren and shadowgraph images from LES data
Elizabeth Luthman
1
· Niccolo Cymbalist
1
· Daniel Lang
1
· Graham Candler
2
· Paul Dimotakis
1
Received: 28 December 2018 / Revised: 6 July 2019 / Accepted: 8 July 2019
© The Author(s) 2019
Abstract
Abstract
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 shad-
owgraph 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.
Graphic abstract
*
Elizabeth Luthman
eluthman@caltech.edu
Extended author information available on the last page of the article
Experiments in Fluids (2019) 60:134
1 3
134
Page 2 of 16
1 Introduction
Schlieren and shadowgraph images are an important com-
ponent of experimental techniques in fluid dynamics. As
computational fluid dynamics evolves as an analysis tool
with more powerful computers and more advanced mod-
els, 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 flow-
field, 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 pres-
sure measurements; ‘eye-norm’ comparisons of experi-
mental 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 sim-
ulation results to experimental data in high-speed flows.
This technique also extracts inferred convection velocity-
field 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 refrac-
tive-index gradients at their interfaces are discernible in
schlieren and shadowgraph systems. The interface geom-
etry 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 cross-
flow fluid.
Cymbalist (
2016
) describes a frame-averaged schlieren-
image 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 pro
-
cessing 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
)]
d
z
; a more compact notation.
Schlieren images are modeled in terms of the first deriva-
tive of the natural log of the IoR field, i.e.,
ln
(
n
)
, while shad-
owgraph 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 sin-
gle 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 gradi-
ent directions for the jet center-plane shock system in their
simulation.
To account for the 3D nature of turbulent flow, rays pass-
ing through near-unity IoR flows, like air, are often approxi
-
mated 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
(1)
휕
2
x
휕
z
2
=
1
n
휕
n
휕
x
,
휕
2
y
휕
z
2
=
1
n
휕
n
휕
y
,
Experiments in Fluids (2019) 60:134
1 3
Page 3 of 16
134
through an IoR field. They simulated schlieren and shadow
-
graph intensity from the concentration of rays or ‘photons’
reaching an image plane with a Gaussian filter for smooth-
ing. Gori and Guardone (
2018
) built on this work to make
use of multi-processor and graphical-processing unit capa-
bilities at the same time, with an extension to non-ideal com-
pressible fluid flows.
Brownlee et al. (
2011
) acknowledged the need to incorpo-
rate 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 opti-
cal 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 image-
formation in simulated optical systems to assess the effect of
optical elements. They applied Snell’s Law at cell bounda-
ries with a 2D computational flow-field acting as a trans-
parency in the test section. In the case of optical elements,
Snell’s Law can provide a reasonable approximation. How
-
ever, Snell’s Law is best suited for interfaces between dis-
tinct media with different indices of refraction. Its applica-
tion to continuously varying IoR fields can yield artifacts
and low fidelity overall.
The following sections describe the experiment, numeri-
cal method used to model the flow field, and the numerical
modeling of the schlieren and shadowgraph image forma-
tion. The resulting simulated schlieren and shadowgraph
images are presented, compared with experimental results,
as well as with images based on other techniques typically
employed.
2
Experimental and numerical turbulent
flow field
The schlieren and shadowgraph simulation technique pro-
posed below can be used with any gas flow field with vari-
able IoR and optical setup. In this case, the technique is
applied to a transverse jet in supersonic flow with applica-
tions 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.
2.1 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, underex-
panded, sonic helium jet is injected into a supersonic nitro-
gen 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 trans-
verse 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 pres-
sure gradients. The distance between the upper and lower
guidewalls is 32.2 mm (
1.26
) at the exit of the converg
-
ing–diverging nozzle. A
30
◦
expansion ramp follows the
quasi-uniform area section and accelerates the flow, isolat-
ing the upstream flow from downstream boundary condi-
tions 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
(2)
J
=
휌
j
u
2
j
휌
∞
u
2
∞
,
Table 1
Flow parameters
T
0,
∞
p
0,
∞
M
∞
u
∞
훾
∞
J
289 K
550 kPa
1.5
430 m/s
1.4
0.9
T
0,j
p
0,j
M
j
u
j
훾
j
295 K
515 kPa
1.0
875 m/s
1.67
Experiments in Fluids (2019) 60:134
1 3
134
Page 4 of 16
replaceable protective 3.2-mm (
0.125
)-thick Pyrex sheets
on both sides of the test-section interior.
2.2
Diagnostics and data processing
A high-speed camera captures time-correlated, high-reso-
lution 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 dura-
tion,
Δ
t
LED
=
0.6
μ
s
) to generate the images in each pair.
The system is capable of capturing up to 1500 image-
pairs per second at HD-level resolution (
1920
×
1080
pixels).
Some 4000 image pairs are captured during a run, in addi-
tion 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.
2.3
Frame‑averaged schlieren image‑correlation
velocimetry
Sets of image pairs are post-processed using the frame-aver
-
aged 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 par
-
tially 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,
̄
I
, 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).
2.4 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.,
(3)
̄
퐂
(
i
,
j
)=
1
N
N
n
=
1
L
−
1
l
=
0
M
−
1
m
=
0
퐀
(
l
,
m
,
n
)
퐁
(
l
−
i
,
m
−
j
,
n
)
.
Fig. 1
Experimental optical setup. Folded-Z schlieren system with
right-propagating light
Experiments in Fluids (2019) 60:134
1 3
Page 5 of 16
134
MacCormack and Candler
1989
; Wright et al.
1996
; Nompe-
lis et al.
2004
; Subbareddy and Candler
2009
). Fluxes are
calculated using the low-dissipation, kinetic-energy-consist-
ent method of Subbareddy and Candler (
2009
) that partitions
fluxes into a non-dissipative, symmetric part and a dissi-
pative flux. The dissipative portion is calculated from the
modified Steger–Warming upwind flux (MacCormack and
Candler
1989
) damped by a shock-sensing factor,
훼
∈[
0,1
]
(Ducros 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 detached-
eddy 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 non-
injecting experiments to provide the guidewall divergence
required to match the near-zero streamwise pressure gradi-
ents in the experiments. In the simulation, this required
0.5
◦
compared to
1.05
◦
in the experiment since the RANS bound-
ary 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 cross-
flow boundary layer (e.g., Fric and Roshko
1994
; Sub-
bareddy 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 Reyn-
olds stresses. Time-averaged statistics compare favorably with
experimental results (Kartha et al.
2015
).
Figure
2
shows an example of the resulting flow field. Q-cri-
terion 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.
3
Simulated shadowgraph and schlieren
images
The method for simulating schlieren and shadowgraph images
presented here can be applied to any optical system and vari-
able-IoR gas-flow experiment. However, the optical elements
and their arrangement, and the IoR field in question are spe-
cific 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 mod-
eling the optical system is discussed, and finally the image-
formation approach.
3.1
Test section ray tracing
The interaction between a collimated light beam and the test-
section flow field is modeled using geometrical optics, which
assume an infinitesimal wavelength,
휆
. This assumption is rea-
sonable 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 equa-
tion is approximated by
(4)
휕
휕
s
n
휕
퐱
휕
s
=∇
n
,
Fig. 2
Half-volume results of jet-injection LES plotting isosurfaces
of Q-criterion colored by streamwise velocity. Freestream flows from
lower right to upper left. The center mid-span on the right plots con-
tours of
∇
휌
showing shock structures and eddy interfaces in this
region. Uniform gray surfaces indicate test section boundaries
Experiments in Fluids (2019) 60:134
1 3
134
Page 6 of 16
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, par
-
axial 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 Car
-
tesian coordinates, Eq. (
4
) is numerically integrated in
s
.
Brownlee et al. (
2011
) define a direction vector,
퐯
, of the
ray as
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 simu-
lated flow field using the following equation:
where
휌
i
is the partial density of the
i
th gas, and
K
i
is its
Dale–Gladstone constant, allowing for full modeling of the
composition and density-dependence of the local IoR.
3.2
Initial conditions and light‑source extent
For each ray traced through the test section, an initial posi-
tion,
퐱
, 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 Car
-
tesian coordinates used in the analysis are then
x
=
stream-
wise,
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 negligi-
ble effect on rays passing through them, while vertical gra-
dients in the floor and ceiling boundary layers are so strong
that, even with the reduced resolution, rays are found to be
(5)
퐯
=
n
휕
퐱
휕
s
,
(6)
휕
퐱
휕
s
=
퐯
n
;
휕
퐯
휕
s
=∇
n
.
(7)
n
=
1
+Σ
i
K
i
휌
i
,
deflected outside the finite lens aperture, as detailed below.
In this case, a greater boundary-layer resolution in the ray-
tracing 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 shadow
-
graph ray-tracing modeling for image formation.
In terms of resolution, we are limited by that of the under
-
lying 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 defini-
tion of the shadow image since the additional rays are trave-
ling 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 discern-
ible 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 pass
-
ing 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 unfold-
ing 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 describ-
ing a typical ray. Here,
휃
and
휙
are angles of elevation and
azimuth, respectively, and
Δ
x
and
Δ
y
are offsets from the
Fig. 3
Point sources on the
LED light source. Red dots
show sources contributing full
intensity, orange show sources
contributing 50% intensity and
green show sources contributing
25% intensity
Experiments in Fluids (2019) 60:134
1 3
Page 7 of 16
134
optical axis, as indicated. From these, the direction vector
is given by
with the new initial conditions,
퐱
=[
x
+Δ
x
,
y
+Δ
y
,0
]
and
퐯
=[
v
x
,
v
y
,
v
z
]
.
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) con-
tribute 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.
3.3 Propagator
On exiting the test section, each ray has a coordinate in rela-
tion 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
:
(8a)
v
x
=
sin
휙
cos
휃
(8b)
v
y
=
sin
휃
(8c)
v
z
=
cos
휙
cos
휃
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 redi-
rected 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 proper
-
ties 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 gener
-
ated thousands of images per second. A less than 50% knife
(9)
⎡
⎢
⎢
⎢
⎣
x
�
휙
�
y
�
휃
�
⎤
⎥
⎥
⎥
⎦
=
⎡
⎢
⎢
⎢
⎣
1
d
00
0100
001
d
0001
⎤
⎥
⎥
⎥
⎦
⎡
⎢
⎢
⎢
⎣
x
휙
y
휃
⎤
⎥
⎥
⎥
⎦
.
(10)
⎡
⎢
⎢
⎢
⎣
x
�
휙
�
y
�
휃
�
⎤
⎥
⎥
⎥
⎦
=
⎡
⎢
⎢
⎢
⎣
1000
0
n
i
∕
n
f
00
0010
000
n
i
∕
n
f
⎤
⎥
⎥
⎥
⎦
⎡
⎢
⎢
⎢
⎣
x
휙
y
휃
⎤
⎥
⎥
⎥
⎦
.
(11)
⎡
⎢
⎢
⎢
⎣
x
�
휙
�
y
�
휃
�
⎤
⎥
⎥
⎥
⎦
=
⎡
⎢
⎢
⎢
⎣
1000
0100
−
1
∕
f
010
00
−
1
∕
f
1
⎤
⎥
⎥
⎥
⎦
⎡
⎢
⎢
⎢
⎣
x
휙
y
휃
⎤
⎥
⎥
⎥
⎦
.
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