Validation: Greyvenstein and Laurie (1994)

Verify the mean-flow solver against Example 3 of Greyvenstein and Laurie (1994), a segregated CFD pipe-network benchmark.

Source notebook: examples/validation/greyvenstein_laurie_network.ipynb. This page is executed from it at build time.

G. P. Greyvenstein and D. P. Laurie, “A segregated CFD approach to pipe network analysis,” Int. J. Numer. Methods Eng. 37, 3685-3705 (1994).

Example 3 (their Fig. 9, Tables III-IV) is the paper’s only compressible case and the natural benchmark for Nefes: a 29-pipe compressed-air distribution network with a constant Darcy friction factor f = 0.03, isothermal at T = 288.15\,\mathrm{K} (15 ^\circC), fed by two 6-bar supplies and drained by thirteen 3-bar demands. Every terminal pressure is prescribed, so it is a pure Dirichlet problem: the mass flows are the answer.

Why this is the right case for Nefes. The paper’s method solves a lumped Darcy-Weisbach pressure-drop law \Delta p = K\,\rho\,Q^2 on each pipe with K = fL/D. Nefes’s pipe element is exactly that law: its docstring calls it the “DUCT (+) LOSS unification (Greyvenstein-Laurie)”: one length-bearing element that drops total pressure by \Delta p_t = K\,\tfrac12\rho u^2, K = fL/D, on the mean flow and carries its length for the acoustic phase. So reproducing Example 3 needs no new physics: just wiring.

What we expect. Nefes solves the fully-coupled, fully-compressible problem (static/total split, energy in the loop, exact complex-step Jacobian), whereas the paper uses a low-Mach, static-\approx-total, frozen-temperature model. At this operating point the flow is firmly low-Mach (M \lesssim 0.05), so the two should agree to 3-4 significant figures, with the small residual being physically real and in Nefes’s favour.

import numpy as np
import plotly.graph_objects as go

import nefes
from nefes.elements import catalog as cat
from nefes.plotting import use_nefes_theme, COLORWAY

use_nefes_theme()  # the bundled Nefes plotly template
'nefes'

1. The network data (Table III)

Each row is one pipe: its upstream node, downstream node, diameter D [m] and length L [m]. The friction factor is fixed at f = 0.03 and the flow area of a pipe is \pi D^2/4 (circular).

R_AIR, GAMMA, T_ISO = 287.0, 1.4, 288.15
FRICTION = 0.03
P_SUPPLY, P_DEMAND = 6.0e5, 3.0e5  # Pa

# Table III: element number -> (upstream node, downstream node, D [m], L [m]).
PIPES = {
    1: (1, 2, 0.019, 200),
    2: (2, 3, 0.01588, 400),
    3: (3, 5, 0.01588, 400),
    4: (5, 6, 0.010, 100),
    5: (6, 7, 0.010, 100),
    6: (6, 8, 0.010, 100),
    7: (8, 9, 0.010, 100),
    8: (8, 10, 0.010, 100),
    9: (5, 11, 0.01588, 400),
    10: (11, 12, 0.010, 100),
    11: (11, 13, 0.01588, 400),
    12: (14, 13, 0.019, 200),
    13: (13, 15, 0.01588, 400),
    14: (15, 16, 0.010, 100),
    15: (17, 15, 0.01588, 400),
    16: (17, 18, 0.010, 100),
    17: (18, 19, 0.010, 100),
    18: (18, 20, 0.010, 100),
    19: (21, 17, 0.01588, 400),
    20: (21, 22, 0.010, 100),
    21: (23, 21, 0.01588, 400),
    22: (23, 24, 0.010, 100),
    23: (24, 25, 0.010, 100),
    24: (26, 23, 0.01588, 400),
    25: (26, 27, 0.010, 100),
    26: (2, 26, 0.01588, 400),
    27: (3, 4, 0.010, 100),
    28: (5, 28, 0.010, 100),
    29: (24, 29, 0.010, 100),
}
SUPPLY_NODES = (1, 14)  # 6.0 bar terminals
DEMAND_NODES = (4, 7, 9, 10, 12, 16, 19, 20, 22, 25, 27, 28, 29)  # 3.0 bar terminals

