Diode laser absorption tomography using data compression techniques

choruspillowUrban and Civil

Nov 29, 2013 (4 years and 7 months ago)


Diode laser absorption tomography using data compression

Chad Lindstrom
, Chung-Jen Tam
, Ryan Givens
, Doug Davis
, and Skip Williams

Air Force Research Laboratory, Propulsion Directorate, 1950 Fifth Street, WPAFB, OH 45433;
Taitech, Inc., 1430 Oak Court, Suite 301, Beavercreek OH 45430;
Air Force Institute of Technology, 2950 Hobson Way, WPAFB, OH 45433
Tunable diode laser absorption spectroscopy (TDLAS) shows promise for in situ monitoring in high-speed flows.
However, the dynamic nature of typical flows of supersonic combustors, gas turbine engines and augmenters can also
lead to inhomogenities that cannot be captured by a single line-of-sight TDLAS measurement. Instead, multiple
measurements varied over several spatial locations need to be made. In the current study, shock train structure in the
isolator section of the Research Cell 18 supersonic combustion facility at Wright-Patterson AFB is measured. Although
only two view angles are available for measurement, multiple absorption features along with a priori computational fluid
dynamics (CFD) simulations enable estimates of two dimensional flow features to be formed. Vector quantization/k-
means data clustering is used to identify key flow features from the temporal history of the raw sinograms. Through the
use of multiple absorption features that are measured nearly simultaneously, an approximate two-dimensional image can
be formed. This image can be further refined through the use of an optimal set of basis functions that can be derived
from a set of CFD simulations that describes the flow shapes.
Keywords: Absorption Spectroscopy, Tomography, Diode Laser, TDLAS, Supersonic Flows, Supersonic Combustion

In the most advanced space launch system in use today, the space shuttle, about half the launch weight is the liquid
oxygen and solid oxidizer it must carry along its way to orbit. Because jet engines need oxygen to burn fuel, and the
upper atmosphere does not have enough of it to sustain combustion, building an aircraft that can soar from a runway to
outer space requires rocket propulsion. However, propulsion systems such as a supersonic combustion ramjet, known as
a scramjet, which can scoop oxygen from the atmosphere during ascent are beginning to open up new opportunities for
access to space. The key engine components of a high speed (> Mach 4), supersonic combustion ramjet (scramjet)
propulsion system are an inlet, isolator, combustor, exhaust nozzle
1, 2
. During operation, the vehicle forebody and inlet
provide the air compression for the engine. As a result of compression, the air undergoes a reduction in Mach number
and an increase in pressure and temperature as it passes through shock waves at the forebody and inlet. The isolator in a
scramjet is a critical component. It allows a supersonic flow to adjust to a static back-pressure higher than the inlet static
pressure. The combustor accepts the airflow and fuel is added to the air to promote efficient fuel–air mixing at several
points along its length. The expansion system, consisting of the exhaust nozzle and vehicle aftbody, controls the
expansion of the highpressure, high-temperature gas mixture to produce net thrust.
In a conventional ramjet, the incoming supersonic airflow is slowed to subsonic speeds by multiple shock waves, created
by back-pressuring the engine. As fuel is added to the subsonic airflow, the mixture combusts, and exhaust gases
accelerate through a narrow throat, or mechanical choke, to supersonic speeds. By contrast, the airflow in a dual-mode
scramjet remains supersonic throughout the combustion process and does not require a choking mechanism
. Rather,
when the combustion process begins to separate the boundary layer, a precombustion shock forms in the isolator. This
thermal choking phenomenon is of great importance in an isolator in dual-mode ramjet/scramjet combustor. The isolator
also enables the combustor to achieve the required heat release and handle the induced rise in combustor pressure
without creating a condition called inlet unstart, in which shock waves prevent airflow from entering the isolator. Inlet
unstart occurs if the heat released by the combustion process becomes greater than a critical value, the Mach number in
Computational Imaging VI, edited by Charles A. Bouman, Eric L. Miller, Ilya Pollak,
Proc. of SPIE-IS&T Electronic Imaging, SPIE Vol. 6814, 68140W, © 2008 SPIE-IS&T ∙ 0277-786X/08/$18
SPIE-IS&T/ Vol. 6814 68140W-1

the engine falls to unity and the flow becomes choked
. This choked flow may, in turn, cause a normal shock or a
pseudo-shock to form at the engine inlet. This process, known as unstart, causes large amounts of drag and radically
reduces performance of the engine at high flight Mach numbers.
Therefore dynamics of the precombustion shock train such as its location, amplitude and frequency characteristics of its
instability, and the variation of the upstream and downstream static pressure ratios and Mach numbers figure
prominently in any overall engine design. Presently, computational fluid dynamics (CFD) modeling is unable to fully
predict theses variables. Ultimately, CFD is used for the design of hardware that can be tested in ground facilities. In
most ground test facilities, wall measurements (e.g., pressure, temperature, and heat flux) dominate the instrumentation
suite. If in-stream information (usually pitot pressure) is available, it is usually sparse and is generally available only at
the inflow and outflow planes of the test article. While valuable for various analyses, wall measurements provide little
or no detailed descriptions of the local state properties within the device. The proposed effort is intended to address
some of these deficiencies using laser-based instrumentation.
Specifically, this paper discusses the development of in-situ performance diagnostics using diode laser sensors for
ground test applications where the determination of species concentration, temperature, and velocities in supersonic
flows is required. In order to accomplish this task, a combination of well-controlled experiments and computational
simulations are used to verify the results. In the experiments, tunable diode laser absorption spectroscopy (TDLAS) is
used to characterize the shock structure of a vitiated supersonic flow in the Research Cell 18 supersonic combustion
facility at WPAFB. TDLAS is a relatively mature technique for making in-laboratory measurements of gas species
concentrations such as water (H
O) (the focus in the current study) and static temperature in homogenous environments.
However, application of this technique for monitoring high-speed flows such as afterburner exhaust plumes of jet
engines, supersonic combustors, and gas turbines is still in its infancy. Here the dynamic nature of the flow and spatial
variations present challenges. For example, the supersonic flow introduces scintillation effects because of the high-
speed turbulence as discussed in one of our previous works
. The dynamic nature of the flow means that measurements
must be made quickly and repeatedly in order to understand the flow environment. Spatial inhomogenities must be
addressed as well because TDLAS is a line of sight technique (LOS). Most attempts at measurements in supersonic
flows have focused to a large extent either on a single path averaged measurement or the use of a single line of sight
scanned across a supersonic flow
In this study, 16 LOS measurements are used to acquire time varying spatially resolved data to enable both spatial and
temporal variations in the flow to be observed. Specifically, the temporal and spatial variations of a back pressure
induced shock train structure in the isolator section of the WPAFB supersonic combustion test facility are quantitatively
observed. Although two view angles significantly limits our ability to reconstruct the flow field in the spanwise
direction, three principle reasons motivated this choice to test our system capabilities. First, the shock structure is
largely quasi-one dimensional so that the principle features are probed along the axial direction with excellent spatial
resolution. Second, spanwise information can be obtained from computational fluid dynamics simulations that can be
used as basis functions to fit to the measured data. Finally, this problem constitutes an important fundamental problem
in the development of supersonic combustion engines as detailed below as well as having impacts on the development of
an in-flight mass flux sensor based on TDLAS.
In tunable diode laser absorption (TDLAS), a single-mode diode laser beam is projected by a send fiber optic assembly
across the gas flow and is collected by a receive fiber assembly with the intensity monitored by a photodiode. The laser
frequency is scanned over the entire absorption feature to enable the baseline transmitted signal without absorbance, I
to be determined. The Beer-Lambert law can be expressed as:


