MetaEnergy
MetaEnergy

MetaFlow's Transient Flow Experiment on PSIG-1010

MetaFlow already calculates transient gas flow in production. On PSIG-1010, the 62-hour field case of Osiadacz and Chaczykowski, we evaluated its method-of-characteristics (MOC) solver alongside four other transient solvers we have built. The production MOC tracks delivery pressure with an RMS error of 31.4 kPa; the selected thermal DGSEM configuration has RMS errors of 28.4 kPa and 0.140 K in one run and is a candidate for future production integration.

MetaFlow already calculates transient gas flow in production. Its steady-state network solver establishes the operating condition of a pipeline network, and its isothermal method-of-characteristics (MOC) solver calculates how pressure and flow then evolve as demands and boundary pressures change. When demand at the far end of a long, large-diameter line rises, the line first draws on the gas stored inside it, the linepack, so the delivery pressure falls before the supply catches up.

Over the past two months we ran a transient flow experiment on one of the most complete public records of such behaviour: the PSIG-1010 field case of Andrzej J. Osiadacz and Maciej Chaczykowski. We evaluated the production MOC alongside four other transient solvers we have built, including a thermal research prototype developed for MetaFlow that uses a discontinuous Galerkin spectral element method (DGSEM) and also calculates the gas temperature. Our aim was to reproduce the authors' comparison with our own solvers and to identify which experimental approach to consider for future production integration. We scored each solver by its root-mean-square (RMS) error against the measured delivery pressure and, for thermal solvers, temperature. Four configurations with the same inputs anchor the comparison:

ConfigurationPressure RMS errorTemperature RMS error
Production isothermal MOC31.4 kPaFixed temperature (isothermal)
Isothermal DGSEM (product code, not enabled)32.5 kPaFixed temperature (isothermal)
Thermal DGSEM with steady heat exchange25.4 kPa0.905 K
Selected thermal DGSEM with wall and soil storage28.4 kPa0.140 K

These results use reconstructed, time-aligned inputs and case-specific settings, as described in Data and methods.

The production MOC tracks the measured delivery pressure through both demand events of the record with an RMS error of 31.4 kPa, on a line whose delivery pressure averages 6,181 kPa. The selected thermal DGSEM configuration adds heat storage in the pipe wall and the surrounding soil; in one run it tracks the delivery pressure with an RMS error of 28.4 kPa and the delivery gas temperature with 0.140 K. That combination, pressure and temperature from a single run, is why we are considering the approach for future integration into MetaFlow's production solver suite.

A field record worth building on

Osiadacz and Chaczykowski presented Verification of Transient Gas Flow Simulation Model as paper PSIG-1010 at the 2010 annual meeting of the Pipeline Simulation Interest Group [1]. For about 62 hours of operation, it publishes both the boundary records that drive a simulation and the delivery measurements that test it, together with the results of the authors' own non-isothermal models, one with a steady and one with a transient heat-transfer model. Chaczykowski's journal article of the same year examines the effect of the thermal model on transient flow [2]. Few published field comparisons are this complete, and it has been a reference throughout MetaFlow's development.

The line, drawn in Figure 1, is a 107 km section of the Yamal-Europe pipeline between Lwówek, the inlet of this section, and the Mallnow compressor station at the delivery end: DN1400-class pipe with a 1.3836 m inner diameter, buried about 1.5 m deep. Over the evaluation window the gas enters at a mean 7,153 kPa and 20.4 °C and leaves at a mean 6,181 kPa and 10.7 °C, at a mean mass flow of 612 kg/s, cooling by 9.7 K along the line. A pressure wave crosses the line in about 4.6 minutes; the gas takes about 3.8 hours.

