Towards robust BOS measurements for axisymmetric flows

Abstract

Accurately reconstructing a radial refractive index (n) field is challenging in axisymmetric background oriented Schlieren (BOS) measurement. In this study, we systematically investigated several widely adopted inversion algorithms in BOS applications. To quantitatively assess the performance of each algorithm, a synthetic experiment mimicking a helium jet discharged into ambient air was established to provide the reference. Relying on the necessity to solve a Poisson equation for a line-of-sight projected variable, tested algorithms were categorized into two groups: direct and indirect. In the direct approach, the algorithm is applied directly to the light deflection angle ($\varepsilon$) to reconstruct the radial $\delta$ field, defined as $\delta =(n-n_{0})/n_{0}$ where $n_{0}$ is the reference refractive index. In the indirect group, the Poisson equation is solved first. Then, an inversion algorithm is subsequently applied to the projected ${{\overline{\delta }}}$ to obtain $\delta$ in the radial plane. The two approaches were compared with the synthetic experiment both using the adaptive Fourier–Hankel methods (AFH). The comparison showed that at the cost of introducing the additional step of solving the Poisson equation, the indirect approach performed more accurately when noises were present in the $\varepsilon$ measurements. To identify the proper inversion algorithm suitable for the indirect approach, we further compared four types of algorithms in the synthetic experiments including AFH, onion peeling (OP), three-point Abel (TPA), and filtered back projection tomography (FBPT). The results showed that TPA had the best performance in terms of the reconstruction accuracy with noisy $\varepsilon$ data. Finally, experiments on axisymmetric helium jets were conducted to confirm the effectiveness of the proposed TPA algorithm in the indirect approach.

Graphic Abstract

1 Introduction

Intensive efforts have been made recently to apply BOS method to measure three-dimensional (3D) refractive index n, or density $\rho$ fields, of non-reactive and reactive flows (Venkatakrishnan and Meier 2004; Goldhahn and Seume 2007; Atcheson et al. 2008; Nicolas et al. 2016; Hayasaka et al. 2016; Grauer et al. 2018). The tomographic reconstruction of the 3D n field requires multiple 2D projected views from different viewing angles. These projected views can be obtained by triggering different cameras simultaneously (Nicolas et al. 2016) or by one camera with a fiber bundle enabling multiple views at the expense of reduced resolution (Liu et al. 2019), or by rotating one camera around the probed volume assuming steady or periodic flows (Ota et al. 2011; Lang et al. 2017). Achieving sufficient spatial resolution and accuracy of the 3D n field is time-consuming and computationally expensive.

Axisymmetric flows consist of a particular type of 3D flows widely encountered in laboratories, especially for laminar jet experiments. The tomographic $\rho$ field of a jet can be reconstructed using a single view angle owing to the symmetry of the flow. Such investigations on axisymmetric flows were conducted shortly after the introduction of the BOS technique (Dalziel et al. 2000; Richard and Raffel 2001; Meier 2002). Abel inversion was first applied directly to the deviation angle of the BOS measurements (Wong 2001), and the accuracy was also explored (Kirmse 2003), although the two references are not easily accessible. A special recursive ring method was designed by discretizing the n field into a set of rings, assuming a constant n within a ring (Kirmse 2003; Klinge et al. 2003). A sufficient number of rings is required to allow accurate reconstruction.

Instead of processing the deflection angle data directly to obtain an axisymmetric density field, a Poisson equation of the line-of-sight integrated ${{\overline{\delta }}}$ was solved first, and the solution field was fed into FBPT algorithm to reconstruct the tomographic field (Venkatakrishnan and Meier 2004; Venkatakrishnan 2005). The preference to FBPT over Abel inversion was attributed to the noise susceptibility and singularity issues of the Abel inversion. However, the introduction of the Poisson equation was argued to be unnecessary and redundant since displacement information can be fed into the FBPT algorithm directly (Goldhahn and Seume 2007). To avoid the computationally expensive FBPT methods, a Fourier-based method, such as AFH method, was proposed (Ma et al. 2008) as a way of accounting for the noise susceptibility and minimizing truncation errors without an additional filter function. Recently, the AFH method was introduced to process the axisymmetric BOS data (Tan et al. 2015), which removes the requirement to solve the Poisson equation. Albeit all these studies dedicated to the axisymmetric BOS technique, a systematic study on the necessity to solve the Poisson equation is much needed. Moreover, detailed comparisons on the performance of different inversion algorithms, especially with the presence of the noise, are still missing in BOS applications.

In this paper, a synthetic experiment mimicking the helium jet discharged into the ambient air was designed to evaluate the necessity to solve the Poisson equation and the performance of various inversion algorithms for the axisymmetric BOS measurement. Laboratory measurements were also conducted over a steady helium jet to verify the selected algorithm.

2 Methodology

2.1 Principle of BOS

A typical BOS setup consists of a camera with a lens with a long focal length (Goldhahn and Seume 2007; Gojani et al. 2013), a background with a suitable pattern compatible with the correlation algorithm (Atcheson et al. 2009; Ota et al. 2011), and a light source for the background illumination.

Fig. 1

A typical BOS setup with an axisymmetric helium jet as the Schlieren object

