Source code for idaes.models_extra.power_generation.unit_models.helm.turbine_outlet

#################################################################################
# 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).
#
# Copyright (c) 2018-2023 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.
#################################################################################
"""
Steam turbine outlet stage model.  This model is based on:

Liese, (2014). "Modeling of a Steam Turbine Including Partial Arc Admission
    for Use in a Process Simulation Software Environment." Journal of Engineering
    for Gas Turbines and Power. v136.
"""
# TODO: Missing docstrings
# pylint: disable=missing-class-docstring

__Author__ = "John Eslick"

from pyomo.environ import Var, sqrt, value, units as pyunits
from idaes.models_extra.power_generation.unit_models.helm.turbine import (
    HelmIsentropicTurbineData,
)
from idaes.core import declare_process_block_class
from idaes.core.util import from_json, to_json, StoreSpec
from idaes.core.solvers import get_solver
import idaes.core.util.scaling as iscale

import idaes.logger as idaeslog

_log = idaeslog.getLogger(__name__)


[docs]@declare_process_block_class( "HelmTurbineOutletStage", doc="Outlet stage steam turbine model", ) class HelmTurbineOutletStageData(HelmIsentropicTurbineData): # Same settings as the default pressure changer, but force to expander with # isentropic efficiency CONFIG = HelmIsentropicTurbineData.CONFIG()
[docs] def build(self): super().build() self.flow_coeff = Var( initialize=0.0333, doc="Turbine flow coefficient [kg*C^0.5/s/Pa]", units=pyunits.kg * pyunits.K**0.5 / pyunits.s / pyunits.Pa, ) self.eff_dry = Var(initialize=0.87, doc="Turbine dry isentropic efficiency") self.design_exhaust_flow_vol = Var( initialize=6000.0, doc="Design exit volumetirc flowrate [m^3/s]", units=pyunits.m**3 / pyunits.s, ) self.efficiency_mech = Var(initialize=1.0, doc="Turbine mechanical efficiency") self.efficiency_isentropic.unfix() self.eff_dry.fix() self.design_exhaust_flow_vol.fix() self.flow_coeff.fix() self.efficiency_mech.fix() self.tel_c0 = Var( initialize=0.0064 * 1e6, units=pyunits.J / pyunits.mol, doc="c0 in tel = c0 + c1*fr + c2*fr**2 + ... + c5*fr**5 (fr is ratio" " of exhaust volumetric flow to design exhaust volumetric flow)", ) self.tel_c1 = Var( initialize=-0.0328 * 1e6, units=pyunits.J / pyunits.mol, doc="c1 in tel = c0 + c1*fr + c2*fr**2 + ... + c5*fr**5 (fr is ratio" " of exhaust volumetric flow to design exhaust volumetric flow)", ) self.tel_c2 = Var( initialize=0.0638 * 1e6, units=pyunits.J / pyunits.mol, doc="c2 in tel = c0 + c1*fr + c2*fr**2 + ... + c5*fr**5 (fr is ratio" " of exhaust volumetric flow to design exhaust volumetric flow)", ) self.tel_c3 = Var( initialize=-0.0542 * 1e6, units=pyunits.J / pyunits.mol, doc="c3 in tel = c0 + c1*fr + c2*fr**2 + ... + c5*fr**5 (fr is ratio" " of exhaust volumetric flow to design exhaust volumetric flow)", ) self.tel_c4 = Var( initialize=0.022 * 1e6, units=pyunits.J / pyunits.mol, doc="c4 in tel = c0 + c1*fr + c2*fr**2 + ... + c5*fr**5 (fr is ratio" " of exhaust volumetric flow to design exhaust volumetric flow)", ) self.tel_c5 = Var( initialize=-0.0035 * 1e6, units=pyunits.J / pyunits.mol, doc="c5 in tel = c0 + c1*fr + c2*fr**2 + ... + c5*fr**5 (fr is ratio" " of exhaust volumetric flow to design exhaust volumetric flow)", ) self.tel_c0.fix() self.tel_c1.fix() self.tel_c2.fix() self.tel_c3.fix() self.tel_c4.fix() self.tel_c5.fix() @self.Expression(self.flowsheet().time, doc="Total exhaust loss curve") def tel(b, t): f = b.control_volume.properties_out[t].flow_vol / b.design_exhaust_flow_vol return ( +self.tel_c5 * f**5 + self.tel_c4 * f**4 + self.tel_c3 * f**3 + self.tel_c2 * f**2 + self.tel_c1 * f + self.tel_c0 ) @self.Constraint(self.flowsheet().time, doc="Stodola eq. choked flow") def stodola_equation(b, t): flow = b.control_volume.properties_in[t].flow_mol mw = b.control_volume.properties_in[t].mw Tin = b.control_volume.properties_in[t].temperature Pin = b.control_volume.properties_in[t].pressure Pr = b.ratioP[t] cf = b.flow_coeff return flow**2 * mw**2 * (Tin) == (cf**2 * Pin**2 * (1 - Pr**2)) @self.Constraint(self.flowsheet().time, doc="Efficiency correlation") def efficiency_correlation(b, t): x = b.control_volume.properties_out[t].vapor_frac eff = b.efficiency_isentropic[t] dh_isen = b.delta_enth_isentropic[t] tel = b.tel[t] return eff == b.eff_dry * x * (1 - 0.65 * (1 - x)) * (1 + tel / dh_isen) @self.Expression(self.flowsheet().time, doc="Thermodynamic power [J/s]") def power_thermo(b, t): return b.control_volume.work[t] @self.Expression(self.flowsheet().time, doc="Shaft power [J/s]") def power_shaft(b, t): return b.power_thermo[t] * b.efficiency_mech
[docs] def initialize_build( self, outlvl=idaeslog.NOTSET, solver=None, optarg=None, calculate_cf=True, ): """ Initialize the outlet turbine stage model. This deactivates the specialized constraints, then does the isentropic turbine initialization, then reactivates the constraints and solves. Args: outlvl : sets output level of initialization routine solver (str): Solver to use for initialization optarg (dict): Solver arguments dictionary """ init_log = idaeslog.getInitLogger(self.name, outlvl, tag="unit") solve_log = idaeslog.getSolveLogger(self.name, outlvl, tag="unit") sp = StoreSpec.value_isfixed_isactive(only_fixed=True) istate = to_json(self, return_dict=True, wts=sp) # sp is what to save to make sure state after init is same as the start # saves value, fixed, and active state, doesn't load originally free # values, this makes sure original problem spec is same but initializes # the values of free vars for t in self.flowsheet().time: if self.outlet.pressure[t].fixed: self.ratioP[t] = value(self.outlet.pressure[t] / self.inlet.pressure[t]) self.deltaP[t] = value(self.outlet.pressure[t] - self.inlet.pressure[t]) # Deactivate special constraints self.stodola_equation.deactivate() self.efficiency_correlation.deactivate() self.efficiency_isentropic.fix() self.deltaP.unfix() self.ratioP.unfix() self.inlet.fix() self.outlet.unfix() super().initialize_build(outlvl=outlvl, solver=solver, optarg=optarg) for t in self.flowsheet().time: mw = self.control_volume.properties_in[t].mw Tin = self.control_volume.properties_in[t].temperature Pin = self.control_volume.properties_in[t].pressure Pr = self.ratioP[t] if not calculate_cf: cf = self.flow_coeff self.inlet.flow_mol[t].fix( value(cf * Pin * sqrt(1 - Pr**2) / mw / sqrt(Tin)) ) super().initialize_build(outlvl=outlvl, solver=solver, optarg=optarg) self.control_volume.properties_out[:].pressure.fix() # Free eff_isen and activate special constraints self.efficiency_isentropic.unfix() self.outlet.pressure.fix() if calculate_cf: self.flow_coeff.unfix() self.inlet.flow_mol.unfix() self.inlet.flow_mol[0].fix() flow = self.control_volume.properties_in[0].flow_mol mw = self.control_volume.properties_in[0].mw Tin = self.control_volume.properties_in[0].temperature Pin = self.control_volume.properties_in[0].pressure Pr = self.ratioP[0] self.flow_coeff.value = value(flow * mw * sqrt(Tin / (1 - Pr**2)) / Pin) else: self.inlet.flow_mol.unfix() self.stodola_equation.activate() self.efficiency_correlation.activate() # Create solver slvr = get_solver(solver, optarg) with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: res = slvr.solve(self, tee=slc.tee) init_log.info( "Initialization Complete (Outlet Stage): {}".format(idaeslog.condition(res)) ) # reload original spec if calculate_cf: cf = value(self.flow_coeff) from_json(self, sd=istate, wts=sp) if calculate_cf: # cf was probably fixed, so will have to set the value again here # if you ask for it to be calculated. self.flow_coeff = cf
def calculate_scaling_factors(self): super().calculate_scaling_factors() for t, c in self.stodola_equation.items(): s = ( iscale.get_scaling_factor( self.control_volume.properties_in[t].flow_mol, default=1, warning=True, ) ** 2 ) iscale.constraint_scaling_transform(c, s, overwrite=False)