The PSIG-1010 line: prescribed inputs, calculated outputs and how a transient calculation starts Conceptual schematic. A 107 km buried DN1400-class pipeline (1.3836 m inner diameter, about 1.5 m deep) runs from the Lwówek sending node to the Mallnow compressor station. Prescribed inputs, the measured records given to the models, enter at both ends: at the inlet the pressure (mean 7,153 kPa) and, for thermal models only, the gas temperature (mean 20.4 °C); at delivery the mass flow (mean 612 kg/s). Calculated outputs leave the delivery end: the pressure and, for thermal models only, the gas temperature, scored against the measured records, whose window means are 6,181 kPa and 10.7 °C. The models calculate pressure and flow along the line; thermal models also calculate temperature. Inset, conceptual: from steady state to transient. 1, a steady solution at the starting boundary values gives the initial operating condition. 2, that state is settled on the solver’s own discrete equations with the inputs held at their starting values. 3, the transient calculation from hour 9 to hour 71 advances pressure and flow, and in thermal models temperature, as the boundary records change. In MetaFlow’s production path the steady-state network solver supplies the starting state and the transient solver settles it on its own equations; the thermal prototype builds and settles its own. Scoring uses 241 points every 15 minutes over hours 11 to 71 of the published record; hours 9 to 11 are not scored. CONCEPTUAL SCHEMATIC One 107 km line, driven from both ends Prescribed input measured record Calculated output compared with measurements At the inlet Lwówek, sending node Prescribed input MEASURED MEAN Pressure 7,153 kPa Gas temperature 20.4 °C thermal models only At delivery Mallnow compressor station Prescribed input MEASURED MEAN Mass flow 612 kg/s buried about 1.5 m gas flow The models calculate pressure and flow along the line; thermal models also calculate temperature 107 km At delivery Mallnow compressor station Calculated output MEASURED MEAN scored against the measured records Pressure 6,181 kPa Gas temperature 10.7 °C thermal models only DN1400-class pipe · 1.3836 m inner diameter A pressure wave crosses the line in about 4.6 minutes; the gas takes about 3.8 hours FROM STEADY STATE TO TRANSIENT (CONCEPTUAL) Steady starting state a steady solution at the starting boundary values: the initial operating condition Settled start steady on the solver’s own equations, inputs held at their starting values Transient calculation hours 9 to 71: pressure and flow, and in thermal models temperature, advance as the boundary records change In MetaFlow’s production path, the steady-state network solver supplies the starting state and the transient solver settles it on its own equations; the thermal prototype builds and settles its own. Evaluation window: hours 11 to 71, 241 points every 15 minutes Not scored hours 9 to 11 11 41 71 Hours on the published time axis
Figure 1. Conceptual schematic of the PSIG-1010 field case of Osiadacz and Chaczykowski (2010). Prescribed inputs enter at both ends: the measured inlet pressure, the inlet gas temperature (thermal models only) and the delivery flow. The models calculate the delivery pressure and, in thermal models, the delivery temperature, which are compared with the measurements; values are measured means over the evaluation window. The inset shows how a transient calculation starts from a steady state.

At the inlet we prescribe the measured pressure and, for thermal models, the measured gas temperature; at the delivery end we prescribe the measured mass flow. The records are published as figures, which we digitized, and we score 241 points at 15-minute spacing from hour 11 to hour 71 of the paper's time axis. The RMS error is the root mean square of model minus measured over those points, the bias is its mean, and the shape error, or bias-removed RMS, is what remains once the mean offset is removed. It shows how well a model follows the rises and falls.

From steady state to transient flow

MetaFlow's production calculation runs in two stages. For prescribed boundary pressures and demands, the steady-state network solver computes the pressure at every node and the flow in every pipe; that solution is the starting state of every transient run, from which the MOC takes over. Both production solvers handle isothermal gas networks of horizontal pipes without control devices.

The starting state is built from the boundary records themselves. The pressure and demand schedules are evaluated at the start time and the steady solver runs on those values, so a schedule whose first value differs from the values stored with the network model does not send a spurious pressure wave down the pipes. Each gas pipe is filled along the pressure profile of steady isothermal flow. That profile is steady for the continuous equations but not exactly for the discrete ones, so before the first step the transient solver adjusts it until one of its own time steps, with the inputs held at their starting values, leaves it unchanged. The run therefore starts in equilibrium with its own equations; the inset in Figure 1 sketches the sequence.

From there the MOC advances pressure and mass flux along the characteristic lines on which pressure waves travel. Friction is treated implicitly, so within a step it can slow the flow but cannot reverse it on its own. The time step is limited by the wave speed and shortened to land exactly on the timestamps of the boundary records, and a step that would leave the method's validity range stops the run with a reported error rather than a doubtful answer.

The steady solver is the starting point of the transient calculation, not one of the transient solvers compared below. The thermal prototype in this experiment builds its own starting state, independently of the production code; we describe it with the selected approach below.

Five transient solvers, three numerical families

A transient pipeline model makes two separate choices. The physical model decides which conservation laws are solved and how heat is treated; the numerical method decides how those equations are discretized in space and time. Figure 2 keeps them on separate axes: five transient solvers in three numerical families, each run with one or more physical models.

We used three physical models:

  • Isothermal: mass and momentum at a fixed gas temperature, here 15.5 °C, with the CNGA real-gas compressibility correlation in the MOC, the isothermal DGSEM and the isothermal implicit FV; it predicts pressure and flow, but not temperature.
  • Thermal with steady heat exchange: adds the total energy equation, with pressure and caloric properties from the AGA8-DC92 equation of state [3] for the paper's gas composition, and a steady heat exchange with the ground at 277 K, set by the overall heat-transfer coefficient to the ground, referred to the inner pipe surface. We call it the heat-transfer coefficient U; our untuned value is 2.09 W/m²K.
  • Thermal with wall and soil storage: replaces the steady exchange with a radial column of wall and soil layers at every node, which takes up heat and gives it back as conditions change.