When light rays pass through transparent inhomogeneous media, the ray deflects. Considering the schematic of a typical BOS setup shown in Fig. 1, the deflection angle $\varepsilon$, e.g., in y-direction, can be calculated via a simple geometric relation, assuming small deflection angles as:

$$\begin{aligned} \varepsilon _y \approx \frac{{{\varDelta }y}}{Z_d}, \end{aligned}$$

(1)

where $Z_d$ is the distance from the center of the Schlieren object to the background and ${{\varDelta }y}$ is the y component of the displacement vector ${{\varvec{{\varDelta }}}}=({{\varDelta }x},{{\varDelta }y})$ in the image plane. The relation between $\varepsilon _y$ and gradient of the refractive index n is given by the following formula:

$$\begin{aligned} \varepsilon _y = \frac{1}{n_0}\int {\frac{\partial {n}}{\partial {y}}}dz, \end{aligned}$$

(2)

where $n_0$ is the refractive index of the ambient gas, z is the line-of-sight direction of the ray, and y is the direction perpendicular to the line-of-sight direction. For two-dimensional flows, Eq. (2) can be simplified by assuming a constant ${\partial {n}}/{\partial {y}}$ over the thickness W of the Schlieren object thus,

$$\begin{aligned} \varepsilon _y = \frac{W}{n_0}\frac{\partial {n}}{\partial {y}}. \end{aligned}$$

(3)

Combining Eqs. (1) and (3) and considering the equivalent equations in x-direction, an equation set can be obtained, which is usually solved in the form of a Poisson equation with proper boundary conditions (Xiong et al. 2020). For axisymmetric flows, ${\partial {n}}/{\partial {y}}$ is not constant along z in Eq. (2) and the thickness of the Schlieren object becomes a function of x. A solution strategy utilizing the symmetry of the flow is required.

2.2 Axisymmetric BOS

Several solutions strategies exist to process axisymmetric BOS data. Depending on the necessity to solve the Poisson equation, these methods can be categorized into direct and indirect groups. In the former group, $\varepsilon$ is directly linked to $\delta$, or n, in the radial plane. In the later one, the line-of-sight integrated $\delta$ field (${{\overline{\delta }}}$) needs to be solved first via the Poisson equation, and then an inversion algorithm is further applied to the field to obtain the $\delta$ field in the radial plane. For convenience, the difference between the direct and indirect approach is illustrated using the Abel transform.

Abel transform describes the projection of a spherically or axially symmetric function f onto a plane. To apply forward and inverse Abel transforms, the field function f(r) has to drop to zero faster than 1/r. In this study, instead of projecting n(r), a modified refractive index field $\delta (r)$ that approaches zero in the far-field is defined in Eq. (4) (shown schematically in Fig. 1):

$$\begin{aligned} \delta =\frac{n}{n_0}-1. \end{aligned}$$

(4)

The relation between $\varepsilon _y$ and $\delta$ can be derived based on Eq. (2) as:

$$\begin{aligned} \varepsilon _y= & {} \frac{1}{n_0}\int {\frac{\partial {n}}{\partial {y}}}dz= \int {\frac{\partial \delta }{\partial {y}}}dz=\frac{\partial }{\partial {y}} \left( \int {\delta }dz\right) \nonumber \\= & {} \frac{\partial {{\bar{\delta }}}}{\partial {y}}, \end{aligned}$$

(5)

where ${{\bar{\delta }}}$ is the line integrated $\delta$ in the z-direction. The forward Abel transform, connecting ${{\bar{\delta }}}$ and $\delta$ together, takes the following form:

$$\begin{aligned} {\bar{\delta }}(y)=2\int _y^\infty \frac{\delta (r)r}{\sqrt{r^2-y^2}}dr, \end{aligned}$$

(6)

and accordingly the inverse Abel transform can be given by:

$$\begin{aligned} \delta (r)=-\frac{1}{\pi }\int _r^\infty \frac{d{\bar{\delta }}(y)}{dy} \frac{dy}{\sqrt{y^2-r^2}}. \end{aligned}$$

(7)

Considering the Eq. (5), it gives:

$$\begin{aligned} \delta (r)= -\frac{1}{\pi }\int _r^\infty {\varepsilon _y} \frac{dy}{\sqrt{y^2-r^2}}. \end{aligned}$$

(8)

Thus $\delta (r)$ correlates directly with the deflection angle $\varepsilon _y$, which could be calculated based on Eq. (1).

The direct approach utilizes Eq. (8) in which $\varepsilon$ and $\delta$ are connected directly without any intermediate step. As a contrast, indirect methods compute $\delta$ via two steps. Firstly, the line-of-sight integrated value ${{\bar{\delta }}}$ is calculated via solving a Poisson equation, which can be derived from Eq. (5) by taking the partial derivative in y-direction and consider the same equation in the x-direction by summing up the two as:

$$\begin{aligned} \frac{\partial ^2 {{\bar{\delta }}}}{\partial x^2}+\frac{\partial ^2 {{\bar{\delta }}}}{\partial y^2} = \left( \frac{\partial \varepsilon _x}{\partial x}+\frac{\partial \varepsilon _y}{\partial y}\right) . \end{aligned}$$

(9)