where α
is the absorption coefficient at a particular frequency, ν, and position along the path, s. For uniform paths, this
expression simplifies to the expression on the right hand side of Equation (1) where the absorption coefficient has been
explicitly written in terms of the cross section, σ
, the concentration of the absorbing species in the quantum-state(s)
probed, N
, and the path length, L. Unlike conventional X-ray based
( )
σα −⎯⎯⎯⎯ →⎯

−= exp)(exp
N Constant
SPIE-IS&T/ Vol. 6814 68140W-2
2 0.08-
fl Post-shock
Pre-shock (lower PT)
Wavenumber (cm)
SPIE-IS&T/ Vol. 6814 68140W-3

tomography which involves ionizing radiation, TDLAS uses single transitions between distinct quantum states and is
therefore sensitive to pressure, temperature, and flow velocity through the spectroscopic lineshape (voigt) lineshape in
addition to the species density of the medium. This sensitivity to the measurement environment can be seen in figure 2
which shows two measurements conducted during the course of the experiment. In the first case, a measurement has
been made before the shock front where the density, pressure, and temperature are lower in the flow. As can be readily
seen the area and the peak width are much less than the measurement made in the flow after the shock. The peak area is
directly proportional to the temperature and the species concentration as can be inferred from equation 1. The width of
the peak is proportional to two factors. The first factor is Doppler broadening which occurs because the thermal
distribution of the velocity of the molecules in the gas. This is independent of the states involved in the transition and
takes the form of a Gaussian lineshape with a half width at half maximum of:

where ν is the line center frequency in wavenumbers, T is the temperature of the gas in Kelvin, and m is the molar mass
in grams
. Collisional broadening is dependent on the quantum mechanical states involved as well as the types of gas
molecules involved in the collisions. So, a different Lorentzian broadening parameter for each transition must be
determined for each type of collision partner over the relevant range of conditions for the gas mixture being probed. A
molar based addition of the broadening coefficients then yields the total collisional broadening coefficient for that gas
mixture. Because the total lineshape involves both a Doppler broadened component as well as a Lorentzian component,
the composite lineshape is a voigt profile
. A voigt profile is the convolution of a Lorentzian lineshape with a Gaussian
and can be computed through the standard Humlíček algorithm
. Because multiple absorption features are often present
in the spectra such as figure 1, it is necessary to use multiple voigt peaks to fit the spectrum and obtain accurate
integrated peak areas. As will be discussed later, large data sets are generated because the submillisecond residence
times of the supersonic flow necessitates kHz scanning frequencies. In order to observe the these dynamics, an
automated peak fitting routine has been developed using Igor Pro to enable fitting of approximately 300,000 spectra in 4
hours using a standard PC. This routine currently uses a constrained free fit of the voigt parameters to determine the
areas of the absorption features.
Although it is possible to determine the temperature directly from the linewidth of the curve fit to the absorption feature,
it is not always possible because the broadening parameters may be unknown or a tomographic method is used that only
determines the frequency integrated absorption coefficient. In this case, it is necessary to use the measured areas of the

(a) (b)
Fig. 1. The left figure (a) shows an absorption spectrum for the line with a ground state energy of 1045 cm
for both for the
flow before the shock front and after the shock front. P1 and P2 represent two different absorption features in the
spectrum which was acquired in each case with ~100 averages for each spectrum. The right figure (b) shows a
Boltzmann plot for one line of sight for data cluster 10 taken after the shock front. The reported errors only take into
account the error due to peak fitting.

absorption features such as those in figure 1 to determine the static temperature and molecular density. Temperature
determination is done by plotting the logarithm of the area of the absorption feature normalized by its cross-section and
degeneracy factor versus the initial state energy divided by k
. The resulting plot, called a Boltzmann plot, should
consist of points lying on a line whose temperature is given by the inverse of the slope if the temperature and density are
uniform along the line of sight. In addition, the intercept of the line is directly related to the molecular density of
absorbing molecules. If the temperature or density is non-uniform along the line of sight, then the resulting curve will
no longer be linear and is an indication of these variations. The right figure in figure 1 shows the Boltzmann plot for one
line of sight in cluster 11 that will be used as an example in the latter half of the paper. It is nearly linear indicating that
the temperature gradients along the line of sight are not large.
Tomography is well-known for its use in medical imaging problems such as X-Ray CAT scans and magnetic resonance
imaging (MRI). These applications use filtered backprojection, but another approach, Algebraic Reconstruction (ART),
is used here because of its suitability
. Figure 2(a) shows a diagram of tomographic reconstruction with the lines of sight
p(s,θ), that are a function of the indicated angle, θ, and the perpendicular distance from the origin, s. The line integrals
for TDLAS can be determined from Equation 1:


where p(s,θ)
is the frequency integrated line integral on the i
written in a parameterized form and α(x,y) is the
frequency integrated absorption coefficient as a function of position in the x-y plane. These line integrals can be made
discrete as shown in Figure 3(b). A system of linear equations results between the projections and the absorption


Equation 4 can be solved by least squares methods for α. The number of projections needed scales with the number of
pixels so M~ 10,000 for a 100x100 region. Clearly, the problem with conventional tomography is that a large number of
beam paths are required to obtain sufficient resolution. However, such an approach does not take advantage of a priori
knowledge that is known about the flow. The key to reducing the number of beam paths is to expand the absorption
function in Equation 4 in terms of a reduced number of basis functions that span over many pixels
. The result can then
be placed into Equation 3 and now only n<<N terms are involved:


where q
are the expansion coefficients for
, φ
(x,y) are the basis functions, and r
is the i
line integral of φ
Equation 5 can also be written in matrix form that is more suitable for computations as:
where R is the m x n projection matrix that is computed by calculating the m path integrals for each of the n basis
functions. This matrix equation is then solved through least-squares methods to determine the expansion coefficient
vector q that can then be used to determine the absorption map for a particular transition. An important side aspect of
determining R is that its condition number will tell us the stability of the resulting linear system of equations for a
particular choice of basis functions and beam paths. Typically it is our goal to achieve conditions numbers below 10
when inverting equation 6 so as not to amplify errors significantly. Thus, suitable choice of basis functions can
dramatically reduce the number of projections required. A critical aspect is the sufficient representation of the flow by
the basis functions.
θ +




,,2,1 ,





),( φαφαα
SPIE-IS&T/ Vol. 6814 68140W-4


Wjj = area ABC/52
(a) (b)
SPIE-IS&T/ Vol. 6814 68140W-5

Creation of the basis functions can be accomplished through several means. Perhaps the easiest to understand is an ad
hoc determination of the basis functions for the flow field based on intuition and experience. For example, in the case of
the McKenna burner described in a previous study a simple three-ring model is found to be sufficient. Of course, the
weakness is that it will be difficult in some cases to use intuition because of complexity of the flow field, although it
should never be underestimated. In contrast, more complex cases require that a formalized procedure exists. In this
case computational fluid dynamics (CFD) simulations or high resolution optical diagnostic data such as Planar Laser
Induced Fluorescence or Rayleigh scattering that give information on both the relative temperature and species
concentration are used to initially generate a training set of absorption images or maps. It is also possible to enhance this
training set through the use of intuitive scaling parameters to take into account simple experimental variations. For
example, the one dimensional measurements made later provide a clear indication of how the two-dimensional CFD
images should be scaled to form a training set. The next step is then to use data clustering/vector quantization to
subdivide the training set into distinct features
. For example, in the case of the isolator, a straightforward approach is
to use the uniform region before the shock, the shock train region, and the post-shock region. The mean value images
can be determined for each region and the projection matrix calculated for each mean image. Then as TDLAS data is
collected, the residuals from the least squares fits to the projection matrices for the mean images are used to determine
the region of the flow that the data is being collected in. Notice that classification of the flow fields can be rather
arbitrary; for example it may be useful to compare CFD simulations with different turbulence models, i.e., RANS vs.
LES derived basis sets could be compared to determine which computation best represented the data In addition, a
formalized technique called the LBG (Linde, Buzo, and Gray) algorithm (also known as the k-means method) can be
used if intuition is not sufficient
. Variations within each data cluster can be further modeled either through the addition
of other basis functions that can allow for such things as baseline adjustment (i.e. a uniform basis function) or through
the use of principle components analysis as detailed in a previous publication

Here we briefly give an overview of the four major aspects needed to understand the experiment. These areas consist of
the supersonic combustion facility, experimental layout, control/data acquisition electronics, and spectroscopic
transitions. We focus first on the facility and then move to the latter three areas.

3.1 Facility
AFRL/RZA owns and operates two direct-connect supersonic combustion facilities. Figure 3 shows a schematic of the
Research Cell 18 where the experiments discussed here were performed. This facility was designed to allow basic
studies of supersonic reacting flows using conventional and non-intrusive diagnostic techniques. It consists of a natural

Figure 2. (a) The indicated geometry for tomography where multiple projections p(s,θ) allow reconstruction of the
shaded region (b) Illustration of breaking the region into a square grid

gas fueled vitiator, transition flange, interchangeable facility nozzle (Mach-1.8 and 2.2 currently available), modular
isolator, modular combustor, instrument section and calorimeter. The rig is mounted to a thrust stand capable of
measuring thrust of up to 2000 lbf. The test rig is provided with continuous flow of up to 30 lbm/sec at 750 psia at a
temperature of up to 1660 °R, and exhausts into a 3.5 psia continuous flow exhauster system. The facility can simulate
flight conditions from Mach 3.5 to Mach 5 at flight dynamic pressures of up to 2000 psf. For the current study, the
combustor was removed and replaced with a large valve that was used to simulate the pressure rise associated with
combustion. Figure 3 also shows a side view of the I-1 isolator. It is fitted with 1” high x 3.4” wide optical access ports
designed to hold 1.7 cm (0.68 in) in thick quartz windows or steel inserts. The quartz windows are wedged at 5º on the
outer surface to prevent etaloning. The nozzle and isolator are instrumented with static pressure taps and thermocouple
with a frequency response of approximately 10 Hz. No difference in isolator performance has been noted between the
using the quartz windows or the thermal barrier coated (TBC) steel inserts despite the differences in surface roughness
between the quartz and the TBC. A key enabling aspect provided by this facility is that the vitiation which is used to
increase the flow enthalpy to the correct level provides the water vapor being probed in this experiment. Although it is
possible to measure oxygen, it is difficult to achieve the excellent signal to noise ratio that can be obtained with water
absorption features. This facility has been described in detail elsewhere
5, 14, 15

3.2 Experimental Layout and Control Electronics
A top view of the experimental layout is shown in figure 4. Three lasers all operating in a different spectral region are
time-multiplexed with a single period lasting 1 millisecond and combined using a Newport 1x4 combiner. The
combined beam is then split using a 10/90 Newport splitter with the 10% beam going to a 50/50 splitter whose output is
routed to either a SiO
etalon with a 2 GHz free spectral range that monitors the frequency spectrum of the lasers or to a
ThorLabs photodetector that is used as a reference. The remaining 90% of the beam goes another 50/50 splitter that then
feeds into the FC/APC bulkhead connectors located on each purge box. These then connect to two 1x 8 Newport fiber
optic splitters that are then collimated using ThorLabs collimators. The collimators are mounted on a precision drilled
baseplate to enable tight beam placement. However, because of the physical limitations of the optical mounts the ray
paths are slightly angled when passing through the probe region as indicated in figure 4. During setup, clipping on the
upstream window edge required that the suboptimal pattern shown in 4(b) be used. The straight-across pattern was
chosen because of the quasi-one dimensional shape of the shock structure in the mid-height measurement plane. In
addition, it allows for relatively straight-forward interpretation of the data. Top to bottom ray paths were not available at
the current time of the study, however, they may become the focus of further studies if more detailed information on the
spanwise structure of the shock is required. Two additional diagonal rays that are 19.3° from the spanwise direction
were also measured at the same time and whose primary purpose is to provide velocity measurements from the Doppler
shift of the peak position.
The raw transmitted laser power is measured simultaneously on all 16 paths through the use of in-house custom designed
InGaAs photodiode arrays. These custom designed arrays amplify the signal from 2mm diameter FGA21 photodiodes
from ThorLabs using a Texas Instruments OPA380 transimpedance amplifier. The signal is then digitized into 14 bits
using a NI PXI-6133 (32 MSamples onboard memory) S Series Multifunction Digital Acquisition (DAQ) Device
sampling at a rate of 2.5 MSamples/s. Because all the channels are all simultaneously sampled at the same time
(18x2.5x2 = 90 MB/s), the PXI bus in the PXI chassis cannot sustain continuous data transfer to the hard drive in the

Figure 3. (a) Schematic of the supersonic combustion test facility (b) Side view of the I-1 isolator with 3.4” wide
optical access ports.
SPIE-IS&T/ Vol. 6814 68140W-6
]rJJ p.
Axial distance from upstream edge ofwindow (cm)
SPIE-IS&T/ Vol. 6814 68140W-7

