Benchmarking background oriented schlieren against interferometric measurement using open source tools

D. Kaganovich; L. A. Johnson; A. A. Mamonau; B. Hafizi

doi:10.1364/AO.406301

1. INTRODUCTION

Many applications of miniature gas jets require characterization of the gas or plasma density distribution. To this end, much research has been performed in a wide range of experimental studies spanning from small space thrusters to laser-generated plasma targets. In laser wakefield accelerators, for example, precise control over the gas and plasma density and its profile is critical to effective acceleration [1]. The flagship gas density measurement diagnostic is interferometry [2]. A typical setup of such a measurement involves a probe laser, beam splitters, mirrors, and imaging optics. In complex experimental configurations, when in situ gas or plasma density measurements are required, it might be quite difficult or even impossible to set up such an interferometer.

In this work, we benchmark a simple and robust diagnostic, background oriented schlieren (BOS), against interferometry. We made the benchmarking process more challenging by performing measurements of a flowing gas using an under-expanded jet nozzle. As a result, the gas density has the pattern associated with the well-known shock diamonds [3]. Shock diamonds have been the objects of previous BOS measurements [4,5]. In this work, we show a simple BOS calibration, provide an analytical model for refractive index calculation, translate this model into a numerical recipe, and compare the results with the interferometry measurements. Description of the principles of the BOS diagnostic can be found in [6].

2. EXPERIMENT

The experimental setup is depicted in Fig. 1. A 2 mm diameter, axisymmetric, under-expanded gas nozzle was placed in one arm of a Mach–Zehnder interferometer. The interferometer used a 532 nm, CW laser with a beam expander to cover the area of interest above the nozzle. The nozzle was imaged onto a CCD camera using a 55 mm macro lens. The camera was a Point Gray model FL3-GE-28S4M-C, with ${{1928}} \times {{1448}}$ pixels sensor. The magnification of the system was 0.54 (6.88 µm of real object per one pixel of the camera). The interference fringes were aligned perpendicular to the gas flow direction.

The BOS diagnostic was positioned perpendicular to the interferometer probe beam, as shown in Fig. 1. The back-illuminated random background was printed on a transparency using a commercial laser printer at 1200 dots per inch (DPI) or about 21 µm per dot resolution. The background was generated using a simple python code (see Appendix A). The random pattern of the background was 400 by 400 elements in a 2 by 2 cm area (50 µm per element). The background was illuminated by a 5 W flashlight powered with a DC power supply. The background was imaged by a camera that is identical to the interferometer CCD camera using a 50 mm objective lens with a 10 mm extension tube for close-up imaging. In additional experiments, a commercial telecentric lens was used for imaging (see below). Magnification of the BOS diagnostic varied between 0.4 and 0.5.

The distance between the gas nozzle and the background was varied between 10 and 20 mm. The nozzle was connected to a compressed air reservoir at 55 PSI. The reservoir was large enough to provide constant air flow for 5 s before the pressure dropped down, and the reservoir had to be refilled. In the middle of this 5 s window, both CCD cameras acquired synchronously triggered images with the same exposure time of 100 µs. Immediately after that, the gas was shut off, and both cameras acquired another image. This second set of images was used as the references for the corresponding diagnostics.

Fig. 1. Experimental setup. Blue components are BOS, gray, Mach–Zehnder interferometer with green laser beam. Pink is the gas jet nozzle.

Download Full Size | PDF

3. OPTICAL CALIBRATION AND SENSITIVITY OF BOS

The optical setup for the BOS diagnostic is depicted in Fig. 2. The most important parameter of the diagnostic is the distance ${Z_3}$ between the background and the nozzle. This distance defines the sensitivity—a magnitude of the shift of the background pixels caused by the angular deflection of the beam in the gas flow. In the paraxial approximation, the sensitivity is linearly proportional to ${Z_3}$, as shown below.

Fig. 2. Ray tracing diagram for the background pixels shift calculation. Cross section of small gas flow element can be approximated as a lens with very long focal length F. Gas flows in the X direction. Black rays represent undisturbed propagation (no gas flow), blue rays travel through and are deflected by the flowing gas.

Download Full Size | PDF