The three numerical families:

  • Method of characteristics (MOC): follows pressure waves along their characteristic lines with an explicit time step, run here on 500 segments of 214 m.
  • Implicit finite volume (FV): balances the fluxes of mass, momentum and, in thermal runs, energy from cell to cell, and allows large implicit time steps.
  • DGSEM: splits the pipe into a few long elements, each carrying a high-order polynomial, coupled to its neighbours through fluxes at the element faces [4].

The five solvers and the role each played, with the paper's two models as the published reference:

  • Isothermal MOC: MetaFlow's production transient solver.
  • Isothermal implicit FV and isothermal DGSEM: in MetaFlow's product code but not enabled; here they cross-check the MOC on the same isothermal model, including the small convective term that the MOC omits.
  • Thermal FV: development code, and our target-blind reference; we froze its settings before we looked at the measured delivery curves and scored it once.
  • Thermal DGSEM: the research prototype, a standalone solver for one horizontal pipe that is not wired into any MetaFlow route. We built it after we had seen the measurements, so we call its unadjusted configuration our untuned baseline rather than target-blind. With wall and soil storage it is the configuration we selected for further development, an engineering choice made during about 240 exploratory runs on this record: with the delivery-flow record left on the paper's time axis, it lowered both errors against the untuned baseline and reached a temperature error of 0.140 K.
Physical models and numerical families of the five transient solvers used on PSIG-1010 Conceptual schematic. Rows are physical models: isothermal (mass and momentum at a fixed 15.5 °C), thermal with steady heat exchange (adds the energy equation with AGA8 gas properties and steady heat exchange with the ground), and thermal with wall and soil heat storage. Columns are numerical families: method of characteristics (500 segments), implicit finite volume (64 or 500 cells) and DGSEM spectral elements (16 or 32 nodes). Isothermal: the method of characteristics is MetaFlow’s production transient solver, highlighted in teal, in MetaFlow today; finite volume and DGSEM are cross-checks in product code, not enabled. Thermal with steady heat exchange: finite volume is development code, the target-blind reference; DGSEM is the research prototype, the untuned baseline. Thermal with storage: finite volume is development code where the storage model was built first; DGSEM is the research prototype in its selected configuration, marked with a dashed outline. The method of characteristics was not run with the thermal models. Production means shipped for isothermal gas networks of horizontal pipes without control devices; the research prototype is a standalone solver, not wired into any MetaFlow route. CONCEPTUAL SCHEMATIC Physics and numerics are two separate choices Five transient solvers in three numerical families, and the role of each on PSIG-1010 NUMERICAL FAMILY: HOW THE EQUATIONS ARE SOLVED Method of characteristics Follows pressure waves along their characteristic lines 500 segments Implicit finite volume Balances fluxes cell by cell, with implicit time steps 64 or 500 cells DGSEM spectral elements High-order polynomials inside a few long elements 16 or 32 nodes PHYSICAL MODEL: WHICH EQUATIONS ARE SOLVED Isothermal Mass and momentum; the gas is held at a fixed 15.5 °C Production solver In MetaFlow today Cross-check In product code, not enabled Cross-check In product code, not enabled Thermal, steady heat exchange Adds the energy equation (AGA8 gas properties) and steady heat exchange with the ground Not run on this case Development code Target-blind reference Research prototype Untuned baseline Thermal with wall and soil storage Adds the heat stored in the coatings, the steel and the surrounding soil Not run on this case Development code Storage built here first Research prototype Selected configuration Production: shipped for isothermal gas networks of horizontal pipes without control devices. Cross-check: in product code, not enabled. Development code: on development branches. Research prototype: a standalone solver, not wired into any MetaFlow route.
Figure 2. Conceptual overview: five transient solvers in three numerical families (columns), run with up to three physical models (rows). MetaFlow's production MOC is highlighted; every other filled cell is labelled with its role and status, and the selected configuration is the thermal DGSEM with wall and soil storage.

Delivery pressure over 60 hours

The record contains two demand events, each opened and closed by a step change in delivery flow. In the first, hours 21 to 25, the measured delivery pressure drops some 220 kPa from its level at hour 18 to about 5,980 kPa. The second is deeper: over hours 45 to 49 it falls about 250 kPa from its level at hour 44.5 to about 5,960 kPa, the lowest of the window, before recovering.

We first compared like with like, with every input exactly as plotted in the paper. Our untuned thermal DGSEM baseline then scores 34.9 kPa RMS against 33.7 and 34.8 kPa for the paper's steady and transient models, and our target-blind thermal FV scores 34.8 kPa. An independent implementation of comparable physics lands where the authors' models do. In the target-blind study, at 64 cells with AGA8-DC92 in both arms, the energy equation with the measured inlet temperature and steady heat exchange to the ground lowered the pressure shape error from 38.3 kPa, with the gas held at 15.5 °C, to 34.3 kPa.