LEAVES = set(SUPPLY_NODES) | set(DEMAND_NODES)
PAPER_NODES = sorted({n for e in PIPES.values() for n in e[:2]})
INTERIOR = [n for n in PAPER_NODES if n not in LEAVES]  # solved junctions
print(
    f"{len(PIPES)} pipes, {len(INTERIOR)} interior junctions, "
    f"{len(SUPPLY_NODES)} supplies + {len(DEMAND_NODES)} demands = {len(LEAVES)} terminals"
)
29 pipes, 14 interior junctions, 2 supplies + 13 demands = 15 terminals

The published results we will check against: Table IV. Note the column the paper heads “Density (kg/m^3)” is actually the element velocity (m/s): every entry equals V = \dot m/(\overline{\rho}\,A) to the printed digits, while the true density \rho = p/RT is different. We confirm that below.

# Table IV: element mass flow [g/s] (+ in the upstream->downstream sense).
PUB_MDOT_GPS = {
    1: 16.461,
    2: 8.033,
    3: 3.596,
    4: 3.375,
    5: 1.782,
    6: 1.593,
    7: 0.797,
    8: 0.797,
    9: -3.596,
    10: 4.437,
    11: -8.033,
    12: 16.461,
    13: 8.429,
    14: 4.135,
    15: -4.293,
    16: 2.800,
    17: 1.400,
    18: 1.400,
    19: -1.493,
    20: 2.986,
    21: 1.493,
    22: 2.800,
    23: 1.400,
    24: 4.293,
    25: 4.135,
    26: 8.429,
    27: 4.437,
    28: 3.816,
    29: 1.400,
}
# Table IV: the "Density" column, actually the element velocity [m/s].
PUB_VELOCITY_MPS = {
    1: 8.56,
    2: 7.19,
    3: 3.77,
    4: 10.07,
    5: 6.05,
    6: 5.37,
    7: 2.78,
    8: 2.78,
    9: -3.77,
    10: 13.14,
    11: -7.19,
    12: 8.56,
    13: 7.65,
    14: 12.47,
    15: -4.73,
    16: 8.77,
    17: 4.81,
    18: 4.81,
    19: -1.75,
    20: 9.60,
    21: 1.75,
    22: 8.77,
    23: 4.81,
    24: 4.73,
    25: 12.47,
    26: 7.65,
    27: 13.14,
    28: 11.72,
    29: 4.81,
}
# Table IV: node pressure [bar].
PUB_PRESSURE_BAR = {
    1: 6.0,
    2: 5.2151,
    3: 4.1131,
    4: 3.0,
    5: 3.8546,
    6: 3.2057,
    7: 3.0,
    8: 3.0423,
    9: 3.0,
    10: 3.0,
    11: 4.1131,
    12: 3.0,
    13: 5.2151,
    14: 6.0,
    15: 3.9848,
    16: 3.0,
    17: 3.5975,
    18: 3.1286,
    19: 3.0,
    20: 3.0,
    21: 3.5478,
    22: 3.0,
    23: 3.5975,
    24: 3.1286,
    25: 3.0,
    26: 3.9848,
    27: 3.0,
    28: 3.0,
    29: 3.0,
}

2. Building the Nefes network

Nefes models each element as a graph node, wired by directed (tail, head, area) edges. So the mapping is:

  • each paper pipe \rightarrow a two-port pipe(L, D, f) element (constant area \pi D^2/4);
  • each interior paper node \rightarrow a static-pressure junction (mass balance + equal static pressure across its ports);
  • each supply leaf \rightarrow a total_pressure_inlet at 6 bar;
  • each demand leaf \rightarrow a pressure_outlet at 3 bar.

