© V. Ramnath, published by EDP Sciences, 2017

1 Introduction

Traditionally in the field of pressure metrology uncertainty quantification (UQ) was performed with the use of the Guide to the Uncertainty in Measurement (GUM) however with the introduction of the GUM Supplement 1 (GS1) the use of Monte Carlo simulations has become an accepted practice for uncertainty analysis in metrology for mathematical models in which the underlying assumptions of the GUM are not valid. Consequently the use of quantile functions (QFs) was developed as a means to easily summarize and report on uncertainty numerical results that were based on Monte Carlo simulations.

Part of the challenge with the application of the use of quantile functions to modelling and summarizing Monte Carlo based uncertainty results is that it is a relatively recent development that is not widely known to metrology researchers and practitioners not directly involved in studying and investigating uncertainty quantification. An additional factor that partially impedes further widespread adoption of quantile functions is that the methodology used to implement the technique requires relatively advanced mathematical/statistical operations and associated software implementation routines when compared to the GUM.

The usage and implementation of these tools and routines may be perceived at least by some non-specialists to be either utilized from commercial software packages such as Matlab, or alternately specialist research codes developed in traditional computer languages such as Fortran or C++. As a result through a combination of some of these perceived factors the utilization of the GS1 and QF methodologies may enjoy less widespread adoption by some metrologists, although the use of the GS1 and QF methodologies has the potential to offer considerable benefits in scientific uncertainty analysis and the associated quality engineering activities in the field of metrology such as the consistent implementation and reporting of uncertainty results.

In this paper, we firstly consider the case of a piston–cylinder gas operated pressure balance where the effective area is modelled in terms of a combination of explicit/implicit and linear/non-linear models that are modelled, refined and solved with a combination of GUM and GS1 methodologies to produce uncertainty analysis data. Our approach incorporates some new developments in the area of pressure metrology through the refinement of the effective area estimates by using a modified weighted arithmetic average approach and through a parameter optimization of slip coefficients in the extended Navier–Stokes equations for transitional flow regimes. Then having developed and performed the uncertainty analysis to generate the associated GUM/GS1 based uncertainty data we then secondly demonstrate how quantile functions may be applied to analyse and compare uncertainty analysis results using some practical implementation simplifications to the underlying algorithms with open source software packages and tools.

2 Mathematical models

The first model that we consider for a pressure balance is the simplest and is based on the assumption that the piston and cylinder have perfectly straight profiles such that the effective area may be calculated as(1)where ⟨r⟩ and ⟨R⟩ are effective radii values for the piston and cylinder respectively. In this arithmetic model ⟨r⟩ and ⟨R⟩ are usually in pressure metrology practice simply estimated as Graybill–Deal estimators, i.e. weighted arithmetic averages for a sequence of values {x₁, …, x_n} such that the expected value y and the variance u²(y) of the expected value y are given by

(2a)(2b) where u(x_i) are the corresponding standard uncertainties. In practice the dimensional measurements of the piston/cylinder geometry are usually averaged diameter measurements along set distances along the length of the piston/cylinder for a sequence of angles reported in a calibration certificate in which case only a single data-set of the radii data is available and the standard weighted mean approach is an appropriate approximation. Other strategies for cylindrical surface measurements include that of generatrix, birdcage and helical line approaches as discussed by Galovska et al. [1]; however, these approaches are in most cases unnecessary for pressure metrology since the pressure balance interface gap is very small when compared to the engagement length.

The main issues with the use of a Graybill–Deal estimator which finds widespread usage by many metrologists are that the estimators variance is underestimated, there are large mean square errors, and that the intervals formed by this estimator have a poor coverage probability for the mean. These limitations of Graybill–Deal estimators in the case of a single measurand are only applicable in the context of the simplest pressure balance model to the extent that the perfectly straight piston/cylinder geometry assumption is valid, however in a practical context it is well known that there are generally fluctuations in the piston/cylinder radii profiles along the pressure balance's engagement length and as a result the assumption of single effective values of ⟨r⟩ and ⟨R⟩ is only an approximation. Whilst more refined estimators have been developed the validity of the use of these alternative estimators assumes a large number of independent laboratory measurements for the underlying data, and as a result the use of Graybill–Deal estimators may not necessarily be considered a reasonable practical first approximation for ⟨r⟩ and ⟨R⟩, and we comment that when access to independent repeat dimensional measurements is available that alternative higher accuracy formulations may be used as discussed by Zhang [2] and other researchers.

Regardless of the level of refinement for the simplest pressure balance model calculations for ⟨r⟩ and ⟨R⟩ the use of the GUM will result in a Gaussian probability density function (PDF) for the measurand which in the case of a univariate model Y = f(X₁, …, X_N), with N random variables where y is the estimate of the actual value Y and x_i are estimates of X_i consistent with statistical sampling from the corresponding probability density functions g_i(ξ_i), i = 1, …, N of X_i, is then of the form

(3a)(3b)(3c) where u_c(y) is the combined standard uncertainty, c_i = ∂f/∂x_i are the sensitivity coefficients of the measurand function and is an estimate for the covariance between two quantities q and r where and with similar expressions for and s²(r) as originally discussed in the GUM [3].

The expanded uncertainty U = ku_c(y) is then obtained by multiplying the combined uncertainty u_c(y) by a coverage factor k_p for a specified probability level p and the general steps are to first calculate the effective degrees of freedom ν_eff with the aid of the Welch–Satterthwaite formula for uncorrelated quantities as(4)where u_i(y) = (∂f/∂x_i)u(x_i) such that where ν_i are the associated degrees of freedom for input x_i which is ν = n − 1 for a single quantity with n observations, ν = n − m for a least-squares fit of m parameters with n data points, and ν ≈ 1/2[Δu/u]⁻² where [Δu/u] is the relative uncertainty for x_i in other cases where further statistical information is not available.

