Unit models and property package for natural gas pipelines (#667)

* function for tracking cost and tests

* __init__

* functions to load piecewise constant inputs into model and tests

* natural gas property package for pipeline

* compressor unit model and tests

* pipeline unit model and tests

* node unit model and tests

* test with pipeline and compressor

* __init__

* update import

* add IDAES boilerplate

* add IDAES boilerplate

* boilerplate

* remove commented code for calculating speed of sound in unit model

* remove redundant instantiation and method calls of IncidenceGraphInterface

* __init__ in gas_distribution/

* add dens_nominal parameter to ParameterBlock

* tests for some previously untested properties

* import ConfigurationError

* use compare_expressions to make sure tracking expressions are correct

* move paper citation to top of file

* typos

* remove todo

* add paper citation and equation number

* fix typo

* update comment

* paper citation and reference to it

* clarify comment

* clean up steady state pipeline tests

* remove commented pipeline.flow_mass variables

* remove comments

* remove unused dulmage mendelsohn calls

* remove commented pipeline.flow_mass

* remove two commented lines

* remove commented code

* remove comment

* deactivate momentum balance at corner condition for given discretization

* test dynami degrees of freedom with different discretizations

* remove now redundant deactivation of momentum balance constraint

* docstrings

* add equality constraints between node and supply temperature and pressure

* update tests for new constraints

* add comment about fixing temperature at supply node

* main runs unittest.main

* tests for some flowsheet configurations

* try to get zero dof

* tests constructing three different small flowsheets

* remove dynamic optimization example from pipeline unit/property branch

Co-authored-by: Ludovico Bianchi <[email protected]>
Co-authored-by: Andrew Lee <[email protected]>

idaes/gas_distribution/properties/tests/test_natural_gas.py

@@ -0,0 +1,160 @@
+#################################################################################
+# The Institute for the Design of Advanced Energy Systems Integrated Platform
+# Framework (IDAES IP) was produced under the DOE Institute for the
+# Design of Advanced Energy Systems (IDAES), and is copyright (c) 2018-2021
+# by the software owners: The Regents of the University of California, through
+# Lawrence Berkeley National Laboratory,  National Technology & Engineering
+# Solutions of Sandia, LLC, Carnegie Mellon University, West Virginia University
+# Research Corporation, et al.  All rights reserved.
+#
+# Please see the files COPYRIGHT.md and LICENSE.md for full copyright and
+# license information.
+#################################################################################
+import pyomo.common.unittest as unittest
+import pytest
+import pyomo.environ as pyo
+from pyomo.util.check_units import (
+    assert_units_equivalent,
+    assert_units_consistent,
+)
+from pyomo.util.calc_var_value import calculate_variable_from_constraint
+
+from idaes.gas_distribution.properties.natural_gas import (
+    NaturalGasParameterBlock,
+)
+
+"""
+Test for simple ideal methane property package created with generic
+property framework.
+"""
+
+@pytest.mark.unit
+class TestNaturalGasPropertyPackage(unittest.TestCase):
+
+    def test_build(self):
+        m = pyo.ConcreteModel()
+        m.properties = NaturalGasParameterBlock()
+        m.state = m.properties.build_state_block(
+            # As usual, if defined_state is True, we need to specify
+            # every mole fraction
+            default={"defined_state": False},
+        )
+        state = m.state
+
+        # The following tests are very specific to this property package.
+        # They test units and values, which may be different for other
+        # property packages that are valid for the pipeline models.
+
+        assert_units_equivalent(state.pressure.get_units(), pyo.units.bar)
+        assert_units_equivalent(
+            state.flow_mol.get_units(), pyo.units.kmol/pyo.units.hr
+        )
+        assert_units_equivalent(state.temperature.get_units(), pyo.units.K)
+        component_list = m.state.config.parameters.component_list
+        j = next(iter(component_list))
+        self.assertIs(state.mole_frac_comp[j].get_units(), None)
+
+        mw = m.state.mw
+        assert_units_equivalent(
+            pyo.units.get_units(mw), pyo.units.kg/pyo.units.kmol
+        )
+        self.assertEqual(pyo.value(mw), 18.0)
+
+        flow_mass = m.state.flow_mass
+        assert_units_equivalent(
+            pyo.units.get_units(flow_mass), pyo.units.kg/pyo.units.hr
+        )
+        m.state.flow_mol.set_value(100.0*pyo.units.kmol/pyo.units.hr)
+        self.assertAlmostEqual(pyo.value(flow_mass), 1800.0)
+
+        kJkmolK = pyo.units.kJ/pyo.units.kmol/pyo.units.K
+        cp_mol_comp = m.state.cp_mol_comp[j]
+        assert_units_equivalent(
+            pyo.units.get_units(cp_mol_comp), kJkmolK
+        )
+        self.assertAlmostEqual(pyo.value(cp_mol_comp), 2.34*18.0)
+
+        cp_mol = m.state.cp_mol
+        assert_units_equivalent(pyo.units.get_units(cp_mol), kJkmolK)
+        self.assertAlmostEqual(pyo.value(cp_mol), 2.34*18.0)
+
+        kJkgK = pyo.units.kJ/pyo.units.kg/pyo.units.K
+        cp_mass = m.state.cp_mass
+        assert_units_equivalent(cp_mass, kJkgK)
+        self.assertEqual(pyo.value(cp_mass), 2.34)
+
+        cv_mol_comp = m.state.cv_mol_comp[j]
+        assert_units_equivalent(pyo.units.get_units(cv_mol_comp), kJkmolK)
+        self.assertAlmostEqual(pyo.value(cv_mol_comp), 1.85*18.0)
+
+        cv_mol = m.state.cv_mol
+        assert_units_equivalent(pyo.units.get_units(cv_mol_comp), kJkmolK)
+        self.assertAlmostEqual(pyo.value(cv_mol), 1.85*18.0)
+
+        cv_mass = m.state.cv_mass
+        assert_units_equivalent(pyo.units.get_units(cv_mass), kJkgK)
+        self.assertAlmostEqual(pyo.value(cv_mass), 1.85)
+
+        heat_capacity_ratio = m.state.heat_capacity_ratio
+        assert_units_equivalent(
+            pyo.units.get_units(heat_capacity_ratio), pyo.units.dimensionless
+        )
+        self.assertAlmostEqual(pyo.value(heat_capacity_ratio), 2.34/1.85)
+
+        speed_of_sound = m.state.speed_of_sound
+        speed_of_sound_eq = m.state.speed_of_sound_eq
+        assert_units_equivalent(speed_of_sound, pyo.units.m/pyo.units.s)
+        assert_units_consistent(speed_of_sound_eq)
+        m.state.temperature.set_value(293.15*pyo.units.K)
+        calculate_variable_from_constraint(speed_of_sound, speed_of_sound_eq)
+        # 370.2 = (0.8*2.34/1.85*8314*293.15/18.0)**0.5
+        self.assertAlmostEqual(pyo.value(speed_of_sound), 370.1628645)
+
+        dens_mol = m.state.dens_mol
+        dens_mol_eq = m.state.dens_mol_eq
+        assert_units_equivalent(dens_mol, pyo.units.kmol/pyo.units.m**3)
+        assert_units_consistent(dens_mol_eq)
+        m.state.temperature.set_value(293.15*pyo.units.K)
+        m.state.pressure.set_value(50*pyo.units.bar)
+        calculate_variable_from_constraint(dens_mol, dens_mol_eq)
+        # 2.654 = 50*1e2/0.8/8.314/293.15
+        self.assertAlmostEqual(pyo.value(dens_mol), 2.5642238411)
+
+        component_list = m.state.config.parameters.component_list
+        j = next(iter(component_list))
+        dens_mol_comp = m.state.dens_mol_comp[j]
+        assert_units_equivalent(
+            pyo.units.get_units(dens_mol_comp),
+            pyo.units.kmol/pyo.units.m**3
+        )
+        self.assertAlmostEqual(pyo.value(dens_mol_comp), 2.5642238411)
+
+        # flow_mol was set to 100 kmol/hr earlier.
+        m.state.mole_frac_comp[j].set_value(1.0)
+        self.assertEqual(
+            pyo.value(m.state.flow_mol_comp[j]),
+            100.0,
+        )
+        assert_units_equivalent(
+            pyo.units.get_units(m.state.flow_mol_comp[j]),
+            pyo.units.kmol/pyo.units.hr,
+        )
+
+    def test_nominal_density(self):
+        m = pyo.ConcreteModel()
+        m.properties = NaturalGasParameterBlock()
+        self.assertEqual(m.properties.dens_nominal.value, 0.72)
+        assert_units_equivalent(
+            m.properties.dens_nominal.get_units(),
+            pyo.units.kg/pyo.units.m**3,
+        )
+
+    def test_compressibility(self):
+        m = pyo.ConcreteModel()
+        m.properties = NaturalGasParameterBlock()
+        m.state = m.properties.build_state_block()
+        self.assertEqual(m.state.compressibility.value, 0.80)
+
+
+if __name__ == "__main__":
+    unittest.main()

idaes/gas_distribution/unit_models/compressor.py

@@ -0,0 +1,176 @@
+#################################################################################
+# The Institute for the Design of Advanced Energy Systems Integrated Platform
+# Framework (IDAES IP) was produced under the DOE Institute for the
+# Design of Advanced Energy Systems (IDAES), and is copyright (c) 2018-2021
+# by the software owners: The Regents of the University of California, through
+# Lawrence Berkeley National Laboratory,  National Technology & Engineering
+# Solutions of Sandia, LLC, Carnegie Mellon University, West Virginia University
+# Research Corporation, et al.  All rights reserved.
+#
+# Please see the files COPYRIGHT.md and LICENSE.md for full copyright and
+# license information.
+#################################################################################
+from pyomo.common.config import ConfigValue
+from pyomo.core.base.constraint import Constraint
+from pyomo.core.base.var import Var
+from pyomo.core.base.units_container import units as pyunits
+
+from idaes.core.unit_model import UnitModelBlockData
+from idaes.core.process_block import declare_process_block_class
+from idaes.core.control_volume_base import MaterialBalanceType
+from idaes.core.util.config import is_physical_parameter_block
+
+"""
+A simple unit model to calculate power required to achieve isothermal
+gas compression. 
+
+Data sources:
+    [1] Stochastic Optimal Control Model for Natural Gas Network
+        Operations. V. Zavala, 2014, Comp. Chem. Eng.
+"""
+
+
+@declare_process_block_class("IsothermalCompressor")
+class IsothermalCompressorData(UnitModelBlockData):
+
+    CONFIG = UnitModelBlockData.CONFIG()
+
+    CONFIG.declare(
+        "property_package",
+        ConfigValue(default=None, domain=is_physical_parameter_block),
+    )
+
+    def build(self):
+        super(IsothermalCompressorData, self).build()
+
+        time = self.flowsheet().time
+        property_package = self.config.property_package
+
+        if len(property_package.phase_list) != 1:
+            raise ValueError(
+                "%s can only be constructed with a "
+                "single-phase property package.\n"
+                "Got phases %s."
+                % (self.__class__, [p for p in property_package.phase_list])
+            )
+        self.phase = next(iter(property_package.phase_list))
+        if self.phase != "Vap":
+            raise ValueError(
+                "%s can only be constructed with a single phase, \"Vap\"."
+                "Got phase %s."
+                % (self.__class__, self.phase)
+            )
+
+        inlet_config = {"defined_state": True}
+        self.inlet_state = property_package.build_state_block(
+            time,
+            default=inlet_config,
+            )
+        self.outlet_state = property_package.build_state_block(time)
+
+        # A little annoying that add_port assumes the state block is indexed
+        # by exactly time. If we supported unindexed state blocks here,
+        # this entire unit model could be agnostic of time.
+        # (This is not strictly true because add_state_material_balances
+        # writes time-indexed balance equations, but we could get around this
+        # by writing our own balance equations, which is easy enough.)
+        self.add_port(
+            name="inlet_port",
+            block=self.inlet_state,
+            doc="The inlet to the compressor",
+        )
+        self.add_port(
+            name="outlet_port",
+            block=self.outlet_state,
+            doc="The outlet from the compressor",
+        )
+
+        self.add_state_material_balances(
+            # TODO: Does componentPhase vs. componentTotal matter here?
+            MaterialBalanceType.componentPhase,
+            self.inlet_state,
+            self.outlet_state,
+        )
+
+        self.add_state_isothermal_equation(self.inlet_state, self.outlet_state)
+
+        # Add boost pressure and pressure equation
+        self.boost_pressure = Var(
+            time,
+            initialize=0.0,
+            units=pyunits.bar,
+            # TODO: Should this variable not be initialized to its bound?
+            bounds=(0.0, None),
+            doc="Increase in pressure provided by this compressor",
+        )
+        self.add_pressure_change_equation(self.inlet_state, self.outlet_state)
+
+        # Add compression coefficient 
+        self.beta = Var(
+            time,
+            initialize=0.28, # Approximately (1.4 - 1.0)/1.4
+            units=pyunits.dimensionless,
+            bounds=(0.0, None),
+            doc="Compression coefficient, beta.",
+        )
+        self.add_beta_equation(self.inlet_state)
+
+        # Add power requirement
+        self.power = Var(
+            time,
+            initialize=0.0,
+            units=pyunits.kW,
+            bounds=(0.0, None),
+            doc="Work required to achieve the desired compression",
+        )
+        self.add_power_equation(self.inlet_state)
+
+    def add_state_isothermal_equation(self, state1, state2):
+        time = self.flowsheet().time
+        # NOTE: We assume our states are only indexed by time and that
+        # temperature is unindexed.
+        def isothermal_rule(b, t):
+            return state1[t].temperature == state2[t].temperature
+        self.state_isothermal_eqn = Constraint(time, rule=isothermal_rule)
+
+    def add_pressure_change_equation(self, inlet, outlet):
+        time = self.flowsheet().time
+        def pressure_rule(b, t):
+            # Here we are using some quantities from the state blocks
+            # and some from the unit model, so we must convert units.
+            return (
+                pyunits.convert(inlet[t].pressure, pyunits.bar)
+                + pyunits.convert(self.boost_pressure[t], pyunits.bar)
+                == pyunits.convert(outlet[t].pressure, pyunits.bar)
+            )
+        self.pressure_change_eqn = Constraint(time, rule=pressure_rule)
+
+    def add_beta_equation(self, inlet_state):
+        time = self.flowsheet().time
+        p = next(iter(self.config.property_package.phase_list))
+        def beta_rule(b, t):
+            # We know that we are constructing this property package with a
+            # single phase, but we use heat_capacity_ratio_phase nonetheless
+            # for compatibility with "generic" property packages.
+            gamma = inlet_state[t].heat_capacity_ratio_phase[p]
+            return b.beta[t] == (gamma - 1.0) / gamma
+        self.beta_eqn = Constraint(time, rule=beta_rule)
+
+    def add_power_equation(self, inlet_state):
+        """
+        This is Equation 2.10 in [1]
+        """
+        time = self.flowsheet().time
+        def power_rule(b, t):
+            cp_mass = inlet_state[t].cp_mol / inlet_state[t].mw
+            power_expr = (
+                inlet_state[t].temperature
+                * cp_mass
+                * inlet_state[t].flow_mass
+                * (((pyunits.convert(inlet_state[t].pressure, pyunits.bar)
+                    + pyunits.convert(self.boost_pressure[t], pyunits.bar))
+                    / pyunits.convert(inlet_state[t].pressure, pyunits.bar)
+                )**b.beta[t] - 1.0)
+            )
+            return b.power[t] == pyunits.convert(power_expr, pyunits.kW)
+        self.power_eqn = Constraint(time, rule=power_rule)