Using Modelica to investigate the dynamic behaviour of the German national standard for high pressure natural gas flow metering

This paper presents a computational model written in Modelica for the high pressure piston prover (HPPP) used as the national primary standard for high pressure natural gas flow metering in Germany. With a piston prover the gas flow rate is determined by measuring the time a piston needs to displace a certain volume of gas in a cylinder. Fluctuating piston velocity during measurement can be a significant source of uncertainty if not considered in an appropriate way. The model was built to investigate measures for the reduction of this uncertainty. Validation shows a good compliance of the piston velocity in the model with measured data for certain volume flow rates. Reduction of the piston weight, variation of the start valve switching time and integration of a flow straightener were found to reduce the piston velocity fluctuations in the model significantly. The fast and cost effective generation of those results shows the strength of the used modelling approach.


Introduction
Natural gas is one of the worlds most frequently used energy carriers and transported over large distances in high pres sure pipelines or as liquefied natural gas using ships. For the trade with natural gas the uncertainty of the gas flow meters is of major importance. The uncertainty depends on the calibration chain and increases with every calibration step. Therefore the German national metrological institute

Measurement Science and Technology
Using Modelica to investigate the dynamic behaviour of the German national standard for high pressure natural gas flow metering M von der Heyde 1 , G Schmitz 1 and B Mickan 2 PhysikalischTechnische Bundesanstalt uses a high pressure piston prover as a national primary standard for the calibra tion of high pres sure natural gas flow meters. Other primary standards for high pressure natural gas metering are for example PVTtstandards or gravimetric standards. PVTt standards measure the pres sure and temperature in a known volume to determine the mass that has flown though a transfer standard in the only inlet during a certain time. Gravimetric standards weigh the gas. International comparisons of the national standards ensure the consistency and accuracy of the measurements.

Modelica
Modelica is a nonproprietary, equationbased and object oriented modelling language [1]. It can be used to assess dynamic physical systems, described as a set of ordinary differ ential, algebraic and discrete equations. Modelica is therefore not limited to any physical or technical domain. The first specification was published in 1997. The nonprofit Modelica association is responsible for the development and publication of the Modelica language as well as the Modelica standard library. The Modelica standard library is a standard ized and free library, currently consisting of 1360 models and 1280 functions from a wide range of engineering domains. Modelica models are highly reusable because they are equationbased and objectoriented. This allows us to build on previous work, which is essential when developing com plex models. For the model presented in this paper it was pos sible to use several components from the Modelica Standard Library containing equations from various physical domains, such as tribology, thermodynamics and fluid flow. Parts of the model, such as the medium description, are replaceable due to the objectorientation of Modelica which allows for easy modification.
As Modelica is equationbased, physical equations can be inserted in the model directly minimizing errors during the modelling process. The conversion from the physical equa tions to the solving algorithm is automatically done by the Compiler.
Several simulation tools based on Modelica are available. One of them is Dymola ® [2], which was used to develop and evaluate the model described in this paper. It is possible in Dymola ® to build up models directly in the Modelica Language and connect those models or any other model graphically.