In general, any medium with a refractive index gradient can be represented as a collection of lens elements (spherical or cylindrical). Spatial variation of the refractive index leads to deflection of the rays in proportion to gradients of the refractive index [7]. The deflection angle $\theta = {h_1}/F$, where $F$ is the focal length of this element, and ${h_1}$ is point of deflection (somewhere in the gas jet). The deflection angle (see Fig. 2) is

(1)$$\theta = \frac{{{h_1}}}{{F(x,y)}} = \frac{{\Delta {h_1}(x,y)}}{{{Z_3}}} = \frac{{\Delta {h_2}(x,y)}}{{M{Z_3}}},$$

where $(x,y)$ are coordinates at the gas jet plane, and $M = \frac{{{Z_2}}}{{{Z_1}}}$ is the optical magnification created by the imaging lens $f$. The focal length $F$ is assumed to be much longer, while the diameter of the gas flow is much smaller than any other longitudinal (${{Z}}$ direction) dimensions in Fig. 2.

The BOS CCD camera must be focused on the background, leaving the nozzle out of focus. For regular (entocentric) imaging optics, objects closer to the lens appear bigger on the camera, and the corresponding optical calibration must be done in the nozzle plane. In general, the magnification of the out-of-focus objects can be calculated, but such calculations require detailed knowledge of the imaging optics. This task might be hard in the case of a compound lens with unknown locations of the principal points. The calibration of the nozzle plane magnification can be done experimentally, for example, by counting the number of pixels in the nozzle diameter. Under optimal conditions of operation (resolution versus sensitivity trade-off [8]), the nozzle image appears blurred on the CCD camera [Fig. 3(a)]. This problem can be resolved by placing a periodic pattern (like transparent ruler) in the nozzle plane. The spatial period of the pattern provides the necessary calibration.

A telecentric lens can provide a simpler, robust, but more expensive solution for the calibration problem. Telecentric lenses were used in BOS measurements [9] and do not require any additional calibration. All objects, in or out of focus, appear with the same magnification on the CCD camera. More expensive commercial telecentric lenses provide an adjustable magnification option. The major disadvantages of the telecentric lenses are a short working distance and limited field of view. The price and size of the lens increase drastically when these parameters are increased. In the present experimental setup, a lower price range lens (Edmund Optics 88385) was sufficient for imaging the gas flow.

4. DETERMINATION OF REFRACTIVE INDEX IN AXISYMMETRIC FLOW

In the geometric optics limit, we derive the relationship between the gas jet density and the pixel shift at the CCD plane. The position ${\bf{s}} = x{\hat{\bf x}} + y{\hat{\bf y}} + z{\hat{\bf z}}$ (coordinates shown in Fig. 2) and wavenumber ${\bf{k}}$ of a ray of light are given by [7,10]

(2)$$\frac{{d{\bf{s}}}}{{dt}} = {\nabla _k}\omega ({\bf{s}},{\bf{k}}) = \frac{{\bf{k}}}{k}\frac{c}{n},$$

(3)$$\frac{{d{\bf{k}}}}{{dt}} = - {\nabla _s}\omega ({\bf{s}},{\bf{k}}) = \frac{{kc}}{{{n^2}}}\nabla n,$$

where $\omega ({\bf{s}},{\bf{k}}) = | {\bf{k}} |c/n({\bf{s}})$ is the local optical frequency (a conserved quantity), $n = n({\bf{s}})$ is the spatially dependent refractive index, $k = | {\bf{k}} |$, and $c$ is the speed of light. We assume there is a uniform background refractive index ${n_0}$ and that the spatially dependent part of the refractive index is a cylindrically symmetric perturbation ${n_1}(x,r)$, where $r = \sqrt {{y^2} + {z^2}}$. The radial extent of ${n_1}$ is ${r_0}$, specifically, ${n_1}(x,r \gt {r_0}) = 0$. The ray trajectories are solved perturbatively with ${n_1}$ being the small parameter: ${\bf{s}}(t) = {{\bf{s}}_0}(t) + {{\bf{s}}_1}(t) + \cdots$, ${\bf{k}}(t) = {{\bf{k}}_0}(t) + {{\bf{k}}_1}(t) + \cdots$, and $n({\bf{s}}) = {n_0} + {n_1}({\bf{s}}) + \cdots$. In the absence of the refractive index inhomogeneity, the ray’s dynamics satisfies