With the time-aligned inputs of the opening table, Figure 3 shows the production MOC and the selected configuration against the measured delivery pressure over the whole window. Both track the same rises and dips: the production MOC with an RMS error of 31.4 kPa, the selected configuration with 28.4 kPa.

Delivery pressure, hours 11 to 71: production MOC and the selected thermal DGSEM configuration Line chart of delivery pressure between about 5,940 and 6,310 kPa over paper hours 11 to 71. The measured trace, MetaFlow's production MOC and the selected thermal DGSEM configuration follow the same rises and the two deep dips, with the largest differences around the four step changes in delivery flow. Delivery pressure, hours 11 to 71 Pressure RMS error Measured (digitized) MetaFlow production MOC 31.4 kPa MetaFlow thermal DGSEM, selected configuration 28.4 kPa Delivery pressure (kPa) 5,900 6,000 6,100 6,200 6,300 11 21 31 41 51 61 71 Hours on the paper’s time axis (evaluation window 11 to 71 h) Measured record digitized from Osiadacz and Chaczykowski (2010).
Figure 3. Delivery pressure over hours 11 to 71 against the digitized measurement: MetaFlow's production MOC, 31.4 kPa RMS, and the selected thermal DGSEM configuration, 28.4 kPa RMS, on 241 scoring points. Configuration as described in Data and methods. Measured record digitized from Osiadacz and Chaczykowski (2010).

Most of what remains sits around the four step changes in delivery flow. Away from them, the selected configuration follows the delivery pressure with a shape error of 4.5 kPa; within three hours of the steps a shape error of about 46 kPa remains, and these windows, a third of the evaluation window, hold 98% of the remaining squared difference. At the steps, the measured response differs from every model we ran and from the paper's two models alike.

Delivery gas temperature

Temperature is the more demanding quantity, and only the thermal solvers calculate it; the isothermal solvers hold the gas at 15.5 °C. Over the whole window the measured delivery temperature moves within about one kelvin, so a model has to get both the timing and the size of small swings right. For scale, a flat line at the measured mean scores 0.177 K RMS.

Figure 4 shows the selected configuration's delivery temperature. Its RMS error is 0.140 K, with a mean offset of only +0.02 K and a shape error of 0.138 K; through the second demand event the measured temperature swings 0.94 K and the model 0.98 K. The heat-storage subsection below explains what damps that response.

Delivery gas temperature, hours 11 to 71: the selected thermal DGSEM configuration Line chart of delivery gas temperature over paper hours 11 to 71. The measured trace moves between about 10.1 and 11.0 °C; the selected MetaFlow thermal DGSEM configuration follows its level and both dips, including the swing of about one kelvin through the second demand event, with departures of up to about 0.4 K, mostly between hours 35 and 51. Delivery gas temperature, hours 11 to 71 RMS error Shape error Measured (digitized) MetaFlow thermal DGSEM, selected configuration 0.140 K 0.138 K Isothermal solvers hold the gas at 15.5 °C and do not calculate delivery temperature. Delivery gas temperature (°C) 10.0 10.5 11.0 11 21 31 41 51 61 71 Hours on the paper’s time axis (evaluation window 11 to 71 h) Measured record digitized from Osiadacz and Chaczykowski (2010).
Figure 4. Delivery gas temperature over hours 11 to 71 against the digitized measurement: the selected thermal DGSEM configuration, 0.140 K RMS with a shape error of 0.138 K, on 241 scoring points. The isothermal solvers hold the gas at 15.5 °C and do not calculate it. Configuration as described in Data and methods. Measured record digitized from Osiadacz and Chaczykowski (2010).

Four configurations side by side

Figure 5 compares the four configurations introduced above, all with the same inputs; Data and methods describes their case-specific settings.

The production MOC and the isothermal DGSEM do not calculate temperature and score 31.4 and 32.5 kPa, two isothermal schemes about 1 kPa apart. The thermal DGSEM with steady heat exchange, our untuned baseline with the same inputs, has the lowest pressure error of the four, 25.4 kPa, and a temperature error of 0.905 K. The selected configuration has a pressure error 3 kPa higher than that variant's and a temperature error more than six times smaller, 0.140 K.

