Quantifying Numerical Uncertainty in Background-Oriented Schlieren

Our study presents and evaluates a method for computing the numerical uncertainty in Background Oriented Schlieren (BOS). We use Richardson extrapolation to assess the uncertainty of numerical integration of density gradients, based on residuals of density results across two grid levels. By integrating this numerical uncertainty with the existing random uncertainty, we obtain the final uncertainty of the density field. We assess the method’s effectiveness using synthetic fields with artificial noise. Our error analysis shows that the sharpness of the density gradient significantly affects bias error and the prediction of numerical uncertainty. The prediction of numerical uncertainty corresponds to variations in bias error, particularly when the noise level and wavelength of the flow field are altered. By accounting for the numerical uncertainty, our method achieves up to 91% accuracy in predicting total uncertainty, as measured against the root-mean-square of the total error. We further demonstrate the utility of our methodology by applying it to experimental BOS images. Our proposed approach offers a more accurate understanding of uncertainty estimation in the BOS technique, with implications for future experiments.


Introduction
Background Oriented Schlieren (BOS) is an image-based density measurement technique based on the apparent distortion of a target pattern viewed through a medium with varying density (Meier 2002).The pattern distortion is due to the refraction of light rays traversing the medium and can be estimated using cross-correlation, tracking, or optical flow algorithms (Meier 2002).The density gradient can be calculated from the displacement field of the pattern based on the optical layout, and numerical integration of the gradient field provides a spatially resolved estimate of the density (Meier 2002).Due to its simple setup and ease of use, it is currently the most widely used method in quantitative Schlieren diagnostics (Raffel 2015).
The BOS measurement chain components can introduce or amplify errors.As shown in Fig. 1 Schematic of the BOS measurement chain.The red arrows indicate the processes that can introduce systematic errors, random error arises from the estimation of the pattern displacement.In contrast, systematic errors can be introduced during the pattern displacement estimation and the calculation and integration of the density gradient.For greater confidence in experimental results, knowledge of the possible systematic error range is important.Therefore, quantifying the uncertainty of BOS measurements is necessary.
A method was recently proposed to estimate the aposteriori, instantaneous, and spatially resolved density uncertainty for BOS (Rajendran et al. 2020c).The method estimates the uncertainty in the pattern displacement field from cross-correlation using the methods developed for Particle Image Velocimetry (PIV) (Sciacchitano 2019), which is then propagated through the measurement chain to estimate the random uncertainty of the measured density.Since then, methods have been proposed to estimate displacement uncertainties for the dot tracking method (Rajendran et al. 2019(Rajendran et al. , 2020b)), and to mitigate the noise sensitivity of the integration procedure by using an uncertainty-based weighted least squares (WLS) approach (Rajendran et al. 2020a).However, the existing uncertainty quantification methods do not account for the systematic uncertainty in the density field, which, based on previous results, could be larger than the random uncertainty, even when the bias error was smaller than the random error in the estimated displacement field (Rajendran et al. 2020c).This systematic error in the measured density is partly due to the discretization errors arising from the numerical integration of density gradients.This work aims to provide uncertainty on this systematic error, referred to as numerical uncertainty.
Numerical uncertainty due to discretization errors is a widely recognized issue in CFD studies and is commonly estimated using Richardson extrapolation (Oberkampf and Roy 2013).This work proposes a method to estimate the numerical uncertainty introduced by the density gradient integration in BOS measurements using Richardson extrapolation.The conventional method for CFD studies is modified to account for the limitations of experimental setups.An approach is also introduced to combine this systematic uncertainty with the random uncertainty to provide an instantaneous, spatially-resolved total uncertainty of the density estimates.The method is validated using synthetic density fields applied to synthetic BOS images and then tested against experimental data.

