Modelling a small steam generation with a turbine

The process to model

Let us say, we have 10 MW of heat available at a sufficiently high temperature level to produce high pressure steam. The task is to sketch a steam boiler system consisting of a boiler feed water source, a steam drum, a boiler and a superheater, as well as a condensing turbine with downstream total condensation. One interesting question is, how much mechanical power can be harvested for a given steam pressure, turbine efficiency, and achievable condenser temperature.

../_images/steam_system.svg

Stream c1 is slightly subcooled boiler feed water, mixed with saturated steam and water in the steam drum, from which a condensate flow c2 is extracted to the boiler and partially evaporated into return flow c3/s3. A small bleed, approximately 1 % of the boiler feed water flow, is removed as stream c4 to avoid accumulation of impurities that otherwise would potentially scale or corrode the boiler coils. The produced saturated steam exits as stream s5 into the steam superheater, continuing as superheated stream s6 into the turbine. The turbine causes partial condensation into stream c7/s7, before the entire stream is condensed into the saturated condensate stream s8.

We could, but do not consider the return of s8 towards the boiler feed water, as this normally involves buffer tanks and initial heat exchangers, exploiting low level waste heat to raise the BFW temperature.

May the system be specified by

  • constant and specified pressure level of 100 bar upstream of the turbine

  • no pressure drop in condenser

  • given BFW feed temperature of 250 °C

  • The blowdown c4 is 1% of the feed flow c1

  • The boiler circulation yields 12 % evaporation in stream c3/s3

  • Total duty in boiler and superheater is 10 MW

  • The turbine has an isentropic efficiency of 80 %

  • There is 10 % condensation in the turbine

  • Given temperature of 50 °C at condenser outlet c8

As such, we expect the duty distribution to yield the desired condensation in the turbine, and the boiler duty to determine the BFW flow. The condenser temperature determines the turbine discharge pressure. This is hence a typical process that is still simple (and simplified), but includes a large recycle and multiple implicit specifications.

Thermodynamic models definition

The process includes steam and condensate, for which MaterialDefinition objects need to be created. For this purpose, we first create the ThermoFrame objects as follows:

 1from simu import (InitialState, MaterialDefinition,
 2                  SpeciesDefinition, ThermoFrame)
 3from simu.app import RegThermoFactory, ThermoStructure, predefined_parameters
 4
 5
 6def _create_frames() -> (ThermoFrame, ThermoFrame):
 7    factory = RegThermoFactory()
 8    species_def = {"H2O": SpeciesDefinition("H2O")}
 9
10    struct = ThermoStructure.from_predefined("IAPWS-Liquid")
11    cond = factory.create_frame(species_def, struct + "GenericProperties")
12
13    struct = ThermoStructure.from_predefined("IAPWS-Gas")
14    steam = factory.create_frame(species_def, struct + "GenericProperties")
15    return cond, steam

Here, we make use of some convenience functionality:

  • The ThermoStructure facilitates handling of thermodynamic model structures, but also offers static methods to create pre-defined ones, such as for the IAPWS equation of state.

  • Likewise, predefined_parameters is a ThermoParameterStore that already includes the parameters that ship with SiMu.

  • A group of ThermoContribution classes may be called Augmenters, defining sets of add-on properties. A pre-defined one is GenericProperties, calculating among other properties enthalpy and mass flow.

In below code, the function to create the frames is called in line 23.

17                     t: float, p: float, n: float) -> MaterialDefinition:
18    initial_state = InitialState.from_si(t, p, [n])
19    return MaterialDefinition(frame, initial_state, predefined_parameters)
20
21
22_f_cond, _f_steam = _create_frames()
23
24hp_steam = _create_material(_f_steam, 600, 100e5, 100)
25hp_condensate = _create_material(_f_cond, 600, 100e5, 100)
26lp_steam = _create_material(_f_steam, 373, 1e5, 100)
27lp_condensate = _create_material(_f_cond, 373, 1e5, 100)

The objects are then used as arguments to the _create_material function that attaches the model parameters and individual initial states to form the MaterialDefinition objects. As we deal with HP steam at 100 bar as well as condensation at subatmospheric pressure, it is advisable to either define individual MaterialDefinition objects or provide better initialisation in other ways.