By applying proper boundary conditions, ${{\bar{\delta }}}$ can be calculated. More details to solve this Poisson equation in BOS measurements can be found in the literature, e.g., (Xiong et al. 2020). The inversion of ${{\bar{\delta }}}$ can be achieved, for example, by inverse Abel transform as shown in Eq. (7). A schematic summarizing the roadmap of the solution methodology within the BOS measurements is shown in Fig. 2, where $\otimes$ denotes the multiple-pass cross-correlation operation, ${{\varvec{\nabla }}} = (\partial /\partial x, \partial /\partial y)$ is the gradient operator and ${{\varvec{\nabla }}}\cdot {{\varvec{\nabla }}} = (\partial ^2 /\partial x^2+\partial ^2 /\partial y^2)$ is the laplacian operator, ${{\varvec{\varepsilon }}}=(\varepsilon _x,\varepsilon _y)$, K is the Gladstone-Dale constant, ${\mathcal {D}}$ is the inversion algorithm in the direct approach and ${\mathcal {I}}$ is the operator for the indirect approach.

Fig. 2

The roadmap for the solution methodology in axisymmetric BOS measurements

3 Experimental setup

3.1 Synthetic experiments

Synthetic experiments are widely used to test the accuracy of the inversion algorithm for the tomography but typically limited to 1D conditions (Dasch 1992; Zhang et al. 2018). In the jet experiments, density gradients occur not only in the radial direction but also in the stream-wise direction (x-direction as shown in Fig. 1). In this regard, a synthetic 2D displacement field, mimicking a helium jet discharged into air ambient, was designed with the non-dimensional displacement field as:

$$\begin{aligned} {\widehat{{\varDelta }x}_a}= & {} -\frac{8}{15}{\widehat{Z_d}}\delta _c\left( 1-{\hat{y}}^2\right) ^{5/2}, \end{aligned}$$

(10)

$$\begin{aligned} {\widehat{{\varDelta }y}_a}= & {} \frac{8}{3}{\widehat{Z_d}}\delta _c({\hat{x}}-2){\hat{y}}(1-{\hat{y}}^2)^{3/2}, \end{aligned}$$

(11)

where the subscript a denotes analytical expression, the normalized variables are defined as ${\hat{x}} = x/R$, ${\hat{y}} = y/R$, ${\hat{r}} = r/R$, ${\widehat{Z_d}} = Z_d/R$ and $\delta _{c} = \left( n_{{j}} - n_{0} \right) /n_{0}$ with the jet nozzle radius R, jet refractive index $n_{{j}}$ and ambient air refractive index $n_{0}$. The selected expressions for ${\widehat{{\varDelta }x}_a}$ can allow analytical solutions for ${\hat{\delta }}_a$ and $\hat{{{\overline{\delta }}}}_a$ as:

$$\begin{aligned} {\hat{\delta }}_ {{a}}= & {} \delta _{c} \left( 1 - {\hat{r}}^{2} \right) ^{2} \left( 1 - \frac{{\hat{x}}}{2} \right) , \end{aligned}$$

(12)

$$\begin{aligned} \hat{{\overline{\delta }}}_ {{a}}= & {} - \frac{8}{15} \delta _{c} ( {\hat{x}} - 2 ) \left( 1 - {\hat{y}}^{2} \right) ^{5 / 2} . \end{aligned}$$

(13)

Note that ${\widehat{Z_d}}$ was set to 2 (other values can be selected also) to obtain Eq. (13) by solving the Poisson equation as:

$$\begin{aligned} \frac{\partial ^{2} \hat{{\overline{\delta }}}}{\partial {\hat{x}}^{2}} + \frac{\partial ^{2} \hat{{\overline{\delta }}}}{\partial {\hat{y}}^{2}} = \frac{1}{{\hat{Z}}_{d}} \left( \frac{\partial {\widehat{{\varDelta }x}}}{\partial {\hat{x}}} + \frac{\partial {\widehat{{\varDelta }y}}}{\partial {\hat{y}}} \right) . \end{aligned}$$

(14)

The ${\hat{\cdot }}$ sign and the subscript a are omitted in the following derivations to alleviate the notation. To exhibit the synthetic experiments intuitively, 2D ${{\varvec{{\varDelta }}}}$ field, ${{\overline{\delta }}}$ field, together with the associated 3D $\delta$ field are shown in Fig. 3. Inside ten iso-surfaces of $\delta$ were presented in 3D space, projected ${\varDelta }$ and ${{\overline{\delta }}}$ fields were positioned at $Z=2$, the center of the jet was located at $y=0, z=0$. $n_{0}=1.000258$ and $n_{{j}}=1.000030809$ (refractive index of helium) were used. In practical BOS measurements, noises inevitably appear in the ${{\varvec{{\varDelta }}}}$ data from various sources. If a high spatial resolution is desired for the BOS measurements, the value of $Z_d$ determining the out-of-focus effect should be minimized together with the measurement sensitivity (Goldhahn and Seume 2007). Thus the resulting ${{\varvec{{\varDelta }}}}$ field can be only with a small signal-to-noise ratio. To take the noise effect into account, white noise term was added to the R.H.S. of Eqs. (10) and (11) as:

$$\begin{aligned} {{\textit{\textbf{w}}}} = \mathrm{SNR}\times {{\varvec{{\mathcal {R}}}}}\times C. \end{aligned}$$

(15)