Numerical uncertainty estimation for BOS measurements with Richardson extrapolation
In BOS, the density-gradient integration is usually performed using a Poisson solver (Venkatakrishnan and Meier 2004),  = ( 2 ) −1 ( • ) = (  ) −1 (  ), (1) where  is the density field, ∇ is the density gradient field, and  is the gradient operator used for discretizing the derivative, which is also the source of numerical uncertainty considered in this work.Richardson extrapolation estimates the truncation error  ℎ in the numerical solution  ℎ obtained on a grid with spacing h as: where  is the actual field,  ℎ is the solution obtained using the same numerical scheme on a coarse grid with spacing rh.Here, r is the grid spacing ratio (set to 2 in the present study) and p represents the order of accuracy.For  ℎ to be accurate, the solutions  ℎ and  ℎ need to be in the asymptotic range, which requires the error to decrease with mesh refinement at the order of accuracy (Oberkampf and Roy 2013).An uncertainty bound (U), also referred to as the Grid Convergence Index (GCI) (Roache 1994), can be determined from the estimated error as: where   is the Factor of Safety, a scaling factor that provides an uncertainty bound for a predefined coverage (Roache 1997).The value of   is chosen empirically to provide a 95% confidence interval for a range of differential equations, discretization procedures, and boundary conditions commonly encountered in CFD problems (Roache 1994).In this study, we propose to use the absolute error from Equation ( 2) to estimate the standard numerical uncertainty defined as the standard deviation of the numerical error distribution.
As indicated in Equation ( 2), the order of accuracy p of the integration process is needed to apply Richardson extrapolation for numerical uncertainty estimation.The order of accuracy can be observed from the variation of the errors with two grid spacing levels: where   is the order of accuracy based on the errors.However, determining   is not feasible for experimental BOS because the true solution is unknown in experiments.The order of accuracy can also be observed from the residuals between three grid spacing levels: where   represents the order of accuracy based on the residuals.For   to be accurate, the solutions  ℎ and   2 ℎ on the coarse grids must be smooth and representative of the base solution  ℎ (Oberkampf and Roy 2013), which may not be possible with the resolution of experimental BOS data.Therefore, we propose to use the order of accuracy of the discretization scheme known as the formal order (   ) to approximate the actual order of accuracy of the numerical integration.In the present study, the order of accuracy is assumed to be 2 since the second-order central (SOC) difference scheme was used for the numerical discretization.The measurement noise in BOS leads to random uncertainty in the obtained density field, and it is desirable to combine the numerical uncertainty with the random uncertainty to provide a total uncertainty estimate.Consider the error of a given BOS measurement to be the sum of two random variables, each drawn from the bias and random error distributions of all possible experimental measurements.The variance of the total error distribution is the sum of the variances of the bias and random error distributions.By interpreting the numerical uncertainty as the standard deviation of the bias error distribution and the random uncertainty as the standard deviation of the random error distribution, the standard total uncertainty defined as the standard deviation of the total error distribution can be expressed as:

Synthetic fields
The proposed method is tested on the synthetic fields created from a hyperbolic tangent scalar field described as: where  represents the scalar field and  represents the wavelength.Such a field is chosen to reflect a real flow field in which BOS might be used, such as that involving a shock.
The field   is discretely sampled on a regular Cartesian grid with 81 × 81 points and grid spacing h.The range of the domain is -0.5 to 0.5 for both X and Y dimensions.For all plots and results, the wavelength normalized with respect to grid size ( n ) is used.The synthetic scalar and gradient fields with   = 0.2 are shown in Fig. 2. The density gradients are integrated using the Poisson solver with the SOC scheme and Dirichlet boundary conditions.The boundary conditions of zero uncertainty and the value of the analytical field are specified only at the left edge.This is done because, in most experiments, boundary condition at only one edge is known (free stream region).The error in the integrated field ( ℎ ) is determined as the deviation from the analytical solution .To evaluate the method's robustness in noise and validate the combination framework for the total uncertainty, the gradient fields are corrupted with noise drawn from a zero-mean normal distribution of a prescribed noise level relative to the maximum gradient in the field.Several (1000) realizations of the corrupted field are generated and integrated.The bias uncertainty is estimated for each realization, and the total uncertainty is calculated using Equation ( 6) with   estimated by propagating the uncertainty in the gradient fields to the density fields as in (Rajendran et al. 2020c).The bias error is calculated as the mean value of the error ( ℎ -), across all frames for each grid point.The random component is calculated as the standard deviation of the error across all trials, and the total error is the root mean square of the error across all trials.
Fig. 2 The synthetic field f and its gradient field with a n = 0.2

