Heat pump flowsheet

[alga821]

Hi,

I modeled a supercritical CO2 heat pump (see figure) in IDAES and it seems to work, i.e., I get an optimal solution which is consistent with my implementation in Excel. However, I still have 1 degree of freedom which I assume/hope is the heat source. However, when I specify its value (0 DoF), the model never converges. I would like to fix this so that I can specify the heat source, and the model will determine the flowrate of the working fluid. Right now, this 1 DoF entails different optimal solutions (all of them correct) with different heat source/working fluid flowrate values. Any help will be much appreciated. You may find the code below.

T = Set(initialize=[0])
m = ConcreteModel(name="Detailed heat pump")
m.fs = FlowsheetBlock(
    dynamic=False,
    time=T,)
m.fs.thermo_params_co2_stream = swco2.SWCO2ParameterBlock()
m.fs.H201 = Heater(
    property_package=m.fs.thermo_params_co2_stream,
    has_phase_equilibrium=False)
m.fs.C201 = PressureChanger(
    dynamic=False,
    property_package=m.fs.thermo_params_co2_stream,
    thermodynamic_assumption=ThermodynamicAssumption.isentropic,
    compressor=True)
m.fs.H202 = Heater(
    property_package=m.fs.thermo_params_co2_stream,
    has_phase_equilibrium=False)
m.fs.T201 = Turbine(
    property_package=m.fs.thermo_params_co2_stream,)
m.fs.s01 = Arc(source=m.fs.H201.outlet, destination=m.fs.C201.inlet)
m.fs.s02 = Arc(source=m.fs.C201.outlet, destination=m.fs.H202.inlet)
m.fs.s03 = Arc(source=m.fs.H202.outlet, destination=m.fs.T201.inlet)
TransformationFactory("network.expand_arcs").apply_to(m)
# Parameters
m.fs.T_sink = Param(default=273.15+120, units=pyunits.K)
m.fs.eff_comp = Param(default=0.8)
m.fs.eff_turb = Param(default=0.9)
m.fs.DT_high = Param(default=10)
m.fs.DT_low = Param(default=5)
# Variables
m.fs.T_source = Var(T, bounds=(273.15+20, 273.15+80), units=pyunits.K)
m.fs.T1 = Var(T, bounds=(200, 600), initialize=350, units=pyunits.K)
m.fs.T2 = Var(T, bounds=(200, 600), initialize=510, units=pyunits.K)
m.fs.T3 = Var(T, bounds=(200, 600), initialize=400, units=pyunits.K)
m.fs.T4 = Var(T, bounds=(200, 600), initialize=280, units=pyunits.K)
m.fs.T4_bis = Var(T, bounds=(200, 600), initialize=280, units=pyunits.K)
m.fs.p_lp = Var(T, bounds=(1e5, 34.85e5), initialize=30e5, units=pyunits.Pa)
m.fs.p_hp = Var(T, bounds=(10e5, 200e5), initialize=100e5, units=pyunits.Pa)
m.fs.fluid = Var(T, bounds=(2, 150), initialize=10, units=pyunits.mol/pyunits.s)
m.fs.Q_source = Var(T, bounds=(100, 100e3), initialize=50e3, units=pyunits.W)
# Equations
def pres_1(b, t):
    return b.control_volume.properties_in[t].pressure == m.fs.p_lp[t]
m.fs.H201.inlet_pres = Constraint(T, rule=pres_1)
def temp_1(b, t):
    return b.control_volume.properties_out[t].temperature == m.fs.T1[t]
m.fs.H201.outlet_temp = Constraint(T, rule=temp_1)
def temp_2(b, t):
    return m.fs.C201.control_volume.properties_out[t].temperature == m.fs.T2[t]
m.fs.C201.outlet_temp = Constraint(T, rule=temp_2)
def temp_3(b, t):
    return m.fs.H202.control_volume.properties_out[t].temperature == m.fs.T3[t]