(4)$$\frac{{d{{\bf{s}}_0}}}{{dt}} = \frac{{{{\bf{k}}_0}}}{{{k_0}}}\frac{c}{{{n_0}}},$$

(5)$$\frac{{d{{\bf{k}}_0}}}{{dt}} = 0,$$

and the lowest order correction is

(6)$$\frac{{d{{\bf{s}}_1}}}{{dt}} = \frac{{{{\bf{k}}_1}}}{{{k_0}}}\frac{c}{{{n_0}}} - \frac{{{{\bf{k}}_0}}}{{{k_0}}}\frac{c}{{{n_0}}}\frac{{{n_1}}}{{{n_0}}} - \frac{{{{\bf{k}}_0}}}{{{k_0}}}\frac{c}{{{n_0}}}\frac{{{{\bf{k}}_0} \cdot {{\bf{k}}_1}}}{{k_0^2}},$$

(7)$$\frac{{d{{\bf{k}}_1}}}{{dt}} = - \frac{{{k_0}c}}{{{n_0}}}{\left. {\nabla {n_1}} \right|_{{\bf{s}} = {{\bf{s}}_0}}}.$$

Fig. 3. (a) Raw image of the background with the nozzle visible on the left side, (b) quiver plot of the shift vector field, and (c) vertical component of the shift. The dashed rectangle in (a) represents the analyzed area.

Download Full Size | PDF

Consider the horizontal principle ray (shown as black in Fig. 2) coming from a point, denoted as $0{\hat{\bf x}} + {h_1}{\hat{\bf y}} - {Z_3}{\hat{\bf z}}$, on the background screen. The unperturbed ray has the wavenumber and trajectory given by ${{\bf{k}}_0} = {k_0}{\hat{\bf x}}$ and ${{\bf{s}}_0}(z) = 0{\hat{\bf x}} + y{\hat{\bf y}} + z{\hat{\bf z}}$, where $y = {h_1}$, and we have changed the independent variable from time $t$ to the unperturbed ray’s $z$ position $z = ct/{n_0}$. The transverse deflection relative to ${{\bf{k}}_0}$ is governed by

(8)$$\frac{{d{{\bf{s}}_{1 \bot}}}}{{dz}} = \frac{{{{\bf{k}}_{1 \bot}}}}{{{k_0}}},$$

(9)$$\frac{{d{{\bf{k}}_{1 \bot}}}}{{dz}} = - {k_0}{\left({{\hat{\bf x}}\frac{{\partial {n_1}}}{{\partial x}} + {\hat{\bf y}}\frac{{\partial {n_1}}}{{\partial y}}} \right)_{{\bf{s}} = {{\bf{s}}_0}}},$$

where ${{\bf{s}}_{1 \bot}} = {\hat{\bf r}} \cdot {{\bf{s}}_1}$, ${{\bf{k}}_{1 \bot}} = {\hat{\bf r}} \cdot {\bf{k}}$, and ${\hat{\bf r}} = {\hat{\bf x}} + {\hat{\bf y}}$. Prior to the gas jet, ${{\bf{s}}_{1 \bot}}(z \lt - {r_0}) = 0$ and ${{\bf{k}}_{1 \bot}}(z \lt - {r_0}) = 0$. After the ray has passed through the gas jet, the transverse wavenumber will reach its final value ${{\bf{k}}_{1 \bot}}(z \gt {r_0})$. This defines the ray’s angular deflection $\boldsymbol{\theta} = {{\bf{k}}_{1 \bot}}(z \gt {r_0})/{k_0}$, where $\boldsymbol{\theta} = {\theta _x}{\hat{\bf x}} + {\theta _y}{\hat{\bf y}}$:

(10)$$\begin{split}\boldsymbol{\theta} (x,y) =\! \int_{- \infty}^\infty\! {\text{d}}z^\prime \left({\frac{{\partial {n_1}(x,\sqrt {{{z^\prime}^2} + {y^2}})}}{{\partial x}}{\hat{\bf x}} + \frac{{\partial {n_1}(x,\sqrt {{{z^\prime}^2} + {y^2}})}}{{\partial y}}{\hat{\bf y}}} \right)\!.\end{split}$$

The integral is along a cord intersecting the gas jet. With a change of variables in integration from $z$ to $r$, then Eq. (10) becomes the Abel transform for ${n_1}$:

(11)$$\begin{split}\boldsymbol{\theta} (x,y) = \int_{\left| y \right|}^\infty \frac{{r{\text{d}}r}}{{\sqrt {{r^2} - {y^2}}}}\left({\frac{{\partial {n_1}(x,r)}}{{\partial x}}{\hat{\bf x}} + \frac{y}{r}\frac{{\partial {n_1}(x,r)}}{{\partial r}}{\hat{\bf y}}} \right).\end{split}$$

The inverse Abel transform can be carried out for both directions, but the $y$ component of the deflection is of particular interest because it is the result of radial gradients in the density. The radial gradient is given by

(12)$$\frac{1}{r}\frac{{\partial n}}{{\partial r}} = - \frac{1}{\pi}\int_r^\infty \frac{{{\text{d}}y}}{{\sqrt {{y^2} - {r^2}}}}\frac{\partial}{{\partial y}}\left({\frac{{{\theta _y}(x,y)}}{y}} \right).$$

The angular deflection in the $y$ direction ${\theta _y}$ is not directly measured, but Eq. (1) shows the relationship between an angular deflection $\theta$ and the transverse displacement at the camera $\Delta {h_2}$. Additionally, the coordinates $x$, $y$, and $r$ are in the object plane of the gas density perturbation, not the plane of the camera where the displacement is measured. The measured displacement $\Delta {h_2}$ is a function of camera plane coordinate ${h_2}(x,y) = My$. Therefore, the Abel transform in the plane of the camera is

(13)$$\begin{split}\frac{1}{r}\frac{{\partial n}}{{\partial r}} = - \frac{M}{{\pi\! {Z_3}}}\int_{\textit{Mr}}^\infty \frac{{{\text{d}}{h_2}}}{{\sqrt {h_2^2 - {M^2}{r^2}}}}\frac{\partial}{{\partial {h_2}}}\left({\frac{{\Delta {h_2}(x,{h_2})}}{{{h_2}}}} \right).\end{split}$$

The gas density approaches the ambient density ${n_0}$ as $r \to \infty$ and allows for integration of the above expression:

(14)$$\begin{split}n(r,x)& = {n_0} + \frac{M}{{\pi\! {Z_3}}}\int_r^\infty {\text{d}}r^\prime r^\prime \int_{{Mr^\prime}}^\infty \frac{{{\text{d}}{h_2}}}{{\sqrt {h_2^2 - {M^2}{{r^\prime}^2}}}}\frac{\partial}{{\partial {h_2}}}\\&\quad\times\left({\frac{{\Delta {h_2}(x,{h_2})}}{{{h_2}}}} \right).\\[-1.5pc]\end{split}$$

The first integration (over ${h_2}$) in Eq. (14) is the inverse Abel transformation from projection coordinate $y$ to radius $r$. The second integration calculates the refractive index from its gradient. As explained in Section 3, pixels of shift in the BOS camera image must be translated into real length units of the gas jet plane. In this case, the magnification $M$ is already accounted for in calibration of $h$ and $r$ values and must be omitted from the numerical calculation of the Eq. (14) integrals presented below.

5. NUMERICAL ANALYSIS OF BOS DATA

A vector field representing the pixels’ shift was constructed using the Python implementation of OpenPIV [11] from the image pair acquired in the experiment. The images were saved in 16 bits tagged image file format (TIFF) files and imported into OpenPIV using the openpiv.tools.imread function. The full frame image was cropped to the region of interest (ROI) shown approximately by the dashed rectangle in Fig. 3(a).

After selecting the ROI, the initial cross correlation was done using the openpiv.process.extended_search_area_piv function. The size of the interrogation window was between 8 and 12 pixels, and the overlap was 2 pixels smaller. The quiver plot coordinates were extracted using the openpiv.process.get_coordinates function, the outliers were removed with the openpiv.validation.sig2noise_val function, and the missing vectors were substituted by the openpiv.filters.replace_outliers function. The final vector shift map [Fig. 3(b)] was scaled from pixels of the camera to millimeters of the gas jet plane using the openpiv.scaling.uniform function. This procedure is similar to an example in the OpenPIV tutorial [12]. This part of the analysis is the most time consuming and might take a few minutes on a desktop computer.

Fig. 4. Cross section of the refractive index change. Coordinate system follows notations of Fig. 2.