MetaFlow's transient flow experiment on PSIG-1010: delivery pressure and temperature errors of the four selected configurations, hours 11 to 71 Horizontal bars of RMS error against the digitized measurement over paper hours 11 to 71 (241 points), all four with the same time-aligned inputs: Production isothermal MOC 31.4 kPa, fixed temperature (isothermal); Isothermal DGSEM (product code, not enabled) 32.5 kPa, fixed temperature (isothermal); Thermal DGSEM with steady heat exchange 25.4 kPa, 0.905 K (lowest pressure error); Selected thermal DGSEM with wall and soil storage 28.4 kPa, 0.140 K. The production MOC and the selected configuration are highlighted. Delivery pressure and temperature errors, hours 11 to 71 Pressure RMS error, kPa Temperature RMS error, K Production isothermal MOC Production isothermal MOC: pressure RMS error 31.4 kPa 31.4 fixed temperature (isothermal) Isothermal DGSEM (product code, not enabled) Isothermal DGSEM (product code, not enabled): pressure RMS error 32.5 kPa 32.5 fixed temperature (isothermal) Thermal DGSEM with steady heat exchange lowest pressure error Thermal DGSEM with steady heat exchange: pressure RMS error 25.4 kPa 25.4 Thermal DGSEM with steady heat exchange: temperature RMS error 0.905 K 0.905 Selected thermal DGSEM with wall and soil storage Selected thermal DGSEM with wall and soil storage: pressure RMS error 28.4 kPa 28.4 Selected thermal DGSEM with wall and soil storage: temperature RMS error 0.140 K 0.140 0 10 20 30 0 0.25 0.5 0.75 1.0 Scores against the delivery record digitized from Osiadacz and Chaczykowski (2010).
Figure 5. Pressure and temperature RMS errors over hours 11 to 71 for the four configurations of the opening table, all with the same inputs. MetaFlow's production MOC and the selected thermal DGSEM configuration are highlighted; the isothermal solvers do not calculate temperature. Configuration as described in Data and methods; scores against the delivery record digitized from Osiadacz and Chaczykowski (2010).

The two finite-volume solvers were run only with the inputs as plotted, so they appear in the numerical checks below rather than here; we do not estimate their scores with time-aligned inputs.

Inside the selected transient-flow approach

The selected configuration couples the thermal DGSEM gas model to heat storage in the pipe wall and the soil. Its calculation works as follows.

Starting state. The prototype marches a steady solution along the pipe from the inlet state and delivery flow applied at the start of the run, then runs its full transient scheme for 5,000 seconds with every boundary value held fixed, so the calculation starts settled on its own discrete equations. Each wall and soil column then starts at its steady temperature profile.

Boundary conditions. At the inlet, the measured pressure and gas temperature enter through a ghost state, a fictitious neighbour built from those values, and the flux at the inlet face. At the delivery end, at every stage of every time step, the solver adjusts the pressure of the outlet ghost state until the flux leaving the pipe matches the measured mass flow. Two conditions at the inflow and one at the outflow are what subsonic pipe flow requires; the delivery pressure and temperature we score are read from the outlet state.

Mass, momentum and energy. The prototype solves the one-dimensional balances of mass, momentum and total energy of the gas along a horizontal line. Wall friction acts in the momentum balance and, because the wall does not move, does no work on the gas. In the energy balance it only turns kinetic into internal energy.

Spatial representation. The pipe is divided into four elements of 26.75 km. In each, the solution is a polynomial of degree three held at four Legendre-Gauss-Lobatto nodes, written P3/E4 (polynomial degree 3, four elements): 16 nodal values for the whole 107 km. Within an element the fluxes are evaluated in a split form that keeps the aliasing errors of a high-order polynomial under control [5]; it needs only the values at the nodes, so no equation-of-state evaluation is made between them. Neighbouring elements exchange mass, momentum and energy through an HLL (Harten-Lax-van Leer) flux at their shared faces. The flow is smooth, with a Mach number of about 0.02, and the prototype uses no limiter or artificial viscosity.

Time advancement. The solution advances with the three-stage strong-stability-preserving Runge-Kutta method (SSP-RK3). The step follows a Courant-Friedrichs-Lewy (CFL) stability rule: 0.9/(2p + 1), with p = 3 the polynomial degree, of the time the fastest wave needs to cross an element, here at most about 8 seconds. The step is shortened when necessary to land exactly on every boundary-record time and every 5-minute output time.

Recorded runtime. The selected 62-hour run took about three minutes of recorded wall time on a shared machine.

Friction. Friction follows the Colebrook law with the authors' roughness of 0.01 mm; at a Reynolds number of about 4.5 × 10⁷ the flow is fully turbulent.

Real-gas properties. Pressure, temperature and every caloric property come from the AGA8-DC92 equation of state [3] with the paper's gas composition. Because the enthalpy carries its full dependence on pressure, Joule-Thomson cooling, the temperature drop of a real gas as it expands, is part of the solution rather than a separate correction.

Energy and heat exchange. Heat leaves the gas toward the ground at 277 K. With steady heat exchange it crosses the gas film, coating, steel, polyethylene and soil at once, through the heat-transfer coefficient U; with storage it passes through the wall and soil column described next. The coefficient sets the level of the delivery temperature rather than the shape of its response.