The conventional Welch–Satterthwaite formula that estimates the effective degrees of freedom for a GUM based uncertainty analysis where there is an absence of correlation between the input random variables, such as for example piston/cylinder dimensional measurements, is usually adequate for most pressure balance effective area calculations that use gases such as nitrogen as a working fluid at low/medium pressures however in pressure balances that use fluids such as some silicon oils at medium/high pressures where correlation effects may be present in certain model inputs, such as for example the oil temperature and viscosity where correlations are implicitly present through the particular fluid's thermodynamic equation of state, a generalized effective degrees of freedom ν_gen developed by Willink [4] may be substituted for ν_eff such that

(5a) (5b)

Once the effective degrees of freedom are approximated the next step is to obtain the corresponding t_p-factor for a desired level of confidence p from tables of a Student's t-distribution where if ν_eff is not an integer one may either interpolate or truncate ν_eff to the next lowest integer and then take the coverage factor as k_p ≡ t_p(ν_eff) so that the expanded uncertainty is U = k_pu_c. An alternative to the use of tables when calculating the coverage factor is a direct numerical solution of the integral equation in order to determine the unknown t_p such that

(6a)(6b) where Γ(t) is the Gamma function as used in the construction of the Student's t-distribution.

The second model for the pressure balance that we consider for the effective area A₀ incorporates a continuum flow assumption within the interface gap following the well known Dadson model [5] such that the effective area is(7)where the pressure utilizes an explicit equation(8)where p₁ is the inlet/source pressure, p₂ the outlet/sink pressure, r(z) the radius of the inner piston, R(z) the radius of the outer cylinder, u(z) = r(z) − r₀ the fluctuation in the piston radial profile, U(z) = R(z) − R₀ the fluctuation in the cylinder radial profile where r₀ and R₀ are the piston and cylinder radii at the inlet respectively, h(z) = R(z) − r(z) is the interface gap, and p_c(z) denotes the pressure profile along the engagement length that is based on a continuum flow model.

Although this model utilizes an explicit mathematical equation for the effective area due to the fact that it is in the form of an integral equation it cannot be easily reduced to an equation that is directly amenable to a traditional GUM based uncertainty analysis that utilizes sensitivity coefficients. For this reason the uncertainty for a pressure balance's effective area that utilizes the Dadson model is computed with the use of the GS1 approach [6] in the case of a univariate model Y = f(X) as conceptually illustrated in Figure 1 and which is applicable for a pressure balance model that determines a single value of the effective area.

When performing a Monte Carlo simulation random numbers ρ are first sampled from a rectangular distribution ρ ∈ R[0, 1] which may then be used to generate sampled values z for other PDF distributions such as from the rectangular distribution R[a, b] which is just z = a + (b − a)ρ, from a Gaussian distribution N(μ, σ²) with a mean μ and variance σ² which is just z = μ + σρ, or for a sampled value ξ from an arbitrary PDF with a continuous distribution function G_X which is obtained from solving G_X(ξ) = ρ.

Further technical details in terms of nomenclature, and recommended statistical algorithms to implement the GS1 may be consulted in the official documentation however we briefly comment that when calculating coverage intervals for a required coverage probability p the conventional approach is to first construct a variable α such that 0 ≤ α ≤ (1 − p). In the event that the PDF of the measurand g_Y(η) where η is a random variable of possible values of Y is non-symmetric one then finds the value of α such that is a minimum where G_Y(η) is the discrete approximation to the measurand's distribution function so that the end points of a 100p% coverage interval for Y are and .

Alternately a choice of α = 1/2(1 − p) will give a coverage interval for a probabilistically symmetric 100p% coverage interval, or if the PDF g_Y(η) is asymmetric then the shortest 100p% coverage interval is given with a α value that satisfies , and we thus remark that there is a certain element of subjectivity when analysing and reporting summaries of Monte Carlo based GS1 uncertainty quantification results.

The third model that we consider is that for mixed continuum/transitional flows that is specified in terms of a non-linear integral equation that was developed by Pitakarnnop and Wongpithayadisai [7] which is an improvement of the original Dadson model and is of the form(9)

This is an implicit equation for the pressure such that

(10a)(10b)(10c)(10d)(10e) where p_t(z) is the pressure profile corresponding to mixed continuum/transitional flow regime, G_p the reduced flow rate, δ the rarefaction parameter, α the tangential momentum accommodation coefficient which is usually set to α = 1 under the assumption of a diffuse reflection for the interaction at the gas/solid interface, and Kn is the Knudsen number at a particular distance along the engagement length L where A_α = 1.1466 and A_β = 0.6470 for flow regimes such that Kn ≤ 0.4 as discussed by Hadjiconstantinou [8].

When calculating the Knudsen number by convention we use a mean free path length λ as indicated above from standard rarefied gas dynamics practice where k_B = 1.3806488 × 10⁻²³ J K⁻¹ is the Boltzmann constant, T is the thermodynamic temperature of the pressure balance gas which by convention is taken at a temperature of t = 20 °C, m = M/N_A is the mass of the gas molecule which for nitrogen gas has a molar mass of M = 28.01348 × 10⁻³ kg mol⁻¹ where N_A = 6.022140857 × 10²³ mol⁻¹ is the Avogadro constant, and for simplicity the gas viscosity is calculated in terms of a Sutherland model as a first approximation of the form(11)which in the case of nitrogen gas have values of T₀ = 273 K, S_μ = 107 K, and μ₀ = 1.663 × 10⁻⁵ N s m⁻² respectively as discussed in White [9].

Since the Pitakarnnop and Wongpithayadisai equation for the mixed continuum/transitional pressure profile is in terms of a non-linear integral it is seen that a non-linear implicit mathematical model occurs for the pressure balance effective area and as a result we have to utilize a GS1 approach to work out the corresponding effective area for the continuum/transitional flow model. The calculation of the effective area is based on a free-body analysis of the constituent forces acting on the piston such that the effective area A_eff is calculated as where the force components are specified as , , , and as discussed by Sabuga [10], i.e. the various flow regime models are used to calculate the constituent force terms such as the bottom/upwards pressure force F_b from the inlet/source pressure p₁, the top/downwards pressure force F_t from the outlet/sink pressure p₂, the retarding pressure frictional force F_f, and the vertical component of the hydrostatic pressure for the particular pressure/flow conditions along the engagement length. Once all the various force terms are calculated from a choice of flow regime model they are then implemented in a free-body diagram analysis using vector superposition of forces acting on the piston to calculate A_eff, and as a result in order to specify A_eff requires a corresponding solution of pressure profile for the choice of flow regime model.