m.fs.H202.outlet_temp = Constraint(T, rule=temp_3)
def temp_4(b, t):
    return b.control_volume.properties_in[t].temperature == m.fs.T4[t]
m.fs.H201.inlet_temp = Constraint(T, rule=temp_4)
def ratio_eqn(b, t):
    return m.fs.p_hp[t]/m.fs.p_lp[t]
m.fs.ratio = Expression(T, rule=ratio_eqn)
def checking(b, t):
    return m.fs.T201.control_volume.properties_out[t].temperature == m.fs.T4_bis[t]
m.fs.check = Constraint(T, rule=checking)
def closed_cycle(b, t):
    return m.fs.T4[t] == m.fs.T4_bis[t]
m.fs.closed_hp = Constraint(T, rule=closed_cycle)
def working_fluid_eqn(b, t):
    return m.fs.C201.control_volume.properties_in[t].flow_mol == m.fs.fluid[t]
m.fs.fluid_val = Constraint(T, rule=working_fluid_eqn)
def heat_source_eqn(b, t):
    return m.fs.H201.control_volume.heat[t] == m.fs.Q_source[t]
m.fs.heat = Constraint(T, rule=heat_source_eqn)
def t1_eqn(b, t):
    return m.fs.T1[t] == m.fs.T_source[t]-m.fs.DT_low
m.fs.t1_cons = Constraint(T, rule=t1_eqn)
def t3_eqn(b, t):
    return m.fs.T3[t] == m.fs.T_sink+m.fs.DT_high
m.fs.t3_cons = Constraint(T, rule=t3_eqn)
def C201_outlet_pres_eqn(b, t):
    return m.fs.C201.deltaP[t] == m.fs.p_hp[t]-m.fs.p_lp[t]
m.fs.C201.outlet_pres = Constraint(T, rule=C201_outlet_pres_eqn)
def T201_outlet_pres_eqn(b, t):
    return m.fs.T201.deltaP[t] == m.fs.p_lp[t]-m.fs.p_hp[t]
m.fs.T201.outlet_pres = Constraint(T, rule=T201_outlet_pres_eqn)
for t in T:
    m.fs.C201.efficiency_isentropic.fix(m.fs.eff_comp)
    m.fs.T201.efficiency_isentropic.fix(m.fs.eff_turb)
    m.fs.T_source.fix((273.15+80)*pyunits.K)
    m.fs.p_lp.fix(34.85e5*pyunits.Pa)
    m.fs.T4.fix((282.57)*pyunits.K)
    # m.fs.Q_source.fix(50e3*pyunits.W)
# Performance assessment
def COP_Carnot_eqn(b, t):
    Tlow = m.fs.T1[t]
    Thigh = m.fs.T2[t]
    return Thigh/(Thigh-Tlow)
m.fs.COP_max = Expression(T, rule=COP_Carnot_eqn)
def Q_hp_eqn(b, t):
    return -m.fs.H202.heat_duty[t]
m.fs.Q_hp = Expression(T, rule=Q_hp_eqn)
def w_c_eqn(b, t):
    return m.fs.C201.work_mechanical[t]
m.fs.w_c = Expression(T, rule=w_c_eqn)
def w_t_eqn(b, t):
    return m.fs.T201.work_mechanical[t]
m.fs.w_t = Expression(T, rule=w_t_eqn)
def P_hp_eqn(b, t):
    return m.fs.w_c[t]+m.fs.w_t[t]
m.fs.P_hp = Expression(T, rule=P_hp_eqn)
def COP_hp_eqn(b, t):
    return m.fs.Q_hp[t]/m.fs.P_hp[t]
m.fs.COP = Expression(T, rule=COP_hp_eqn)
def eff_hp_eqn(b, t):
    return m.fs.COP[t]/m.fs.COP_max[t]
m.fs.eff_hp = Expression(T, rule=eff_hp_eqn)
print('Remaining DOF: ', degrees_of_freedom(m))
# Solution
solverpath = 'C:/Workdir/Programs/Anaconda3/Lib/site-packages/idaes/core/solvers/ipopt'
solver = SolverFactory('ipopt', executable=solverpath)
results = solver.solve(m, tee=True)