current configuration. Instead a short segment of data (0.2 s) is captured, processed, and saved to the hard drive. This
process takes approximately 5 s and then another 0.2 s of data is captured, processed, and saved. Although the shock is
unsteady it has been found that pushing the shock slowly from the rear of the isolator to the front enabled the key shock
features to be measured. Because of the unsteady nature of the shock, a NI Intel Core Duo 2.0 GHz Processor with 2GB
of DDR2 host RAM is used to average 200 spectra, determine peak areas through trapezoidal integration, and then create
a preliminary image of the flow field. This provided feedback as to the success of the data runs. The data
acquisition/processing software was developed in LabWindows to allow a scalable architecture for developing TDLAS
based tomography. It achieves this because 8 additional channels can be simply added by inserting an additional PXI-
6133 DAQ board, changing software settings, and interfacing directly to an in-house detector array. Of course, data
acquisition rates will vary based on the channel count however at least 2 seconds of continuous data can be acquired
because of the large on-board memory. Currently, it is expandable to up to 80 channels. It is anticipated that in the near
future significant speed up in the data acquisition process can be achieved through the use of either a RAID0 disk array
or DDR RAM based hard drive using a PCI ExpressCard slot as well as software upgrades. The next level of COTS
based solutions is PCI Express which may be able to allow up to 80 channels to disk at digitizer rates of 2.5 MS/s.
A final aspect of the experimental setup is the metal purge boxes. They serve multiple purposes. The sidewalls of the
isolator can reach 600-700 K during a test run and the purge boxes shield the detector arrays from the heat. In addition,

