A primary standard for the volume flow rate of natural gas under high pressure based on laser Doppler velocimetry

The optical primary standard is connected to the German gas transmission network and it uses high-pressure natural gas in a pressure range between 1.6 MPa and 5.5 MPa and flow rates up to 1600 m³ h⁻¹. The facility entrance includes a filter to skim particles like flash rust out of the gas flow and to preserve the installations from deposits. The pressure and temperature regulation is provided by the pigsar™ system to ensure appropriate conditions for all flow rates taking gas from the upstream pipe network 'Gescher' and finally delivering into the local pipe network 'Glückauf' (figure 2). A test section with a high efficiency flow straightener can be used to prove equipment under test; length compensators (L) allow its unproblematic installation. The test section is followed by two static mixers from Sulzer, Winterthur in Switzerland [2]. Between them is a seeding device that injects small particles in the range of 1 µm used as scatterers for the optical measurement. The mixers provide a good homogeneity of the particle distribution in the flow. Behind the seeding and mixing area a Zanker-type flow conditioner diminishes the angular momentum of the flow and a diffuser enlarges the pipe diameter from DN 200 to DN 500. The connected settling pipe is a 10 m long pipe of DN 500 with the largest cross section and an inner diameter of 423 mm. The mean flow velocity in the pipe is typically 0.5 m s⁻¹ and the turbulence intensity is expected to be as high as 30% at the first section of the pipe, but could not be measured. The pipe with a length of 24D is to reduce the pre-existing turbulence from the flow conditioner and forms a pipe flow halfway to equilibrium conditions at the entrance to the optical test section. Here, the flow exhibits the wanted rotational symmetric shape of the flow profile. At the entrance, woven wire screens mainly reduce streamwise velocity fluctuations with little effect on the change in flow direction, and reduce the turbulent boundary layer. Vortices with the size of the wires behind the screen superpose with their adjacent vortices, dissolve and homogenize the momentum. This suppresses the formation of vortices larger than the mesh size, see [3]. The boundary layer also is affected by this homogenization. The vortices from the mesh transport momentum from the core flow into the boundary layer, making it thin. A normal, new turbulent boundary layer with vortices of the order of the boundary layer thickness will form after a short settling length, see [4]. A short stilling section after the mesh allows the decaying of the small-scale turbulence of the flow that is induced by the mesh.

Zoom In Zoom Out Reset image size

The LDV nozzle module with optical access to the test section and a connected nozzle bank with critical nozzles form the end of the test facility. The nozzle bank with eight nozzles sets the flow rate and serves as a transfer standard to the volumetric flow standard High-Pressure Piston Prover (HPPP, see figure 23) and also stabilizes the flow during the optical measurement.

Optical methods for measuring the flow velocity need a sufficiently large number of particles with a known size distribution. Because of the high pressure in the pipe flow, a Laskin nozzle aerosol generator for high-pressure seeding has been developed. The use of a slit-type nozzle increases the dispersed phase production and reduces the volume of compressed air required, see [5]. The nozzle with its jet impactor plate operates below the liquid surface in the pressurized vessel. The nozzle is constructed for a pressure difference of up to 2.0 MPa between the nozzle and the seeding fluid. The vessel has two glass sight funnels to check the liquid level and to see if the system is operating properly. In order to reduce—in the small internal space of the vessel—the effect of splashing and to separate large particles, a perforated baffle was placed in front of the output [5].

For the production of particles, an overpressure of the carrier gas to spray the seeding fluid is necessary. The real gas from the facility is taken as the carrier gas. Laskin-type aerosol generators show a linear production rate with pressure while maintaining a homogeneous particle size distribution. For a slit-type nozzle with a 100 µm slit and a 1 l quantity of seeding fluid in the vessel, operation has shown that about 0.1 MPa gauge pressure is sufficient for the particle rate needed, which is equivalent to about 2 ml h⁻¹ of fluid.

The working pressure at the facility is operated from 1.6 MPa to 5.5 MPa. To obtain a gauge pressure of 0.1 MPa for the aerosol generator, a pressure difference regulation valve was inserted in the pressure supply. In case of operation needed at very high pressures, a reciprocating gas compressor driven by the high-pressure gas supply can be switched on. The aerosol generator is designed for operating at up to 10 MPa; the vessel conforms with the European pressure equipment directive (PED) and is approved for application in real gas facilities. The nature of the seeding fluid used is di-ethyl-hexyl-sebacate (DEHS) with a density near that of water (figure 3).

Zoom In Zoom Out Reset image size