[dallan-keylogic]

A few things:

# Variables
m.fs.T_source = Var(T, bounds=(273.15+20, 273.15+80), units=pyunits.K)
m.fs.T1 = Var(T, bounds=(200, 600), initialize=350, units=pyunits.K)
m.fs.T2 = Var(T, bounds=(200, 600), initialize=510, units=pyunits.K)
m.fs.T3 = Var(T, bounds=(200, 600), initialize=400, units=pyunits.K)
m.fs.T4 = Var(T, bounds=(200, 600), initialize=280, units=pyunits.K)
m.fs.T4_bis = Var(T, bounds=(200, 600), initialize=280, units=pyunits.K)
m.fs.p_lp = Var(T, bounds=(1e5, 34.85e5), initialize=30e5, units=pyunits.Pa)
m.fs.p_hp = Var(T, bounds=(10e5, 200e5), initialize=100e5, units=pyunits.Pa)
m.fs.fluid = Var(T, bounds=(2, 150), initialize=10, units=pyunits.mol/pyunits.s)
m.fs.Q_source = Var(T, bounds=(100, 100e3), initialize=50e3, units=pyunits.W)
# Equations
def pres_1(b, t):
    return b.control_volume.properties_in[t].pressure == m.fs.p_lp[t]
m.fs.H201.inlet_pres = Constraint(T, rule=pres_1)
def temp_1(b, t):
    return b.control_volume.properties_out[t].temperature == m.fs.T1[t]
m.fs.H201.outlet_temp = Constraint(T, rule=temp_1)
def temp_2(b, t):
    return m.fs.C201.control_volume.properties_out[t].temperature == m.fs.T2[t]
m.fs.C201.outlet_temp = Constraint(T, rule=temp_2)
def temp_3(b, t):
    return m.fs.H202.control_volume.properties_out[t].temperature == m.fs.T3[t]
m.fs.H202.outlet_temp = Constraint(T, rule=temp_3)
def temp_4(b, t):
    return b.control_volume.properties_in[t].temperature == m.fs.T4[t]
m.fs.H201.inlet_temp = Constraint(T, rule=temp_4)

Instead of using Vars and Constraints for these relationships, I'd suggest using Pyomo References when you want to create an alias for a unit model variable on the flowsheet. Edit: I realized that temperature is not a state variable, and therefore exists only as an Expression in the property package.

Could you please post the IPOPT log?

Have you tried troubleshooting with the Diagnostics Toolbox and SVD Toolbox yet? For the SVD toolbox, once you run svd_tbx.display_underdetermined_variables_and_constraints(), print svd_tbx.s to see the singular values (there's a PR for them to be printed automatically, but it hasn't been merged yet). Any singular value below 1e-8 is concerning, and something below 1e-12 usually means the matrix is numerically singular.

Instead of manually writing constraints between m.fs.T201.outlet and m.fs.H201.inlet, I'd add an arc, then deactivate the constraint corresponding to flow_mol on the arc_expanded block. That will write a constraint for enthalpy equality instead of temperature equality. In general I'd be wary of specifying temperature instead of enthalpy. We're operating very close to the two-phase window by necessity, and temperature isn't going to convey complete information about the stream. Additionally, there is nonsmoothness when passing between the one-phase and two-phase region, and I don't believe any of the Helmholtz property packages try any tricks to smooth the derivatives. You might try bounds on temperature using b.control_volume.properties_out[0].temperature_sat as a reference to guarantee that streams are the phase you're expecting.

Edit: I don't know how b.control_volume.properties_out[0].temperature_sat behaves for supercritical fluids. Also, be careful to stay away from the critical point. Equations of state tend to struggle near it, and large changes in physical properties in response to small changes in temperature and pressure are an inherent property of the working fluid there.

[alga821]

