Experimental and numerical investigations on the shape and roughness of cylindrical critical flow venturi nozzles (CFVN)

This paper describes a part of a research project aimed at reducing the uncertainty of the determination of the CFVN discharge coefficient. In the current state of the work, there are two parallel studies: an experimental and a numerical. Experiments have been carried out to determine the impact of roughness on the CD and to demonstrate the need for dimension parameters, such as diameter and roughness coefficient characterisation. Numerical simulations are used in the same CFVN geometry with a smooth surface to determine the main parameters which can influence the numerical prediction of the discharge coefficient and the flow behaviour. The results show that an increase of roughness may lead to a significant decrease of the discharge coefficient, particularly in the high Reynolds number range.


Introduction
Natural gas has a growing role as a primary energy source. In the last 10 years, the production of natural gas has increased worldwide by more than 20%. In earlier procedures, flow meters were calibrated by air and this method has now been extended to natural gas applications. The range of gas flow rate meters should be between 10 to 200 m 3 /h, with pressures of up to 70 bar. Obviously, the Reynolds number range for these applications exceeds the upper limit given by the ISO 9300 standard [1]. Therefore, this paper is part of a project whose aim is to extend the Reynolds number limit to over 1×10 7 . Transport and distribution operations of natural gas outside the range require uncertainties below 0.1%.
The goal of the project is to reduce the uncertainty of primary and secondary standards based on CFVNs by investigation and analysis of the impact of the nozzle wall roughness on the discharge coefficient of the flow. This study will be concerned with experiments on roughness and will also involve numerical simulations. Currently, an increasing number of researchers are taking notice of the effect of roughness on the discharge coefficient. Major research organizations are also realizing that roughness has a certain influence on the discharge coefficient. However, most of the research is currently being done with toroidal geometry while the second standard, the cylindrical standard, lacks investigation. Further investigation of roughness effects on the discharge coefficient are required to improve the ISO 9300 reference equations.

General consideration 2.1 On the discharge coefficient of CFVNs
Assuming a critical condition at the throat means sonic conditions, (1) The mass flow through a convergent-divergent nozzle is approximately constant, subject to viscous effects. Its theoretical value is computed from the expression (1). The real flow in the nozzle throat is not isentropic due to viscosity and is not one-dimensional, due to the nozzle shape as the wall has a curvature radius. Therefore, formula (1) must be completed by the introduction of two coefficients, as in (2).
The inclusion of the so-called discharge coefficient CD (2) corrects the effect of geometry and the effect of viscosity on the flow [2], and the non-uniformity of the flow in the throat.
where C* represents the critical flow function. It is usually defined as the ratio between the actual flow rate and the theoretical flow rate (3).
This paper focuses on the cylindrical nozzle throat [1] whose geometry is illustrated in Figure 1. The discharge coefficient is commonly expressed using three constants: As the discharge coefficient is partially influenced by the gas viscosity, it clearly depends on the Reynolds number ReD in the nozzle. An appropriate definition was made in (5).
where μ 0 is the molecular viscosity evaluated at the upstream stagnation conditions. Table 1 summarizes some of the proposed coefficient values for the coefficients "a", "b" and "n" for equation (4). Usually the "n" coefficient is taken as 0.5 for laminar development of the boundary layer and between 0.2 and 0.13 for turbulent, as is the case in the equation for the development of the boundary layer on a flat plate. For the cylindrical nozzle shape, extensive research was carried out during the early 1970s and, after this period, research began to orient towards the toroidal shape. One reason for this was because of the higher CD coefficients. In France, the cylindrical nozzle shape is still in use, although research and knowledge dating from the 1970s is difficult to access. This highlights that similar research still has to be carried out for a cylindrical reference equation for the CD equation calculations.

