The first process model

Recap

A brief description on how models are created is given in the Hello World paragraph of the getting started section. The example was

from simu import Model, NumericHandler, Quantity


class Square(Model):
    """A model of a square"""

    def interface(self):
        self.parameters.define("length", 10, "m")
        self.properties.declare("area", "m^2")

    def define(self):
        self.properties["area"] = self.parameters["length"] ** 2


def main():
    numeric = NumericHandler(Square.top())
    func = numeric.function
    print(func.arg_structure)
    print(func.result_structure)

    args = numeric.arguments
    print(args)

    args[NumericHandler.MODEL_PARAMS]["length"] = Quantity(20, "cm")
    result = func(args)

    print(f"{result[NumericHandler.MODEL_PROPS]['area']:.3fP~}")


if __name__ == '__main__':
    main()

To recap, the model declares an interface telling other (parent) models that it has a parameter called length and calculates a property called area. Then, the definition implements the relationship.

The attentive reader might at this point have realised that no thermodynamic model was inv olved. The entire example is somehow remote to classical process engineering unless one calculates the cross-section of a square duct. Let us get this rectified in this section, and bring in our ideal gas model from the previous section.

To do so, we first need to create a simu.MaterialDefinition object, and by this follow the proper way to build up a simulation, and in practice, you might soon build up a repository of materials required for your field of application.

Creating a material definition

We can reuse much of the previous code, resulting into a simu.ThermoFrame object for a pure methane ideal gas. Only now, we store configuration data in yaml files. Let us have one file with some chemical species and their formulae (species_db.yml):

Methane: CH4
Hydrogen: H2
Nitrogen: N2
Pentamminerhodium chloride: "[RhCl(NH3)5]:2+"

Further, we store the thermodynamic model structures in another file (thermo_model_structures.yml):

simple_ideal_gas:
    state: GibbsState
    contributions:
      - H0S0ReferenceState
      - LinearHeatCapacity
      - IdealMix
      - GibbsIdealGas
simple_ideal_solution:
    state: GibbsState
    contributions:
      - H0S0ReferenceState
      - LinearHeatCapacity
      - IdealMix
      - ConstantGibbsVolume

And finally, the actual thermodynamic parameters (ideal_gas_param.yml):

H0S0ReferenceState:
  T_ref: 25 degC
  dh_form:
    Methane: -74.87 kJ/mol
  p_ref: 1 bar,
  s_0:
    Methane: 188.66 J/K/mol
LinearHeatCapacity:
  cp_a:
    Methane: 35.69 J/K/mol
  cp_b:
    Methane: 50 mJ/K**2/mol

Let’s first import some classes and read in those files:

from pathlib import Path
from yaml import safe_load

from simu import (
    InitialState, MaterialDefinition, ThermoParameterStore,
    StringDictThermoSource, SpeciesDB)
from simu.app import RegThermoFactory

CURRENT_DIR = Path(__file__).parent

# load species database
with open(CURRENT_DIR  / "species_db.yml") as file:
    species = SpeciesDB(safe_load(file))

# load model structure database
with open(CURRENT_DIR  / "thermo_model_structures.yml") as file:
    model_structures = safe_load(file)

# load thermodynamic parameter database
with open(CURRENT_DIR  / "ideal_gas_param.yml") as file:
    parameter_source = StringDictThermoSource(safe_load(file))

Now, as before, we can create the factory and from there the frame for our thermodynamic model:

factory = RegThermoFactory()
frame = factory.create_frame(species.get_sub_db(["Methane"]),
                             model_structures["simple_ideal_gas"])

Next, we go for the material definition, requiring also the initial state and a simu.ThermoParameterStore object:

initial_state = InitialState.from_si(400, 2e5, [1.0])
store = ThermoParameterStore()
ch4_ideal = MaterialDefinition(frame, initial_state, store)

This parameter store can be shared among multiple – normally all – material definitions, and thus holds a global set of thermodynamic parameters. Multiple stores are only required if two materials containing the same chemical species need to receive distinct values for the same thermodynamic parameter.

So far, we did not provide the parameter values, but the material definition has already told the store which parameters are required. We can query the super-set of the parameter names and units of all missing parameters required:

>>> pprint(missing_symbols)
{'H0S0ReferenceState': {'T_ref': 'K ',
                        'dh_form': {'Methane': 'J / mol '},
                        'p_ref': 'Pa ',
                        's_0': {'Methane': 'J / K / mol '}},
 'LinearHeatCapacity': {'cp_a': {'Methane': 'J / K / mol '},
                        'cp_b': {'Methane': 'J / K ** 2 / mol '}}}

Finally, we provide the already read parameters to the store:

store.add_source("my_source", parameter_source)

This time, there are no more missing symbols, and the print statement prints an empty dictionary:

>>> pprint(store.get_missing_symbols())
{}

In real applications, storing the meta-data and parameters in yml files is not the most stupid idea, but you might connect to any other file format or database of your choice, as long as the source can provide the nested dictionary of properties as requested by the store.

Note

As multiple parameter sources can be stacked in one store, we recommend to assign one source per bibliographic source of parameters. The models can then easily be queried for the names of the used sources and by that keep these sources traceable.

Using a material in a model