High pressure piston prover
The HPPP is operated and owned by the German national metrological institute PhysikalischTechnische Bundesanstalt (PTB) and currently installed on the calibration site for gas flow meters pigsar TM in Dorsten, Germany.
The HPPP can be operated with inlet pressures up to 90 bar and flow rates up to 480 m h 3 1 . The central element of the HPPP is a piston in a cylinder. Indicators are mounted on the cylinder to signal the piston position. Pressure and temperature are measured downstream of the cylinder. Some other components are needed for the calibration of transfer standards and have to be considered here as their dynamic behaviour has an effect on the piston motion, such as the transfer standards themselves, several valves and a nozzle bank.
The whole calibration setup is shown in figure 1. The closing of start valve 2 commences the runningin phase (start valve 1 is already open as it is only needed to prevent movement of the piston in between calibration runs due to the pressure drop across start valve 2). The motion of the piston is indicated by the piston position indicator a 1 . The measurement phase starts as the piston passes indicators a b c , , 2 2 2 and ends as the piston passes the indicators a b c , , 4 4 4 . The volume flow rate is determined as stated in equation (1) from the volume in between the indicators V PP and the time span t PP ∆ as given by the piston position indicator signals. It is therefore traceable to standards of length and time.
Turbine meters (TM) are currently used as transfer stand ards. TM measure the mass flow rate using the rotational speed of a turbine inserted in the fluid flow. The rotational speed of the turbine is metered by counting magnetically induced discrete signals. The volume flow rate V TM can be determined using a relationship between the number of sig nals per time period indicated by the TM and the volume flow rate, known from previous calibration of the TM. Two TM are connected in a row to minimize random measuring errors. The pressure at the TM is measured at their refer ence point and the temperature 2 diameters downstream of the TM.
The nozzle bank is used to set the flow rate. The critical nozzles are not necessary for the operation of the HPPP but provide the advantage of decoupling the calibration setup from pressure fluctuations downstream of the nozzle bank. It consists of several critical flow nozzles in parallel connection. The pressure downstream of the nozzles is always low enough to ensure critical flow in the nozzles.
The 4way valve is needed to revert the gas flow direction and move the piston back to its starting position after each calibration run. The check valve is used to prevent gas from flowing past the piston during the start of the reverse move ment and a safety valve is included to prevent high forces on the piston at the end of the piston reverse movement.
The HPPP is further described in detail in the [3,4] as well as the calibration facility pigsar TM [5,6].

Calibration uncertainty
The result of a calibration run is the relative deviation f of the corrected volume flow rate indicated by the TM V c TM and the corrected volume flow rate indicated by the HPPP V c PP as shown in equation (2). Several corrections are used in equation (2) to improve the calibration accuracy.
1. The volume flow rate indicated by the turbine meters V TM is corrected as shown in equation (3) to prevent an error caused by the discrete nature of the TM signals. t PP ∆ is the duration of the measurement phase as determined from the piston position indicator signals. t TM ∆ is the time span from the first TM signal after the start of the measurement phase to the first TM signal after the end of the measurement phase.
2. The temporal mean density over the measurement phase at the piston prover PP ρ and at the TM TM ρ , both deter mined from measured pressure and temperature, are used to take the density changes along the gas flow direction into account as shown in the first term of equation (4).
3. The temporal change of the stored mass in between the cylinder and the TM during the measurement phase is taken into account as shown in the second term of equa tion (4), with V encl being the enclosed volume, S ρ the spatial mean density in the enclosed volume at the start of the measurement phase and E ρ the spatial mean density in the enclosed volume at the end of the measurement phase.
Several possible errors in the calibration process lead to the measurement uncertainty of the calibrated TM. These are 1. uncertainty in the determination of the volume in between piston position indicators, 2. uncertainty in the determination of the mean density, 3. repeatability of the TM measurement, 4. leakage between piston and cylinder, 5. dynamic error of the TM, 6. uncertainty in the determination of the stored mass in the enclosed volume.
The dynamic error of the TM is a consequence of the incor rect measurement of fast fluctuating volume flow rates due to turbine inertia. This error can be diminished using a math ematical correction method [7], but the correction method as well has uncertainties.
The resulting uncertainty of the calibrated transfer stand ards is 0.06% [4,6]. The last two listed errors combined lead to an uncertainty of 0.035% [7]. They are of dynamic nature and a consequence of piston velocity fluctuations. Figure 2 shows the measured data for such piston velocity fluctuations in the measurement phase. The model is built to assess the dynamics of the HPPP and to find ways to reduce the piston velocity fluctuation and therefore also the measurement uncertainty of the HPPP.