The aerosol generator is supplied with a differential pressure regulator valve that maintains the set differential pressure between the inlet pressure and the outlet pressure of the aerosol generator when the working pressure of the pipe facility is changed and also functions independently of the input pressure of the high pressure supply. The additional gates in figure 4 are used for refilling the seeding fluid as well as venting (safety shut-off valve SAV) and pressurizing the generator and for safety purposes (return valve RV).

Zoom In Zoom Out Reset image size

A nozzle contraction ratio c = 18 was designed to cover the flow rate range of all G1000 working standards used by pigsar™ (100 m³ h⁻¹ to 1600 m³ h⁻¹, see figure 23) in such a way that the maximum velocity at the nozzle outlet does not exceed approximately 50 m s⁻¹ for the selected flow rates. This ensures optimal conditions for LDV signal acquisition. The nozzle was designed as a compatible insert of the LDV nozzle module. The module is fitted with two high-pressure glass windows to provide the optical access required for LDV velocity measurements at the nozzle outlet surface. The flow temperature is monitored with a PT100 inside the wall of the nozzle's throat. Additionally, the temperature balance at all points of temperature measurement upstream and downstream of the optical-based system (see figure 2) is kept within 0.2 K to ensure an uncertainty of temperature measurement below 0.1 K. The static pressure is read at the inlet of the nozzle and through four static pressure holes of 1 mm in diameter at the nozzle exit.

A diffuser is installed downstream of the nozzle (to the right in figure 5) to ensure further downstream flow stability. With optical access, it is possible not only to operate a standard LDV to measure the centre line velocity, but also to use the optical boundary layer sensor to investigate small profiles.

Zoom In Zoom Out Reset image size

A bank of eight critical nozzles was installed in the pipework covering a flow range from 25 m³ h⁻¹ to 1600 m³ h⁻¹. The position of the nozzles was selected so that they can be operated in series either with the piston prover, the new optical standard or the two new meter runs. AGA-8 calculation algorithms are implemented in the data acquisition system in order to calculate the critical flow factor c^*, the thermodynamic property essential to critical nozzle mass flow measurements. Temperature and pressure are measured at each nozzle, and gas composition through an online process gas chromatograph. The function of the critical nozzles is manifold: firstly, they simply stabilize flow, which is important for LDV measurements. Secondly, the critical nozzles represent secondary standards of pigsar™. They make it possible to compare the two primary standards and check the reproducibility of the working and primary standards.

The measurement of flow rate in a turbulent pipe flow by integrating the flow velocity over the pipe area with the requirement of a low uncertainty is not compatible with short measuring times because of the level of turbulence. Hence a nozzle is a way of shortening measuring times. A nozzle contracts the flow by converting the pressure into kinetic energy. The right shape of its contour contributes to the same amount of pressure conversion for all streamlines and if the energy increase is homogeneous in the main flow, the velocity differences in the nozzle exit will be small. A higher contraction causes a more planar velocity profile.

The computation of a large number of nozzle shapes undertaken by Börger [4] resulted in a set of profile parameters for different Reynolds numbers and contraction ratios that have been adapted to the LDV nozzle.

For the optimization of the nozzle contour a potential flow inside the nozzle has been assumed. The calculus uses a singularity model. The walls have been replaced by vortex layers; at the inlet is a source disc and at the outlet of the nozzle is a sink disc. The wall friction was taken into account by boundary layer calculations using the method of Bradshaw et al [6]. A profile is regarded as optimized when—in the exit of the nozzle—there is a plane velocity distribution within 0.1% and when—in the boundary layers of the contour—the friction coefficients are larger than 0.002 to suppress flow separation [6].

A nozzle with a diameter of 100 mm (see table 1) has been manufactured and the nozzle exit plane roundness has been calibrated by PTB Department 5.3.

The optical flow rate measurement is traced back to a precise velocity measurement at the LDV nozzle outlet of the LDV nozzle module (figure 5). The velocity u of the gas flow is measured by a dual scatter LDV based on the differential Doppler technique and is determined by

$$\begin{equation} u=d\cdot f_{{\rm D}} \end{equation} \tag{ 2 }$$

with f_D as the Doppler frequency of the measuring signal generated by the scattered light of tracer particles embedded in the flow and d as the fringe spacing in the LDV measuring volume.

The measuring volume is given by the intersection volume of two laser beams crossing each other at an angle of ε. The fringe spacing is ideally given by

$$\begin{equation} d=\frac{\lambda}{2\sin \varepsilon_{\frac{1}{2}}}. \end{equation} \tag{ 3 }$$