Since the continuum/transitional flow regime model is non-linear integral equation one may either the method of successive approximations or alternately a direct Nystrom solution. In this paper we opt for the method of successive approximations by solving the integral equation from the general non-linear form

(12a)(12b) where φ_i(x) are successive approximations for i = 1, 2,…, N and where we simply use the original continuum flow Dadson model as an initial starting solution so that convergence is reached relatively quickly. We utilize this practical approach since we wish to reduce the numerical solution time for the Pitakarnnop and Wongpithayadisai model in order to maximize the potential number of Monte Carlo simulation events for the subsequent GS1 uncertainty analysis of model 3.

Since the form of the Pitakarnnop and Wongpithayadisai equation is a type of hybrid Volterra/Fredholm integral equation it will exhibit numerical instability if the outlet sink pressure p₂ is too low as this will result in the Knudsen number being too large. In order to stay within the validity limits of the continuum/transitional flow model which has an upper Knudsen number limit of Kn ≤ 0.4 we restrict the outlet/sink pressure to p₂ < 25 kPa and use an inlet/source pressure of p₁ = 185 kPa.

We have considered three measurement models for the pressure balance's effective area as firstly for a purely geometrical arithmetic approximation, secondly as A₀ for a continuum flow approximation using the Dadson model, and thirdly as for a continuum/transitional flow approximation using the approach of Pitakarnnop and Wongpithayadisai. Of the three models considered only the first model is amenable to a GUM based uncertainty analysis whilst the other two models must use incorporate Monte Carlo simulations using the GS1 approach in order to estimate the uncertainty.

The application of quantile functions (QFs) thus offers a relatively simple conceptual approach that avoids certain subjective aspects when Monte Carlo simulation results are utilized in order to analyse and compare measurement results and uncertainties as recently developed by Harris et al. [11] where the intrinsic features are that a quantile function Q(ρ) satisfies the relationship(13)where g(η) is the PDF of the measurand Y and G(η) the corresponding distribution function such that g(η) = G′(η) as previously discussed, and where the random variable η is determined as

(14) (15)

Based on the properties of quantile functions four useful consequences reported by Harris et al. that arise are that:

• Q(ρ) is a random draw from the PDF of Y if ρ is a random draw from the rectangular distribution R[0, 1].

• Coverage intervals for a probability p are [Q(α), Q(p + α)] where α is a parameter such that 0 ≤ α ≤ (1 − p) and the corresponding probabilistic symmetric coverage interval is [Q(p₁), Q(p₂)] with p₁ = 1/2(1 − p) and p₂ = 1/2(1 + p).

• The expectation E(Y) and variance V(Y) are calculated as and respectively.

• Any function Q(ρ) that is increasing on the interval is a legitimate quantile function for some corresponding PDF.

To exploit the above properties in order to avoid the subjectivity of analysing/summarizing Monte Carlo simulation data for uncertainty analysis requires either a choice of quantile family, or alternately the use and application of quality computer software for validated and verified (V&V) numerical analysis results from a quality engineering perspective. Of the various numerical algorithms involved it may be observed that the main essential underlying requirement is the use of a good quality spline interpolation routine.

Earlier work undertaken by Willink [12] investigated the use of extended lambda distributions (ELDs) as one particular class of quantile functions which we will utilize in the form

(16a) (16b)

The main observations from a choice of ELD as a QF are that the dimensionless parameters a and b determine the shape of the PDF whilst the parameters c and d have the same dimensions as the random variable η and influence the scale and location of the distribution. In addition the associated PDF to the QF is of the form(17)where the constraints on the parameters a, b, c, d to ensure statistical consistency are that

(18a)(18b) which result in the following limits of g(η) as

(19a) (19b) (19c)

The mean μ and variance σ² for the above choice of ELD are then

(20a) (20b)

Different possible approaches were originally investigated by Willink to solve for the QF parameters a, b, c, d and in this paper we opt for the approach of fitting the parameters using four quantiles for i = 1, 2, 3, 4 where in general p₁ < p₂ < p₃ < p₄. Based on numerical experiments Willink deduced that in the special case where

(21)(22) that a unique solution exists where

(23)(24) and the parameter b is numerically solved from the non-linear equation

(25) (26)

Once the value of b is determined then the remaining parameters may be calculated as(27)(28)(29)

Since the above equations are valid for any appropriate choice of quantiles a particular choice of(30)will be used as per the earlier recommendations of Willink in the subsequent numerical simulations. We comment that the parameters when determined from the above choice of approach are by definition a “best-fit” curve fit however the specification of the parameter values by themselves is incomplete since estimates of the associated parameters should also technically be specified for completeness.

Due to the fact that a GS1 approach is not valid for a regression analysis as discussed by Elster and Toman [13] the parameter uncertainties corresponding to the above estimates may be determined either by a Markov chain Monte Carlo (MCMC) approach as discussed by Forbes [14] or alternately by a conventional Levenberg–Marquardt approach as discussed by Press et al. [15]. Following the discussion by Saunders [16] the parameters a = [a₁, …, a_N] for a weighted least squares y(x) curve fit with M data points may be obtained by minimizing the merit function

(31) (32)

Since the earlier work by Willink has already implemented an alternative approach to determine the best-fit parameter a = (a, b, c, d) for the QF system under consideration we do not need to minimize the merit function since we already have a solution as previously outlined and can then use the original merit function to estimate the covariance matrix C by calculating the inverse of the matrix α for the parameters such that