Description of the model
A graphical representation of the model is shown in figure 3.
Several general assumptions were used in the model.
1. Pressure losses are proportional to the dynamic pressure, 2. the gas flow is one dimensional, 3. the system is adiabatic, 4. potential energy of the gas can be neglected, 5. no leakage occurs, 6. the heat transfer in the gas can be neglected in comparison to convective energy transport.
The gas at the inlet to the HPPP is assumed to have a con stant temperature and pressure. This is consistent with meas urement data and is modelled using a supply volume of infinite size. Equations (5) and (6) set these boundary conditions with T IN being the inlet temperature and p IN the inlet pressure.
The nozzle bank sets the other boundary condition. The used nozzles comply with ISO 9300 [8]. The nozzles are model led using a constant critical volume flow rate V N through all nozzles as shown in equation (7).
The valves have a linear opening function Y(t) and the mass flow ṁ V is proportional to the pressure drop across the valve p V ∆ as shown in equation (8) with ṁ n V and p n V ∆ being the nominal mass flow and pressure drop.
The start valves are used in the model to eliminate the influ ence of guessed initial conditions on the piston movement, as the stationary flow condition at the beginning of the running in phase is not known.
The HPPP is used to meter natural gas flow. Due to the high pressure and the high precision a real gas model is nec essary. A Modelica implementation of GERG 2008 with a constant gas composition out of 10 molecules is used. GERG 2008 derives the equation of state for natural gas from the free energy. It is described in detail in the references [9,10].
The enclosed gas volumes in the measuring cylinder change with piston movement. They can store mass m and internal energy mu as stated in equations (9) and (10). The volumes have i inlets or outlets. h is the specific enthalpy, v the mean velocity in a cross area A, p the static pressure and V the Volume. The pressure losses at inlets and outlets p ∆ are considered using constant coefficients A ζ as shown in equa tion (11) with ρ being the mean density. No gradient for the thermodynamic state and no fluid friction is considered in the volumes.
The position of the piston is determined from the equa tion of motion (12) with F F,P being the friction force on the piston, p P ∆ the pressures difference across the piston, A P the piston cross area, m P the piston weight and a P the piston acceleration.
The friction force on the piston is described as the sum of velocity independent Coulomb friction F C , velocity propor tional friction F v P P and Stribeck friction F e kv S P − as stated in equation (13).
The coulomb friction F C is modelled as a function of the piston position s P . This function is determined by measuring the power consumption of a linear motor moving the piston slowly through the cylinder. Measured data is only available for 80% of the cylinder length. After that the coulomb friction is assumed constant.
The velocity proportional friction F P was determined from measuring the pressure difference across the piston for different piston velocities and by linear interpolation of the measured data points.
The pipes can store mass m, internal energy mu and momentum mv as stated in equations (14)-(16). A finite volume discretisation in the direction of fluid flow is used. Each finite volume reaches from cross area i to cross area i + 1. In equations (14)-(16) ṁ is the mass flow, h the specific The turbine meters are modelled using a constant pressure loss coefficient TM ζ as stated in equation (17) with p TM ∆ being the pressure loss, ρ the spatial mean density and v A the mean velocity in the cross area A.
The pressure loss coefficient TM ζ is taken from measure ments. The relation between the indicated volume flow rate V TM ind and the true volume flow rate V TM in the TM can be model led as shown in equation (18) [7]. The time constant τ is modelled as stated in equation (19) using a linear relation between the time constant and the initial mass flow ṁ IC as approximately found in measurement.

Verification
As measure for the verification and validation, the relative deviation of the piston velocity from it's mean velocity in the measurement phase v ∆ , is used. v ∆ represents the piston velocity fluctuations and is calculated as shown in equa tion (20) using the piston velocity v P and the mean piston velocity v P determined from the distance l ∆ and the duration of the measurement phase t ∆ . For an easy comparison of dif ferent volume flow rates a normalized time t n is used in the figures. It is determined as in equation (21) from the time t, the start time of the measurement phase t s and the duration of the measurement phase t ∆ . v For the timeintegration the Solver Dassl included in Dymola ® is used with a relative tolerance of 10 −6 [2].   A further decrease of the relative tolerance does not change the model trajectory as shown in figure 4. No major change in the trajectories is detected when using other high order vari able step solvers implemented in Dymola ® .
Due to calculation time constraints it is not possible to use a high number of finite pipe volumes in conjunction with real gas behaviour. Here 4 discrete volumes in the first pipe and 2 volumes in the 2nd and 3rd pipe are used.
The verification of the model shows an increasing frequency and decreasing deflection of the relative piston velocity devia tion for increasing inlet pressures as shown in figure 5.
Due to a shorter duration of the runningin phase, the piston velocity fluctuation resulting from piston acceleration remains active during the measuring phase for higher volume flow rates as shown in figure 6.