For a well-aligned optical set-up and Gaussian beams crossing each other at their beam waists, the fringe spacing is constant within the measuring volume. This is not the case when thick windows for optical access have to be used. Thus, the LDV has to be calibrated including the high-pressure glass window. The calibration is done by measuring the local fringe spacing and determining its deviations over the length of the measuring volume. The deviation of the fringe spacing within the measuring volume determines the uncertainty of the LDV velocity measurement and also gives the impression of a non-existing turbulence intensity of the flow. The LDV with a wavelength of 852 nm and a power of 26 mW was used in a forward scattering mode. The focal length of 600 mm generates a fringe distance of 10.2 µm, a diameter of 163 µm and a length of 4 mm of the measurement volume (figure 7). This allows the LDV to measure in the centre of the flow as near as 7 mm from the nozzle exit plane.

The calibration method is to generate precise particle velocities and to measure the Doppler frequency in order to calculate the fringe spacing at different positions of particle transitions in the measuring volume. The most precise method is to use single particles on a rotating glass wheel with its well-known diameter. The wheel is driven at a constant rotation speed by a stepping motor to ensure a very smooth and accurate circumferential velocity. After calibration of the LDV, the LDV and its power supply are installed in a container that protects them from explosion risk in the zoned area 2. The container with the LDV inside is installed at the window of the nozzle module (figure 6) and is flushed with nitrogen before operation. The velocity signals of the LDV are sent via a 30 m long fibre optical cable to the signal acquisition system.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

This LDV system has been used in backscatter mode in order to monitor the centre line velocity of the top-hat profile of the nozzle and is designed to have access to the centre of the nozzle flow only.

The boundary layer LDV has to measure the shape of the boundary layer at the nozzle exit. The spatial resolution of conventional laser Doppler systems is limited by the dimensions of the measurement volume. In applications with high velocity gradients such as in boundary layers, the integration of the local velocity over the extension of the measurement volume cannot be neglected. The method used locates the particle trajectory inside the measurement volume [7].

Instead of analysing the scattered light from two beams (dual backscatter mode), only the scattered light of the illuminating beam is analysed with light from a reference beam (reference beam technique). The reference beam is sent parallel and 2 mm downstream of the nozzle exit plane, whereas the illuminating beam intersects the reference beam near the wall of the nozzle exit. Inside the reference beam leaving the nozzle area, a fibre bundle is placed to receive the light of the reference beam as well as scattered light from the illuminating beam (figures 8 and 9).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Each fibre of the bundle 'sees' a slightly tilted virtual fringe system and fringe distance d_i in the illuminating beam, and the phase differences Δφ between pairs of signals are used to calculate the position z of traversing particles in the nozzle exit plane. Three of the fibres have been utilized in order to increase the phase resolution by a factor of two and for a better signal validation and noise reduction [8].

With the distance r of the measuring volume to the receiver and the distance D between the receivers, the position z to the crossover point of the beams is

$$\begin{equation} z\approx \frac{\lambda r}{4\pi D\sin \varepsilon_{\frac{1}{2}}}\Delta \varphi =\frac{\bar{{d}}r}{2\pi D}. \end{equation} \tag{ 4 }$$

The fringe distance d'(z) of the multiple reference beam LDVs shows a small gradient that can be corrected by

$$\begin{equation} d'(z)\approx \frac{r}{r+z}\cdot d_{i} . \end{equation} \tag{ 5 }$$

Because the fringe distances d_i of the fibres are slightly different from each other, they have been synchronously calibrated with a calibrated standard LDV in the centre of the nozzle (figure 10).

Zoom In Zoom Out Reset image size

To ensure traceability to the units 'metre' and 'second' for the measurement of centre line velocity, the LDV was properly calibrated by measuring the local fringe spacing and determining its deviations inside the region of the measuring volume where the intensity of light is larger than 1/e² of the maximum value.

The calibration method used to measure the local fringe distance is to generate precise particle velocities and to measure the Doppler frequency in order to calculate the fringe spacing at different positions of particle transitions along the measuring volume. We used single particles on a rotating glass wheel with its well-known diameter driven by a stepping motor to ensure a very constant and accurate speed. For the calibration of the LDV the high-pressure glass window, which is used at the test section, was included. The equipment and a detailed uncertainty budget for the particle velocity on the glass wheel are given in [9].

Calibration of the LDV is done according to the relevant calibration and measurement capability (CMC) of PTB (see Key Comparison Database [KCDB], 9.7.1, DE39) with an expanded uncertainty of U = 0.11% including a term of reproducibility checked by four repeated calibrations in a time scale of several weeks.