Important

The objects created from line 25 downwards can live in the global space and be imported as needed when defining the process (parts). Alternatively, we could have stored them in a data structure, but the important aspect is that no copies shall or need to be created for each application.

This is not only due to performance reasons, but also to share a common database of thermodynamic parameters, which in other cases is essential to perform consistent fits of thermodynamic parameters.

Process model development

When building up a process model in \(\color{#0044aa}\mathtt{SigmaMu}\), a good approach is to think hierarchically and in terms of reusable and encapsulated model parts. There is no considerable overhead in doing so, as the numerical process model representation is not impacted by the chosen level of granularity.

A model part does not need to represent an entire unit operation, but it can also represent a reusable aspect of it. The hierarchical modelling approach greatly reduces the development effort when scaling up the scope of a process model.

Gas liquid equilibrium

The first module to introduce is one to constrain a steam and a condensate phase to be at equilibrium.

 1from simu import AModel
 2
 3
 4class VLE(AModel):
 5    def interface(self):
 6        self.md("gas")
 7        self.md("liquid")
 8
 9    def define(self):
10        g, l = self.m["gas"], self.m["liquid"]
11        self.ra("T_eq", l["T"] - g["T"], "K")
12        self.ra("p_eq", l["p"] - g["p"], "bar")
13
14        for s in g["mu"].keys() & l["mu"].keys():
15            self.ra(f"mu_eq_{s}", g["mu"][s] - l["mu"][s], "kJ/mol")

As the interface, the module expects two material connections, namely gas and liquid, both of them being arbitrary materials. In the definition part, the model defines \(2+n_c\) constraints, where \(n_c\) is the number of species being in common in both phases, in our case one, being \(\mathrm{H_2O}\).

\begin{align*} T_l &= T_g\\ p_l &= p_g\\ \mu_{l,i} &= \mu_{g,i}\quad \text{for}\ i \in \mathcal S_l \cap \mathcal S_g \end{align*}

Material balance

In an equally generic manner, we can define a node that assures species conservation over a node with variable number of inlets and outlets.

As for the equilibrium module, this is an example of generic validity, and, as we will see, does not require equal species sets or any other constraints on the connected materials. The following constructor does not contain \(\color{#0044aa}\mathtt{SigmaMu}\) specific code:

1from simu import AModel
2
3
4class SpeciesBalance(AModel):
5    def __init__(self, num_in: int = 1, num_out: int = 1, tol_unit = "kmol/h"):
6        self._num_in, self._num_out = num_in, num_out
7        self._tol_unit = tol_unit
8        super().__init__()

All it does is to store attributes for later use, namely the number of inlet and outlet ports to be defined, as well as the tolerance unit for the molar flows. Regarding the latter, kmol/h is a suitable choice for industrially sized processes, but for modelling of laboratory-scale equipment, mol/h or mmol/s might be a better choice.

The following generators are used to produce the port names as in_01, in_02 etc. and likewise out_01, out_02 etc:

24    def in_names(self):
25        for i in range(self._num_in):
26            yield f"in_{i + 1:02d}"
27
28    def out_names(self):
29        for i in range(self._num_out):
30            yield f"out_{i + 1:02d}"

Hence nothing magic is going on here either. Next, we define the interface of the model:

10    def interface(self):
11        for name in self.in_names():
12            self.md(name)
13        for name in self.out_names():
14            self.md(name)

This interface simply defines the ports for all input and output ports. Note again that \(\color{#0044aa}\mathtt{SigmaMu}\) does not support the notion of port direction inherently, as it often does not apply, but is rather an implementation detail of the individual models.

Keeping the best for last, the actual definition of the model:

16    def define(self):
17        inlets = [self.m[name] for name in self.in_names()]
18        outlets = [self.m[name] for name in self.out_names()]
19
20        d_n = sum(m["n"] for m in inlets) -  sum(m["n"] for m in outlets)
21        for species, dn_i in d_n.items():
22            self.ra(species, dn_i, self._tol_unit)

While line 17 and 18 simply collect the connected material objects, the actual calculation happens in line 20: The mole vectors m["n"] are of type QuantityDict and thus define the + operator and can be used as arguments to the sum function. When adding two QuantityDict objects, non-existing keys are treated as zero, such that the result object contains the union of both species sets. d_n is therefore the residual of molar flows for all involved species.

Lines 21-22 then define the constraints to force these residuals to zero, applying the before mentioned tolerance unit.

Boiler

At this point, we can start to define simple unit operations, one of them being the boiler. For simplicity, the model is implemented specifically for the intended purpose without generalisation. As such, we specify zero pressure drop and the fraction of water being evaporated, while the duty is not locked, but only reported as a property. The class head and interface then looks as follows:

 1from simu import AModel
 2
 3from models.species_balance import SpeciesBalance
 4from models.vle import VLE
 5
 6
 7class Boiler(AModel):
 8    def interface(self):
 9        self.md("feed")
10        self.md("steam_out")
11        self.md("cond_out")
12
13        self.pad("vap_frac", 12, "%")
14        self.prd("duty", "MW")

Here, the (condensate) feed and steam and condensate outlet are defined as material ports. Further, the parameter vap_frac is declared to represent the vapour fraction, and the duty parameter to export the unit’s heat duty. The model definition is

16    def define(self):
17        f, c, s = self.m["feed"], self.m["cond_out"], self.m["steam_out"]
18        self.ra("p_f", f["p"] - s["p"], "bar")  # zero pressure drop
19
20        with self.ha("n_balance", SpeciesBalance, num_out=2) as unit:
21            unit.mcm(in_01=f, out_01=s, out_02=c)
22
23        with self.ha("vle", VLE) as unit:
24            unit.mcm(gas=s, liquid=c)
25
26        # vapour fraction from boiler
27        self.ra("vap_frac", self.pa["vap_frac"] * f["N"] - s["N"], "kmol/h")
28        self.pr["duty"] = c["H"] + s["H"] - f["H"]

Here, line 18 constraints the steam outlet pressure to be the steam pressure. No pressure constraint is explicitly given for the condensate outlet, as this is covered by the VLE sub-model. First the above developed MaterialBalance model is utilised in lines 20-21 to conserve the total amount of water, then the VLE model, included in lines 23-24, assures phase equilibrium.

Line 27 specifies the vapour fraction, whereas a little trick is applied. Two alternatives to define these residuals are

\[\frac{\dot n_s}{\dot n_f} - \beta^\mathrm{spec} = 0 \quad\text{or}\quad \dot n_s - \beta^\mathrm{spec}\cdot \dot n_f = 0\]

While the first form may be more intuitive, the second form is linear in \(\dot n_i\), which might have a slight positive effect on solver convergence, and is as such a good habbit. Line 28 defines the duty as the difference in enthalpies.

Superheater

As the superheater only deals with one-phase flow, its definition is ridiculously simple and needs no further explanation beyond what is already valid for the boiler above:

31class SuperHeater(AModel):
32    def interface(self):
33        self.md("inlet")
34        self.md("outlet")
35        self.prd("duty", "MW")
36
37    def define(self):
38        i, o = self.m["inlet"], self.m["outlet"]
39        self.ra("p", i["p"] - o["p"], "bar")  # zero pressure drop
40
41        with self.ha("n-balance", SpeciesBalance) as unit:
42            unit.mcm(in_01=i, out_01=o)
43
44        self.pr["duty"] = o["H"] - i["H"]

Total Condenser

While the superheater was simple, we might learn a thing about thermodynamics on this one, as the task is to receive a partially condensed flow, and bring it down to its exact boiling point. We could be pragmatic and define steam and condensate flows as both inlets and outlets, and so specify the vapour fraction to be small. If there were inert gases present and vacuum ejectors installed, the approach were entirely suitable, but here we stubbornly obliged ourselves to implement total condensation.

To do this, we need a trial gas phase that does not enter the material balance, but stands in phase equilibrium with the actual outlet stream. The size of the trial gas phase is arbitrary, represented by the _trial_phase_size parameter in below definition:

47class TotalCondenser(AModel):
48    def interface(self):
49        self.md("steam_in")
50        self.md("cond_in")
51        self.md("outlet")
52        self.pad("temperature", 50, "degC")
53        self.prd("duty", "MW")
54
55        self.pad("_trial_phase_size", 1, "mol/s")

Further, a parameter is defined for the desired temperature, which will via the boiling point fix the pressure level. We also announce to calculate the duty of the unit.

57    def define(self):
58        i_s, i_c, o = self.m["steam_in"], self.m["cond_in"], self.m["outlet"]
59        t = self.mcf("steam_trial", i_s.definition)
60
61        with self.ha("n-balance", SpeciesBalance, num_in=2) as unit:
62            unit.mcm(in_01=i_s, in_02=i_c, out_01=o)
63
64        with self.ha("vle", VLE) as unit:
65            unit.mcm(gas=t, liquid=o)
66        self.ra("trial_n", t["N"] - self.pa["_trial_phase_size"], "mol/s")
67
68        self.ra("pressure", i_s["p"] - o["p"], "bar")
69        self.ra("temperature", o["T"] - self.pa["temperature"], "K")
70
71        self.pr["duty"] = i_s["H"] + i_c["H"] - o["H"]

Here, line 59 creates the trial gas phase. Instead of giving a concrete MaterialDefinition object from among the ones defined in above section on thermodynamic models, we just use the definition given for the inlet steam, thereby making our model capable of dealing with other feeds, such as steam with minor amounts of inert gases.

In line 64-66, the outlet steam is brought to equilibrium with our trial phase, and the trial phase size is set. This will provide in total 4 equations, one more than the degrees of freedom introduced with the trial phase. The result is the link between temperature and pressure by the boiling point.

The remaining part contains nothing new. The temperature is set to the specified value, and pressure drop is defined as zero. This means however that the upstream unit (the turbine) cannot also specify its discharge pressure, as it is now fixed to the boiling point at 50 °C.

Steam Turbine

This is the last thermodynamic lecture for today – covering entropy and isentropic efficiency. The isentropic expansion is the process described by reducing the pressure of the flow at constant entropy:

\[\dot H^\mathrm{out}_\mathrm{rev} = \dot H(S^\mathrm{in}, p^\mathrm{out}, \dot{\vec n})\]

Thanks to the canonical formulation of our thermodynamic models, there is not more to it, and we do not need to start juggling adiabatic exponents and ideal gas approximations. The reversible (maximal reachable) work extraction is then simply

\[\dot W_\mathrm{rev} = \dot H^\mathrm{in} - \dot H^\mathrm{out}_\mathrm{rev}\]

From here, the actual outlet enthalpy computes by removing the actual work, defined by the isentropic efficiency \(\eta_s\), from the inlet enthalpy:

\[\dot H^\mathrm{out} = H^\mathrm{in} - \dot W \quad\text{with}\quad \dot W = \eta_s\,\dot W_\mathrm{rev}\]

The only complication here is the condensing operation, requiring to calculate all output enthalpies (reversible and actual) as the sum of liquid and gas enthalpy. Implementation-wise, we can first implement a model part that partially condensates a gaseous feed:

 1from simu import AModel
 2from models.vle import VLE
 3from models.species_balance import SpeciesBalance
 4
 5
 6class PartialCondensation(AModel):
 7    def interface(self):
 8        self.md("feed")
 9        self.md("steam_out")
10        self.md("cond_out")
11        self.prd("duty", "MW")
12        self.prd("vapour_fraction", "%")
13
14    def define(self):
15        f = self.m["feed"]
16        s = self.m["steam_out"]
17        c = self.m["cond_out"]
18
19        with self.ha("vle", VLE) as unit:
20            unit.mcm(gas=s, liquid=c)
21        with self.ha("n-balance", SpeciesBalance, num_out=2) as unit:
22            unit.mcm(in_01=f, out_01=s, out_02=c)
23
24        self.pr["duty"] = f["H"] - c["H"] - s["H"]
25        self.pr["vapour_fraction"] = s["N"] / f["N"]

The model PartialCondensation just establishes the species balance and phase equilibrium in the outlet. No calorimetric or pressure-related constraints are yet defined, but for convenience, we calculate the duty and the outlet vapour fraction. These are to be further exported for the actual process, but for curiosity, we can for instance monitor the vapour fraction of the isentropic expansion.

We can now use this model twice in the turbine model: once for the isentropic expansion, and once for the real expansion. For this, the interface is defined as follows:

28class Turbine(AModel):
29    def interface(self):
30        self.md("inlet")
31        self.md("steam_out")
32        self.md("cond_out")
33        self.pad("isentropic_efficiency", 80, "%")
34        self.pad("vapour_fraction", 90, "%")
35        self.prd("duty", "MW")

Here, the isentropic efficiency and the desired vapour fraction are defined as input parameters, while duty is further exported to be constrained in the context of the overall process.

In the first part of the define method, the streams are referenced and defined:

37    def define(self):
38        f = self.m["inlet"]
39        o_c = self.m["cond_out"]
40        o_s = self.m["steam_out"]
41        i_c = self.mcf("cond_iso", o_c.definition)
42        i_s = self.mcf("steam_iso", o_s.definition)

Like for the total condenser, we create materials locally (line 41, 42), here to represent the isentropic expansion that encapsulated within the turbine and not visible to the parent model. The material definitions are obtained from the materials connected to the actual output streams.

The next part instantiates our PartialCondensation module for the real and the isentropic process:

44        with self.ha("isentropic", PartialCondensation) as unit:
45            unit.mcm(feed=f, steam_out=i_s, cond_out=i_c)
46        iso_duty = unit.pr["duty"]
47        with self.ha("actual", PartialCondensation) as unit:
48            unit.mcm(feed=f, steam_out=o_s, cond_out=o_c)
49        duty, v_frac = unit.pr["duty"], unit.pr["vapour_fraction"]

Here, we store the duties and the actual vapour fraction in local variables. The model is completed by defining the missing constraints:

51        # set vapour fraction
52        self.ra("vapour_frac", v_frac - self.pa["vapour_fraction"], "%")
53
54        # constraints on expansion
55        self.ra("discharge_pressure", o_s["p"] - i_s["p"], "bar")
56        self.ra("entropy_balance", f["S"] - i_s["S"] - i_c["S"], "kW/K")
57
58        # duties and efficiencies
59        eta = self.pa["isentropic_efficiency"]
60        self.ra("duty", duty - eta * iso_duty, "MW")
61        self.pr["duty"] = duty
  • Line 52 sets the vapour fraction in the outlet.

  • Line 55 constraints that both the isentropic and the real process expand to the same pressure.

  • Line 56 conserves the entropy for the isentropic process, and

  • Line 60 specifies the duty, exported as a parameter in line 61, of the real process to \(\eta_s\,\dot W_\mathrm{rev}\).

Steam drum

The steam drum is modelled as an adiabatic vessel, and its complexity is merely given by the number of connected streams. There is the fresh boiler feed water inlet bfw of under-saturated liquid, a blow-down outlet of saturated condensate blow_down, and the outlet of saturated steam steam.

Additionally, saturated condensate boiler_feed exits to the boiler(s) and returns as inputs into the drum in two phases as boiler_return_cond and boiler_return_steam.

We also define the drum pressure and the amount of blow-down as parameters. The interface and the first part of the definition looks as follows:

 1from simu import AModel
 2from models.species_balance import SpeciesBalance
 3from models.vle import VLE
 4
 5
 6class SteamDrum(AModel):
 7    def interface(self):
 8        self.md("bfw")
 9        self.md("blow_down")
10        self.md("boiler_feed")
11        self.md("boiler_return_cond")
12        self.md("boiler_return_steam")
13        self.md("steam")
14
15        self.pad("pressure", 100, "bar")
16        self.pad("blow_down", 1, "%")
17
18
19    def define(self):
20        bfw = self.m["bfw"]
21        bd = self.m["blow_down"]
22        bf = self.m["boiler_feed"]
23        bc = self.m["boiler_return_cond"]
24        bs = self.m["boiler_return_steam"]
25        s = self.m["steam"]
26
27        # pressure and blowdown flow spec
28        self.ra("p_steam", self.pa["pressure"] - s["p"], "bar")
29        self.ra("blow_down", self.pa["blow_down"] * bfw["N"] - bd["N"], "mol/s")
30
31        # species balances
32        with self.ha("n_balance", SpeciesBalance, num_in=3, num_out=3) as unit:
33            unit.mcm(in_01=bfw, in_02=bc, in_03=bs,
34                     out_01=bd, out_02=bf, out_03=s)
35
36        # enthalpy balance
37        res = bfw["H"] + bc["H"] + bs["H"] - bf["H"] - bd["H"] - s["H"]
38        self.ra("h_balance", res, "MW")

Here, lines 28, 29 specify the pressure and the blow-down flow relative to the total BFW flow. Lines 32-34 define the species balance, and lines 37-38 the enthalpy balance. The remaining constrains are to set the steam in equilibrium with both condensate flows- the blow-down and the circulation flow. We can easily reuse our VLE module for one of them, e.g. the circulation flow:

40        # phase equilibria
41        with self.ha("vle", VLE) as unit:
42            unit.mcm(gas=s, liquid=bf)

Next, we meet however a thermodynamic challenge. Applying the same construct to the blow-down would re-impose the constraint that temperature is given by equilibrium for a set pressure. There are (for later) thermodynamic sound ways to escape this dilemma, but for the sake of simplicity, we assume for now a pure component flow, for which it is valid to equalise simply temperature and pressure:

44        self.ra("t_bd", bd["T"] - s["T"], "K")
45        self.ra("p_bd", bd["p"] - s["p"], "bar")

If we were dealing with a component mixture, additional constraints would be required for all but one species, and there are several ways to do this.

Overall process model

The overall model of a larger process might well be further divided into process sections. A typical feature is however the occurance of more material (streams). While we assigned the materials to individual local variables in the sub-models, at some point one can with advantage chose a different approach:

 1from simu import AModel
 2
 3from thermo import hp_condensate, hp_steam, lp_condensate, lp_steam
 4from models.heat_exchanger import Boiler, SuperHeater, TotalCondenser
 5from models.drum import SteamDrum
 6from models.turbine import Turbine
 7
 8
 9class SteamGeneration(AModel):
10    def interface(self):
11        self.pad("T_bfw", 250, "degC")
12        self.pad("blow_down", 1, "%")
13        self.pad("total_duty", 10, "MW")
14
15    def define(self):
16        # [(definition, prefix, [<stream numbers>], ...]
17        streams = [(hp_steam, "s", [3, 5, 6]),
18                   (hp_condensate, "c", [1, 2, 3, 4]),
19                   (lp_steam, "s", [7]),
20                   (lp_condensate, "c", [7, 8])]
21
22        for definition, prefix, numbers in streams:
23            for i in numbers:
24                self.mcf(f"{prefix}{i}", definition)
25        m = self.materials

For convenience, let us revisit the process flowsheet

../_images/steam_system.svg

Here, we also import for the first time the material definitions that we initially defined. The interface now cannot define any material ports, as these would have to be connected and thus require another parent model. It really does make sense to define the battery limits still as ports (c1, c4 and c8) and wrap a slim layer around. In our case, the outer layer would only need to specify temperature and pressure of the boiler feed water, as the flow demand is already determined by the process. However, we leave this generalization as an exercise to the attentive user.

Lines 17-20 define a data structure on which materials to instantiate, while lines 22-24 perform the instantiation. Line 25 simply shortcuts self.materials to m, so we can from there conveniently address materials by for instance m["c1"].

The following code mostly inserts and connects the previously developed sub-models:

27        with self.ha("boiler", Boiler) as unit:
28            unit.mcm(feed=m["c2"], cond_out=m["c3"], steam_out=m["s3"])
29        q_boiler = unit.pr["duty"]
30
31        with self.ha("drum", SteamDrum) as unit:
32            unit.mcm(bfw=m["c1"], blow_down=m["c4"], steam=m["s5"],
33                     boiler_feed=m["c2"], boiler_return_cond=m["c3"],
34                     boiler_return_steam=m["s3"])
35
36        with self.ha("superheater", SuperHeater) as unit:
37            unit.mcm(inlet=m["s5"], outlet=m["s6"])
38        q_superheater = unit.pr["duty"]
39
40        with self.ha("turbine", Turbine) as unit:
41            unit.mcm(inlet=m["s6"], steam_out=m["s7"], cond_out=m["c7"])
42
43        with self.ha("condenser", TotalCondenser) as unit:
44            unit.mcm(steam_in=m["s7"], cond_in=m["c7"], outlet=m["c8"])
45
46        # distribute duty for given vapour fraction
47        res = q_superheater + q_boiler - self.pa["total_duty"]
48        self.ra("total_duty", res, "MW")

Note that we rescure the duties of boiler and super-heater in lines 29 and 38 for the given specification of the total duty (line 47-48).

The remaining code is, as the necessity having been mentioned before, specifying the BFW inlet conditions:

51        self.ra("T_bfw", self.pa["T_bfw"] - m["c1"]["T"], "K")
52        self.ra("p_c1", m["s5"]["p"] - m["c1"]["p"], "bar")

Running the model

The code to run the model itself is rather generic and as follows:

 1from pprint import pprint
 2from simu import NumericHandler, SimulationSolver, Quantity
 3from process import SteamGeneration
 4
 5
 6def main():
 7    numeric = NumericHandler(SteamGeneration.top(), port_properties=False)
 8    solver = SimulationSolver(numeric)
 9    result = solver.solve()
10    print_result(result)

Line 10 evaluates all properties of the process model, including the residuals and the numerical vectors for bounds and residuals, useful for debugging and analysis. The physically relevant entries are props["thermo_props"] and props["model_props"].

To print the results a bit nicer, we here coded a function that filters the results by excluding keys of minor interest and already formatting the printing of Quality instances, and another one to print first the stream results and then the model properties:

18def print_result(result):
19    props = result.properties
20    print("Stream properties:")
21    excluded = "A G Mw S U V m mu mw n w x".split()
22    thermo_props = filter_results(props["thermo_props"], excluded)
23    pprint(thermo_props)
24    print()
25    print("Model properties:")
26    pprint(filter_results(props["model_props"]))
27
28
29def filter_results(results, excluded_keys=None):
30    excluded_keys = [] if excluded_keys is None else excluded_keys
31    try:
32        items = results.items()
33    except AttributeError:
34        return f"{results:.5g~}"
35    return {key: filter_results(value, excluded_keys)
36            for key, value in items if key not in excluded_keys}

Interesting results there are:

  • Given the 10 MW totally available duty, The turbine power is 4.40 MW.

  • The boiler duty is 7.25 MW, while 2.75 MW are required in the superheater to avoid too much condensation.

  • The demand on BFW is 16 t/h, while the boiler water circulation flow is 165 t/h, as expected a bit more than 10 times the feed flow.

And more …

This is where the fun really starts, as it is possible to investigate the effect of variations. Mostly this is left as an exercise, but just as an example, we can check the impact of improving our cooling towers and being able to cool the condensate to 35 °C:

12    temp = Quantity(35, "degC")
13    solver.model_parameters["model_params"]["condenser"]["temperature"] = temp
14    result = solver.solve()
15    print_result(result)

Now:

  • With an increase by 372 kW, The turbine now gives 4.72 MW of power.

  • The boiler duty must however be reduced to 7 MW, giving 3 MW to the superheater, limited by the amount of condensation.

  • The demand on BFW is reduced to 15.5 t/h.

The required temperature of superheated steam was 490 °C before, but with better cooling downstream, 521 °C are now required.

Note

Constrained by the margins on materials (creep strength of the steel used in the blades and the superheater tubes), and practically, such configuration is not found in real installations. The turbine would typically be split into several stages. Re-heating of steam can then happen on a lower temperature level and as such normally utilizing more of the available heat in a given process.

Though entirely feasible to model in \(\color{#0044aa}\mathtt{SigmaMu}\), this is also left as an exercise to the ambigious user.

Pre-built models

Todo

Define some of the above developed models in the simu.app.models module, so they do not need to be reinvented. This is definitely VLE, SpeciesBalance, and Turbine, maybe some HX. Maybe make this an own chapter in the tutorial.