Application of quantile functions for the analysis and comparison of gas pressure balance uncertainties

Traditionally in the field of pressure metrology uncertainty quantification 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 was developed as a means to easily summarize and report on uncertainty numerical results that were based onMonte Carlo simulations. In this paper, we considered the case of a piston–cylinder operated pressure balance where the effective area is modelled in terms of a combination of explicit/implicit and linear/non-linear models, and how quantile functions may be applied to analyse results and compare uncertainties from a mixture of GUM and GS1 methodologies.


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.

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 A jj 0 may be calculated as 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 fx 1 , ..., x n g such that the expected value y and the variance u 2 (y) of the expected value y are given by 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 1 , ..., 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 (j i ), i = 1, ..., N of X i , is then of the form sðq; rÞ ¼ 1 nðn À 1Þ where u c (y) is the combined standard uncertainty, c i = ∂f/ ∂x i are the sensitivity coefficients of the measurand function and sðq; rÞ is an estimate for the covariance between two quantities q and r where q ¼ ð1=nÞ P n k¼1 q k and s 2 ðqÞ ¼ ð1=ðn À 1ÞÞ P n j¼1 ðq j À qÞ 2 with similar expressions for r and s 2 (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 n eff with the aid of the Welch-Satterthwaite formula for uncorrelated quantities as where u i (y) = (∂f/∂x i )u(x i ) such that n eff P N i¼1 n i where n i are the associated degrees of freedom for input x i which is n = n À 1 for a single quantity with n observations, n = n À m for a least-squares fit of m parameters with n data points, and n % 1/2[Du/u] À2 where [Du/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 n gen developed by Willink [4] may be substituted for n eff such that n gen ¼ max f1; n eff 0 g; ð5aÞ 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 n eff is not an integer one may either interpolate or truncate n eff to the next lowest integer and then take the coverage factor as k p ≡ t p (n eff ) so that the expanded uncertainty is U = k p u 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 ∫ t p Àt p fðu; n eff Þ du ¼ p; where G(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 0 incorporates a continuum flow assumption within the interface gap following the well known Dadson model [5] such that the effective area is where the pressure utilizes an explicit equation where p 1 is the inlet/source pressure, p 2 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 0 the fluctuation in the piston radial profile, U(z) = R(z) À R 0 the fluctuation in the cylinder radial profile where r 0 and R 0 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 r are first sampled from a rectangular distribution r ∈ 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)r, from a Gaussian distribution N(m, s 2 ) with a mean m and variance s 2 which is just z = m + sr, or for a sampled value j from an arbitrary PDF with a continuous distribution function G X which is obtained from solving G X (j) = r.
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 a such that 0 a (1 À p). In the event that the PDF of the measurand g Y (h) where h is a random variable of possible values of Y is non-symmetric one then finds the value of a such that ½G À1 Y ðp þ aÞ À G À1 Y ðaÞ is a minimum where G Y (h) is the discrete approximation to the measurand's distribution function so that the end points of a 100p% coverage interval for Y are G À1 Y ðaÞ and G À1 Y ðp þ aÞ. Alternately a choice of a = 1/2(1 À p) will give a coverage interval for a probabilistically symmetric 100p% coverage interval, or if the PDF g Y (h) is asymmetric then the shortest 100p% coverage interval is given with a a value that satisfies g Y ðG À1 Y ðaÞÞ ¼ g Y ðG À1 Y ðp þ aÞÞ, 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 This is an implicit equation for the pressure such that Kn where p t (z) is the pressure profile corresponding to mixed continuum/transitional flow regime, G p the reduced flow rate, d the rarefaction parameter, a the tangential momentum accommodation coefficient which is usually set to a = 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 a = 1.1466 and A b = 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 l as indicated above from standard rarefied gas dynamics practice where k B = 1.3806488 Â 10 À23 J K À1 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 À3 kg mol À1 where N A = 6.022140857 Â 10 23 mol À1 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 which in the case of nitrogen gas have values of T 0 = 273 K, S m = 107 K, and m 0 = 1.663 Â 10 À5 N s m À2 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 A Ã 0 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 [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 1 , the top/downwards pressure force F t from the outlet/sink pressure p 2 , the retarding pressure frictional force F f , and the vertical component of the hydrostatic pressure F v 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 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 2 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 2 < 25 kPa and use an inlet/source pressure of p 1 = 185 kPa.
We have considered three measurement models for the pressure balance's effective area as firstly A jj 0 for a purely geometrical arithmetic approximation, secondly as A 0 for a continuum flow approximation using the Dadson model, and thirdly as A Ã 0 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(r) satisfies the relationship where g(h) is the PDF of the measurand Y and G(h) the corresponding distribution function such that g(h) = G 0 (h) as previously discussed, and where the random variable h is determined as Based on the properties of quantile functions four useful consequences reported by Harris et al. that arise are that: -Q(r) is a random draw from the PDF of Y if r is a random draw from the rectangular distribution R[0, 1].
-Coverage intervals for a probability p are [Q(a), Q (p + a)] where a is a parameter such that 0 a (1 À p) and the corresponding probabilistic symmetric coverage interval is [Q(p 1 ), Q(p 2 )] with p 1 = 1/2(1 À p) and -Any function Q(r) that is increasing on the interval r∈I ¼ ½0; 1 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 h ¼ QðrÞ; ð16aÞ 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 h and influence the scale and location of the distribution. In addition the associated PDF to the QF is of the form where the constraints on the parameters a, b, c, d to ensure statistical consistency are that which result in the following limits of g(h) as The mean m and variance s 2 for the above choice of ELD are then 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 u p i ¼ Qðp i Þ for i = 1, 2, 3, 4 where in general p 1 < p 2 < p 3 < p 4 . Based on numerical experiments Willink deduced that in the special case where that a unique solution exists where and the parameter b is numerically solved from the nonlinear equation Once the value of b is determined then the remaining parameters may be calculated as Since the above equations are valid for any appropriate choice of quantiles a particular choice of ½p 1 ; p 2 ; p 3 ; p 4 ¼ ½0:025; 0:25; 0:75; 0:975; ð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 1 , ..., a N ] for a weighted least squares y(x) curve fit with M data points may be obtained by minimizing the merit function 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 a for the parameters such that a k;ℓ ¼ 1 2 where the diagonal matrix elements of C will indicate the variance of the parameters, i.e. C i,i = Var(a i ) and the offdiagonal matrix elements will indicate the respective covariances as C i,j = Cov(a i , a j ) respectively.

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 mm and u(R) = ±0.050 mm for 0 z L. Making the additional assumption that dimensional measurements along the engagement length L are approximately uðLÞ ¼ 1=2 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð0:1Þ 2 þ ð2LÞ 2 q Â 10 À6 m we may reasonably take an upper bound for an engagement length of L = 30 mm so that u(z) = ±0.058 mm 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(m;s 2 ) as j = m + 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 A jj 0 ¼ fð⟨r⟩; ⟨R⟩Þ ¼ ðp=2Þ½⟨r⟩ 2 þ ⟨R⟩ 2 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 The main issues with the use of Graybill-Deal estimators, i.e. the use of weighted means for data sets fx 1 , ..., x n g with an estimatex ¼ ½ P i w i x i =½ P i w i and associated variance for the estimate VarðxÞ ¼ 1=½ P i w i with a weighting factor w i ¼ s À2 i as commonly used by many metrologists is that if the actual variance s 2 i is approximated as s 2 i % s 2 i with the sample variances s 2 i then there is an under-estimate of the variance of the weighted mean VarðxÞ, if s 2 i 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 w i ¼ ½y DL þ s 2 i À1 instead of w i ¼ 1=½s 2 i . The optimum value of y DL is first calculated as where p is the number of independent measurements which have associated data fx i , s i g for i = 1, ..., p andx GD is the Graybill-Deal estimate calculated as Then the Graybill-Deal estimate and associated variance are modified to calculate the corresponding DerSimonian and Laird estimate and variance as 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 A jj 0 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 where u is the tangential fluid velocity, u w is the tangential wall velocity, n is the normal vector at the point of contact between the fluid and the wall, a is a tangential momentum accommodation coefficient, l is the mean free path of the fluid molecules, and A a and A b are constants that characterize the second order velocity slip model. Some typical values for A a and A b 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 . We can utilize this information to estimate the corresponding uncertainties for A a and A b by examining the behaviour of G p (d) 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 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 halfchannel 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) 2 = 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 a and A b in the extended Navier-  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 2 % 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 1 ∼ 100 kPa and an outlet/sink pressure of p 2 ≲ 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 a and A b 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 2 = 0.5 kPa for an interface gap of h = 1 mm in order to illustrate the general calculation approach to optimize A a and A b and quantify their optimal fit parameter uncertainties (Fig. 4).
where it is assumed that the tangential momentum accommodation coefficient is a = 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 w i ¼ 1 for M data points the merit function is then of the form so that 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 x 2 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 x 2 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 1 which will be in the continuum flow regime, and a small section near the outlet pressure p 2 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 w i for the inverse Knudsen numbers d i when constructing the merit function 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 F v ¼ 2pr 0 ∫ L 0 p½r 0 ðzÞ dz and corresponding frictional drag F f ¼ À2pr 0 ∫ L 0 ðh=2Þ½p 0 ðzÞ dz 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 w i ∝d i by for example setting w i ¼ Ad B i where A is a constant of proportionality and B is a constant power to refine the magnitudes, then the x 2 merit function optimization to determine the "best" values of the parameters A a and A b 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 w i 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 2 = 16.91 kPa for a nominal pressure balance interface gap of h % 1 mm, 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 2 = 16.91 kPa are applicable.
Under these circumstances the corresponding variances and covariances in the parameters A a = 1.1466 and A b = 0.6470 for the Hadjiconstantinou fit may be approximately estimated as C ¼ 1 a 11 a 22 À a 12 a 21 a 21 Àa 12 Àa 21 a 11 !
; ð47Þ The above partial derivatives may be approximated through the use of central differences for a multi-variable function u(x, y) as as discussed by Tannehill et al. [23]. The parameters in the Pitakarnnop and Wongpithayadisai model then have the following approximate uncertainties uðA a Þ % 1:224 Â 10 À4 ; ð51Þ CovðA a ; A b Þ % 8:996 Â 10 À14 : We now have sufficient information to generate the corresponding PDFs g(h) for h ¼ A 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 GðhÞ ¼ ∫ h À∞ gðjÞ dj for easy visualization comparison purposes (Fig. 6).
Noting that h = Q(r) and r = G(h) we can simply use the previous cumulative distribution functions to determine the values u p i ¼ Qðp i Þ for i ∈ [1,2,3,4]. When this is implemented we then have three sets of non-linear equations specified as ½p b 1 À ð1 À p 1 Þ b =½p b 2 À ð1 À p 2 Þ b ¼ j 0 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 This process is illustrated in Figure 7 for the three models with f i (b i ) for i = f1, 2, 3g which can be used to estimate b 1 , b 2 and b 3 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 m for each of the models using the determined quantile function parameters such that where the associated probability density functions are illustrated in Figure 8 which are obtained with 0.001 r 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(h) 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.

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.

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. inferred uncertainties u(A a ) and u(A b ) 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 x 2 . 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 l is approximated as where m 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, 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 À23 J K À1 , u(k B ) = 0.00000079 Â 10 À23 K À1 , N A = 6.022140857 Â 10 23 mol À1 , and u(N A ) = 0.000000074 Â 10 23 mol À1 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 À3 kg and a viscosity at a temperature of 300 K of m = 17.8771 Â 10 À6 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 À3 and u(m) % 0.0893855 Â 10 À6 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 À5 )p] Pa and that u(T) = 0.050 K respectively. Under these assumptions we obtain the uncertainties listed in Table C.   V. Ramnath: Int. J. Metrol. Qual. Eng. 8, 4 (2017) 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).
Once the uncertainties u(d i ) have been determined then the weighting factors may be calculated as where u Q i signifies the uncertainty in the Poiseuille flow rate, previously assumed as Q(d i ) = 10 À5 , which is fixed as previously discussed and u d i signifies the uncertainty for the inverse Knudsen number which is calculated as u d i ¼ exp ðm log ðd i Þ þ cÞ using the straight line fit parameters.