The description of the local mean velocity in turbulent boundary layers as a function of the wall distance y is always related to logarithmic dependences (see for this especially the universal logarithmic law of the wall, e.g. [10]). The reduction of large sets of single data should therefore have a relation to logarithmic or exponential functions but should additionally have good flexibility to fit to experimental data and should have a reasonably small number of parameters. The choice in our case was a modified hyperbolic tangent function as given in

$$\begin{equation} u(y)=U_{0} \cdot \tanh \left[ {k\cdot \left( {y-y_{{\rm w}}} \right)^{b}} \right]. \end{equation} \tag{ 8 }$$

The parameters U₀ (velocity at the outer edge of the boundary layer) and y_w (position of the wall in the coordinate system of the optical boundary layer LDV) are general values for each velocity profile measured. The parameters k and b in equation (8) are determined by means of a non-linear least-squares fit to adapt the specific shape of the velocity profile⁵. Figure 15 shows one very typical example of such an approximation.

With the functional representation of the mean velocity according to equation (8) we obtained a good base for the determination of further important data as is illustrated in figure 15. The first additional data were the local fluctuations in the velocity component s_u (turbulent intensity in the flow direction) at a specific measurement location in the wall direction. The second were in a similar way the local fluctuations in the location y_meas for a specific level of mean velocity u/U₀ in the flow profile. This second value s_y does not yet have an established use in the literature. It is quite a new situation generated by the features of the boundary layer sensor and we have to emphasize that this additional information was the main basis for the successful treatment of all data.

The first important conclusion drawn from the data s_y was the determination of the limit location up to where the boundary layer of the nozzle flow was influenced by the free stream vortex. As one can see in figure 17, the dependence of the local fluctuation s_y shows a good self-similarity for all profiles (by means of normalization with the boundary layer thickness δ) and all profiles have a significant local minimum of s_y at the same (normalized) location. We assume a no-slip condition at the wall for the nozzle flow and therefore the turbulent fluctuation has to decrease continuously to zero towards the wall. Hence, this local minimum indicates the limit position of the free stream vortex influence, which causes this deviation from the expected behaviour of nozzle flow. With this information (available separately for each profile) we were able to determine the range of the mean velocity profile, which we can use for further calculations assuming a conventional wall bounded nozzle flow (blue dotted part of the mean velocity profile in figure 15). The upper limit of the mean velocity profile usable for the subsequent calculations was set at the location of u/U₀ = 0.995 because the hyperbolic tangent function does not provide a limited boundary layer u/U₀ = 1 at y = δ as is necessary for our calculation approach (see the following section).

The turbulent boundary layers have been the subject of a very large number of experimental investigations and detailed analyses in the past which have led to well-established half-empirical models⁶. We will focus here only on those approaches which were finally important for our solution because it is impossible to give a comprehensive overview of everything (for this please refer to the standard literature, e.g. [10]).

For all turbulent boundary layers with no-slip condition at the wall, general behaviour of self-similarity was found if the velocity profile is normalized with the appropriate characteristic values:

$$\begin{equation} \frac{u}{u_{\tau}}=u^{+}=G\left( {\frac{y\cdot u_{\tau}}{\nu};\frac{y}{\delta}} \right)=G\left( {y^{+};\frac{y}{\delta}} \right), \end{equation} \tag{ 9 }$$

where u_τ is the so-called wall shear rate (or wall shear velocity), δ is the boundary layer thickness (location y = δ where u/U₀ = 1) and ν is the molecular kinematic viscosity of the fluid).

Equation (9) splits the boundary layer into an inner part G_i(y⁺) and an outer part G_o(y/δ). The inner part is described by the well-known logarithmic law of the wall and the outer part is described by the velocity defect law. The logarithmic law of the wall has a universal character for turbulent shear layers with no-slip condition but the velocity defect law is specific for a certain type of flow configuration. Nevertheless, for a given flow configuration also the outer part follows the rules of self-similarity.

The most common starting point for the description of mean velocity in the boundary layer is the basic relationship between the total shear stress τ_total, the mean velocity gradient and the turbulent shear stress 〈u'v'〉:

$$\begin{equation} \tau_{{\rm total}} =\tau_{{\rm lam}} +\tau_{{\rm turb}} =\rho \left[ {\upsilon \frac{{\rm d}u}{{\rm d}y}-\left\langle {{u}'{v}'} \right\rangle} \right]. \end{equation} \tag{ 10 }$$