Shape of the cylindrical nozzle
As seen in Figure 1, the cylindrical nozzles are manufactured with respect to the shape of the inner contour, the convergent inlet with a curvature radius of rc = d (d the diameter of nozzle throat), the cylindrical throat of a constant diameter d over the same length and, finally, the conical diffusor with a length of 7•d and an angle of 7 degrees.  Figure 1) in the nozzle throat in the cylindrical region (HEXAGON™ measurement) For this study, nozzles with three different diameters have been investigated. Table 2 gives the nominal throat diameter and information on each cylindrical nozzle investigated. For each diameter size, three different roughness levels were considered. In order to investigate the influence of different surface qualities on the main flow rate characteristics, the nozzle designs allow the inner flow channels to have average roughness heights Ra between 0.4μm and 1.2μm.

Standard facilities used for the flow rate measurements
At this time, CFVN must be calibrated before it can be used. To ensure reliable and accurate results, the primary volume and flow rate standards of PTB-PIGSAR TM and CESAME-EXADEBIT were applied. In cooperation with, and under the supervision of, PTB, PIGSAR TM is responsible for maintaining and disseminating the German reference unit of volume for natural gas under high-pressure conditions. The PIGSAR TM test facility is described in detail in references [7] and [8] and some of its characteristics are summarized in Table  3. It uses natural gas at pressures between 16 and 50 bar. Working fluid Natural gas with uncertainty of C* estimated at 0.065%, (k = 2) and molar mass uncertainty estimated at 0.1% (k = 2)

PTB-PIGSAR TM
The traceability is based on the primary method of the geometrically-measured volume of a high-pressure piston prover and gas density measurements made with a buoyancy balance. The natural gas used comes from the Groningen region in the North Sea. The critical flow factor, C*, is calculated using AGA8-DC92 state equations.
The uncertainty of C* is also calculated with the GERG2004 state equations. The uncertainty for the calibration value of a sonic nozzle at PIGSAR TM can be claimed as 0.15% (0.13% to 0.16% depending on the flow rate) (k=2). The test facility was designed for measurements up to 90bar as pressure and 480 m 3 /h as a flow rate. The HPPP setup is pictured in Figure 2 with the line employed presently. The HPPP is constituted of a piston inside a cylinder.
Many measurement points are located within the cylinder with an appropriate probe for each one to record the position of the piston and a sensor for temperature and pressure. The primary test facility HPPP calibrates turbine meters as transfer standards. Two of these turbine meters are used in order to minimise measurement errors on the test facility. Then these turbine meters can be used to calibrate flow meters on the test facility. Some of the measurements presented in this paper were obtained, with natural gas, on this test facility for higher pressure, others were obtained on another PTB test facility, with air as a working fluid and at lower pressures of up to 16bar.

CESAME-EXADEBIT
The CESAME-EXADEBIT test facility principle is described in references [9] [15] and pictured on the left of Figure 3. The main characteristics are summarized in Table 4. For the calibration of sonic nozzles, the technique is based on the pVT,t method under dry air, whereby a perfectly known volume capacity is filled by registering the total pressure and total temperature of the flow upstream of the test nozzle and the downstream pressure inside the capacity. With the knowledge of the filling time, the mass flow rate of the nozzle can be easily calculated. The working pressure range goes up to 60 bar. The test facility operates under water in order to accelerate the stabilisation process and limit vibrations. The uncertainty for the calibration value of a sonic nozzle at CESAME-EXADEBIT can be claimed as less than 0.11% on A•CD value for pressure up to 60 bar (k=2).For each facility, the nozzles in this study were tested within the mass flow ranges presented in Figure 4. The range of the mass flow rate corresponds to the range of Reynolds numbers which are related to pressure levels. In Figure 4, the blue curves represent the range covered by CESAME-EXADEBIT and the red show the range covered by PTB and PIGSAR TM . In blue, the tested range of experiments was larger due to the higher availability of both the primary (indicated by the blank squares in Figure 4) and secondary test facility (coloured squares).

Uncertainty budget
For CESAME-EXADEBIT's primary test facility, the uncertainty calculation comes from the requirement of the Guide for the expression of uncertainty in measurement (GUM) [JCGM 100:2008] and the associated law of propagation of uncertainties. It is calculated with the composed uncertainties detailed in equation (7) from the formulation of equation (6) about the measurand ACD.