(a) (b)
Fig. 4. (a) A top view of the experimental setup with an overlay depicting the shock train structure in the isolator is shown.
The measurement plane occurs at mid-height in the isolator. For clarity, only 8 beams (red lines) going from the
bottom purge box to the top are shown. However, another 8 lines of sight (LOS) from the top to the bottom box
are also used. The indicated flow direction is what will be used in subsequent figures. (b) The actual LOS in
robe re
ion after correctin
for refraction at the windows as well an absor
tion ma

enerated fro

nitrogen gas is bubbled through liquid nitrogen providing a cool, dry nitrogen gas flow used to displace the lab air during
the test run with a ~ 100:1 reduction in humidity along the beam paths. The purge gas displaced essentially all of the
water outside of the flow environment inside the isolator. This proved to be important, because, as can be seen in figure
4, relatively long path lengths inside the purge boxes are required to get the proper beam spacing before entering the
isolator. The purge boxes are also mounted on a sliding base plate that moves with the test rig as it thermally expands to
ensure that the same region of the flow is being monitored.
3.3 Spectroscopic Transitions
Water vapor has several strong rovibrational spectroscopic transitions from the visible to the mid-IR. The availability of
compact, spectroscopic quality lasers at telecom wavelengths, 1250-1650 nm, have made the 2ν
, 2ν
, and ν
+ ν

vibrational bands particularly attractive for diagnostic applications as evidenced by several recent studies. The ν
corresponds to the symmetric stretch vibration and the ν
mode is the asymmetric stretch vibration in water. Therefore,
the transitions accessible in the spectral region are vibrational overtone transitions and combination bands which are
typically much weaker than the fundamentals in the mid-IR region. However, the abundance of water vapor as a
primary combustion product and the relatively large linestrengths for these transitions enable these transitions to be used
with good sensitivity for many applications. Single mode fiber-coupled, DFB lasers are used to generate near-infrared
radiation in the laboratory to measure H
O line intensities near 1348-1397 nm corresponding to the spectral region
ranging from 7155 to 7417 cm
which is of interest for the in situ monitoring of H
O in a variety of environments to
include combustion studies and atmospheric monitoring. Principally, four transitions of the ν