Heat storage in the pipe wall and the soil

Why should a 19.22 mm steel wall matter for a pipe 1.38 m across? Per metre of pipe, the steel holds about 1.4 times the heat capacity of the gas inside it and exchanges heat with the gas in about 5 minutes.

A steady heat-exchange law has no store. Its only path is a direct loss to the ground, with a time constant of about 7 hours for this line at the untuned coefficient. The demand events last a few hours, so in a steady-exchange model the gas compresses and expands with little heat exchange and its temperature swings too far. With storage, the steel returns heat within minutes as the gas cools and takes it up as the gas warms, pulling the gas toward isothermal behaviour on the time scale of the events, which matches the damped swings of the measured record.

Each node of the thermal DGSEM carries its own radial column of wall and soil layers (Figure 6). The column takes one implicit step per gas step and receives exactly the heat the gas gives up. With the inputs as plotted, adding storage reduced the temperature shape error from 0.511 K to 0.112 K, without additional parameter fitting in this paired comparison. The same pull toward isothermal behaviour softens the pressure response: in the same comparison it adds about 2.7 kPa to the pressure RMS at a fixed heat-transfer coefficient.

Wall and soil heat storage around the buried pipe Conceptual schematic, not a measured temperature field and not to scale. In the storage model, a wedge of the pipe cross-section at every node runs from the gas outward: a gas film, an internal coating 0.5 mm, steel 19.22 mm, external polyethylene 3 mm, then soil shells that are thin near the pipe and widen outward, ending at the outer edge of the soil column, sized from the 1.5 m burial depth, where the ground is held at 277 K. Heat moves both ways between the gas and the steel, and slowly outward through the soil. The steel holds about 1.4 times the heat capacity of the gas in the same length of pipe. Time scales on a logarithmic axis: gas and steel exchange heat in about 5 minutes; demand events on the line last a few hours; without storage the gas exchanges heat with the ground over about 7 hours. CONCEPTUAL SCHEMATIC, NOT A MEASURED TEMPERATURE FIELD Heat stored in the pipe wall and the soil The layers the storage model places around the gas at every node (not to scale) Gas Soil shells thin near the pipe, wider outward slow loss to the ground Gas film Coating 0.5 mm Steel 19.22 mm Polyethylene 3 mm The steel holds about 1.4 times the heat capacity of the gas in the same length of pipe fast two-way exchange Ground held at 277 K at the outer edge of the soil column, sized from the 1.5 m burial depth WHY STORAGE MATTERS: THREE TIME SCALES Gas and steel exchange heat about 5 min Demand events on the line a few hours Gas to ground, without storage about 7 h 1 min 10 min 1 h 4 h 12 h
Figure 6. In the storage model, every node carries its own radial column of wall and soil layers that take up heat and give it back to the gas. The steel exchanges heat with the gas within minutes while demand events last hours, so storage damps the temperature swings that a steady heat-exchange model overstates.

Checking the numerics

Figure 7 collects the checks on these computations, all with the inputs as plotted. Panel (a) compares discretizations of the same physics (the MOC omits the small convective term). The production MOC, the isothermal DGSEM and the isothermal implicit FV score 40.9, 41.7 and 42.5 kPa RMS. Nearly all of that spread is a constant offset that each scheme inherits from its own discrete steady starting state; their bias-removed pressure RMS scores span 0.096 kPa. Close scores do not mean identical curves. The three isothermal pressure curves differ by about 2 to 4 kPa RMS, almost entirely as that offset. On 16 nodes the thermal DGSEM baseline matches the target-blind thermal FV scores within 0.11 kPa and 0.02 K, and their pressure curves differ by 0.52 kPa RMS.

Panel (b) is a mesh-refinement study on a logarithmic axis. For the isothermal DGSEM, the inter-grid RMS differences in delivery pressure from its 64-node run are 2.271, 0.173 and 0.081 Pa on 12, 24 and 32 nodes. For the thermal DGSEM baseline, with steady heat exchange, the inter-grid RMS differences from its 128-node run are 62 Pa on 16 nodes and 9 Pa on 64. These are differences between grids, not errors against the measurement, which stay at 41.65 kPa and 34.89 to 34.90 kPa on every grid. Doubling the number of elements of the storage model changes its score by less than 0.001 kPa.

In the thermal DGSEM the gas mass and energy balances close at every time step to within about one part in 10¹⁵. Our untuned storage model, built first in our thermal FV solver, scores 0.648 K there on 64 cells against 0.660 K in the thermal DGSEM, with shape errors of 0.114 and 0.112 K. The production MOC passes code-level tests against exact solutions: a closed-pipe resonance reaches its first trough at 29.982 s against the analytical 30 s, and a pressure pulse arrives at 22.4955 s against 22.5 s. These are checks of consistency and resolution, not of field accuracy: the tested mesh differences are much smaller than the mismatch with the digitized field record. Every MetaFlow score in this article is recomputed from its recorded run output.