A pipe gets two edges: one into it from the upstream node (port 0) and one out to the downstream node (port 1).

def build_network():
    net = nefes.Network(gas=nefes.perfect_gas(R_AIR, GAMMA))
    node_of = {}
    for n in INTERIOR:
        node_of[n] = net.add(cat.junction(name=f"n{n}"))
    for n in SUPPLY_NODES:
        node_of[n] = net.add(cat.total_pressure_inlet(P_SUPPLY, T_ISO, name=f"supply{n}"))
    for n in DEMAND_NODES:
        node_of[n] = net.add(cat.pressure_outlet(P_DEMAND, T_ISO, name=f"demand{n}"))

    pipe_edges = {}  # element -> (in_edge, out_edge, area)
    for k, (u, d, D, L) in PIPES.items():
        area = np.pi * D * D / 4.0
        pk = net.add(cat.pipe(L, D, FRICTION, name=f"p{k}"))
        e_in = net.connect(node_of[u], pk, area)  # port 0: inflow from upstream
        e_out = net.connect(pk, node_of[d], area)  # port 1: outflow to downstream
        pipe_edges[k] = (e_in, e_out, area)
    return net, node_of, pipe_edges


net, node_of, pipe_edges = build_network()
net
Network · 58 elements, 58 edges | perfect gas (R=287 J/kg/K, gamma=1.4) | p=101325 Pa  ·  T=300 K  ·  mdot=0.066 kg/s (auto)
Elements
# name type
0 n2 JunctionStaticP
1 n3 JunctionStaticP
2 n5 JunctionStaticP
3 n6 JunctionStaticP
4 n8 JunctionStaticP
5 n11 JunctionStaticP
6 n13 JunctionStaticP
7 n15 JunctionStaticP
8 n17 JunctionStaticP
9 n18 JunctionStaticP
10 n21 JunctionStaticP
11 n23 JunctionStaticP
12 n24 JunctionStaticP
13 n26 JunctionStaticP
14 supply1 TotalPressureInlet
15 supply14 TotalPressureInlet
16 demand4 PressureOutlet
17 demand7 PressureOutlet
18 demand9 PressureOutlet
19 demand10 PressureOutlet
... (38 more)
Edges
# connection area [m²] name
0 supply1 → p1 0.0002835 e0
1 p1 → n2 0.0002835 e1
2 n2 → p2 0.0001981 e2
3 p2 → n3 0.0001981 e3
4 n3 → p3 0.0001981 e4
5 p3 → n5 0.0001981 e5
6 n5 → p4 7.854e-05 e6
7 p4 → n6 7.854e-05 e7
8 n6 → p5 7.854e-05 e8
9 p5 → demand7 7.854e-05 e9
10 n6 → p6 7.854e-05 e10
11 p6 → n8 7.854e-05 e11
12 n8 → p7 7.854e-05 e12
13 p7 → demand9 7.854e-05 e13
14 n8 → p8 7.854e-05 e14
15 p8 → demand10 7.854e-05 e15
16 n5 → p9 0.0001981 e16
17 p9 → n11 0.0001981 e17
18 n11 → p10 7.854e-05 e18
19 p10 → demand12 7.854e-05 e19
... (38 more)

Topology

The 29-pipe layout mirrors the paper’s Fig. 9: two 6-bar supplies feeding a trunk, thirteen 3-bar demand branches, and a single closed loop (2 - 26 - 23 - 21 - 17 - 15 - 13 - 11 - 5 - 3 - 2) that makes this a genuine network rather than a tree: historically the hard case for pipe solvers. Each pipe is a two-port pipe element; each interior node is a static-pressure junction; supplies are total-pressure inlets, demands are static-pressure outlets.