band have been used
in this study. The line parameters are taken from the HITRAN database and are listed in Table 1
. The line intensities
at different temperatures are derived using the same procedures and partition function expressions as defined in
HITRAN. Line broadening parameters were taken from Delaye, et al, for CO
, O
, N
, and H
O when determining the
pressure from the line widths

Table 1: Spectral line parameters for transitions used in this study. All transitions originate in the ground (0,0,0) vibrational
state (connoted with “) and terminate (1,0,1) vibrational state (connoted with ‘).
λ (nm) ν
S at 296 K
E" (cm
' J",K

1391.67 7185.6 7.947e-22 1045.06 6,6,0 6,6,1
1392.53 7181.1558 1.505e-20 137.17 2,0,2 3,0,3
1392.81 7179.752 2.299e-22 1216.19 7,6,2 7,6,1
1396.37 7161.4101 1.174E-20 224.84 3,1,3 4,1,4


The first step in the data analysis is to determine the frequency integrated peak areas for each line of sight and absorption
feature. This is done by splitting the time multiplexed data in the three spectral regions and then using the reference path
to determine the absorption spectrum. The absorption features in each spectrum are then fit using a voigt profile as
described in section 2. This is all accomplished within an automated function within Igor Pro. The output of this step
yields a time history of frequency integrated absorption as shown in figure 5. As can be seen in figure 5, the axial
position of the shock front can shift at up to 100 Hz due to instability in how the shock anchors in the duct. It also
indicates that the amplitude of these oscillations are approximately 4 cm. Already this is a significant undertaking,
however, several more steps in the data processing are required before accurate estimates of the shock shape can be
SPIE-IS&T/ Vol. 6814 68140W-8

25x1 0
E 20
= 15
2- 10-
Axial Distance Across Wndow (cm)
25x1 0
Axial Distance Across Window (cm)
Axial Distance Across Window (cm)
Freq Integrated Absorbance (cm1)
Time = 140 ms
Member of Cluster II

SPIE-IS&T/ Vol. 6814 68140W-9

4.1 Data clustering of time varying data
Because of the unsteady nature of the shock train, it is necessary to use data clustering on the collected data to average
the key features that are observed. The use of averaging is important because although the single shot spectra are usually
quite good, turbulence noise does exist in the flow and generates outliers (light and dark blue spots in figure 5) that can
bias individual results. In addition, current state-of-the-art computational fluid dynamics (CFD) simulations typically
only employ steady-state models of the flow which should more closely correspond to averaged data. Due to the
irregular nature, and large amplitude of the oscillations, simple time averaging techniques should not be applied to this
problem because doing so significantly broadens, and hence, misrepresents the time-varying shcok structure. The
approach that is used here is the k-means data clustering which assigns members to clusters based on their distances
from the cluster center using a distance metric such as the standard Euclidean sum of squares
. It does this through an
iterative process that converges to a local minimum. Of course, it is necessary to ensure through varying the starting
conditions as well as the distance metric that it is converging to a correct solution. Conveniently, Igor Pro provides a
software routine that allows different starting conditions, number of classes, and the distance metric to be varied. A key
input to the algorithm is the number of classes. It was found that approximately 20 classes did an excellent job of
describing the major flow features of the data run consisting of 6000 members used in the analysis conducted here.
These flow features should be taken as representative of the underlying shock structure in the flow, however, they are
not eigenmodes or coherent structures and they will depend to some degree on how the k-means algorithm is set up.
Figure 6 shows the 20 classes that have been derived from the 6000 members for the 7185.6 cm
absorption line.
Interestingly, if these classes are plotted versus a weighted average of backpressure measured near the backpressure
valve these classes show the shock position response to backpressure on average. The reason a weighted average must
be used is because the pressure will vary as function time, but the pressure data that was currently available was only