(33)(34) where the diagonal matrix elements of C will indicate the variance of the parameters, i.e. C_i,i = Var(a_i) and the off-diagonal matrix elements will indicate the respective covariances as C_i,j = Cov(a_i, a_j) respectively.

Fig. 1 Illustration of GUM Supplement 1 (GS1) methodology.

3 Numerical simulations

We investigate the application of the use of quantile functions for pressure balance effective areas by considering the geometry data specified in Figure 2 which was previously reported by Ramnath [17]. Due to the fact that all of the three models that we have considered utilize a radial symmetry assumption we assume constant representative values for corresponding piston/cylinder radii along the engagement length by using the standard uncertainties for the external diameter d corresponding to the piston and the internal diameter D corresponding to the cylinder as u(D) = 100 nm and u(d) = 170 nm respectively. Since the radial symmetry assumption specifies that r = 1/2d and R = 1/2D it then follows that u(r) = ±0.085 μm and u(R) = ±0.050 μm for 0 ≤ z ≤ L. Making the additional assumption that dimensional measurements along the engagement length L are approximately we may reasonably take an upper bound for an engagement length of L = 30 mm so that u(z) = ±0.058 μm for all points along the engagement length.

As our estimates for the dimensional uncertainty terms u(r), u(R), and u(z) all emanate from data from experiments and inter-comparisons that utilize the GUM for the uncertainty estimates it follows that statistical sampling from these estimates must follow the associated underlying probability density function distributions. By definition GUM based results are reported in terms of either a Student's t-distribution or a Gaussian distribution, and as a result we opt to statistically sample from the corresponding standard normal distribution N(μ;σ²) as ξ = μ + z where z is a draw from the corresponding standard Gaussian distribution N(0, 1) as discussed in the GS1 documentation.

For the first model we have that so in the absence of correlation, i.e. by making the assumption that cylinder radial measurements do not influence the piston cylinder radial measurements and vice versa by setting r(⟨r⟩, ⟨R⟩) =0 it follows that the combined standard uncertainty for the first model is(35)

The main issues with the use of Graybill–Deal estimators, i.e. the use of weighted means for data sets {x₁, …, x_n} with an estimate and associated variance for the estimate with a weighting factor as commonly used by many metrologists is that if the actual variance is approximated as with the sample variances then there is an under-estimate of the variance of the weighted mean , if is too small then the weighting factors become too large, and that if the number of samples n is too small, e.g. if n ≈ 8 as discussed by Rukhin [18] then the Graybill–Deal estimators significantly under-estimate the true variance.

In our particular case for a piston–cylinder pressure balance it is very unlikely that a large number of independent piston/cylinder radii measurements are available and for practical physical equipment calibration reasons the typical number of radii measurements for the piston with a length of L_p = 60 mm are usually around n_p = 20, and around n_c = 14 for a cylinder with a length of L_c = 33 mm for typical large cross-sectional area pressure balances that use gas as a working fluid in the range 5 ≤ p/[kPa] ≤ 200.

Further developments after the earlier work by Zhang on the use of Graybill–Deal estimators have been undertaken by Rukhin in the context of weighted mean statistics for inter-laboratory studies who investigated maximum likelihood estimates and one particular option that was proposed to mitigate against the issues of Graybill–Deal estimators is to use what is known as the DerSimonian and Laird procedure which modifies the weighting factors through the introduction of a small positive constant y_DL such that the weighting factors are now instead of . The optimum value of y_DL is first calculated as(36)where p is the number of independent measurements which have associated data {x_i, s_i} for i = 1, …, p and is the Graybill–Deal estimate calculated as(37)

Then the Graybill–Deal estimate and associated variance are modified to calculate the corresponding DerSimonian and Laird estimate and variance as

(38) (39)

Due to the fact that we assume an absence of correlation between the physical dimensions of the piston and cylinder effective radii ⟨r⟩ and ⟨R⟩ which is a reasonable assumption since piston radii measurements do not influence cylinder radii measurements and vice versa it follows that the Welch–Satterthwaite formula may be used to calculate the corresponding effective degrees of freedom. The above aspects will be implemented for the first model for results with a standard uncertainty with a confidence level of p = 95.45%, noting that the GUM based approach results for the effective area associated with model 1 will technically follow a Student's t-distribution which is very similar to that of a normal distribution for the associated degrees of freedom as calculated with the Welch–Satterthwaite formula. In most practical cases it is unnecessary to utilize an exact mathematical formula for a non-centred Student's t-distribution since it is relatively straightforward to apply appropriate scaling/shifting of the corresponding symmetric Student's t-distribution, however exact mathematical formulae for the corresponding probability density function (PDF) and cumulative distribution function (CDF) in terms of hyper-geometric functions are included in Appendix A for the sake of completeness.

Following Pitakarnnop and Wongpithayadisai [7] the velocity slip model for two parallel walls is modelled as(40)where u is the tangential fluid velocity, is the tangential wall velocity, n is the normal vector at the point of contact between the fluid and the wall, α is a tangential momentum accommodation coefficient, λ is the mean free path of the fluid molecules, and A_α and A_β are constants that characterize the second order velocity slip model. Some typical values for A_α and A_β that were reported by Pitakarnnop and Wongpithayadisai of various researchers for a variety of fluid flows are indicated below.

In our case we are interested in using the velocity slip model to calculate the pressure profile within the interface of the piston–cylinder pressure balance which is modelled as two parallel walls. The main physical quantity that is utilized to calculate the pressure profile is that of the reduced flow rate which is specified as where δ is the rarefaction parameter defined as . We can utilize this information to estimate the corresponding uncertainties for A_α and A_β by examining the behaviour of G_p(δ) for parallel walls that is predicted by the second order velocity slip model with results predicted by the linearised Boltzmann equation using the well known Bhatnagar–Gross–Krook (BGK) approximation for the collision integral terms for isothermal flows.

For Poiseuille flow between two parallel plates assuming perfectly diffuse accommodation at the plate surfaces the BGK–Boltzmann equation reduces to the integral equation system