The turbulent shear stress 〈u'v'〉 (correlation between fluctuation of velocity in the main flow direction and direction normal to the wall) is, in most cases, not measured and not known; therefore, we have here the key point for the different approaches to describe the boundary layer by modelling. The oldest and most simple (algebraic) models⁷ are the concepts of eddy viscosity ν_turb introduced by Boussinesq [11] and the mixing length l_m by Prandtl [12] with which equation (10) is reformulated as

$$\begin{eqnarray} &&\fl\tau_{{\rm total}} =\tau_{{\rm lam}} +\tau_{{\rm turb}} =\rho \left[ {\upsilon \frac{{\rm d}u}{{\rm d}y}+\upsilon_{{\rm turb}} \frac{{\rm d}u}{{\rm d}y}} \right]\nonumber\\ &&=\rho \left[ {\nu \frac{{\rm d}u}{{\rm d}y}+l_{{\rm m}}^{2} \left| {\frac{{\rm d}u}{{\rm d}y}} \right|\frac{{\rm d}u}{{\rm d}y}} \right]. \end{eqnarray} \tag{ 11 }$$

A normalization using the molecular kinematic viscosity ν, the shear stress rate u_τ and the wall shear stress $\tau_{{\rm w}} =\rho u_{\tau}^{2}$ yields

$$\begin{eqnarray} &&\fl\frac{\tau_{{\rm total}}}{\tau_{{\rm w}}}=\tau^{+}=\left( {1+\upsilon _{{\rm turb}}^{+}} \right)\frac{{\rm d}u^{+}}{{\rm d}y^{+}}=\left( {1+l_{{\rm m}}^{+^{2}} \left| {\frac{{\rm d}u^{+}}{{\rm d}y^{+}}} \right|} \right)\frac{{\rm d}u^{+}}{{\rm d}y^{+}}.\nonumber\\ \end{eqnarray} \tag{ 12 }$$

The concept of mixing length introduced by Prandtl [12] is perhaps the best known because, with the simple modelling of the mixing length using the wall distance and the empirical constant κ (introduced by Karman [13])

$$\begin{equation} l_{{\rm m}} =\kappa y, \end{equation} \tag{ 13 }$$

there is a short and direct way to the basic formulation of the logarithmic law of the wall:

$$\begin{equation} u^{+}=\frac{1}{\kappa}\ln \left( {y^{+}} \right)+C. \end{equation} \tag{ 14 }$$

As mentioned above, the logarithmic law of the wall is valid in the inner part of the boundary layer but not in the outer part. Consequently, for the complete boundary layer a model approach for either the eddy viscosity ν_turb or the mixing length l_m as a function of the overall location of the layer is necessary. Established models for this are given by Cebeti and Smith [14] and by Baldwin and Lomax [15] for the eddy viscosity as well as by Michel [16], Escudier [17] and Szablewski [18] for the mixing length. All of them have in common that in the inner part (small values of wall distance y, smaller than a certain distance y_i/o), equations (13) and (14) are valid. The modelling of the eddy viscosity by [14] as well as [15] in the outer part of the layer is characterized by a significant decrease in the viscosity for larger distances to the wall. The models for the mixing length operate with a constant mixing length l_m,o in the outer part which is in the order 0.04δ to 0.09δ.

As mentioned above, the self-similar outer part of the boundary layer is specific for a certain type of flow configuration. Therefore, these established model approaches to cover the inner as well as the outer part do not give necessarily reliable solutions for all cases. And in fact, the trial to fit our boundary layer profiles with the models using a constant mixing length for the outer part did not lead to a sufficient solution.

It is, nevertheless, necessary to make the choice of an algebraic model and to determine the triple value [u_τ, δ, y_w]⁸ for each of the profiles being measured so that the experimental data of mean velocity are well represented by the model function for local velocity u versus wall distance y.

Any approach which we can choose here will make use, more or less, of some information based on the established half-empirical solutions out of the literature. This generates an additional challenge to our central aim of realizing a real independent optical primary flow rate standard. All information which we introduce or which we make use of has to be evaluated finally with its uncertainty and its impact on the uncertainty of our final measurand. This is not possible for most of the information from the literature. For example, one can look at the van Driest constant A⁺ and the Karman constant κ (see equations (15) and (16)). These constants are determined to fit the related models to experimental data. Although it seems that these constants have a universal character, there is an impressive discussion about their values which can vary more or less in the publications depending on the specific conditions of the experiment or the measurement. To get an impression one can refer to [19] (κ and A⁺) as well as [20] (κ)⁹.