CESAME-
EXADEBIT's primary test facility (pVT,t method) on the left and secondary test facility (right).
The different terms in equation (7) are the variables that were measured with the uncertainty of the associated measurement material. Even if the volume faces variations between the initial state and the final state, these are integrated into the equation as a correction. In equations (6) and (7), "i", "f" and "0" are respectively the initial, final and total conditions. The standard deviation σ is also in the composed uncertainty calculation and depends on the repeatability range. This range cannot be less than three repeatable measurement points.
When the uncertainty of the measurement is increasing, the number of measurement points is increased in order for repeatability to be satisfied. Using the A•CD, instead of the more commonly-used CD, avoids the need to know the exact nozzle diameter precisely. Between the calibration of the nozzle and its utilization in industrial processes, the diameter remains unchanged under the same temperature conditions. In the case where the discharge coefficient is needed, for example, in international material comparisons, it can be easily extracted using knowledge of the precise diameter of the cross-section area.

Experimental results
Only the curves of three nozzles of equivalent throat diameters equal to 5 mm size but with different surface roughness are given in Figure 5. Nozzle 1 had the smoothest roughness under the parameter Ra; the second, Nozzle 2, a medium roughness definition, while the surface of Nozzle 3 was rougher. The different surface qualities of the aforementioned nozzles are indicated by different symbols on the curves: circles, triangles and diamonds for the smoother, the intermediate and the rougher surface respectively. The red curves were obtained by PTB and the blue curves were obtained with the pVT,t method at CESAME-EXADEBIT. Figure 5 also presents the theoretical curves obtained by Vincent [3] and Mickan [6] from the 2018 experiments, and the ISO 9300 curve.
The recent points obtained by PTB at low Reynolds numbers exactly fit the laminar curves established in 2018 [6]. Points obtained with the pVT,t method cannot be certified under 6bar for this nozzle diameter size, and a miscorrespondence of the curves is possible at low Reynolds numbers. Supplementary experiments were performed on the secondary test facility of CESAME-EXADEBIT (M1 facility) with Nozzle 2 at low pressures from P0=0.9 to 4.5bar. The results are in very good agreement with the curves established in 2018 [6] and are inside the range of uncertainty.
As can be seen in Figure 5, the tested nozzles numbered 1 to 3 lead to, for higher Reynolds numbers, curves located between around 0.4 to 0.6% lower than the reference curve of the ISO 9300 standard for cylindrical shape. This disparity is not modified significantly by carrying out supplementary measurements of the diameter throat. This disparity comes mostly from the roughness of the internal nozzle surface but also possibly from flaws in the shape or finally the divergent length. In Figure 5, an inspection of the calibration curves reveals an increase of the CD value for Reynolds numbers between 5×10 4 to 2×10 5 . This is consistent with a mainly laminar development of the boundary layer according to references [10]. This increase is followed by a drop in the CD value between Reynolds numbers from 3×10 5 to 6×10 5 for the three tested nozzles, then a stabilization of the CD value within the higher Reynolds number range. This last part is more likely to normally correspond to a turbulent development of the boundary layer according to references [10]. This drop in CD values could be defined by the so-called transition range of the boundary layer from the laminar development to the turbulent development. The maximum CD value is for Reynolds numbers between 3×10 5 and 4×10 5 . Then values are approximately constant beyond Reynolds numbers of 1×10 6 .
Overall, comparing the results from PTB for higher Reynolds numbers shows that the differences correspond to 0.2% on the CD values from the pVT,t method. This is superior to the recommended 0.1% difference for comparing the two facilities. The measurement can still be considered comparable for this study. For this nozzle size, PTB has always obtained higher values than CESAME-EXADEBIT as can be assessed by the difference obtained by Mickan in 2016 and 2018 [6] for air and natural gas (NG). This gap can be attributed to the difference of the technologies involved and in the nature of the facilities, as well as the uncertainty at high pressure and the C* coefficient equation type. For these pressure levels, PTB works with the secondary test facility and natural gas (NG) and the primary method pVT,t uses air. Finally, it can be seen that the magnitude of the CD drop increases with the increase of roughness criteria, with 0.36% for Nozzle 1 to 0.79% for Nozzle 3. Those three curves have the maximum and minimum CD for transition in the same range of Reynolds numbers between 2×10 5 and 3×10 5 for the maximum and at 6×10 5 for the end of the transition. For this chart, an equation for the laminar part can be calculated and the coefficient can be compared in Table 1 with a=0.998, b=5.96 and n=0.5. These coefficients are highly comparable to those in the laminar part of the curve, especially to Vincent [3] and Mazen [4]. The behaviour of these curves for higher Reynolds numbers could be mainly attributed to the change in the transition point, i.e. the point at which the boundary layer transitions from laminar to turbulent. When roughness increases, the transition point occurs slightly earlier than in smoother nozzles. As a consequence, the boundary layer thickens in the throat, creating a mass flow deficit and, subsequently, a smaller CD coefficient.