Download Full Size | PDF

Fig. 5. (a) Cross sections of the gas density measured from BOS and (b) interferometry in the same shot of the gas jet. (c) Lineouts of the data in (a) and (b) along white lines. (d) Similar lineouts for a different shot with telecentric lens used in BOS. Green dotted line in (c) and (d) shows the ambient gas density.

Download Full Size | PDF

The vertical ($y$ direction in notations of Fig. 2) component of the shift vectors created a schlieren image of the gas flow [Fig. 3(c)]. The gas flow from under-expanded nozzle produced a periodic structure that is characteristic of shock diamonds, and the gradient of the refractive index was encoded in this image. Recall that the gas flow was along the positive $x$ -axis. The Fig. 3(c) image is anti-symmetric with respect to the horizontal axis that passes through the center of the gas flow. The inverse Abel transformation in Eq. (14) represents the vertical shift component divided by the corresponding vertical projection coordinate $\frac{{\Delta {h_2}(x,{h_2})}}{{{h_2}}}$. This ratio is expected to be a symmetric function of ${h_2}$. To ensure this, it is possible to introduce a new anti-symmetric vertical coordinate $y$ in the range $[- {y_{{\max}}}/2,{y_{{\max}}}/2]$ with the line $y = 0$ at the symmetry axis (in $[x,y]$ plane). The final symmetrical function that is suitable for Abel inversion has the appearance of $\frac{{\Delta {h_2}(y) - \Delta {h_2}(- y)}}{{2y}}$. In case the new $y$ coordinate is exactly zero at the axis (odd number of pixels in the vertical dimension of the ROI), the neighbor value is employed.

We used the Python PyAbel library [13] to transform the projection coordinates of the refractive index gradient into three-dimensional (3D) cylindrically symmetric coordinates. The symmetry axis in PyAbel must be vertical, and the horizontally symmetrical image must be transposed prior to inversion. We used the basex transformation method [14] and, because the integration in PyAbel uses pixels of the image, applied the experimental scaling factor (corrected to the overlap parameter used in OpenPIV), to keep the numerical values in real physical units. Finally, the second integral of Eq. (14) was calculated using the standard Python NumPy cumtrapz function applied to the result of the Abel transform multiplied by the coordinate $y$ introduced above. The result was divided by factor from Eq. (14) to get values of the refractive index change. The noise was reduced by a Gaussian filter, and the final result is shown in Fig. 4.

6. RESULTS AND DISCUSSION

The radial dependence of the gas density in the flow can be calculated from the refractive index value using the Gladstone–Dale relation [15]:

n(r) - 1 = ({n_0} - 1)\frac{{\eta (r)}}{{{\eta _0}}},

where $\eta (r)$ is the gas number density to be determined, ${\eta _0} = 2.5 \times {10^{19}}\;{\text{c}}{{\text{m}}^{- 3}}$ is the number density of unperturbed gas, and ${n_0} = 1 + 2.3 \times {10^{- 4}}$ is the refractive index of unperturbed air. The value of the refractive index change $n(r) - 1$ calculated from Eq. (14) is shown in Fig. 4. The number gas density in the jet is shown in Fig. 5(a).

Gas density measurements using interferometry were extracted from the raw interferograms using standard techniques described elsewhere [16]. Inverse Abel transformation of the interferometry data employed the same approach as we used in the BOS calculations. The cross section of the gas density extracted from the interferometry is presented in Fig. 5(b). Lineouts of BOS and interferometry density profiles 0.3 mm from the gas flow axis are plotted in Fig. 5(c). Since the CCD cameras (Fig. 1) were not lined up, we shifted the horizontal position of the curves to match the phase of the density modulations in Figs. 5(c) and 5(d). BOS measurement lineout in Fig. 5(d) is an example of using a telecentric lens in the BOS diagnostic. In both cases, agreement between the diagnostics is remarkable.

To extract BOS data, we used only vertical components of the shift vector field. The horizontal component can also be used for the density calculation, but the simple recipe and analysis provided here requires a different approach and will be discussed in a forthcoming paper. In fact, the horizontal component of the shift is analogous to an interference pattern with fringes oriented along the flow. Since the main interest is to obtain the radial dependents of the gas density, such orientation is not optimal.

APPENDIX A: BACKGROUND GENERATION CODE