(41a)(41b)(41c)(41d) as originally discussed by Lo and Loyalka [19] where d = 1/Kn is the reciprocal of the Knudsen number, and T_n(x) is the corresponding Abramowitz function valid for integers n ≥ −1 and real numbers such that x ≥ 0.

The solution of this integral equation whilst computationally demanding is in principle amenable to a solution with an arbitrary accuracy if sufficient computational resources are available. In practical quality engineering terms the required accuracy for this integral equation will be specified through a combination of the desired metrological accuracy of a pressure scale realization and the physical accuracy limits of the dimensional measurements of the pressure balance's piston and cylinder radii (Fig. 3).

The main numerical challenge in solving this integral is due to the presence of singularities. Aspects of this problem were previously addressed by Li et al. [20] who numerically solved the linearised BGK equation for steady Couette flow for Knudsen numbers in the range 0.003 ≤ Kn ≤ 10 for two parallel plates which move in opposite directions however their results are in terms expressed in a slightly different form which is more amendable to the flow rate for a half channel width. In addition their results do not provide explicit expressions for the Chebyshev power series expansions for the Abramowitz functions in the BGK–Boltzmann equation and as a result we cannot directly utilize their results.

Later numerical simulations undertaken by Jiang and Luo [21] also in the context of Couette flow similar to that studied by Li et al. for the flow distribution again reported results of a slightly modified integral equation for a Knudsen number range (using their symbols) of 0.003 ≤ k ≤ 100. Their final results accurate to twelve decimal places were reported for the half-channel velocity and half-channel flow rate.

As a result we also again cannot directly use their simulation results however the key advantage by Jiang and Luo is that they provide a convenient to use formula for the Abramowitz functions which is documented in Appendix B for the sake of completeness. For comparison in order to illustrate the computational complexity to solve the linearised BGK–Boltzmann equations it was observed by Jiang and Luo in their paper that Yap and Sader [22] used a Gauss–Legendre quadrature scheme to study steady and oscillating Couette flows which resulted in a corresponding linear system of (12800)² = 163840000 simultaneous equations which they solved with arithmetic packages that offered resolutions of 30 digits.

Unfortunately even with these extreme computational measures they were still only able to achieve accuracies between 5 and 8 significant digits for Knudsen numbers in the range 0.01 ≤ k ≤ 100. Based on these observations in many cases it is therefore not feasible to directly compute solutions to non-continuum flows in terms of modifications to the Boltzmann equation unless one has access to substantial computational resources.

As a result based on these practical considerations in this paper we opt to simply utilize the existing reported data of Lo and Loyalka in our analysis noting that the original results of Lo and Loyalka for the flow rate Q(d) of two parallel plates corresponds to the particular situation of a piston cylinder interface gap with stationary walls. This particular dataset which is known to be accurate to five decimal places is sufficiently accurate for our present purposes in this paper and may be used to determine the corresponding fitness and accuracies of the parameters A_α and A_β in the extended Navier–Stokes equation system developed by Pitakarnnop and Wongpithayadisai in model 2.

The original model of Pitakarnnop and Wongpithayadisai was valid in a Knudsen number range 0 ≤ Kn ≤ 0.4 since it used values of A_α = 1.1466 and A_β = 0.6470 that were previously obtained by Hadjiconstantinou. Some typical values of the Knudsen number Kn in terms of variation of pressure for a characteristic length of L_C = 1 ×10⁻⁶ m, i.e. when the interface gap is approximately h ≈ 1 μm are indicated in Table 1.

From the data in Table 1 it is seen that the current upper validity of the Hadjiconstantinou corresponds to a outlet/sink pressure of approximately p₂ ≈ 16.91 kPa and as a result the Pitakarnnop and Wongpithayadisai model is then valid for mainly gauge pressure scale realizations.

Although the variation in a pressure balance's effective area as computed with gauge pressure and absolute pressure measurements is usually minimal in many cases an absolute pressure measurement with an inlet/source pressure of p₁ ∼ 100 kPa and an outlet/sink pressure of p₂ ≲ 10 Pa.

In practice the outlet/sink pressure with vacuum pumps is usually pumped down to smaller pressures but in principle the desired physical characteristic when using pressure balances in absolute mode is to have the outlet pressure as small as possible. As a result it will be beneficial to optimize the parameter values A_α and A_β to encompass a wider Knudsen number range to accommodate a lower outlet pressure and we opt for a compromise of an upper Knudsen number of Kn = 13 which roughly corresponds to an outlet pressure of p₂ = 0.5 kPa for an interface gap of h = 1 μm in order to illustrate the general calculation approach to optimize A_α and A_β and quantify their optimal fit parameter uncertainties (Fig. 4).

Utilizing the regression model of Saunders with a = [A_α, A_β]^T it follows that(42)where it is assumed that the tangential momentum accommodation coefficient is α = 1. Since we do not have specific knowledge of the uncertainties of d and Q(d) with the exception of the dataset by Lo and Loyalka we assume constant uncertainties so that unity weighting factors apply. Based on this approach where we set for M data points the merit function is then of the form so that(43)

The above merit function is then optimized in terms of the parameters A_α and A_β for the specified Knudsen number range by solving the simultaneous system ∂χ²/∂A_α = 0 and ∂²χ²/∂A_β = 0 where(44)and(45)

Results for the unweighted optimization for a Knudsen number range of 0.002 ≤ Kn ≤ 13.0 are illustrated in Figure 5 based on the previously reported Lo and Loyalka data for 0.002 ≤ Kn ≤ 100 from which it is seen that the unweighted optimized fit more closely matches the linearised BGK–Boltzmann solution at higher values of d and diverges somewhat at lower values of d, whilst the Hadjiconstantinou fit whilst not the best fit at higher values of d generally performs reasonably well over the entire range of d values. We comment that the data that is plotted in Figure 5 is only a qualitative indication since the quantitative contributions in the χ² merit function are cumulatively larger in magnitude at larger inverse Knudsen number d values, i.e. the overall best fit is achieved when the bulk of the data more closely matches the expected Poiseuille flow rate Q(d) data at medium and larger d values and that a few outlying deviations at lower values of d does not significantly alter the magnitude of the χ² merit function.