Numerical results
For a better understanding of the phenomena occurring inside the CFVN, this project includes numerical simulations of the flow field. As a prerequisite for simulating complex wall surface effects on CD evolution in CFVNs, an efficient numerical strategy is defined. The present study is the first report of some tests carried out to check the sensitivity of CD evaluation to parameters such as mesh refinement and time convergence. This section will cover the definition of the numerical strategy, description of the flow and comparison with hydraulically smooth results.

Numerical strategy
Compressible Navier-Stokes equations in axisymmetric formulation are solved with a finite volume approach. The numerical code used in this project is based on the C++ toolbox OpenFOAM. The transient compressible solver density rhoCentralFoam has been chosen for its ability to capture shocks, in particular through the use of centralupwind schemes of Kurganov and Tadmor [11] and [12]. An implicit second-order backward time scheme is used for time discretization, but only a limited CFL constraint can be imposed due to the severity of flow conditions. Simulations required around 0.02ms (6×10 7 integration steps with time steps of order 10 -10 s) to reach a relevant convergence of the solution (CD converged at less than 0.02%). The domain is relatively long so, in order to have the whole flow, relatively long simulation times are required compared to the characteristic time of the flow. The Gauss linear discretisation scheme is used for the gradient schemes, which requires the interpolation of values from cell centres to face centres. A second-order flux reconstruction is done with a Van Leer limiter. Viscous terms are integrated by Gauss integration with linear interpolation. The computational domain is built according to the norm [1] requirements for the cylindrical nozzle shape as shown in Figure 6. The domain topology of the mesh was organized in a different number of sub-domain areas to enable modifications of the near-wall treatment and the core flow grid distribution separately.
For the post-process applications, it has been proven that boundary layer thickness estimation along the wall is made easier by maintaining the thickness of the near-wall region with quasi-orthogonal mesh lines ( Figure 6). The grid was successively refined until a converged estimation of the solution was reached; then, this solution was compared with the experimental results.
The main parameter for this estimation was the CD (at ReD = 1505000 with a d=10mm nozzle throat) coefficient as can be compared in Table 5. The final version was built with 33 860 cells. The addition of more cells only modified the CD estimation by less than 0.01%. The diffuser half angle was 3.5°, the divergent length was 7•d and the diameter at the throat is 10 mm. The domain is decomposed and parallelized onto 60 processors to decrease the time calculation and accelerate the convergence.
The working fluid was assumed to be air (γ=1.4) entering in the nozzle through the left boundary, with the total temperature given at T0 = 300 K at different inlet pressures, P0, from 3 bar up to 60 bar, which corresponds to the ReD = 4.5×10 5 up to 8.9×10 6 for this diameter throat. The perfect gas assumption was used for the simulations. The Prandtl number was 0.72 and the viscosity was evaluated with Sutherland law.
The code was initialized by imposing the upstream total pressure of the divergent nozzle region and the atmospheric condition inside the rest of the nozzle. This technique provides an acceleration of the convergence by comparison with a simple initialization of flow at rest and upstream conditions imposed at the inlet section only. The calculation was run within the tested pressure levels at the inflow for a constant outflow pressure level prescribed at the outflow boundary.
The computational grid was refined in the normal direction of the wall and close to the nozzle surface. The mesh was generated with an elliptic smoothing method such that quasi-orthogonal mesh lines are generated near the wall and the cell-to-cell stretching factor is limited to 1.2% towards the solid wall. The resolution was sufficient enough to provide 30 cells within the boundary layer. The position of the first point close to the wall was set to respect the criterion of Δy + < 1.
Simulations have been carried out in both laminar and turbulent regimes with the Spalart Allmaras model. This model employs one-equation, based on a modified turbulence viscosity, that solves a modelled transport equation for the kinematic eddy turbulent viscosity by using a viscosity like variable [13]. Fig. 6.