The Python code below generates a ${{50}} \times {{50}}$ pixels random background pattern and saves it to a background.png file. The final image file must be scaled appropriately prior to printing.

#!/usr/bin/env python

# coding: utf-8

import numpy as np

import matplotlib.pyplot as plt

# number of pixels

row = 50

col = 50

# figure size on the screen

plt.rcParams[“figure.figsize”] = [5,5]

# random pattern generation

image = np.random.randint(2, size = (row, col))

# show result and save image to file

plt.imshow(image, ’gray’)

plt.axis(’off’)

plt.savefig(’background.png’, bbox_inches = ’tight’, pad_inches = 0)

Funding

U.S. Naval Research Laboratory.

Disclosures

The authors declare no conflicts of interest.

REFERENCES

1. S. Bulanov, N. Naumova, F. Pegoraro, and J. Sakai, “Particle injection into the wave acceleration phase due to nonlinear wake wave breaking,” Phys. Rev. E 58, R5257–R5260 (1998). [CrossRef]

2. E. P. Goodwin and J. C. Wyant, Field Guide to Interferometric Optical Testing (SPIE, 2006), Vol. FG10.

3. M. L. Norman, K. H. A. Winkler, L. Smarr, and M. D. Smith, “Structure and dynamics of supersonic jets,” Astron. Astrophys. 113, 285–302 (1982).

4. N. Van Hinsberg and T. Rösgen, “Density measurements using near-field background-oriented Schlieren,” Exp. Fluids 55, 1720 (2014). [CrossRef]

5. F. Nicolas, D. Donjat, O. Léon, G. Le Besnerais, F. Champagnat, and F. Micheli, “3D reconstruction of a compressible flow by synchronized multi-camera BOS,” Exp. Fluids 58, 46 (2017). [CrossRef]

6. M. Raffel, “Background-oriented Schlieren (BOS) techniques,” Exp. Fluids 56, 60 (2015). [CrossRef]

7. L. D. Landau and E. M. Lifshitz, The Classical Theory of Fields, Vol. 2 of Course of Theoretical Physics Series (Butterworth-Heinemann, 1980).

8. A. B. Gojani, B. Kamishi, and S. Obayashi, “Measurement sensitivity and resolution for background oriented Schlieren during image recording,” J. Visualization 16, 201–207 (2013). [CrossRef]

9. M. Ota, F. Leopold, R. Noda, and K. Maeno, “Improvement in spatial resolution of background-oriented Schlieren technique by introducing a telecentric optical system and its application to supersonic flow,” Exp. Fluids 56, 48 (2015). [CrossRef]

10. T. Khan and A. Thomas, “On derivation of the radiative transfer equation and its spherical harmonics approximation for scattering media with spatially varying refractive indices,” Technical Report (Clemson University Mathematical Sciences, 2003), pp. 1–49.

11. https://github.com/OpenPIV/openpiv-python.

12. https://openpiv.readthedocs.io/.

13. https://zenodo.org/record/3370021.

14. D. D. Hickstein, S. T. Gibson, R. Yurchak, D. D. Das, and M. Ryazanov, “A direct comparison of high-speed methods for the numerical Abel transform,” Rev. Sci. Instrum. 90, 065115 (2019). [CrossRef]

15. J. H. Gladstone and T. P. Dale, “Researches on the refraction, dispersion, and sensitiveness of liquids,” Philos. Trans. R. Soc. London 153, 317–343 (1863). [CrossRef]

16. S. Sugawara, S. Nakao, Y. Miyazato, Y. Ishino, and K. Miki, “Three-dimensional reconstruction of a microjet with a mach disk by Mach–Zehnder interferometers,” J. Fluid Mech. 893, A25 (2020). [CrossRef]

Benchmarking background oriented schlieren against interferometric measurement using open source tools

Abstract

1. INTRODUCTION

2. EXPERIMENT

3. OPTICAL CALIBRATION AND SENSITIVITY OF BOS

4. DETERMINATION OF REFRACTIVE INDEX IN AXISYMMETRIC FLOW

5. NUMERICAL ANALYSIS OF BOS DATA

6. RESULTS AND DISCUSSION

APPENDIX A: BACKGROUND GENERATION CODE

Funding

Disclosures

REFERENCES

Cited By

Figures (5)

Equations (15)

Applied Optics