This is the big moment, as we now can use the material definition in an actual process model. The above created simu.MaterialDefinition object can be global for the entire project along with all other material definitions that you might need.

The following model is a hello world example for using such material:

from simu import Model
from .material_factory import ch4_ideal_gas

class Source(Model):
    """A model of a methane source"""

    def interface(self):
        self.parameters.define("T", 25, "degC")
        self.parameters.define("p", 1, "bar")
        self.parameters.define("V", 10, "m^3/hr")

    def define(self):
        src = self.materials.create_flow("source", ch4_ideal_gas)
        self.residuals.add("T", self.parameters["T"] - src["T"], "K")
        self.residuals.add("p", self.parameters["p"] - src["p"], "bar")
        self.residuals.add("V", self.parameters["V"] - src["V"], "m^3/h")

Here we first define the three parameters T, p and V that determine our system. These three parameters also constitute the interface of our model. The definition creates a methane flow simu.Material object from our definition. Finally, lines 14-16 constrain the system to the direct specifications of the parameter variables.

Well, the above syntax is very verbose, but might get into the way with regards to coding efficiency and the ambitions to keep lines short and to the point. For this reason, we define a subclass to simu.Model, namely simu.AModel that does nothing but defining abbreviations. As such, we can reduce the above model to:

from simu import AModel
from .material_factory import ch4_ideal_gas

class Source(AModel):
    """A model of a methane source"""

    def interface(self):
        self.pad("T", 25, "degC")
        self.pad("p", 1, "bar")
        self.pad("V", 10, "m^3/hr")

    def define(self):
        src = self.mcf("source", ch4_ideal_gas)
        self.ra("T", self.pa["T"] - src["T"], "K")
        self.ra("p", self.pa["p"] - src["p"], "bar")
        self.ra("V", self.pa["V"] - src["V"], "m^3/h")

This is as much as we can do without entirely drowning out the pythonic way of coding.

Either way, here we are with a complete process model. By creating a simu.NumericHandler, we obtain a function object representing our model. The function argument and result is a nested structure of quantities:

>>> from simu import NumericHandler
>>> numeric = NumericHandler(Source.top())
>>> args = numeric.arguments
>>> pprint(args)
{'model_params': {'T': <Quantity(25.0, 'degree_Celsius')>,
                  'V': <Quantity(10.0, 'meter ** 3 / hour')>,
                  'p': <Quantity(1.0, 'bar')>},
 'thermo_params': {'default': {'H0S0ReferenceState': {'T_ref': <Quantity(25, 'degree_Celsius')>,
                                                      'dh_form': {'Methane': <Quantity(-74.87, 'kilojoule / mole')>},
                                                      'p_ref': <Quantity(1, 'bar')>,
                                                      's_0': {'Methane': <Quantity(188.66, 'joule / kelvin / mole')>}},
                               'LinearHeatCapacity': {'cp_a': {'Methane': <Quantity(35.69, 'joule / kelvin / mole')>},
                                                      'cp_b': {'Methane': <Quantity(50.0, 'millijoule / kelvin ** 2 / mole')>}}}},
 'vectors': {'states': <Quantity([400, 200000, 1], 'dimensionless')>}}

Firstly, we can recognize the model parameters, the thermodynamic parameters, and the thermodynamic state of our material. The latter is stored in a dimensionless vector for the purpose of numerical solving. Later-on, we show how this vector, and/or individual parameters can be substituted by CasADi symbols and thus become free variables in a calculation.

Further, we can query the result by calling the function with this argument:

>>> pprint(numeric.function(args))
{'residuals': {'T': <Quantity(-101.85..., 'kelvin')>,
               'V': <Quantity(-0.0138..., 'meter ** 3 / second')>,
               'p': <Quantity(-100000.0, 'pascal')>},
 'thermo_props': {'source': {'S': <Quantity(194.09666..., 'watt / kelvin')>,
                             'T': <Quantity(400.0, 'kelvin')>,
                             'T_ref': <Quantity(298.15, 'kelvin')>,
                             'V': <Quantity(0.01662..., 'meter ** 3 / second')>,
                             'mu': {'Methane': <Quantity(-148614.30..., 'joule / mole')>},
                             'n': {'Methane': <Quantity(1.0, 'mole / second')>},
                             'p': <Quantity(200000.0, 'pascal')>,
                             'p_ref': <Quantity(100000.0, 'pascal')>}},
 'vectors': {'bounds': <Quantity([2.e+05 1.e+00 4.e+02], 'dimensionless')>,
             'residuals': <Quantity([-1.018...e+09 -4.986...e+08 -1.00000000e+07], 'dimensionless')>}}

Here we see the residuals as physical quantities, but also converted to a dimensionless vector, representing the quotient of residuals and their tolerances. Thermodynamic properties are included, and model properties would, if there were any.

The volume is calculated to 59.86 m3/hr, but we specified 10 m3/hr, and also the pressure and temperature are not yet as desired. The specifications are only fulfilled once the residuals are brought down to values below their tolerances.

Summary / Outlook

• Based on the previously defined thermodynamic model, we created a material definition object.

• For good house-keeping, we can move most of the static configuration for instance into yaml files.

• Such material definition can be utilized in process models to initiate a state or flow of that material.

• Once instantiated in a model, the properties of the material can be used to derive new properties or to define model constraints (residuals).