Numerical checks: scores of the same physics on different discretizations, and inter-grid RMS differences in delivery pressure beside the error against the measurement (a) Error against the measurement with inputs as plotted and no settings fitted by us. Isothermal: MOC (production), 500 segments RMS 40.875 kPa, bias-removed 38.284 kPa; DGSEM, P3/E8 RMS 41.654 kPa, bias-removed 38.346 kPa; Implicit FV, 500 cells RMS 42.535 kPa, bias-removed 38.380 kPa; bias-removed RMS scores spanning 0.096 kPa, RMS errors spanning 1.660 kPa; the three curves differ from one another by about 2 to 4 kPa RMS. Thermal, steady heat exchange with the ground at U 2.09: FV, 64 cells (target-blind) RMS 34.781 kPa, bias-removed 34.304 kPa, temperature RMS 0.861 K; DGSEM, P3/E4 (untuned baseline) RMS 34.893 kPa, bias-removed 34.396 kPa, temperature RMS 0.880 K; the two pressure curves differ by 0.52 kPa RMS. (b) Log scale. Inter-grid RMS difference of the delivery-pressure trace over all output times: isothermal DGSEM against its P3/E16 run (64 nodes) P2/E4 2.271 Pa, P2/E8 0.173 Pa, P3/E8 0.081 Pa; thermal DGSEM baseline physics (steady heat exchange) against its P3/E32 run (128 nodes) P3/E4 61.7 Pa, P2/E8 57.3 Pa, P3/E8 23.2 Pa, P3/E16 8.9 Pa. Error against the measurement on every grid: isothermal DGSEM 41.65 kPa, thermal DGSEM baseline 34.89 to 34.90 kPa. (a) Same physics, different discretizations Error against the measurement, inputs as plotted, no settings fitted by us shape error (bias-removed RMS) RMS error Temperature RMS ISOTHERMAL, FIXED 15.5 °C GAS MOC (production), 500 segments MOC (production), 500 segments: bias-removed RMS error 38.284 kPa MOC (production), 500 segments: RMS error 40.875 kPa 38.284 40.875 n/a DGSEM, P3/E8 DGSEM, P3/E8: bias-removed RMS error 38.346 kPa DGSEM, P3/E8: RMS error 41.654 kPa 38.346 41.654 n/a Implicit FV, 500 cells Implicit FV, 500 cells: bias-removed RMS error 38.380 kPa Implicit FV, 500 cells: RMS error 42.535 kPa 38.380 42.535 n/a bias-removed RMS scores span 0.096 kPa RMS errors span 1.660 kPa the three curves differ from one another by about 2 to 4 kPa RMS THERMAL, STEADY HEAT EXCHANGE WITH THE GROUND, U 2.09 W/m²K FV, 64 cells (target-blind) FV, 64 cells (target-blind): bias-removed RMS error 34.304 kPa FV, 64 cells (target-blind): RMS error 34.781 kPa 34.304 34.781 0.861 K DGSEM, P3/E4 (untuned baseline) DGSEM, P3/E4 (untuned baseline): bias-removed RMS error 34.396 kPa DGSEM, P3/E4 (untuned baseline): RMS error 34.893 kPa 34.396 34.893 0.880 K RMS scores differ by 0.11 kPa differ by 0.02 K the two pressure curves differ by 0.52 kPa RMS 32 34 36 38 40 42 kPa (b) Mesh refinement Delivery pressure, inputs as plotted, no settings fitted by us; log scale inter-grid RMS difference from the reference run error against the measurement (RMS) ISOTHERMAL DGSEM INTER-GRID RMS DIFFERENCE FROM ITS P3/E16 RUN (64 NODES) P2/E4 · 12 nodes P2/E4 (12 nodes): inter-grid RMS difference from its P3/E16 run 2.271 Pa 2.271 Pa P2/E8 · 24 nodes P2/E8 (24 nodes): inter-grid RMS difference from its P3/E16 run 0.173 Pa 0.173 Pa P3/E8 · 32 nodes P3/E8 (32 nodes): inter-grid RMS difference from its P3/E16 run 0.081 Pa 0.081 Pa THERMAL DGSEM, BASELINE PHYSICS (STEADY HEAT EXCHANGE) INTER-GRID RMS DIFFERENCE FROM ITS P3/E32 RUN (128 NODES) P3/E4 · 16 nodes P3/E4 (16 nodes): inter-grid RMS difference from its P3/E32 run 61.7 Pa 61.7 Pa P2/E8 · 24 nodes P2/E8 (24 nodes): inter-grid RMS difference from its P3/E32 run 57.3 Pa 57.3 Pa P3/E8 · 32 nodes P3/E8 (32 nodes): inter-grid RMS difference from its P3/E32 run 23.2 Pa 23.2 Pa P3/E16 · 64 nodes P3/E16 (64 nodes): inter-grid RMS difference from its P3/E32 run 8.9 Pa 8.9 Pa ERROR AGAINST THE MEASUREMENT (RMS) · EVERY GRID Isothermal DGSEM isothermal DGSEM, error against the measurement on every grid: 41.65 kPa 41.65 kPa Thermal DGSEM baseline thermal DGSEM baseline, error against the measurement on every grid: 34.89 to 34.90 kPa 34.89 to 34.90 kPa 0.01 Pa 0.1 Pa 1 Pa 10 Pa 100 Pa 1 kPa 10 kPa 100 kPa
Figure 7. Numerical checks, inputs as plotted. (a) Scores of the same physics on different discretizations: the isothermal schemes' bias-removed pressure RMS scores span 0.096 kPa, and the thermal FV and thermal DGSEM scores differ by 0.11 kPa and 0.02 K; close scores are not identical curves. (b) Inter-grid RMS differences in delivery pressure against each scheme's reference run (P3/E16 for the isothermal DGSEM, P3/E32 for the thermal baseline), beside the errors against the measurement.