Network topology
# Node coordinates (x, y) laid out to resemble the paper's Fig. 9.
# Reused by the solved-flow figure in section 4 to place the nodes.
POS = {
    1: (0, 2),
    2: (1, 2),
    3: (2, 2),
    5: (3, 2),
    11: (4, 2),
    13: (5, 2),
    14: (6, 2),
    26: (1, 0),
    23: (2, 0),
    21: (3, 0),
    17: (4, 0),
    15: (5, 0),
    4: (2, 3),
    6: (3, 3.2),
    28: (2.4, 2.9),
    7: (2.5, 4.3),
    8: (3.5, 4.3),
    9: (3.2, 5.3),
    10: (3.9, 5.3),
    12: (4, 3),
    27: (0.4, 0),
    24: (2, -1.2),
    25: (1.6, -2.3),
    29: (2.4, -2.3),
    22: (3, -1.2),
    18: (4, -1.2),
    19: (3.6, -2.3),
    20: (4.4, -2.3),
    16: (5.6, 0),
}

3. Solve

A pure pressure-driven mesh with 14 unknown junction pressures and no fixed mass flow anywhere: precisely the regime Nefes’s vanishing-friction homotopy is built for.

sol = net.solve(verbose=1)
assert sol.converged, (sol.residual_norm, sol.print_residuals())
print(f"\nconverged={sol.converged}  iterations={sol.iterations}  residual={sol.residual_norm:.2e}")
kappa=0.1   ->   6 iters, ||R_hat||=2.547e-11, converged=True
kappa=0.01  ->   4 iters, ||R_hat||=1.527e-14, converged=True
kappa=0     ->   2 iters, ||R_hat||=2.403e-12, converged=True

converged=True  iterations=12  residual=2.40e-12
# extract Nefes results, keyed by paper element / node number
mdot = sol.field("mdot")
p = sol.field("p")
rho = sol.field("rho")
mach = sol.field("M")
area = sol.field("area")

fns_mdot = {k: mdot[pipe_edges[k][1]] * 1000.0 for k in PIPES}  # g/s, up->down positive
fns_p = dict(PUB_PRESSURE_BAR)  # terminals are prescribed
for n in INTERIOR:
    e_out = pipe_edges[next(k for k, (u, d, _D, _L) in PIPES.items() if d == n)][1]
    fns_p[n] = p[e_out] / 1e5  # junction static pressure [bar]

max_m = max(abs(fns_mdot[k] - PUB_MDOT_GPS[k]) / abs(PUB_MDOT_GPS[k]) for k in PIPES)
max_p = max(abs(fns_p[n] - PUB_PRESSURE_BAR[n]) for n in INTERIOR)
print(f"max mass-flow error   : {max_m:.3%}")
print(f"max node-pressure error: {max_p*1000:.2f} mbar  ({max_p / np.mean(list(PUB_PRESSURE_BAR.values())):.3%})")
print(f"peak Mach number       : {np.max(np.abs(mach)):.4f}")
max mass-flow error   : 0.086%
max node-pressure error: 1.81 mbar  (0.050%)
peak Mach number       : 0.0457

4. The solved network

Node colour is the solved pressure; pipe thickness is |\dot m|. Hover for the numbers.

fig2 = go.Figure()
m_max = max(abs(v) for v in fns_mdot.values())
for k, (u, d, _D, _L) in PIPES.items():
    w = 1.5 + 7.0 * abs(fns_mdot[k]) / m_max
    fig2.add_trace(
        go.Scatter(
            x=[POS[u][0], POS[d][0]],
            y=[POS[u][1], POS[d][1]],
            mode="lines",
            line=dict(color="#b0b0b0", width=w),
            showlegend=False,
            hoverinfo="text",
            text=f"pipe {k}: {fns_mdot[k]:+.3f} g/s",
        )
    )