(a) (b)
Fig. 5. (a) The figure shows a waterfall plot of the frequency integrated absorbance along each path (except for angled beams)
versus time that depicts the shock front oscillations. The solid black line is meant as guide for the reader as to the shock
front position. Because of the slight diagonal nature of some of beams, the x dimension is the average axial position of
each beam. (b) Two members of cluster 11 at different times are shown to indicate the variation that occurs.
C1 PB=l7.69PSIA
02u--- B1797 PSIA07u-- B1985 PSIA
Freq. mt. Absorp
Freq. mt.
Axial Distance Across Window (cm)
E 20I C8PBl9.93PSIA
I 09u-- B1995 PSIA
E 20
B20°5 PSIA
- 1 CllU--PB=2O.I2PSIA
Axial Distance Across Window (cm)
Axial Distance Across Window (cm)
Axial Distance Across Window (cm)
° PB0
] °
a) _____________ I
2468 2468
Axial Distance Across Window (cm)
Axial Distance Across Window (cm)
Absorption (cm1)
C 01 C
Absorption (cm1
C 01 C
C20*-- B21°6 PSIA
0 0
2468 2468
Axial Distance Across Window (cm)
Axial Distance Across Window (cm)
SPIE-IS&T/ Vol. 6814 68140W-10

Fig.6. The data clusters (C1-C20) derived from a data set with 6000 members versus backpressure derived from the frequency
integrated absorbance of the 1045 cm
line. The axial location of the shock correlates with the backpressure, however,
some discre
ancies do exist and are likel
due to inaccuracies in estimatin
the back
ressure as ex
lained in the text.
I II I01I I I
Axial Distance Across Window (cm)
Axial Distance Across Window (cm)
SPIE-IS&T/ Vol. 6814 68140W-11

sampled at 10 Hz. Because of this mismatch in sampling rates, it is not possible to detect the slight pressure variations
that could be influencing the shock location complete accuracy. This is the reason for example that C12 does not come
before C10 and C11 for example. However, the general trend is quite clear.
4.2 Uniform spanwise estimates of T and P
The next level of analysis assumes no variation in the flow properties in the spanwise direction. The reason for this
analysis is two fold. First, physical intuition as well as CFD shows that the predominant variations in pressure and
temperature occur in the axial direction. Second, this approach yields a simple, intuitive model of the flow. All that is
required is to take the same classes as computed in figure 6 and use the class membership to generate the class averages
for the remaining absorption lines. Then Boltzmann plots similar to figure 1(b) are generated for each line of sight to
determine the temperature and water concentration. Once the temperature is known, the pressure can be determined
from the collision broadening coefficient through a voigt curvefit. The results are shown in figure 7 along with a

comparison with the CFD for a similar condition. The most obvious aspect that appears is that the axial node spacing in
the CFD shock structure is incorrect, i.e., the first shock reflection is quite apparent in the data pressure and temperature
curves but absent in the CFD curves. This is perhaps not to surprising because prior studies that attempted to predict the
shock anchoring point in the isolator encountered difficulty using only the CFD
. A likely explanation for this will
become apparent after looking at the 2D reconstructed images versus the CFD. Another important aspect is that in the
cases of C1 (6.4 psia vs 6.2 psia) and C20 (20.5 psia vs 19.4 psia) the side wall pressure taps are in good agreement with
the pressure determined from the voigt curvefit. Although no measurements of the static temperature are available
because of the difficulty of currently employed apporaches, the averaged static temperature measured using TDLAS is in
agreement the calculated properties of the freestream flow for C1.
4.3 Basis function analysis
Figure 8 shows a pseudo-color plot of the CFD temperature profile within the plane of measurement for a back pressure
of 20.9 psia. The image is typical for this isolator for a variety of backpressure conditions. This figure is shown to
illustrate two important constraints that are present in this problem. The first is that a mirror plane of symmetry exists in
the measurement plane at a spanwise distance of 5 cm. Although this symmetry can be broken by differing forces in the
spanwise direction such as different wall roughness, it is a useful first approximation to make when performing the
inversion. Second, when the shock separates, the thermal boundary region next to the wall reaches a static temperature
very near the total temperature of the gas which is equivalent to the temperature attained if the velocity of the flow is
brought to a stop in an isentropic manner. In addition to these constraints, any solution must obey the equation of state

(a) (b)
Fig. 7 (a) The static temperature and (b) static pressure assuming uniformity along the spanwise direction for C11. This is
constrasted with the CFD solution that has been averaged in the same manner as the TDLAS measurements.
Degrees (K)
Axial Distance (cm)
SPIE-IS&T/ Vol. 6814 68140W-12

of the gas, and, because the water vapor is created upstream from the facility nozzle, the water vapor is a constant molar
fraction of the total density of the gas. These constraints are shown in table 2, and any reconstructed image should
satisfy these constraints. In addition to these constraints, the ideal solution should obey the fluid flow equations (i.e. the
Navier-Stokes equations). Although it is very difficult to derive or even compute a full solution for this problem, an
approximation known as the Reynolds Averaged Navier-Stokes Equations can be used
. In this approximation, a
steady-state solution is computed using a model for the flow turbulence for the given condition. This is shown in figure
8 as computed with CFD++
. However, this model has several shortcomings. Clearly, it will not capture the dynamics
of the flow because it is a steady-state solution. In addition, the solution varies with the turbulence model used.
Improper accounting for shock-turbulent boundary layer interactions could be reason for example that the CFD++
solution overestimates the shock reflection spacing in figure 7. To account for the for inaccuracies in the shape
information presented by CFD it is possible to recreate the correct node separation by photographically scaling the image
by a scale factor, β. This can be straightforwardly implemented in the region where the shock nodes exist because the
window that we are observing the flow with is much smaller than the entire isolator as shown in figure 6 by the black
rectangle. However, in the spanwise direction this could lead to larger than anticipated thermal boundary layers because
as can be seen in figure 6 the thermal boundary layer scales with the axial distance from the shock front. To correct for
this an additional scaling must be performed in this direction. However, a problem exists because the spanwise
dimension is fixed at approximately 10 cm so a simple linear scaling is not possible. Instead a non-linear scaling
procedure must be used. One approach is simply to scale the interior of the image linearly and leave the remaining
boundary layer. However, this will introduce a sudden discontinuity that is unphysical. Instead this transition is
modeled as a soft transition using the following scaling method which is one possible approach:

Fig.8. Plot of the CFD temperature in the isolator at the measurement plane for a backpressure of 20.9 psia. The black
rectangle illustrates the portion of the flow that is actually sampled by the window. The white line is where the
spanwise slice is taken for figure 9.

Table 2: List of simple constraints that are imposed on the reconstruction of the image


Mirror plane at 5 cm in spanwise direction

= T
(Boundary temperature = Total temperature after flow separation)

P = NkT/V (Ideal Gas Equation of State)

= χN (Water is a constant mole fraction)
Core Flow Region
  - scaled
Spanwise Distance (cm)
SPIE-IS&T/ Vol. 6814 68140W-13

{ }
flow core :)('


( 7)
where η is the scaling coefficient, L
is the original boundary layer length, y
is the spanwise center coordinate, y

is the lower limit for the upper boundary layer, y is the spanwise coordinate, L is the boundary layer length, and r is an
exponent that adjusts the rate at which the transistion from the central core flow to the original boundary layer occurs.
Notice that if y = L
then y’ becomes equal to the compressed spanwise coordinate and if y = 0 then y’ = 0 so that the
original spanwise dimensions are retained but yet a much smoother transition is obtained from the core flow to the
boundary layer. An application of this to a slice in the spanwise direction represented by the white line in figure 8 is
shown in figure 9 with η=1.3, r = 1, and L
= 2.6 cm. Clearly this method preserves the relative scaling between the
boundary layer and the core flow but it will change the gradients between them. The next aspect that will be addressed in
the fit of the CFD basis functions to the data is the amplitude of the features of the absorption maps. Relative amplitude
variations is accounted for by solving equation 5 and determining the expansion coefficient. However,, it is also
desirable to take into account a baseline offset. This baseline variation is accomplished through adding a uniform (i.e.
constant value at all locations) basis function to equation 5. Finally, the shock is moving inside the window so the CFD
basis function depends on the shock front position relative to the window. To model this the window in figure 8 is
allowed to slide to various positions in the CFD simulations a function of its upstream edge position, x. Ultimately then
we can write the basis function fit to the flow field as:

is the frequency integrated absorption coefficient for the absorption feature with E’’ as its ground state
energy, x is the axial dimension, y is the spanwise dimension, β is the axial dimension scaling coefficient, η is the
spanwise scaling coefficient, r is the spanwise exponent, and x
is the upstream edge position of the window. Although
equation 8 is a function of r we have found little dependence on r with the optimal value for η changing from 1.34 to 1.4
when r changes from 1 to 2 for the current case. In addition, we define L
as the point at which the temperature has
dropped by 90% from the wall temperature to its minimum. The dashed lines indicate where this would be in figure 9.
So we must optimize equation 8 over all 16 lines of sight as well as the four different absorption features in order to
obtain an optimal shape as a function of β, η, q
, and q
. It is important to realize that in order to do this optimization

Fig.9. Plot of the scaled CFD temperature in the isolator at the white line in figure 6. In this case η=1.3, L
2.6 cm, and r =1.

' I
400 500 6000.51.01.5
Temperature (K)
Pressure (atm)
e (cm)
nwise Distan
anwise Distar
I I I I0
1i 0
Axial Distance Across Window (cm)
Axial Distance Across Window (cm)
Original CFD


400 500 600
Temperature (K)
Pressure (atm)
10- 10
Distance (cri
I\) -
o-r1I I0
6 66
Axial Distance Across Window (cm)Axial Distance Across Window (cm)
Reconstructed Images
SPIE-IS&T/ Vol. 6814 68140W-14

Fig.10. The upper plot shows the original hand aligned CFD images and the lower plot shows the reconstructed images
that have been fit to the experimental data.

first equation 8 must be computed then the projection integrals in equation 5 must be determined numerically which is
computationally intensive. Although it is possible to use optimization codes to carry out the analysis, a direct search
method that examines a grid of search points was employed in the initial analysis primarily because it will yield a
globally convergent solution (here we know the physically significant search regions). Because the primary sensitivity
of the measurements is in the axial flow direction, the optimization was first carried out for β. Here, the β was varied
and the subsequent linear system (i.e. equation 6) was solved for q
and q
. This yielded a value of β = 0.54 for cluster
11 (C11 in figure 6) and suggests that the CFD node spacing is off by almost 50%. For the next step in the analysis, the
axially scaled image with β =0.54 was optimized in the y-direction. This was carried out by adjusting η over a range of
values until an optimal value was found. It was found that the optimal fit was obtained when η = 1.34. The value
obtained here depends critically on L
and typically η will decrease as L
decreases. Although this might suggest that
we do not adjust the spanwise profile a subsequent analysis found that if the no adjustment in the spanwise dimension is
made then the boundary layer temperature will reach unphysically high values (i.e. T
~ 670 K in this case but it yields
>770 K if no adjustment is made). What this suggests is that qualitatively the shock is wider in the spanwise
direction than the computed by the CFD. Of course, because of the LOS used to probe the flow only limited sensitivity
to the spanwise direction exists so this should not be considered a firm confirmation that the shock is wider. Finally, the
images for each absorption feature can be fit using Boltzmann plots to determine an estimate of the temperature and
water concentration in the plane of measurement. Using the constraint of constant molar fraction, it is then possible to
use the ideal gas law to derive an estimate of the pressure. The resulting images can then be compared to either existing
measurements or the original unscaled CFD.
Figure 10 shows the results of the fit and compares the result to the original unscaled CFD that has been shifted so that
the first shock node aligns with the data. Clearly, the second shock node is not properly aligned. The most likely reason
for this can also be seen in the figure. The dark lines illustrate the shock angle. In the CFD image, the shock angle is
less than in the reconstructed image. The shock angle depends on a number of factors such as the turbulence as well as
the wall roughness that must be modeled so this could be indicating problems with these models. Interestingly, despite
the poor spacing between the shock nodes a point by point comparison between the CFD and the reconstructed image
shows for example that the maximum temperature in the first node is only 20 K hotter than the CFD (~ 4% difference).

To our knowledge this represents the first time (millisecond) and spatially (5 mm) resolved in-flow measurement of
static temperature, static pressure, and density of a supersonic isolator shock-train. In addition, this study represents one
of the largest simultaneous tunable diode laser absorption experiment that has been performed at similar supersonic
combustion facilities with 16 paths sampled over 3 spectral regions (48 spectra) at rates greater than 1 kHz. This work
lays the groundwork for understanding the effects of shock train structure on TDLAS based mass-flux sensors as well as
improving fundamental understanding of the dynamic flow physics of these shock structures. This study used a back
pressure valve to generate the shock train, but this approach can easily be applied to the study of precombsution shock
train amplitude and frequency instabilities resulting from combustion. Despite not having the complication the shock
train being produced by combustion, an apparent discrepancy regarding the axial node spacing has been found between
the CFD simulations and the measured data. This discrepancy is clearly established because of the good experimental
spatial resolution in the axial direction. However, a basis function/scaling approach was applied and can provide
additional insight into the shock structure in the spanwise direction as well as enabling comparison between the CFD
simulations and the measured values along the center of the flow. Specifically, it was shown that the CFD node spacing
is off by almost 50% and that qualitatively it is likely that the shock is wider than computed. A point by point
comparison between the CFD and the reconstructed image suggests the maximum temperature in the first node is only
20 K hotter than the CFD (~ 4% difference). In other words, the discrepancy between the CFD and measured values is
much less for the maximum pressure and temperature inside the first shock suggesting that the CFD is reproducing the
magnitudes in the shock train with reasonable accuracy. Also, it is important to note that despite the fact that the original
CFD solution provided a poor comparison with the raw data, the approach outlined here was able to utilize the pressure,
temperature and density contour information provided by the CFD to gain insight into the measured data. Furthermore,
the reconstructed images provided additional feedback in the sense that it provided information regarding how “the CFD
should have looked” to agree with the data.
This approach can be carried over to other parts of the supersonics combustion flowpath that involve complex flow
distributions such as the exhaust plane of the supersonic combustion facility at WPAFB. This region is of significant
SPIE-IS&T/ Vol. 6814 68140W-15

interest because it can allow estimates of the spatially resolved combustion efficiency. In addition, other techniques
other than CFD can be used such as flow sampling or PLIF imaging to help aid in development of the basis functions to
improve the accuracy of the technique. Finally, because it is always desirable to be free of any reliance on basis
functions it should be possible in the future to start the development of a system that could have up to 80 simultaneously
sampled paths built around the core data acquisition system developed during this measurement campaign.

This work was supported by the Air Force Office of Scientific Research (Dr. Michael Berman, Program
Manager) and C.D. Lindstrom was sponsored under a NRC Fellowship. The authors thank William Terry for his
great assistance in hardware design and setup. The authors also thank Professors F. C. Gouldin and W. Bailey
for helpful discussions.

[1] Andreadis, D., "Scramjets integrate air and space," Industrial Physicist, 10 (4), pp. 24-27 (2004).
[2] Curran, E. T., Heiser, W. H., and Pratt, D. T., "Fluid phenomena in scramjet combustion systems," Annual
Review of Fluid Mechanics, 28, pp. 323-360 (1996).
[3] Matsuo, K., Miyazato, Y., and Kim, H.-D., "Shock train and psuedo-shock phenomena in internal gas flows,"
Progress in Aerospace Sciences, 35 (1), pp. 33-100 (1999).
[4] Miyazato, Y., Masuda, M., Matsuo, K., Kashitani, M., and Yamaguchi, Y., "One-dimensional analysis of
thermal choking in case of heat addition in ducts," Journal of Thermal Science, 9 (3), pp. 224-229 (2000).
[5] Williams, S., Barone, D., Barhorst, T., Jackson, K., Lin, K.-C., Masterson, P., Zhao, Q., and Sappey, A., "Diode
Laser Diagnostics of High Speed Flows," 14
AIAA/AHI International Space Planes and Hypersonic Systems
and Technologies Conference AIAA 2006-7999(2006).
[6] Liu, J. T. C., Rieker, G. B., Jeffries, J. B., Gruber, M. R., Carter, C. D., Mathur, T., and Hanson, R. K., "Near-
infrared diode laser absorption diagnostic for temperature and water vapor in a scramjet combustor," Applied
Optics, 44 (31), pp. 6701-6711 (2005).
[7] Bernath, P. F., [Spectra of Atoms and Molecules], Second ed.: Oxford University Press. New York, NY, 2005.
[8] Humlíček, J., "Optimized Computation of the Voigt and Complex Probability Functions," J. Quant. Spectrosc.
Radiat. Transfer, 27 (4), pp. 437-444 (1982).
[9] Kak, A. C. and Slaney, M., [Principles of Computerized Tomographic Imaging]: Society of Industrial and
Applied Mathematics, 2001.
[10] Torniainen, E. D., Hinz, A. K., and Gouldin, F. C., "Tomographic Analysis of Unsteady, Reacting Flows:
Numerical Investigation," AIAA Journal, 36 (7), pp. 1270-1278 (1998).
[11] Jain, A. K., Murty, M. N., and Flynn, P. J., "Data Clustering: A Review," ACM Computing Surveys, 31 (3), pp.
264-322 (1999).
[12] Linde, Y., Buzo, A., and Gray, R. M., "An Algorithm for Vector Quantizer Design," IEEE Transactions on
Communications, 28 (1), pp. 84-95 (1980).
[13] Lindstrom, C., Tam, C.-J., Davis, D., Eklund, D., and Williams, S., "Diode Laser Absorption Tomography of
2D Supersonic Flow," 43rd AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit AIAA-2007-
[14] Gruber, M., Donbar, J., Jackson, K., Mathur, T., Baurle, R., Eklund, D., Smith, C., and . "Newly Developed
Direct-Connect High-Enthalpy Supersonic Combustion Research Facility," Journal of Propulsion and Power,
17, pp. 1296-1304 (2001).
[15] Mathur, T., Gruber, M., Jackson, K., Donbar, J., Donaldson, W., Jackson, T., and Billig, F., "Supersonic
Combustion Experiments with a Cavity-Based Fuel Injector," Journal of Propulsion and Power, 17, pp. 1305-
1312 ( 2001).
[16] Rothman, L. S., Jacquemart, D., Barbeb, A., Benner, D. C., M. Birkd, Browne, L. R., Carleerf, M. R., Jr., C. C.,
Chancea, K., Couderth, L. H., Danai, V., Devic, V. M., Flaudh, J.-M., Gamache, R. R., Goldman, A.,
Hartmannh, J.-M., Jucksl, K. W., Makim, A. G., Mandini, J.-Y., Massien, S. T., Orphalh, J., Perrinh, A.,
Rinslando, C. P., Smitho, M. A. H., Tennysonp, J., Tolchenovp, R. N., Tothe, R. A., Auweraf, J. V., Varanasiq,
P., and Wagner, G., "The HITRAN 2004 molecular spectroscopic database," Journal of Quantitative
Spectroscopy & Radiative Transfer, 96, pp. 139-204 (2005).
SPIE-IS&T/ Vol. 6814 68140W-16

[17] Delaye, C., Hartmann, J.-M., and Taine, J., "Calculated tabulations of H
O line broadening by H
O, N
, O
, and
at high temperature," Applied Optics, 28 (23), pp. 5080-5087 (1989).
[18] Lin, K.-C., Tam, C.-J., Jackson, K. R., Eklund, D. R., and Jackson, T. A., "Characterization of Shock Train
Structures inside Constant-Area Isolators of Model Scramjet Combustors," 44th AIAA Aerospace Sciences
Meeting and Exhibit AIAA 2006-0816(2006).
[19] Metacomp, "http://www.metacomptech.com/index.html,
" 2006.

SPIE-IS&T/ Vol. 6814 68140W-17