Structure of the flow
Visualization of the flow inside the convergent part and cylindrical throat at x/d=1 (streamwise location) with the mesh refinement finer close to the wall. The thick red lines are to extract data along the wall for the boundary layer evaluation. Figure 7 shows an example of the visualization of the Mach number distribution obtained inside the nozzle. Inside the convergent part, the conditions are subsonic and the flow accelerates as the cross-section reduces. At the throat, inside the cylindrical section, the Mach number is close to unity. The flow conditions are critical, as is characteristic of the behaviour of the CFVNs. In the divergent part, the flow is supersonic and, with the increase of the crosssection, the flow accelerates and this leads to complex flow structures associated with nozzle jet separation, including the separation zone and the separation shock structures As illustrated in Figure 8, the change of wall curvature at the beginning of the cylindrical part (from a constant curvature radius to an infinite one) can lead to the formation of a particular flow structure for a low to moderate Reynolds number range (including the transition region). The differential acceleration of the flow leads to a strong deformation of the sonic line and the apparition of a supersonic bulb close to the wall. This supersonic region extends around the beginning of the cylindrical part of the nozzle throat. A recompression occurs at the end of this supersonic region. It is associated with a local sudden thickening of the boundary layer and a subsequent formation of compression wave patterns propagating towards the axis where the flow has reached a supersonic state. Several successive reflections of these compression wave patterns are then observed between the nozzle axis and the cylindrical wall.
Such a flow structure yields strong similarities with shock-train like structures that could be observed in supersonic diffusers but has never been reported for the present nozzle geometry. It is worth noting that a sufficient mesh refinement (not necessarily used in previous numerical studies) is necessary to observe the formation of this  particular structure which would otherwise appear numerically diffused leading to a full subsonic flow up to the end of the cylindrical throat region. The sensitivity of the computed solution to the turbulence closure has been checked by testing k-omega SST, k-epsilon and Spalart Allmaras models. These different models lead to rather similar boundary layer growth and similar flow topology in the inviscid flow region. However, this particular flow topology is not observed for higher Reynolds numbers (i.e. for ReD>5.0×10 6 according to the present tests), as illustrated in Figure 8 (bottom). Further investigation would be required to sustain these numerical observations and identify the geometrical conditions possibly leading to the formation of this shock-train-like pattern. However, it is important to note that its presence is likely to modulate the growth of the boundary layer within the cylindrical throat and significantly impact the evolution of CD.
This weak shock-train may appear even if the intensity of the shocks and triggers is low but sufficient to modulate the boundary conditions seen by the boundary layer, which recreates an acceleration and re-compression effect that modulates the thicknesses of the boundary layer, enough to have a higher CD in the presence of the shocktrain. Much like in a shock network in a supersonic jet, the thickening of the mixture layer is delayed compared to a case where there would be no shock-train. Without a shock-train, the thickness of the boundary layer is only a function of the length of the cylinder, whereas with successive reaccelerations there are changes in the thickness in the cylindrical part. Figure 9 reports the radial distributions of ρ×U. As the density and velocity change strongly, they are made dimensionless by the reference values in the inviscid flow. This is a ratio between the mass flow per unit area and the reference mass flow. The positioning on the abscissa allows a better identification of regions where mass flow deficit can be observed, whatever its relative importance. The CD is higher than unity for an excess of mass flow with respect to the reference (1D) theoretical value. There is a small mass flow deficit in the central part but it does not change significantly with the upstream flow conditions. However, a significant change of the distribution is observed in the remaining part of the flow. It is worth noting that we observe not only an expected mass flux deficit within the boundary layer region, but also an excess of mass flux associated with the inviscid region beyond the boundary layer edge, which is associated with the strong initial streamline curvature. Fig. 9. ρ×U profiles at the entrance of the throat (at 0.01mm or x/d=1) above and at the exit (at 0.02mm or x/d=2) below for four pressure levels P0.