What comes next: considering production integration

On this field case the production MOC tracked the delivery pressure through 60 hours of changing demand with an RMS error of 31.4 kPa, and the selected thermal DGSEM configuration, in one run, had RMS errors of 28.4 kPa in pressure and 0.140 K in temperature.

The selected configuration's results make it a candidate for future integration into MetaFlow's production solver suite. Its combination of transient pressure calculation, temperature modelling and a compact numerical representation, 16 nodes for a 107 km line, gives us a concrete basis for the next stage of development.

We are considering that integration alongside broader network testing, additional field comparisons and the engineering work needed to meet MetaFlow's production standards. Our existing production solvers remain the foundation of the platform while this evaluation continues.

We are grateful to Andrzej J. Osiadacz and Maciej Chaczykowski for publishing a field comparison detailed enough to build on. If you operate or study a line with records like these, we would be glad to hear from you.

Data and methods

Records and scoring. We digitized all records from the published figures of the PSIG-1010 paper [1]: inlet pressure, inlet temperature and delivery flow as inputs, delivery pressure and temperature as targets. We score 241 points every 15 minutes over hours 11 to 71 of the paper's time axis and read the temperature figure on its °F scale.

Time-aligned inputs. In the published figures, at two inlet steps during flat demand, delivery pressure completes half its rise before inlet pressure does, which a passive pipeline driven only by these records cannot do. For results with time-aligned inputs we apply the digitized inlet pressure record, and for the thermal results in this article also the inlet temperature record, 50 minutes earlier than plotted. This lowers the pressure error of every model we reran, the production MOC included, by 9 to 10 kPa. Delivery flow, calculated output curves and measured targets are not shifted. We treat this as a reconciliation of the published graphical records, not as an identified error in the original measurement system; which record carries the offset is not known.

Calibration. Every model uses the authors' pipe roughness of 0.01 mm, which they calibrated on this record. For the selected configuration we set U = 1.78 W/m²K on this record through the soil conductivity. In a split-sample test on this record, values fitted on one half carried over to the other for U and for a large inlet advance; neither half fixes the exact advance.

References

  1. Osiadacz, A. J., and Chaczykowski, M. (2010). Verification of Transient Gas Flow Simulation Model. Paper PSIG-1010, PSIG Annual Meeting, Bonita Springs, Florida, May 2010. Pipeline Simulation Interest Group. Available through OnePetro.
  2. Chaczykowski, M. (2010). Transient flow in natural gas pipeline: the effect of pipeline thermal model. Applied Mathematical Modelling, 34(4), 1051 to 1067. https://doi.org/10.1016/j.apm.2009.07.017
  3. Starling, K. E., and Savidge, J. L. (1994). Compressibility Factors of Natural Gas and Other Related Hydrocarbon Gases. Transmission Measurement Committee Report No. 8, American Gas Association, Arlington, VA.
  4. Kopriva, D. A. (2009). Implementing Spectral Methods for Partial Differential Equations: Algorithms for Scientists and Engineers. Springer, Dordrecht. https://doi.org/10.1007/978-90-481-2261-5
  5. Kennedy, C. A., and Gruber, A. (2008). Reduced aliasing formulations of the convective terms within the Navier-Stokes equations for a compressible fluid. Journal of Computational Physics, 227(3), 1676 to 1700. https://doi.org/10.1016/j.jcp.2007.09.020

About the author

Technology executive for critical infrastructure with 18+ years of experience improving energy operations. Combines practical operator knowledge with SCADA monitoring and control, remote measurements (telemetry), data analysis, forecasting and applied mathematics.

Blog