pv = [fns_p[n] for n in PAPER_NODES]
fig2.add_trace(
    go.Scatter(
        x=[POS[n][0] for n in PAPER_NODES],
        y=[POS[n][1] for n in PAPER_NODES],
        mode="markers+text",
        text=[str(n) for n in PAPER_NODES],
        textposition="top center",
        marker=dict(
            size=18, color=pv, colorscale="Turbo", colorbar=dict(title="p (bar)"), line=dict(width=1, color="white")
        ),
        customdata=pv,
        hovertemplate="node %{text}: %{customdata:.4f} bar<extra></extra>",
        showlegend=False,
    )
)
fig2.update_xaxes(visible=False)
fig2.update_yaxes(visible=False, scaleanchor="x", scaleratio=1)
fig2.update_layout(title="Solved mean flow: node pressure (colour), pipe |mass flow| (width)", height=560)
fig2.show()

5. Nefes vs. the paper: parity

Every point sits on y = x: the mass flows agree to well under a percent and the node pressures to a few millibar.

ks = sorted(PIPES)
ns = sorted(PAPER_NODES)

fig3 = go.Figure()
fig3.add_trace(
    go.Scatter(
        x=[PUB_MDOT_GPS[k] for k in ks],
        y=[fns_mdot[k] for k in ks],
        mode="markers",
        marker=dict(size=9, color=COLORWAY[0]),
        text=[f"pipe {k}" for k in ks],
        hovertemplate="%{text}<br>paper %{x:.3f}<br>Nefes %{y:.3f} g/s<extra></extra>",
        name="pipes",
    )
)
fig3.add_trace(go.Scatter(x=[-9, 18], y=[-9, 18], mode="lines", line=dict(dash="dash", color="#888"), name="y = x"))
fig3.update_layout(
    title="Mass flow: Nefes vs. Table IV", xaxis_title="paper (g/s)", yaxis_title="Nefes (g/s)", height=460
)
fig3.show()

fig4 = go.Figure()
fig4.add_trace(
    go.Scatter(
        x=[PUB_PRESSURE_BAR[n] for n in ns],
        y=[fns_p[n] for n in ns],
        mode="markers",
        marker=dict(size=9, color=COLORWAY[1]),
        text=[f"node {n}" for n in ns],
        hovertemplate="%{text}<br>paper %{x:.4f}<br>Nefes %{y:.4f} bar<extra></extra>",
        name="nodes",
    )
)
fig4.add_trace(go.Scatter(x=[2.9, 6.1], y=[2.9, 6.1], mode="lines", line=dict(dash="dash", color="#888"), name="y = x"))
fig4.update_layout(
    title="Node pressure: Nefes vs. Table IV", xaxis_title="paper (bar)", yaxis_title="Nefes (bar)", height=460
)
fig4.show()

Where the (tiny) differences come from

The residual is not solver error: it is the physics the paper’s low-Mach model drops, and it is in Nefes’s favour. The per-pipe mass-flow error is largest on the highest-flow pipes; the node-pressure error is largest at the fast supply junctions (nodes 2 and 13), where the dynamic head \tfrac12\rho u^2 that separates static from total pressure is biggest. Both are O(M^2) effects, landing in the 3rd-4th significant figure: consistent with the peak Mach number below.

fig5 = go.Figure()
fig5.add_trace(
    go.Bar(x=ks, y=[(fns_mdot[k] - PUB_MDOT_GPS[k]) / abs(PUB_MDOT_GPS[k]) * 100 for k in ks], marker_color=COLORWAY[0])
)
fig5.update_layout(
    title="Per-pipe mass-flow error vs. Table IV", xaxis_title="pipe", yaxis_title="error (%)", height=380
)
fig5.show()

fig6 = go.Figure()
fig6.add_trace(go.Bar(x=ks, y=[abs(mach[pipe_edges[k][1]]) for k in ks], marker_color=COLORWAY[3]))
fig6.update_layout(
    title="Pipe-exit Mach number (peak ~0.05: firmly low-Mach)", xaxis_title="pipe", yaxis_title="M", height=380
)
fig6.show()