In order to interpret the physical significance of the respective fits we note that d = 1/Kn and that the bulk of the flow within the interface gap will be in the transitional flow regime apart from a small section near the inlet pressure p₁ which will be in the continuum flow regime, and a small section near the outlet pressure p₂ which may possibly be in the molecular flow regime respectively. Since a larger portion of the flow will be in the transitional flow regime it will generally be advantageous to utilize a fit that more closely matches the linearised BGK–Boltzmann solution at higher values of d since even if there are any discrepancies at lower values of d, i.e. higher Knudsen numbers near the outlet since the effect of these possible discrepancies is less physically significant due to the fact that there will be substantially smaller viscous drag forces exerted on the surface of the piston near the outlet. These physical observations may be loosely restated as that it is more important to adequately mathematically model the flow effects near the inlet region and in the intermediate regions than it is near the outlet regions based on the physical considerations that viscous drag effects are quantitatively larger and more significant in those regions. We comment that the relative importance of the inlet and intermediate region for the parameter determination may if necessary be incorporated through an appropriate selection of weighting factors for the inverse Knudsen numbers d_i when constructing the merit function χ²([A_α, A_β]^T).

Referring to the data in Table 1 where it is observed that p ∝ d we may then infer that there will be larger viscous forces and corresponding frictional drag forces which are present and which taken together would overall contribute a larger contribution to the resultant force acting on the piston when the pressure is larger. Based on this physical reasoning the force components which are used in the calculation of the pressure balance's effective area are then more significant when the pressure is larger. Since p ∝ d it follows that if the weighting factors are similarly proportional to the pressure which is achieved by setting by for example setting where A is a constant of proportionality and B is a constant power to refine the magnitudes, then the χ² merit function optimization to determine the “best” values of the parameters A_α and A_β in the context of a pressure metrology analysis for a pressure balance's effective area may be determined in a weighted optimization search. We elaborate on this approach to determine the weighting factors in more technical detail in Appendix C for the sake of completeness.

In practical terms the Hadjiconstantinou fit whilst exhibiting slight discrepancies at smaller Knudsen numbers when compared to the optimized fit over a wider Knudsen number range 0.002 ≤ Kn ≤ 13.0 is nevertheless reasonably accurate for most practical pressure metrology purposes if the upper limit of the maximum Knudsen number Kn = 0.4 is taken into account. In our particular case this amounts to strictly limiting the outlet pressure to an approximate value of p₂ = 16.91 kPa for a nominal pressure balance interface gap of h ≈ 1 μm, since the parameter fit of Hadjiconstantinou in the Pitakarnnop and Wongpithayadisai model will significantly diverge if Knudsen numbers larger than Kn = 0.4 corresponding to pressures substantially lower than p₂ = 16.91 kPa are applicable.

Under these circumstances the corresponding variances and covariances in the parameters A_α = 1.1466 and A_β = 0.6470 for the Hadjiconstantinou fit may be approximately estimated as(46)(47)(48)

The above partial derivatives may be approximated through the use of central differences for a multi-variable function u(x, y) as

(49)(50) as discussed by Tannehill et al. [23]. The parameters in the Pitakarnnop and Wongpithayadisai model then have the following approximate uncertainties

(51) (52) (53)

We now have sufficient information to generate the corresponding PDFs g(η) for in model 1, η = A₀ in model 2, and in model 3 respectively. Numerical simulations are performed in the open source software packages Gnu Octave 4.0.0 for the arithmetic and Monte Carlo calculations which are then post-processed using the Python(x, y) 2.7.10.0 distribution and scipy based interpolated univariate spline (IUS) and graphical sub-routines in order to implement the GUM and GS1 approaches.

From these plots the PDF is first normalized such that noting that the effective area η = A_eff/[m²] units are formally expressed in square meters in order to be dimensionally consistent. Then the PDFs are post-processed in order to generate the cumulative distribution function for easy visualization comparison purposes (Fig. 6).

Noting that η = Q(ρ) and ρ = G(η) we can simply use the previous cumulative distribution functions to determine the values for i ∈ [1, 2, 3, 4]. When this is implemented we then have three sets of non-linear equations specified as for each of the three models which must be solved for the parameter b.

Whilst these equations should ideally be numerically solved using a routine such as Newton's method in practical terms since we already know that b < 1.4 we can simply plot the function(54)over a suitable range of possible values such that b_min ≤ b < 1.4 and choose the value that is closest to zero that would formally solve the non-linear equation f(b) = 0.

This process is illustrated in Figure 7 for the three models with f_i(b_i) for i = {1, 2, 3} which can be used to estimate b₁, b₂ and b₃ for each of the respective three models as indicated in Table 2.

Utilizing the data in Table 2 then allows us to calculate the expected value μ for each of the models using the determined quantile function parameters such that

(55)(56)(57) where the associated probability density functions are illustrated in Figure 8 which are obtained with 0.001 ≤ ρ ≤ 0.999. Referring to Figure 8 we observe that the PDFs for each of the three models are completely and accurately specified by just the specification of the four quantile parameters a, b, c and d, and that the expanded uncertainty can now be directly computed with the associated PDF g(η) specified in terms of the QF parameters a, b, c and d for any desired confidence interval.

As a result the provision of a extended lambda distribution based QF parameterization may in principle provide equivalent and additional information over that of a GUM based uncertainty analysis in measurement inter-comparisons and calibration certificates which traditionally only report a measurand's expected value, variance, and expanded uncertainty.

Fig. 2 Interface gap data of pressure balance.

Fig. 3 Poiseuille flow rate data computed by Lo and Loyalka [19] for parallel plates for non-continuum flow regimes.