Experimental BOS images
The proposed method is demonstrated on measurements from a BOS experiment.The experiment is the supersonic flow over a wedge described in (Rajendran et al. 2020c).The BOS images were taken in a supersonic wind tunnel, of a Mach 2.5 flow over an 11.5• wedge with a base of 1 cm and a height of 2.5 cm.One image from the experiment is shown in Fig. 3.The density in the free stream was 0.49 kg/m 3 .The dot pattern consisted of 0.15 mm diameter dots randomly distributed on a transparency with about 25 dots per 32 × 32 pixel window.Further details of the experimental setup are provided in (Rajendran et al. 2020c).The images were processed using a multi-pass window deformation approach for two passes with identical window sizes and window resolution of 32×32 px 2 .The window overlap was set to 0%.The intermediate passes were smoothed but not validated to avoid the vectors in the shock region being identified as outliers (sharp changes are characteristic of this flow).The processing and uncertainty calculation was performed using an open-source code PRANA (https://github.com/aetherlab/prana/).Displacement uncertainty calculation from image matching is used (Rajendran et al. 2020c).After processing, the vectors are corrected to account for camera vibrations by subtracting the mean displacement in free stream region (which is not supposed to show any displacement vectors) from all the displacement vectors.
The depth-averaged density gradient field ∇ is calculated from the displacement field using the BOS optical model as: where Δ ⃗ is the displacement,  is the magnification of the dot pattern,   is the distance between the dot pattern and the mid-point of the density gradient field,  0 is the ambient refractive index,  is the Gladstone-Dale constant (= 0.225 × 10 −3 m 3 /kg for air) (Raffel 2015).The gradients are spatially integrated using the Poisson solver with a Dirichlet boundary condition of 0.49 kg/m 3 at the right edge (free stream) to obtain the projected density field.This is done for all 5000 images of the second pass.The displacement uncertainty is propagated similarly with a Dirichlet boundary condition of 0.0 kg/m 3 at the right edge.For each grid point, the ensemble average of numerical uncertainty is calculated, and the RMS of the random uncertainty is calculated, across all frames.Finally, the total uncertainty is calculated from Equation ( 6).

Analysis with synthetic fields
To determine the order of accuracy of the density gradient integration,   and   were calculated using Equations ( 4) and ( 5), respectively, and compared to the formal order (  ) from the density fields with  n varying from 0.1 to 1.0.As shown in Table 1, the formal order (  ) matched the error-based order (  ) for all cases, while the residual-based order (  ) began to deviate from the other two as the wavelength decreased.
Table 1 The order of accuracy determined using different methods for the numerical integration of the density gradient with a SOC scheme Fig. 4 shows the trend of the normalized (with respect to the maximum value of density) bias and random errors and uncertainties with changing wavelength and noise levels.For a  n of 0.1, the shock is captured by 3 pixels, and for  n = 0.4, the gradient is captured by 9 pixels (5).The relative bias error at  n = 0.1 is 95% lower than that at  n = 0.05.The random error and random uncertainty increase with both sharpness of the gradient and increasing noise level.The numerical uncertainty has a direct relation with noise but an inverse relationship with wavelength.5 shows how a sharper gradient (lower  n ) affects the bias error and numerical uncertainty prediction.
With a sharper gradient of  n = 0.1, the root mean square (RMS) of the total uncertainty is within 21% of the RMS of total error for 1% noise, and within 9% of the total error for 10% noise (Fig. 6).For  n = 0.4, (less sharp gradient), the total uncertainty is within 8% for both the lowest and highest noise levels.Another key observation is the difference in numerical uncertainty as noise level increases for the same wavelength.
For  n = 0.025, the relative increase in numerical uncertainty from 1% to 10% noise level is just 3% while for  n = 0.1 the increase is over 4 times that of the value at 1% noise level.For sharper gradients, as the bias error and numerical uncertainty are already quite high, their increase is relatively smaller.

Analysis with experimental BOS images
The density field obtained from the experimental BOS images and the numerical uncertainty results (normalized with respect to the maximum value in the density field) are shown in Fig. 7 and Fig. 8, respectively.The shock is captured across three pixels as can be seen from the density field.The resulting uncertainty distribution is bimodal.The peaks correspond to the random uncertainty before and after the shock respectively (Rajendran et al. 2020c).In this case the RMS of random uncertainty is greater than that of numerical uncertainty (Fig. 9).This may be because Image Matching is used for calculating the displacement field uncertainty, which reports a higher random uncertainty compared to other methods (Rajendran et al. 2020c) for this experiment.

Discussion and conclusions
In this study, we proposed and applied a method to estimate the numerical uncertainty for BOS density measurement.Based on the differences between the density results obtained on multiple grid levels, the numerical uncertainty is estimated using Richardson extrapolation.Compared to the framework of applying Richardson extrapolation to CFD simulations, several adaptations are made considering the limitations of experimental BOS measurements.First, the formal order of the numerical discretization is used as the order of accuracy of the numerical integration for performing Richardson Fig. 6 The density fields and spatial distributions of bias error and numerical uncertainty of the synthetic field a n =0.1 and n= 0.4 at 5% noise Fig. 7 The ensemble -averaged density field, with the shock indicated Fig. 5 The statistical distributions of the total error and uncertainty from the synthetic fields with a n of 0.1 and (a) 1% noise and (b) 10% noise extrapolation.Previous CFD studies showed that the observed order of accuracy can be lower than the formal order for the simulations of the flow fields governed by hyperbolic equations with strong discontinuities (Oberkampf and Roy 2013), thus the observed order of accuracy using Equation ( 4) or ( 5) are normally preferred.However, determining the order of accuracy from solution errors is not feasible as the true solution is normally unavailable for experimental BOS, and observing the order of accuracy from the solution residuals can fail due to the limited spatial resolution of the experimental data.Since the Poisson equation solved for the density gradient integration is elliptic, the order reduction concern may not apply to BOS.As shown in Table 1, the formal order matched with the order of accuracy observed from the solution errors, while the order of accuracy observed from the solution residuals failed for low-wavelength fields due to the insufficient spatial resolutions of the coarser grid levels, suggesting that the formal order is a suitable choice for the BOS numerical uncertainty estimation.Moreover, instead of determining the numerical uncertainty based on the empirically selected factor of safety (Fs), we propose to use absolute error obtained from Equation (2) to estimate the standard numerical uncertainty, which is validated by the comparison between the bias error with the numerical uncertainty from the synthetic fields as shown in Fig. 4. In addition, we proposed a framework to obtain the total uncertainty of the density measurement by combining the numerical uncertainty with the random uncertainty which is not seen in CFD simulations.
The analysis with the experimental BOS images and synthetic field test cases indicated that the contribution of the numerical uncertainty to the systematic error of the density field varies with the sharpness of the density gradient.Fig. 6 shows that the numerical uncertainty increases with the distance from the right boundary, for which the uncertainty was set to be 0. Additionally, the numerical uncertainty in the regions of the shock is elevated, which is expected due to the sharp density gradient present which cannot be well resolved.
There are several limitations of the proposed numerical uncertainty prediction method.First, the proposed method can be affected by the measurement noise, as seen from the synthetic field results in Fig. 4, where greater noise leads to the overprediction of bias errors.Additionally, shocks and discontinuities, like that present in the experimental images, invalidate the Taylor series approximation and hence also the Richardson Extrapolation (Roache 1994), leading to potential inaccuracy in estimating numerical uncertainty.It is also worth mentioning that the numerical integration is not the only source of bias error in the BOS measurement chain.Other sources, for instance, the error due to the simplified optical model for evaluating the density gradient from light ray displacement, will be investigated in future work.

Declarations Ethical Approval
This declaration is not applicable.

Fig. 1
Fig. 1 Schematic of the BOS measurement chain.The red arrows indicate the processes that can introduce systematic errors

Fig. 3 AFig. 4
Fig.3A sample image taken in the experiment with the region of interest indicated(Rajendran et al. 2020c)

Fig. 8
Fig.8The spatial distributions of the uncertainty and its components from the experimental BOS images

Fig. 9
Fig. 9 The statistical distributions of the total error and uncertainty from the experimental BOS images.The straight lines are the RMS values of the respective quantities