Parameters from flow profile inside cylindrical throat
At the end of the cylindrical portion, when the Reynolds number is increased, we observe both an increase of the near-wall region where the boundary layer thickening causes a mass flux deficit, but also a reduction of the excess of mass flux observed beyond the edge of boundary layer. The transition is a very important phenomenon because it drives the spreading rate of the boundary layer, but the mass deficit in the boundary layer does not appear to be sufficient in itself to explain the indirect effect produced in the rest of the flow. Even if a 1% gain seems large, it only happened in a portion of the flow (because it is not an integral value).

Numerical estimation of CD
The displacement thickness is evaluated by applying a vorticity-based sensor of the boundary layer edge and integration of wall-normal profiles of mass flow deficit. Its evolution along the nozzle is shown in Figure 10. The most salient feature is the sudden jump of displacement thickness associated with the end of the convergent region where a discontinuity of wall curvature is imposed. This jump is associated with the recompression that occurs at the beginning of the cylindrical nozzle throat region.
The simulation is assumed to be hydraulically smooth, provided that no transitional model has yet been implemented inside the code formulation. Only smooth surface nozzles are then compared in Figure 11. The squares around the numerical points of measurement correspond to error boxes for the results of the numerical values. Generally, the numerical error is less than 0.1% of the calculated numerical CD value.
Then it is possible to observe a slight decrease of the CD corresponding to the experiments and an increase for higher Reynolds numbers. Everything seems more sensitive in the throat area than in a toroidal nozzle because of the presence of the cylindrical area. It is worth noting that the numerical and experimental uncertainty range leads to the conclusion that the evolution of the predicted CD value remains globally similar to the experimental evaluation. The present results suggest that the evolution of CD may be driven not only by the relative thickening of the boundary layer but also by a non-negligible modification of the rest of the inviscid flow structure (variable jump of velocity and density).

Conclusion
Surface roughness effects on the CD of CFVNs has been characterized experimentally in a wide range of Reynolds numbers. It has been shown that the CD evolution remains unaffected in a low Reynolds number range but drops significantly for high Reynolds numbers beyond 1.0×10 6 for which it yields a plateau with a lower level of CD for higher roughness levels.
A complementary numerical study for a hydrodynamically smooth surface has been carried out to study the flow structure in the throat region. It is suggested that the highly curved convergent section followed by the cylindrical zone at a constant radius necessarily leads to the formation of a supersonic pocket followed by a small recompression wave close to the wall. However, the evolution of the flow structure in the rest of the internal part of the cylindrical region remains unclear. A rather classical topology is observed for very high Reynolds numbers with a smooth and monotonous flow acceleration up to the end of the cylindrical part of the geometry. However, for lower Reynolds numbers, the formation of a successive reflection of compression waves (pseudo-shock train) within the cylindrical wall region has been observed, which strongly modifies the inviscid flow structure which reciprocally modulates the boundary layer growth significantly within the cylindrical nozzle region. Whatever the exact inviscid flow topology, it is suggested that the global CD evolution, observed as a function of upstream total pressure, could result not only from the change of mass flow deficit within the boundary layer region but also from this variable alteration of the inviscid flow structure.
Further investigation is required to identify more precisely the geometrical and physical conditions which may lead to the original flow topology observed except for high Reynolds numbers. These will possibly confirm its physical relevance and help in determining how wall roughness is likely to modify the whole flow structure.