Fig. 4 Illustration of a comparison of flow rate data predictions between actual results obtained by Lo and Loyalka using the linearised BGK–Boltzmann equation and corresponding results obtained by Hadjiconstantinou (2003) with [Aα = 1.1466, Aβ = 0.6470], Deissler (1964) with [Aα = 1, Aβ = 1.1250], and Cercignani (1964) with [Aα = 1.1466, Aβ = 0.9756] in an extended Navier–Stokes equation formulation for a Knudsen number range 0.005 ≤ Kn ≤ 0.650.

Fig. 5 Illustration of unweighted optimized fit for parameters Aα and Aβ in the Pitakarnnop and Wongpithayadisai model for a Knudsen number range 0.002 ≤ Kn ≤ 13.

Fig. 6 Illustration of cumulative distribution functions associated with respective models.

Fig. 7 Illustration of determining the b parameter for the respective models.

Fig. 8 Representative PDFs for a gas pressure balance operated in gauge model obtained with quantile function parametrizations for 0.001 ≤ ρ ≤ 0.999 using M = 250 Monte Carlo simulation events for model 2 and model 3 respectively to illustrate the ability of quantile functions to model arbitrary PDFs that may exhibit asymmetry and/or skewness.

4 Discussion

In this paper we have considered a pressure balance's effective area A_eff which is specified by an explicit linear model that can be analysed using the GUM, and an explicit non-linear model and an implicit non-linear model which both have to analysed by the GS1 method. All of the three models were numerically solved in order to determine A_eff which produced PDF data for each of the respective models.

The research work that was utilized to solve the models incorporates new developments in the area of pressure metrology through the refinement of the effective area estimates by using a modified weighted arithmetic average approach for model 1, and through a parameter optimization of slip coefficients in the extended Navier–Stokes equations for transitional flow regimes in model 3. Both of these research developments offer the potential of new insights into the benefits and limitations of using each of three pressure balance models in different situations, particularly when results are to bench marked against the conventional Dadson approach of method 2.

Having developed some refinements of the physical models used in pressure metrology we then demonstrated how quantile functions may be applied to analyse and compare uncertainty analysis results. Our research approach in this paper incorporated some practical implementation simplifications to the underlying algorithms for the utilization and application of quantile functions to uncertainty analysis where we adopted open source software packages and tools, and in the process we have presented a practical example of how the GS1 and QF methodologies may be applied by metrologists in other domains and fields.

5 Implications and influences

The mathematical research work reported in this paper through the refinement and uncertainty quantification of the slip coefficient parameters in the extended Navier–Stokes equations will have practical benefits in the field of gas pressure metrology as now metrologists will have the information to further study, investigate, test and calibrate the differing behaviors and characteristics of gas pressure balances when operated in gauge mode.

In addition the reported process to determine the slip parameter uncertainties may in turn be used by computational fluid dynamics practitioners to more accurately model and design instruments which exhibit slip behavior such as spinning rotor gauge based vacuum pressure instruments in the aerospace and space industry sectors.

Finally the statistical research reported in this paper may be used to present the case for the wider adoption of quantile function techniques in uncertainty analysis and the utilization of open source software tools in metrology.

Appendix A Exact mathematical expressions for GUM PDF and CDF

The GUM methodology for uncertainty results for a measurement model may either correspond to a Gaussian distribution or alternately a scaled/shifted Student's t-distribution depending on the magnitude of the effective degrees of freedom. Although the GUM does not specifically provide an explicit analytical expression for a non-centred scaled/shifted Students t-distribution the corresponding probability density function may nevertheless be explicitly mathematically specified in terms of a confluent hyper-geometric function ₁F₁(a;b;z) such thatwhere x is a random variable, μ is the expected value, and σ is a variance as discussed by Lenth [24]. More detailed statistical information may be consulted by referring to work by Hogben et al. [25] and Van Aubel and Gawronski [26].

The generalized hyper-geometric function _pF_q(a;b;z) for vectors a = [a₁, a₂, …, a_p]^T and b = [b₁, b₂, …, b_q]^T has the properties that it is convergent for any ||z|| <∞ or if ||z||<1 and p = q + 1, is divergent for any z ≠ 0 if p > q + 1, and is defined in terms of a series expansion aswhere (x)_n = (Γ(x + n))/(Γ(x)) is the Pochammer symbol. The corresponding cumulative distribution function (CDF) for a non-central Student's t-distribution may be expressed aswhere Γ(z) is the gamma function and I(a;b;z) is the regularized incomplete beta function defined as

The above exact expressions for the GUM based PDF and CDF for the measurand y = f(x₁, …, x_n) whilst mathematically complex are in most practical cases not directly utilized since the underlying PDF may be simply recovered from the Student's t-distribution by using the transformed variable t = (x − μ)/σ and applying appropriate scaling/shifting for qualitative visualization purposes.

Appendix B Approximate mathematical expressions for the Abramowitz functions

The Abramowitz functions are defined by

The Abramowitz functions were earlier computed by Macleod [27] by using a Chebyshev expansion and Li et al. specify a power series expansion for J₁(x) aswhere γ is the Euler–Mascheroni constant.

The above formula for J₁(x) may then using the properties of Abramowitz functions be used to calculate J₋₁(x) through a combination of derivatives and application of the recurrence relationships as indicated below.

• Work out J₁(x) using the power series expansion

• Use finite differences to estimate J₀(x) =− (dJ₁/dx) as

• Use the recurrence formula with n = 3 to work out J₃(x) as

• Use finite differences to work out J₂(x) =− (dJ₃/dx) as

• Use the recurrence formula with n = 2 to work out J₋₁(x) as

Appendix C Technical details for weighting factors in a non-continuum flow optimization

The weighted χ² merit function is defined aswhere M is the number of data points for the chosen Knudsen number flow range, specified as 0.002 ≤ Kn ≤ 13 in the main paper for illustration purposes, Q_i is the Poiseuille flow rate and d_i the corresponding inverse Knudsen number respectively. As the original paper by Lo and Loyalka where it was reported that they solved the linearised BGK–Boltzmann equation to an estimated five decimal places it follows that u(Q_i) = ±10⁻⁵ for i = 1, …, M.