To overcome the latter problem, two approaches will be applied and will be fitted to our experimental data; both approaches will use different sources of information from our experimental data as well as out of the established descriptions of boundary layers. The final evaluation of uncertainty for the discharge coefficient c_D then involves the differences resulting from the two approaches.

4.2.1. Approach 1: self-similarity in the region of defect law.

The first approach sets the self-similarity of the velocity profiles in the defect law region into the focus. We determined the related triples [u_τ, δ, y_w]_i for all individual profiles¹⁰ so that they collapsed early into one normalized profile from a certain common location y ⩾ y_i/o. This is shown as a blue curve in the defect law model in figure 16, the normalized self-similar velocity profile needed to be adapted to the wall (y = 0) and to the core flow.

Zoom In Zoom Out Reset image size

The adaption to the wall has been based on the basic concept of the mixing length in the region of the logarithmic law of the wall with a transition to the region of the defect law. We give here the formulation of Szablewski [18] for the mixing length (equations (15) and (16))^11,12:

$$\begin{equation} \frac{l_{{\rm m,Swabl}}}{\delta}=\kappa \frac{y}{\delta}\cdot \exp \left( {-\frac{y}{y_{{\rm i/o}}}} \right)\tqs {\rm for}\,\,y\le y_{{\rm i/o}} , \end{equation} \tag{ 15 }$$

$$\begin{equation} \frac{l_{{\rm m}}}{\delta}=\frac{l_{{\rm m,Swabl}}}{\delta}\left[ {1-\exp \left( {-\frac{y^{+}}{A^{+}}} \right)} \right]. \end{equation} \tag{ 16 }$$

According to equation (12) there is also the need to define the total shear stress τ_total as a function of wall distance in the region where we want to apply the mixing length model. It can be assumed that near the wall the total shear stress is a linear function with a gradient slightly below zero¹³, see also figure 20.

All together we obtained 65 parameters which had to be optimized by means of least-squares principles to fit our experimental profile data: 21 triples [u_τ, δ, y_w]_i, the common transition point y_i/o/δ between the inner and outer parts of the boundary layer and the gradient of the total shear stress near the wall. The boundary conditions for the optimization are, besides the best collapsing of experimental data to an average profile, the smooth connection of the near wall function (equations (15) and (16)) to the outer part at y_i/o/δ in the local velocity as well as the first derivative of velocity.

As mentioned in section 4.1, the representations of experimental data by the hyperbolic tangent function (equation (8)) were used up to velocities u ⩽ 0.995U₀. To close the last gap to the core flow we applied a simple third-order polynomial as given in equation (17):

$$\begin{eqnarray} &&\fl\frac{U_{0} -u}{u_{\tau}}=a_{0} +a_{1} \frac{y}{\delta}+a_{2} \left( {\frac{y}{\delta}} \right)^{2}+a_{3} \left( {\frac{y}{\delta}} \right)^{3}\quad {\rm for} \ \frac{y}{\delta}\ge 0.75\nonumber\\ \end{eqnarray} \tag{ 17 }$$

with values of a₀ to a₃ so that the adaption is smooth up to the first derivative at y/δ = 0.75 and at y = δ (with additional condition U₀ = u at y = δ).

The result of the whole process is shown in figure 16 using the velocity defect model. It documents the very good agreement of the experimental velocity data with the rules of self-similarity. The near wall region (which is not supported by measured data and which we modelled using the mixing length model) covers about 10% of the boundary layer. Due to the velocity gradient in this part, this implies that about 35% to 38% of the discharge coefficient is directly dependent on the quality of the near wall model. This illustrates the necessity to verify these results by another approach which is explained in the following section.

4.2.2. Approach 2: modelling of the kinematic eddy viscosity.