Validation
The model accuracy is highly relevant due to the low measuring uncertainty of the high pressure piston prover. It depends on the uncertainty of the measured parameters used for the cali bration of the model, the mentioned general assumptions, the simplified mathematical description and the accuracy of the numerical algorithm.
Measured data for the piston velocity is used to validate the model. The piston velocity was measured for volume flow rates up to 100 m h 3 1 − using a laser distance measurement system.
The model validation shows relatively good accordance of the piston velocity fluctuations with measurement data for a volume flow rate of 100 m h 3 1 − , as shown in figure 7. The ground oscillation as well as the superimposed high fre quency oscillation show similar characteristics for a volume flow rate of 100 m h 3 1 − . This is also valid for different inlet pressures.
For lower volume flow rates the model is not able to repro duce the measured piston velocity fluctuations. As an example the piston velocity fluctuation in the model is compared with     − . Measuring uncertainties of the parameters used for model calibration and of the data used for validation might play an important role for low volume flow rates, as the mean velocity is lower. Also the used friction model for the piston is pre sumed to not be accurate enough for low volume flow rates. But as a volume flow rate of 25 m h 3 1 − is at the lower end of the high pressure piston prover operating area, this deficiency can be accepted.

Results
The model is used here to evaluate three different ways to reduce the piston velocity fluctuations in the measuring phase. The maximum deviation of the piston velocity from it's mean velocity v max ∆ is used as a measure for the piston velocity fluctuations. The mean piston velocity deviation would not be an adequate measure, as it does not limit the important piston velocity deviation at the start and end of the measurement phase, whereas the maximum piston velocity deviation does.
The results are shown for a volume flow rate of 100 m h 3 1 − and an inlet pressure of 20 bar. The solver is set according to verification and validation. Figure 9 shows the maximum relative deviation of the piston velocity from it's mean velocity in the measuring phase for different piston weights. A lower piston weight leads to lower maximum piston velocity fluctuation in the model. The real piston weight is 21.7 kg. Reducing the piston weight by 50% would lead to a significant drop of the piston velocity fluctuations. A way to achieve this reduction can be a change of the piston material from aluminum to fiber reinforced polymers.
Another way to reduce the piston velocity fluctuations in the model is shown in figure 10. As can be seen, the switching time of start valve 1 has a strong influence on the maximum deviation of the piston velocity from its mean velocity during the measuring phase below a switching time of 0.4 s. The switching time could be adopted to other volume flow rates using a controller. Figure 11 shows the maximum relative deviation of the piston velocity from it's mean velocity in the measuring phase as a function of the pressure loss coefficient of pipe 1. A sig nificant reduction of the piston velocity fluctuations can be achieved for higher pressure loss coefficients. A possibility to raise the pressure loss coefficient would be the integration of a flow straightener. The error due to a higher pressure loss between the HPPP and the TM can be avoided using correc tion step 2 as described in section 4.

Conclusions
Three independent measures to reduce the piston velocity fluctuations of the high pressure piston prover were identi fied using the developed model. A significant reduction of the maximum piston velocity fluctuation during the measuring phase was found achievable by lowering the piston weight, an appropriate setting of the start valve switching time and the integration of a flow straightener. These measures are expected to reduce the high pressure piston prover measuring uncertainty and can be realized with low effort.
As the computational model presented in this paper is written in Modelica, it was possible to reuse models from the Modelica Standard Library keeping the modelling time and cost low. The results show, that a modelling approach, as chosen in this work, can be very cost and time effective.