We comment that the accuracy of the solution Q(d) should in fact be tested against other numerical solutions in order to infer the actual accuracy, and that this is only technically possible if there are other data sets available for the d and Q(d) solution of the linearised BGK–Boltzmann equation. In practical terms the accuracy of the Poiseuille flow rate Q(d) would be mathematically calculated in a similar manner in which key comparison reference values (KCRV's) are normally calculated, and this is only possible if all the data sets have associated estimates for the corresponding uncertainties. Since no other uncertainty information is currently easily accessible for other Poiseuille flow rate data we opt to assume unity weighting factors in the main paper when performing the optimization to determine the values of the parameters A_α and A_β. Our approach of assuming constant weighting factors to determine the inferred uncertainties u(A_α) and u(A_β) errs on the side of caution due to the lack of additional verified quantitative information on the uncertainties in d and Q(d).

In the event that verified information on the uncertainties arising from the numerical solution of the linearised BGK–Boltzmann equation does become available then one may perform a weighted optimization to minimise the merit function χ². This may be accomplished if we determine the corresponding uncertainties for the inverse Knudsen numbers d_i.

Utilizing conventional practise from the field of rarefied gas dynamics we have that the mean free path of a gas λ is approximated aswhere μ is the gas viscosity, p the gas pressure, T the gas absolute temperature, k_B the Boltzmann constant, and m is the mass of the gas molecules which may be calculated as m = M/N_A where M is the particular gas species molecular weight and N_A is the Avogadro constant. By simple rearrangement of the above equation we immediately have that

In principle the above formula for the inverse Knudsen number is of the form d = f(p, μ, M, k_B, N_A, T) and as a result formal application of the GUM methodology will then specify the corresponding uncertainty for d as

The uncertainties of the physical constants for the Boltzmann constant k_B and the Avogadro number N_A viz. u(k_B) and u(N_A) may be obtained from the CODATA internationally recommended values of the physical constants (see for example http://arxiv.org/pdf/1507.07956 for data supplied by the NIST for the physics and metrology community). Referring to the CODATA values we have k_B = 1.38064852 × 10⁻²³ J K⁻¹, u(k_B) = 0.00000079 × 10⁻²³ K⁻¹, N_A = 6.022140857 × 10²³ mol⁻¹, and u(N_A) = 0.000000074 × 10²³ mol⁻¹ respectively.

The uncertainties for thermo-physical properties in the particular case of gases such as nitrogen and argon may be obtained by Lemmon and Jacobsen [28] who specify for nitrogen gas the values M = 28.01348 × 10⁻³ kg and a viscosity at a temperature of 300 K of μ = 17.8771 × 10⁻⁶ Pa s. Typically in pressure metrology practice measurements are performed or calculated for a temperature of 20 °C =293.15 K and as a result thermo-physical properties at a temperature of 300 K are adequate in most practical circumstances. Based on the discussion by Lemmon and Jacobsen an approximate uncertainty for the viscosity of nitrogen is 0.5%, and as a result we may reasonably specify the thermo-physical property uncertainties as u(M) = 0.00001 × 10⁻³ and u(μ) ≈ 0.0893855 × 10⁻⁶ Pa s.

The uncertainties for the experimental measurements of the pressures and temperatures for each of the inverse Knudsen number values are also necessary and we assume for illustration purposes that the standard uncertainties for the pressure and temperature are u(p) = 1/2[1 + (2.5 × 10⁻⁵)p] Pa and that u(T) = 0.050 K respectively. Under these assumptions we obtain the uncertainties listed in Table C.

If the uncertainty in d is modelled as u(d) ∼ Ad^B then by fitting a straight line through log(d) versus log(u(d)) data points such that log(u(d)) = m log(d) + c it then follows that u(d) = exp(m log(d) + c) for a final fitted approximation (Fig. C).

Fig. C Illustration of approximations to model the inverse Knudsen number u(d).

Once the uncertainties u(d_i) have been determined then the weighting factors may be calculated aswhere signifies the uncertainty in the Poiseuille flow rate, previously assumed as Q(d_i) = 10⁻⁵, which is fixed as previously discussed and signifies the uncertainty for the inverse Knudsen number which is calculated as using the straight line fit parameters.

References

All Tables

All Figures

	Fig. 1 Illustration of GUM Supplement 1 (GS1) methodology.
In the text

	Fig. 2 Interface gap data of pressure balance.
In the text

	Fig. 3 Poiseuille flow rate data computed by Lo and Loyalka [19] for parallel plates for non-continuum flow regimes.
In the text

Fig. 4 Illustration of a comparison of flow rate data predictions between actual results obtained by Lo and Loyalka using the linearised BGK–Boltzmann equation and corresponding results obtained by Hadjiconstantinou (2003) with [Aα = 1.1466, Aβ = 0.6470], Deissler (1964) with [Aα = 1, Aβ = 1.1250], and Cercignani (1964) with [Aα = 1.1466, Aβ = 0.9756] in an extended Navier–Stokes equation formulation for a Knudsen number range 0.005 ≤ Kn ≤ 0.650.

In the text

	Fig. 5 Illustration of unweighted optimized fit for parameters Aα and Aβ in the Pitakarnnop and Wongpithayadisai model for a Knudsen number range 0.002 ≤ Kn ≤ 13.
In the text

	Fig. 6 Illustration of cumulative distribution functions associated with respective models.
In the text

	Fig. 7 Illustration of determining the b parameter for the respective models.
In the text

	Fig. 8 Representative PDFs for a gas pressure balance operated in gauge model obtained with quantile function parametrizations for 0.001 ≤ ρ ≤ 0.999 using M = 250 Monte Carlo simulation events for model 2 and model 3 respectively to illustrate the ability of quantile functions to model arbitrary PDFs that may exhibit asymmetry and/or skewness.
In the text

	Fig. C Illustration of approximations to model the inverse Knudsen number u(d).
In the text