The second approach is based on modelling the eddy viscosity ν_turb in the boundary layer. The concept of eddy viscosity was introduced by Boussinesq [11] and derives from a similarity consideration of viscosity in kinetic gas theory where the molecular kinematic viscosity is derived from the relation of viscosity ν with the average molecule velocity $\bar{{c}}$ and the average free path length λ yielding $\nu =\textstyle{1 \over 3}\bar{{c}}\cdot \lambda$. Boussinesq set the turbulence intensity equivalent to the molecule velocity and the free path length to a certain length scale l_MFW which is a characterizing parameter of the turbulent flow: $\bar{{c}}\Rightarrow \sqrt {\left\langle {u'^{2}} \right\rangle}$; λ ⇒ l_MFW. The length scale l_MFW is in close relationship with the mixing length which was introduced later by Prandtl [12].

It is a new and unique situation that the application of the boundary layer sensor (section 2.6) allows us to derive two additional values out of the velocities usable to establish a quite similar approach: the standard deviation of velocity fluctuation s_u (which is identical to conventional turbulence $\sqrt {\langle {u'^{2}} \rangle} )$ and the standard deviation of location s_y (which has a similar meaning to s_u but for the length), see for this also figure 15. We will now define the turbulent kinematic (eddy) viscosity ν_turb proportional to the product of s_y and s_u for reasons of consistency of physical units:

$$\begin{equation} \nu_{{\rm turb}} =\alpha \cdot s_{y} \cdot s_{u} =l_{{\rm m}}^{2} \cdot \frac{\partial u}{\partial y}. \end{equation} \tag{ 18 }$$

The task is to determine for all profiles the related values of s_y and s_u as well as the proportionality factor α¹⁴. We were starting with the assumptions that s_y and s_u follow the basic rules of self-similarity as well as that the proportionality factor α is a characteristic constant valid in the whole range of boundary layers and for all profiles out of the same flow configuration.

The self-similarity of the normalized values of s_y and s_u is demonstrated in figures 17 and 18. It allows the representation of s_y as well as s_u by two general average functions versus the wall distance valid for all profiles independent of the Reynolds number. This self-similarity was consequently introduced as an additional condition for the least-squares approximation to define the best fitting parameter triples [u_τ, δ, y_w]_i.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The general average functions for s_y and s_u in figures 17 and 18 consequently establish the turbulent kinematic (eddy) viscosity ν_turb (figure 19) according to equation (18). The parameter α has to be determined so that the turbulent shear stress calculated using equation (12) fits the experimental velocity data in an optimal way. As one can see in figure 20, the total shear stress also follows the principles of self-similarity and can be represented by a general function for all profiles.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The finally determined turbulent kinematic viscosity from the experimental data in figure 19 can be compared with the values which are given by the model approach of Cebeti and Smith [14] and Cebeti and Bradshaw [21]. This approach is characterized by a viscosity ν_turb,∞ in the outer part (y > y_i/o) which is damped towards the core flow by a so-called intermittency function:

$$\begin{equation} \nu_{{\rm turb}} =\nu_{{\rm turb},\infty} \left[ {1+5.5\left( {\frac{y}{y_{19}}} \right)^{6}} \right]^{-1}. \end{equation} \tag{ 19 }$$

y₁₉ is the location where (U₀ − u)/u_τ = 0.19 applies according to Gersten [22].

The transition to the inner part of the boundary layer at the transition point y_i/o is given where the viscosity according to equation (19) is larger than the turbulent viscosity derived from the universal logarithmic law of the wall.

We get ν_turb,∞ = 0.017 and y_i/o = 0.03δ for our measurements which are in good agreement with the literature values ν_turb,∞ = 0.014 to 0.016 and y_i/o = 0.032δ to 0.034δ [22] and which therefore confirm the reliability of our approach.

In figure 21 the 21 velocity profiles are shown which were calculated using the general (average) functions for turbulent kinematic (eddy) viscosity (figure 19) and the total shear stress (figure 20) according to equation (12) together with the related values for the triples [u_τ, δ, y_w]_i and the parameter α = 0.07 (see equation (18)). The general behaviour and level of results are very close to the results obtained with approach 1 in section 4.2.1 (figure 16).

Zoom In Zoom Out Reset image size

After the introduction of the calculation of the mean velocity profiles out of the single velocity data measured by means of the LDV-based boundary layer sensor, we could get the discharge coefficient of the nozzle flow as a function of the Reynolds number directly by integrating the profiles given in the velocity defect law model (see figures 16 and 21). The results of this integration are given for all 21 profiles in figure 22 for both the first and second approaches of boundary layer calculation.

Zoom In Zoom Out Reset image size

The single results for the discharge coefficient can be summarized as a function depending on the Reynolds number because of the relation to the displacement thickness δ^* (equation (7)) for turbulent boundary layers which generally follows the functionality δ^* ∼ Re⁻ⁿ with values of n in the range between 0.132 and 0.25. We have in our case a turbulent boundary layer which is close to the situation of a layer at a flat plate¹⁵. Therefore, we adopt the value of n = 0.2 which is commonly found in the literature. The complete function describing the discharge coefficient in figure 22 is given as a set in equations (20) to (22) with related parameters (determined by non-linear least-squares fit) given in table 2.

$$\begin{equation} c_{D} =s_{1} \cdot c_{D,1} +s_{2} \cdot c_{D,2} , \end{equation} \tag{ 20 }$$

$$\begin{equation} c_{{\rm D},i} =1+\frac{b_{i}}{Re_{{\rm D}}^{0.2}}, \end{equation} \tag{ 21 }$$

$$\begin{equation} s_{i} =0.5\left\{{1+\left( {-1} \right)^{i}\tanh \left[ {k_{{\rm Tr}} \log \left( {\frac{Re_{{\rm D}}}{Re_{{\rm D,Tr}}}} \right)} \right]} \right\} \end{equation} \tag{ 22 }$$

with s₁ + s₂ = 1 for all Reynolds numbers Re_D.

It is necessary to explain why we have two different sections of c_D dependence according to equation (21), which are connected by a transition function (equation (22)). As mentioned above at the beginning of section 4, the Reynolds number $Re_{\delta^{\ast}}$—with respect to the displacement thickness—clearly indicates a fully turbulent flow regime in the boundary layer. Therefore, this cannot be the laminar–turbulent transition (this would not be compatible also because of n = 0.5 for laminar flow). The key point for the explanation is the relation to roughness. We have measured the roughness by means of a stylus instrument in some parts of the inner nozzle surface and got an average peak-to-valley height of Rz ≈ 9 µm which is in good agreement with the specification given for the manufacturing of the nozzle¹⁶. With the wall shear rate u_τ and the molecular kinematic viscosity ν we obtain the normalized roughness height $k_{{\rm s}}^{+} ={{\rm Rz}\cdot u_{\tau}}/\nu ={{\rm Rz}\cdot Re_{{\rm D}} \cdot u_{\tau}}/{2\cdot R\cdot U_{0}}$. In figure 22 it is to be seen that the transition from the first to the second part of the c_D function starts at about Re_D = 1.9×10⁶ (value for s₁ = 0.05 and s₂ = 0.95 respectively). With the related values for the wall shear rate of u_τ/U₀ = 0.0295 (which we got from the boundary layer calculation for the profiles at this Reynolds number) and the nozzle radius of R = 50 mm, a final value of $k_{{\rm s}}^{+}\approx {\rm 5}$ results. It is a general statement in fluid dynamics (see, e.g., [10]) that a surface can be assumed to be hydraulically smooth as long as the value of $k_{{\rm s}}^{+}$ is ⩽5. This is equivalent to the thickness of the so-called laminar (viscous) sub-layer of a turbulent shear layer at the wall. If the roughness is larger than this, then the viscous sub-layer vanishes and the wall starts to display the characteristic of a rough surface for the flow. With this we have a reliable explanation of our c_D behaviour versus the Reynolds number, which reflects the transition from a smooth to a rough surface at a high level of consistency with the literature values.

The definition of the generalized c_D function allows us to discuss the uncertainty of the c_D calculation out of our determination of boundary layer profiles according to sections 4.2.1 and 4.2.2. Here we have to consider two parts: firstly, the differences coming out of the two different approaches for boundary layer calculation, and secondly, the deviations of the c_D values for the single profiles from the generalized c_D function.

The differences in the results between the two approaches for one profile do not exceed 0.05%. Hereby, the first approach made use of the mean velocity information and a simple mixing length model for the wall region only. The second approach additionally includes information about the turbulence structure based on the measurements. The results for the turbulence structure (eddy viscosity ν_turb) had been found very close to well-established concepts in the literature and we therefore have a high level of confidence in our results. Hence, the difference which is present in the two approaches can be used here as a measure for the uncertainty caused by the model structure for the algebraic solution of the boundary layer calculation, which is consequently conservatively estimated as 0.05% of the c_D value¹⁷.

Additionally to the differences in the calculation approaches, we can derive from figure 22 that the quality of the profile measurement is different for the single profiles and generates a higher uncertainty in the order of 0.15%. This is visible in figure 22 considering the fact that the discharge coefficient should theoretically follow the straight function versus the Reynolds number as given in equations (20) to (22). For example, the c_D values for Reynolds numbers between Re_D = 1.9 × 10⁶ and Re_D = 8 × 10⁶ for profiles measured at pressures of p = 3.3 MPa and 5.0 MPa are very close together while the values for p = 1.9 MPa differ here. The main reasons are found in the quality of optical access, which changes with operation time as was deduced in section 3 for the determination of the centre line velocity.

In total, we have to consider an expanded uncertainty of 0.154% for the c_D value calculated with equations (20) to (22) to cover 95% of our experimental result. This is the largest contribution to the uncertainty budget as listed in table 3.