Inside $\mathrm{SNR}$ is the signal-to-noise ratio varying from 0 to 20%, ${{\varvec{{\mathcal {R}}}}}$ is a randomly generated matrix vector of the same size as ${{\varvec{{\varDelta }}}}$ with values between –1 and 1, and C is chosen to be the peak value of $|{{\varvec{{\varDelta }}}}|$. Examples of $({{\varDelta }x},{{\varDelta }y})$ fields with the noise addition are shown in Fig. 4. Since the noise matrix is randomly generated in each run, one hundred of runs were computed for each test condition to obtain the statistics such as the variances and standard deviations.

Fig. 3

Direct visualization of the fields of $\delta$, $|{{\varvec{{\varDelta }}}}|$ and ${{\overline{\delta }}}$ in the synthetic experiment

Fig. 4

The fields of ${{\varDelta }x}$ and ${{\varDelta }y}$ in the synthetic experiment with a 20% noise addition

3.2 Helium jet experiment

BOS measurements over an axisymmetric helium jet were conducted to obtain the n field. The helium of 99.9% purity was discharged from the gas cylinder, passed through a mass flow controlling valve, and entered a cylindrical cavity. The flow was then straightened by a honeycomb to homogenize the velocity disturbances before the convergent section upstream of the jet pipe. The jet pipe outlet was round-shaped with $R = 7$ mm. To isolate the helium jet from surrounding flow disturbance, acrylic glass walls of a thickness $H_G = 4$ mm were built around the jet with a width of 480 mm and a height of 670 mm. A loudspeaker, driven by a voltage amplifier, was mounted at the bottom of the cavity to enable acoustic modulation of the jet. A function generator provided the sinusoidal signals to the amplifier at desired frequencies. $Z_d$ was set to 219 mm. A similar setup can be found in our previous study (Xiong et al. 2020).