Hi,
Without fixing the heat source (Q_source) I get the following IPOPT logfile (Optimal solution found) logfile_1DoF.log and the following troubleshooting results:

2 WARNINGS

WARNING: 1 Degree of Freedom
WARNING: Structural singularity found
    Under-Constrained Set: 37 variables, 36 constraints
    Over-Constrained Set: 0 variables, 0 constraints

Constraints and Variables associated with smallest singular values

[2.19241215e-09 9.27579753e-05 1.90954072e-04 2.00573510e-04
2.58259520e-04 4.20392610e-04 3.76274200e-03 3.76280209e-03
7.26264413e-03 7.26265292e-03]

(only showing those with singular value < 1e-8)

Smallest Singular Value 1:

    Variables:

        fs.C201.control_volume.properties_out[0].pressure
        fs.C201.properties_isentropic[0].pressure
        fs.H202.control_volume.properties_out[0].pressure
        fs.T201.control_volume.properties_in[0].pressure
        fs.C201.control_volume.deltaP[0]
        fs.p_hp[0]
        fs.H202.control_volume.properties_in[0].pressure
        fs.T201.control_volume.deltaP[0]

    Constraints:

        fs.T201.state_material_balances[0,CO2]
        fs.T201.control_volume.material_balances[0,CO2]

[alga821]

When fixing the Q_source I get the following IPOPT logfile logfile_0DoF.log and the following troubleshooting:

0 WARNINGS

Constraints and Variables associated with smallest singular values

[3.94281886e-10 2.20157469e-09 9.27579883e-05 1.90954340e-04
2.58259289e-04 4.20392527e-04 3.76274193e-03 3.76279926e-03
7.26264245e-03 7.26265290e-03]

(only showing those with singular value < 1e-8)

Smallest Singular Value 1:

    Variables:

        fs.C201.control_volume.properties_in[0].flow_mol
        fs.H201.control_volume.properties_in[0].flow_mol
        fs.H201.control_volume.properties_out[0].flow_mol
        fs.C201.properties_isentropic[0].flow_mol
        fs.C201.control_volume.properties_out[0].flow_mol
        fs.H202.control_volume.properties_in[0].flow_mol
        fs.H202.control_volume.properties_out[0].flow_mol
        fs.T201.properties_isentropic[0].flow_mol
        fs.T201.control_volume.properties_out[0].flow_mol
        fs.T201.control_volume.properties_in[0].flow_mol
        fs.fluid[0]

    Constraints:

        fs.H201.control_volume.material_balances[0,CO2]

Smallest Singular Value 2:

    Variables:

        fs.C201.control_volume.properties_out[0].pressure
        fs.C201.properties_isentropic[0].pressure
        fs.H202.control_volume.properties_out[0].pressure
        fs.T201.control_volume.properties_in[0].pressure
        fs.C201.control_volume.deltaP[0]
        fs.p_hp[0]
        fs.H202.control_volume.properties_in[0].pressure
        fs.T201.control_volume.deltaP[0]

    Constraints:

        fs.T201.state_material_balances[0,CO2]
        fs.T201.control_volume.material_balances[0,CO2]

[dallan-keylogic]

Ah, I think I have an idea about what might be going wrong---you still need to initialize the unit models. Before you add the flowsheet-level variables and constraints, initialize the enthalpy at the entrance of H201 using the htpx helper function, then initialize each unit model in turn (You'll need to fix certain degrees of freedom for each unit model so that they have zero degrees of freedom with the inlet state fixed. For example, you need to fix the heat duty in the heaters at some initial guess.) and use propagate_state to propagate the unit model's solved outlet state to the next unit model's inlet state.

IPOPT might be confused because the CO2 is liquid at the point the variables are initialized at.

[adowling2]

Last month (?), I created a simple refridgeration cycle simulation/optimization using IDAES. The code is here: https://github.com/dowlinglab/idaes-hvacr-cycles

I just switched the repository to public. Feedback either here in this discussion or on the linked repo are appreciated.