6. The Table IV “Density” column is really velocity

The paper’s element table heads its third column “Density (kg/m^3)”, but the numbers are the element velocity V = \dot m/(\overline{\rho}\,A) evaluated at the mean density \overline{\rho} = \tfrac12(\rho_\mathrm{up}+\rho_\mathrm{down}): exactly the reference velocity Nefes’s pipe kernel forms internally. Below, Nefes’s mean-density velocity lands on the published column (y=x), while the true density \rho = p/RT is a different curve entirely: so the label is a misprint.

fns_V, fns_rho_bar = {}, {}
for k in PIPES:
    e_in, e_out, _a = pipe_edges[k]
    rho_bar = 0.5 * (rho[e_in] + rho[e_out])
    fns_rho_bar[k] = rho_bar
    fns_V[k] = mdot[e_out] / (rho_bar * area[e_out])  # signed, up->down positive

fig7 = go.Figure()
fig7.add_trace(
    go.Scatter(
        x=[PUB_VELOCITY_MPS[k] for k in ks],
        y=[fns_V[k] for k in ks],
        mode="markers",
        marker=dict(size=9, color=COLORWAY[2]),
        name="Nefes velocity  V = mdot/(rho_bar A)",
        text=[f"pipe {k}" for k in ks],
        hovertemplate="%{text}<br>column %{x:.2f}<br>Nefes V %{y:.2f} m/s<extra></extra>",
    )
)
fig7.add_trace(
    go.Scatter(
        x=[PUB_VELOCITY_MPS[k] for k in ks],
        y=[np.sign(PUB_VELOCITY_MPS[k]) * fns_rho_bar[k] for k in ks],
        mode="markers",
        marker=dict(size=9, color=COLORWAY[1], symbol="x"),
        name="true density  rho = p/RT",
        text=[f"pipe {k}" for k in ks],
        hovertemplate="%{text}<br>column %{x:.2f}<br>rho %{y:.2f}<extra></extra>",
    )
)
fig7.add_trace(go.Scatter(x=[-8, 14], y=[-8, 14], mode="lines", line=dict(dash="dash", color="#888"), name="y = x"))
fig7.update_layout(
    title='Paper\'s "Density" column vs. Nefes velocity and Nefes density',
    xaxis_title='Table IV "Density" column',
    yaxis_title="Nefes quantity",
    height=460,
)
fig7.show()

7. Summary

quantity agreement with Table IV
pipe mass flows (29) < 0.1\%
node pressures (14 interior) < 2\ \mathrm{mbar} (< 0.05\%)
convergence ~10 Newton iterations, residual \sim 10^{-15}
peak Mach \approx 0.05

Nefes reproduces the paper’s only compressible benchmark to the printed precision, straight out of the box, using the pipe element that is the paper’s Darcy-Weisbach abstraction. The sub-percent residual is the compressible/dynamic-pressure physics the paper’s low-Mach model neglects: physically real, in Nefes’s favour, and exactly the O(M^2) size the Mach numbers predict. Along the way the run confirms two things the paper’s tables obscure: node 28 is a thirteenth 3-bar demand (easy to miss), and the Table IV “Density” column is actually velocity.

The rigorous check lives in tests/test_greyvenstein_laurie.py.

Export for Nemo

The network is in scope as net (a Network) and its converged mean flow as sol (a Solution). Save either to a UI-readable YAML: sol.to_yaml(path) embeds the mean-flow result fields, net.to_yaml(path) writes the topology only: then open the file in Nemo.

import os, tempfile

_out = os.path.join(tempfile.mkdtemp(), "greyvenstein_laurie_network.yaml")
sol.to_yaml(_out)  # embeds the mean-flow results; use net.to_yaml(_out) for topology only
print("saved case:", _out)
saved case: /tmp/tmp4xonqjnh/greyvenstein_laurie_network.yaml
Back to top