The BOS system consisted of a high-speed camera (Lavision, HighSpeed-Star), an over-driving LED light (HARDSoft, IL-106X), and a background printed on a paper with random dots and sandwiched between two thin acrylic glasses. The dot size in the background pattern is designed to be about 2–3 pixels to minimize the peak-locking effects when processed with the cross-correlation algorithm (Raffel et al. 2018). The camera was equipped with a Nikon micro-lens of $f = 200$ mm, and the camera aperture was set at $f_{\#}$ = 32 to minimize the defocusing effect. This lens of a large focal length is preferred in BOS measurements for the sensitivity consideration (Goldhahn and Seume 2007; Gojani et al. 2013). Continuous LED illumination was adopted here.

Recorded images in BOS were processed with the multiple-pass cross-correlation algorithm in DAVIS 8.4 (Lavision). The large correlation window size was 48 by 48 pixels with an overlap of 50%, and the small window size was 12 by 12 pixels with an overlap of 75%. The pixel dimension in the background plane was approximately 0.13 mm. Since the displacements were mainly below 1 pixel in the current configuration, the capability of the DAVIS software to evaluate the sub-pixel shift is essential. DAVIS 8.4 is documented to capture the sub-pixel displacement accurately down to 0.03 pixel in PIV applications. With a carefully designed background pattern in BOS measurements, peak-locking and out-of-plane errors can be minimized. According to the synthetic experiments, the DAVIS algorithm can calculate sub-pixel displacements with a high level of precision down to 0.01 pixel in this study.

4 Results

4.1 Direct vs Indirect approach

The necessity of solving the Poisson equation needs to be clarified before selecting the optimal inversion algorithm for axisymmetric BOS measurements. The same type of algorithms should be adopted to compare the direct and indirect approaches. In this study, AFH was adopted as the type of algorithm. The AFH algorithm utilizes the relation that the Fourier transform of ${{\overline{\delta }}}$ equals to the zero-order Hankel transform of $\delta$ as (Mook 1983):

$$\begin{aligned} {{\mathcal {H}}}{{\mathcal {T}}}[\delta (r)]={\mathcal {F}}{\mathcal {T}}[{{\overline{\delta }}}(y)], \end{aligned}$$

(16)

where ${\mathcal {H}}{\mathcal {T}}$ and ${\mathcal {F}}{\mathcal {T}}$ are the Hankel and Fourier transform respectively. Applying the inverse Hankel transform on both side gives (Ma et al. 2008):

$$\begin{aligned} \delta (r)= & {} {\mathcal {H}}{\mathcal {T}}^{-1}[{\mathcal {F}} {\mathcal {T}}[{{\bar{\delta }}}(y)]] \nonumber \\= & {} {\frac{1}{2 \pi } \int _{0}^{\infty } {\mathcal {F}}{\mathcal {T}}[{{\overline{\delta }}}(y)] \omega J_{0}(\omega r) \mathrm {d} \omega .} \end{aligned}$$

(17)

The numerical integration of Eq. (17) requires ${{\overline{\delta }}}$ field to be obtained first by solving the Poisson equation. To apply AFH directly, instead of ${{\bar{\delta }}}(r)$, the deflection angle $\varepsilon$ should be used. Following the derivations in Chehouani and Fagrich (2013), the following equation is found:

$$\begin{aligned} {\delta (r)=\frac{1}{ 2 \pi i} \int _{0}^{+\infty } {\mathcal {F}}{\mathcal {T}}[\varepsilon _y] J_{0}(\omega r) \mathrm {d} \omega ,} \end{aligned}$$

(18)

where $J_{0}(\omega r)=\frac{1}{2 \pi } \int _{0}^{2 \pi } \exp (-{\hat{i}} \omega r \mathrm {sin} \theta ) \mathrm {d} \theta$ is the zero order Bessel function. By introducing a factor $\alpha$ with $0 < \alpha \le 1$ for choosing the frequency spacing, following numerical integration of Eqs. (17) and (18) can be obtained as:

$$\begin{aligned} \delta (r_i) = \mathrm {D}_{ij}{{\overline{\delta }}}(y_j)={\mathbb {D}}_{ij}\varepsilon _y(y_j). \end{aligned}$$

(19)

Inside the discretization coefficient matrix is of the form as (Ma et al. 2008; Chehouani and Fagrich 2013; Tan et al. 2015):

$$\begin{aligned} \mathrm {D}_{ij}= & {} -\frac{\alpha }{N}\sum _{j=0}^N\sum _{k=0}^{{N_\alpha }} \mathrm {sin}\left( \frac{\alpha \pi jk}{N}\right) J_0\left( \frac{\alpha \pi ki}{N}\right) , \end{aligned}$$

(20)

$$\begin{aligned} {\mathbb {D}}_{ij}= & {} \frac{\alpha ^2\pi }{2NR}\sum _{j=0}^N \sum _{k=0}^{{N_\alpha }}k\mathrm {cos}\left( \frac{\alpha \pi jk}{N}\right) J_0\left( \frac{\alpha \pi ki}{N}\right) . \end{aligned}$$

(21)

Inside $N=R/{\varDelta }{r}$ is the number of intervals with distance ${\varDelta }r$, and $N_\alpha$ is the closest integer less than or equal to $N/\alpha$. $\alpha$ can take any value between 0 and 1. When choosing a smaller $\alpha$, more discrete frequencies are taken into consideration, and the truncation error of the method is reduced. At the same time, the computation time increases significantly as $\alpha$ is reduced. A value of 0.2 is recommended in the literature (Chehouani and Fagrich 2013). In this paper value of 0.1 was used to decrease the truncation error further.

Applying the AFH methods in both the direct and indirect approaches allows a direct evaluation of the necessity of solving the Poisson equation. Consider first the case with no noise appearing in the synthetic experiment. The 2D error distributions $\epsilon$ associated with the indirect and direct approaches are shown in Fig. 5 with ${\epsilon }_{\delta }$ defined as:

$$\begin{aligned} {\epsilon }_{\delta }=\frac{\delta -\delta _{a}}{\delta _0}. \end{aligned}$$

(22)

Similar definition of $\epsilon _{{{\overline{\delta }}}}$ is achieved by replacing the variable $\delta$ with ${{{\overline{\delta }}}}$ in the equation. In the indirect approach, the first step error $\epsilon _{{{\overline{\delta }}}}$ is from solving the Poisson equation for ${{{\overline{\delta }}}}$ as shown in Fig. 5a. Subsequent inversion results in $\epsilon _{\delta }$, as shown in Fig.  5b. Overall the $\epsilon _{{{\overline{\delta }}}}$ is one order of magnitude smaller compared to the $\epsilon _{\delta }$ indicating the major error comes from the AFH inversion step in the indirect approach. Figure  5c exhibits $\epsilon _{\delta }$ from the direct AFH inversion. The comparison between Fig. 5b and c show that the direct approach is more accurate. This could be explained by the additional error introduced by solving the Poisson equation as shown in Fig. 5a. Thus without the noise contamination, the direct approach should be preferred over the indirect approach owing to the smaller error and also to the lower computational cost without the need to solve the Poisson.

Fig. 5

Error distributions from solving the Poisson equation (a), from the AFH inversion in the indirect approach (b), and from the AFH inversion in the direct approach (c) without the noise addition

In practice, noises are always appearing in the BOS measurement from multiple sources (Xiong et al. 2020), thus it is critical to consider the impact from the noises when evaluating the indirect and direct approaches. In Fig.  6, plots similar to Fig. 5 are listed by assuming $\mathrm{SNR} = 10$% in Eq. 15. Comparing the fields of $\epsilon _{\delta }$ in Fig. 4b and c, a significant reduction in the error can be identified from the indirect approach. The advantage of the indirect approach here relies on the additional step of solving the Poisson equation. Comparing the noisy displacement fields shown in Fig. 4 and the field of $\epsilon _{{\overline{\delta }}}$ shown in Fig. 6a, a much homogenized error field can be immediately identified owing to the solving of the Poisson equation.

Fig. 6

Error distributions from solving the Poisson equation (a), from the AFH inversion in the indirect approach (b), and from the AFH inversion in the direct approach (c) with $\mathrm{SNR}=10$%

To comparing the two approaches quantitatively, the radial distributions of $\overline{\epsilon _{\delta }}(r)$, obtained by averaging $\epsilon _{\delta }(x,r)$ along the x direction, from both approaches with $\mathrm{SNR} = 10$% were shown in Fig. 7a. From the figure, the peak ${\overline{\epsilon }}(r)$ always appears near the central axis. In the indirect approach, ${\overline{\epsilon }}(r)$ is approximately only half of the value from the direct approach. This is consisted with the observations in Fig. 6b and c.

$\mathrm{SNR}$ values are further varied from 0 to 20% to clarify the effect of the level of $\mathrm{SNR}$ over $\epsilon _{\delta }$. A global error-index $\overline{{\overline{\epsilon }}}$, obtained by further averaging ${{\overline{\epsilon }}(r)}$ in the radial direction, was adopted. The performance of $\overline{{\overline{\epsilon }}}$ with respect to $\mathrm{SNR}$ is shown in Fig. 7b. As the $\mathrm{SNR}$ increases, $\overline{{\overline{\epsilon }}}$ from both the direct and indirect approach increases linearly with respect to $\mathrm{SNR}$. The value of $\overline{{\overline{\epsilon }}}$ obtained from the direct approach is about 2.5 times larger than the values from the indirect approach, confirming the superiority of the indirect approach with noisy displacement fields.

Fig. 7

Comparing ${\overline{\epsilon }}(r)$ with $\mathrm{SNR}=10$% in the radial direction (a) and $\overline{{\overline{\epsilon }}}$ with $\mathrm{SNR}$ from 0 to 20% (b) for both the direct and indirect approach

4.2 Optimal algorithms for indirect approach

In the indirect approach, there are several other widely adopted algorithms including OP, TPA, and FBPT to invert ${\overline{\delta }}(r)$ in addition to the indirect AFH. Details on each method and associated parameters can be found in the appendix. To select the most robust algorithm among the four in the indirect approach, associated ${\overline{\epsilon }}(r)$ distributions with ${\mathrm{SNR}} = 10$% are presented in Fig. 8a. Note the r coordinate is non-equally spaced from 0 to 0.1 and 0.9 to 1 to illustrate the comparison better. Overall, all four algorithms perform similarly in regions for $0.1\le r \le 0.9$. Noticeable difference can be identified especially near the central axis ($0<r<0.1$) and lateral boundary ($0.9<r<1$). TPA can be identified clearly to yield the lowest inversion error. This conclusion on the superiority of the TPA algorithm is consisted with previous exploration in a 1D problem (Dasch 1992). The level of $\mathrm{SNR}$ was again varied from 0 to 20% to characterize the impact of the $\mathrm{SNR}$ on the performance of each algorithm, and the behaviors of $\overline{{\overline{\epsilon }}}$ were compared in Fig. 8b. Generally, the superiority of the TPA algorithm can be confirmed throughout all $\mathrm{SNR}$. Nevertheless, other algorithms also perform well in the indirect approach indicating that the inclusion of the Poisson equation in the solution strategy is the essential step compared to the selection of the specific inversion algorithm.

Fig. 8

Comparing ${\overline{\epsilon }}(r)$ with $\mathrm{SNR}=10$% in the radial direction (a) and $\overline{{\overline{\epsilon }}}$ with SNR from 0 to 20% (b) for four different algorithms after solving the Poisson equation in the indirect approach

4.3 Importance of centering

All the analysis in previous sections assumes the flow is strictly axisymmetric. However, this is rarely the case. In practice, due to the imperfection of the setup and alignment, a slight asymmetry is usually present in flows, thus contributing to the inversion error. To evaluate the error introduced by the asymmetry quantitatively, we intentionally shifted the central axis in the analytical expression in Eq. 13 from 0 to 3 pixels. Then the TPA inversion algorithm with the indirect approach was applied, and results are shown in Fig. 9 at $x = 0.05$ for the illustrative purpose. Clearly, the error $\delta -\delta _{a}$ is proportional to the offset value. An auto-correlation algorithm was thus applied to correct such asymmetry. The basic principle behind the correction algorithm was to compare the left and right half profile of $|{{\varvec{{\varDelta }}}}|$ at specific x using the auto-correlation algorithm. Thus the lag between the two profiles could be estimated. Then the profile of ${{\varvec{{\varDelta }}}}$ is shifted with respect to the identified lag at each x. Finally, the left and right profiles are further averaged to achieve the exact symmetry.

Fig. 9

Quantitative comparison of the error $\delta -\delta _{a}$ induced by the asymmetric offset ranging from 0 to 3 pixels for $x=0.5$ in the analytical expression in Eq. 13

4.4 Helium jet experiment

A steady helium jet experiment was established to confirm the selection of the TPA inversion algorithm in the indirect approach. A reference background image was recorded firstly with no jet. After several minutes, a data image with a steady jet was subsequently recorded. The cross-correlation algorithm in DAVIS 8.4 was then applied to the reference/data image set to obtain the ${{\varvec{{\varDelta }}}}$ field as shown in Fig. 10a. Inside the arrows represent only the direction of ${\varvec{{\varDelta }}}$ while $|{\varvec{{\varDelta }}}|$ is denoted by the color. Spatially non-uniform spurious displacements (SDs) can be clearly identified in regions far away from the jet. Similar observations and discussions on the possible origins of the SDs were also proposed by the same group (Xiong et al. 2020). A heuristic approach was adopted here to remove the SDs via reconstructing a non-uniform field of SDs based on their distributions in regions away from the helium jet. Resulting filtered ${{\varvec{{\varDelta }}}}_{f}$ is shown in Fig. 10b and a direct comparison with the raw field ${{\varvec{{\varDelta }}}}_{\mathrm{raw}}$ in Fig. 10a immediately confirms the successful removal of the majority of the SDs. However, as indicated by the boundary contours marked in white lines, a slight asymmetry of the field still occurs due to the imperfections in the experiments. The symmetry correction algorithm discussed in Sect. 4.3 was thus applied, and the obtained centered ${{\varvec{{\varDelta }}}}_{c}$ field is shown in Fig. 10c.

Fig. 10

Two dimensional displacement fields without the correction on SDs (a), with the correction on SDs by subtracting the reconstructed spurious displacement field (b), with a further correction on the central axis location (c)

Obtained ${{\varvec{{\varDelta }}}}_{c}$ field is further utilized in the Poisson equation from the indirect approach for solving the ${{\overline{\delta }}}$ field. Then the TPA algorithm is applied to the ${{\overline{\delta }}}$ field to solve for the $\delta {r,y}$ field in the radial plane as shown in Fig. 11a. A steady jet contour with two diffusion layers can be identified. To confirm the accuracy of the inversion, the quantitative comparison of $\delta$ profiles with $\delta _{\mathrm{air}}$ and $\delta _{\mathrm{helium}}$ were conducted at $y=10$, 30, and 50 mm respectively as shown in Fig. 11b–d. A nearly perfect match between the $\delta$ with the reference values in the helium and airside confirms the effectiveness of the proposed BOS solution strategy for axisymmetric flows.

Fig. 11

The $\delta$ field obtained by the TPA from the indirect approach in the radial plane (a) and quantitative comparisons between the $\delta$ profiles with $\delta _{\mathrm{helium}}$ and $\delta _{\mathrm{air}}$ at $y=10$ mm (b), 30 mm (c) and 50 mm (d)

5 Conclusions

The specific inversion strategy and algorithm for axisymmetric BOS measurements have been systematically explored in this paper. A 3D synthetic experiment was established to mimic the helium jet discharged into the ambient air. With minimal noises appearing in the displacement measurements, solving the Poisson equation introduces additional computational cost and raises the error in reconstructing the radial refractive index field. Thus the direct approach is preferred in such a situation. As the noise within the measurement grows, e.g., with $\mathrm{SNR} > 2\%$, solving the Poisson equation reduces the solution error significantly thus, the indirect approach should be selected. The superiority of solving the Poisson equation with noisy displacement measurements can attribute to the fact that formulating BOS problem in solving a Poisson equation equals to find the solution to the BOS problem in the least-squared sense. Thus the high-frequency spatial noise can be smoothed out inherently. Within the indirect approach for $\mathrm{SNR} > 2\%$, the superiority of the TPA algorithm has been verified compared to OP, FBPT, and AFH algorithms widely used in the BOS community. The importance of accurately centering the displacement data was also emphasized, and an effective asymmetry correction algorithm was proposed. An experiment of a steady helium jet was conducted and confirmed the performance of the TPA algorithm in the indirect approach together with the centering algorithm.

Appendices

Appendices

Three inversion methods within the frame of the indirect approach will be introduced including OP, TPA and FBPT.

1.1 Onion-peeling method

Onion peeling is a classical algorithm to discretize Eq. (6) for axisymmetric flows. Axisymmetric flows are divided into N concentric rings of constant $\delta _j$ with $j=1,2,...,N$ as shown in the Fig. 12 (Dasch 1992):

$$\begin{aligned} \begin{aligned} {\bar{\delta }}(r_i)&= 2\sum _{j=i}^{N}\delta (r_j)\int _{r_j-{\varDelta }r/2}^{r_j+{\varDelta }r/2}\frac{ r}{\sqrt{r^2-r_i^2}}dr \\&= {\varDelta }r\sum _{j=i}^{N}W_{ij}\delta (r_j), \end{aligned} \end{aligned}$$

(23)

where the matrix coefficient $W_{ij}$ has the following form:

$$\begin{aligned} \begin{aligned} W_{i j}&= \sqrt{(2j+1)^2-4i^2}-\sqrt{(2j-1)^2-4i^2}, \\&\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \text {For} \ j>i. \\&= \sqrt{(2j+1)^2-4i^2}, \\&\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \text {For} \ j=i. \\&= 0, \\&\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \text {For} \ j<i. \end{aligned} \end{aligned}$$

(24)

Since this is discretation of the inverse Abel transform, Eq. 23 is recast back to the form of Eq. (19) with ${\mathbb {D}}_{ij}=(W^{-1})_{ij}$. The accuracy of onion peeling can be further improved by applying a regularization technique to the problem. A comprehensible explanation of how to implement Tikhonov regularization can be found in (Daun et al. 2006). The choice of the regularization is not trivial and is often found through trial and error. The following implementation of onion peeling does not include any regularization.

Fig. 12

Schematic of onion peeling algorithm. The axisymmetric flow is split into N concentric ring of constant $\delta _j$. Projection of the flow in the line-of-sight direction results to the distribution of ${{\bar{\delta }}}$

1.2 Three-point Abel

The method discretizes Eq. (8) with a quadratic expansion of ${\bar{\delta }}$ around each projection point. According to the comparison study (Dasch 1992), TPA possesses the best noise performance among all the tested methods under 1D conditions. A detailed derivation can be found in (Dasch 1992). The coefficient matrix ${\mathbb {D}}_{ij}$ of Eq. (7) can be expressed as (Karl Matthew Martin 2002; Dasch 1992) :

$$\begin{aligned} \begin{aligned} {\mathbb {D}}_ {ij}&= 0, \\&\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \text {For} \ j<i-1, \\&= H^0_{i,j+1}-H^1_{i,j+1}, \\&\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \text {For} \ j=i-1, \\&= H^0_{i,j+1}-H^1_{i,j+1}+2H^1_{i,j}, \\&\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \text {For} \ j=i, \\&= H^0_{i,j+1}-H^1_{i,j+1}+2H^1_{i,j}-H^0_{i,j-1}-H^1_{i,j-1}, \\&\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \text {For} \ j\le i+1, \\&= H^0_{i,j+1}-H^1_{i,j+1}+2H^1_{i,j}-2H^1_{i,j-1}, \\&\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \text {For} \ i=0,\ j=1. \\ \end{aligned} \end{aligned}$$

(25)

The coefficient matrices $H^0_{i,j}$ and $H^1_{i,j}$ take the following form:

$$\begin{aligned}&\begin{aligned} H^0_{ij}&= 0, \\&\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \text {For} \ j=i=0 \ \text {or} j<i, \\&= \frac{1}{2\pi }\ln \left[ \frac{\sqrt{(2j+1)^2 -4i^2}+2j+1}{2j}\right] , \\&\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \text {For} \ j=i\ne 0, \\&= \frac{1}{2\pi }\ln \left[ \frac{\sqrt{(2j+1)^2-4i^2} +2j+1}{\sqrt{(2j-1)^2-4i^2}+2j-1}\right] , \\&\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \text {For} \ j>i. \\ \end{aligned} \end{aligned}$$

(26)

$$\begin{aligned}&\begin{aligned} H^1_{ij}&= 0, \\&\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \text {For} \ j<i, \\&= \frac{1}{2\pi }\sqrt{(2j+1)^2-4i^2}-2jH^0_{ij}, \\&\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \text {For} \ j=i\ne 0, \\&= \frac{1}{2\pi }\left[ \sqrt{(2j+1)^2-4i^2}- \sqrt{(2j-1)^2-4i^2}\right] -2jH^0_{ij}, \\&\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \text {For} \ j>i. \\ \end{aligned} \end{aligned}$$

(27)

1.3 Filtered back projection tomography

A schematic of the projection of an arbitrarily distributed field $\delta (x,y)$ onto a plane (${{\bar{\delta }}}(\theta ,l)$) is shown in Fig. 13. The projection data ${{\bar{\delta }}}(\theta ,l)$ is obtained via the integration along a path with angle $\theta$ and a distance l from the center along the direction s. The relation between ${{\bar{\delta }}}(\theta ,l)$ and $\delta (x,y)$ is given by:

$$\begin{aligned}&{{\bar{\delta }}}(\theta ,l) \nonumber \\&\quad =\int _{-\infty }^{\infty }\int _{-\infty }^{\infty } \delta (x,y){\varLambda }(x\mathrm {cos}\theta +y\mathrm {sin}\theta -l)dxdy, \end{aligned}$$

(28)

where ${\varLambda }$ is the Dirac delta function. Equation (28) is also known as the Radon transform. The inverse Radon transform calculates $\delta (x,y)$ from the projection data ${{\bar{\delta }}}(\theta ,l)$. FBPT is a typical method to perform such an inversion. The method is based on the projection slice theorem, which states that the slice of angle $\theta$ in the frequency space of the 2D Fourier transform of $\delta (x,y)$ is equal to the 1D Fourier transform of the projection profile at angle $\theta$, ${{\bar{\delta }}}(\theta )$ as:

$$\begin{aligned} \begin{aligned} G(\varrho ,\theta )=F(\varrho \mathrm {cos}\theta , \varrho \mathrm {sin}\theta ). \end{aligned} \end{aligned}$$

(29)

$G(\varrho ,\theta )= FT_{1D} \left[ {{\bar{\delta }}}(l,\theta ) \right]$, $F(u,v)=FT_{2D} \left[ \delta (x,y)\right]$, where $FT_{1D}$ and $FT_{2D}$ are the 1D and 2D fast Fourier transforms respectively. By expressing the inverse 2D Fourier transform in cylindrical coordinates, one can show that:

$$\begin{aligned} \begin{aligned} \delta (x,y)&= \int _{-\infty }^{\infty }\int _{-\infty }^{\infty } F(u,v)e^{2\pi (ux+vy)}dudv \\&= \int _{0}^{2\pi }\int _{0}^{\infty }F(\varrho \mathrm {cos}\theta , \varrho \mathrm {sin}\theta )e^{2\pi l}\varrho d{\varrho }d\theta \\&= \int _{0}^{\pi } \left[ \int _{-\infty }^{\infty }| \varrho |G(\varrho ,\theta )e^{j2\pi \varrho l} d{\varrho }\right] d\theta . \end{aligned} \end{aligned}$$

(30)

Provided with the projection field $G(\varrho ,\theta )$, field $\delta (x,y)$ thus can be evaluated. To correct for the uneven sampling, various filters, such as the Shepp–Logan, Ram-Lak or the Hann filter have been devised. The choice of filter affects the trade-off between contrast and noise. In this study, the FBPT algorithm with the Shepp–Logan filter was adopted using the MATLAB function iradon.

Fig. 13

Schematic of filtered back projection technique