Heat pump design and part load efficiency evaluation

This notebook presents a design of a water-to-water heat pump with internal heat exchanger and its offdesign performance. The compressor model uses manufacturer data for the performance maps and for heat exchangers UA based scaling approaches are taken. The model is implemented using TESPy.

1. The compressor map is built on frequency dependent EN 12900 efficiency maps implementing a PolynomialCompressor subclass that interpolates between the maps for different frequencies.

2. A cycle design is carried out selecting an actual compressor.

3. The design simulation is executed with the compressor map data and fixed pinches on the heat exchangers to dimension them (design UA values).

4. The part load analysis is carried out by running the model at different compressor frequencies as well as heat source and heat sink temperatures. For the compressor performance the respective performance map data are evaluated, for the heat exchangers the a Reynolds exponent scaling approach for UA is implemented as available from the MovingBoundaryHeatExchanger class [3].

Note: There are a variety of simplifying assumptions here, e.g. non-existant pressure drops, scaling of UA applied to the complete side although multiple phases are present, guessed alpha ratio values, etc.. The purpose of this is to show the process of moldeing in tespy and the rough concepts of creating the design and then running offdesign analysis.

Flowsheet and assumptions

The following assumptions are taken

• Refrigerant: R134a

• The heat source is water at 12 °C, the heat sink is water at 40 °C in design conditions. Evaporator and condenser both are designed to a minimum pinch of 5 K.

• The condenser outlet is saturated liquid and the evaporator outlet is superheated in all operating modes. The IHX drives suction side superheating and liquid subcooling.

• The motor is modeled with a constant 2.5 % loss and 5 % of the compression power is dissipated into the working fluid (effectively reduces isentropic efficiency).

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib.patches import Patch
from matplotlib.path import Path
from tespy.components import (
    CycleCloser, Source, Sink, Valve, Motor, PowerSource,
    MovingBoundaryHeatExchanger, PolynomialCompressor
)
from tespy.components.component import component_registry
from tespy.components.displacementmachinery.polynomial_compressor import (
    calc_EN12900, generate_eta_polys_from_data
)
from tespy.connections import Connection, PowerConnection
from tespy.models import ModelTemplate
from tespy.networks import Network
from tespy.tools.data_containers import GroupedComponentProperties as dc_gcp
from tespy.tools.data_containers import SimpleDataContainer as dc_simple
from tespy.tools.fluid_properties import isentropic
from tespy.tools.fluid_properties.functions import T_bubble_p, T_dew_p

wf = "R134a"

# categorical color assignment used in all figures
COLORS = {
    "hot": "#eb6834", "cold": "#2a78d6",
    "20": "#2a78d6", "25": "#eb6834", "30": "#1baf7a",
}

Compressor map

The BITZER SP-111-1 catalogue [4] provides EN 12900 performance tables for a frequency controlled semi-hermetic screw compressor (20 to 90 Hz). The tables provide cooling capacity $\dot Q_0$ and electrical power $P_e$ over evaporation and condensation temperature at 30, 50 and 90 Hz based on 10 K suction superheat without liquid subcooling.

Data

nan = float("nan")
DATA = {
    30: {
        "t_evap": [10, 7.5, 5, 0, -5, -10],
        "t_cond": [30, 40, 50],
        "Q": [
            [194300, 176900, 160700, 131800, 107200, 86200],
            [175900, 159600, 144500, 117600, 94700, 75200],
            [153000, 138300, 124600, 100300, 79700, 62300]
        ],
        "P": [
            [24.8, 25.7, 26.3, 27.1, 27.3, 27.3],
            [34.1, 34.5, 34.7, 34.9, 34.8, 34.7],
            [45.2, 45.3, 45.3, 45.2, 45.1, 45.3]
        ],
    },
    50: {
        "t_evap": [10, 7.5, 5, 0, -5, -10],
        "t_cond": [30, 40, 50, 60],
        "Q": [
            [325500, 296900, 270300, 222800, 182100, 147500],
            [296400, 269700, 245000, 200800, 163000, 130900],
            [262200, 237900, 215300, 175200, 141100, 112200],
            [223800, 202300, 182300, 147000, nan, nan]
        ],
        "P": [
            [41.2, 42.5, 43.5, 44.6, 44.9, 44.6],
            [56.4, 56.9, 57.2, 57.1, 56.6, 55.8],
            [73.7, 73.6, 73.3, 72.4, 71.3, 70.3],
            [93.2, 92.6, 91.9, 90.4, nan, nan]
        ],
    },
    90: {
        "t_evap": [10, 7.5, 5, 0, -5, -10],
        "t_cond": [30, 40, 50, 60],
        "Q": [
            [564100, 515300, 469900, 388700, 318900, 259500],
            [515700, 470200, 428100, 352600, 288100, 233200],
            [461600, 420100, 381600, 312900, 254200, 204500],
            [402700, 365500, 331100, 269800, nan, nan]
        ],
        "P": [
            [76.3, 78.7, 80.4, 82.3, 82.5, 81.6],
            [104.1, 104.9, 105.2, 104.8, 103.3, 101.2],
            [134.9, 134.4, 133.6, 131.2, 128.2, 125.2],
            [167.8, 166.0, 164.1, 160.0, nan, nan]
        ],
    },
}
T_SH, T_SC = 10.0, 0.0
DISPLACEMENT = 144.0          # m³/h ...
FREQ_DISPLACEMENT = 20.0      # ... at this frequency
ETA_MOTOR = 0.975             # assumed motor efficiency

TESPy specific implementation

TESPy transposes the data to represent polynomials for isentropic and volumetric efficiency following Cecchinato et al. [5]. In contrast to cooling capacity and power polynomials, the efficiency maps transfer to cycles with different superheat and subcooling.

The catalogue power is the total electrical power of the semi-hermetic machine, i.e. it includes the motor losses. Since the plant model uses an explicit Motor component, an assumed motor loss of 2.5 % is removed from the power data before fitting — the efficiency maps then represent the shaft power, and the motor restores the electrical power in the model without counting its losses twice.

eta_s_polys, eta_vol_polys = {}, {}
for frequency, data in DATA.items():
    cooling = pd.DataFrame(
        data["Q"], index=data["t_cond"], columns=data["t_evap"], dtype=float
    )
    power = pd.DataFrame(
        data["P"], index=data["t_cond"], columns=data["t_evap"], dtype=float
    ) * 1e3 * ETA_MOTOR  # electrical to shaft power
    eta_s_polys[frequency], eta_vol_polys[frequency] = generate_eta_polys_from_data(
        power, cooling, wf, {
            "T_sh": T_SH, "T_sc": T_SC, "frequency_poly": float(frequency),
            "displacement": DISPLACEMENT,
            "frequency_displacement": FREQ_DISPLACEMENT,
        }
    )
    print(
        f"{frequency} Hz: eta_s(2 °C, 57 °C) = "
        f"{calc_EN12900(eta_s_polys[frequency], 2, 57):.3f}, "
        f"eta_vol(2 °C, 57 °C) = {calc_EN12900(eta_vol_polys[frequency], 2, 57):.3f}"
    )

30 Hz: eta_s(2 °C, 57 °C) = 0.505, eta_vol(2 °C, 57 °C) = 0.903
50 Hz: eta_s(2 °C, 57 °C) = 0.570, eta_vol(2 °C, 57 °C) = 0.915
90 Hz: eta_s(2 °C, 57 °C) = 0.576, eta_vol(2 °C, 57 °C) = 0.920

Data visualization

Below we plot the isentropic efficiency maps at 30 Hz and 90 Hz. The gray markers show the original points, the shaded areas mark regions outside of the operating envelope.

# application envelopes in the displayed window, per frequency
ENVELOPE = {
    30: [(-10, 30), (10, 30), (10, 50), (-10, 50)],
    50: [(-10, 30), (10, 30), (10, 60), (0, 60), (-10, 50)],
    90: [(-10, 30), (10, 30), (10, 60), (0, 60), (-10, 50)],
}


def shade_outside_envelope(ax, T_evap, T_cond, envelope):
    # light transparent red over the region outside the envelope polygon
    X, Y = np.meshgrid(
        np.linspace(T_evap.min(), T_evap.max(), 300),
        np.linspace(T_cond.min(), T_cond.max(), 300),
    )
    path = Path(envelope + [envelope[0]])
    points = np.column_stack([X.flatten(), Y.flatten()])
    outside = ~path.contains_points(points, radius=1e-9)
    ax.contourf(
        X, Y, outside.reshape(X.shape).astype(float),
        levels=[0.5, 1.5], colors=["#e34948"], alpha=0.25, zorder=2
    )
    return Patch(facecolor="#e34948", alpha=0.25, label="outside envelope")


T_evap, T_cond = np.meshgrid(np.linspace(-10, 10, 60), np.linspace(30, 60, 60))

eta_fields = {
    f: calc_EN12900(eta_s_polys[f], T_evap, T_cond) for f in (30, 90)
}
levels = np.linspace(
    min(field.min() for field in eta_fields.values()),
    max(field.max() for field in eta_fields.values()),
    13
)

fig, axes = plt.subplots(1, 3, figsize=(13, 4.5))
for ax, (frequency, field) in zip(axes, eta_fields.items()):
    cs = ax.contourf(T_evap, T_cond, field, levels=levels, cmap="Blues")
    cl = ax.contour(
        T_evap, T_cond, field, levels=levels, colors="white", linewidths=0.6
    )
    ax.clabel(cl, fontsize=8, fmt="%.2f")
    data = DATA[frequency]
    points = [
        (te, tc)
        for i, tc in enumerate(data["t_cond"]) for j, te in enumerate(data["t_evap"])
        if not np.isnan(data["Q"][i][j])
    ]
    ax.plot(*zip(*points), "o", ms=4, color="#777", label="catalogue data")
    patch = shade_outside_envelope(ax, T_evap, T_cond, ENVELOPE[frequency])
    handles, labels = ax.get_legend_handles_labels()
    ax.legend(handles=handles + [patch], loc="lower right", fontsize=8)
    ax.set_xlabel("evaporation temperature in °C")
    ax.set_title(rf"$\eta_s$ at {frequency} Hz")
axes[0].set_ylabel("condensation temperature in °C")

# speed dependence at a typical operating point of our application
frequencies = sorted(eta_s_polys)
f_fine = np.linspace(20, 90, 100)
for polys, label, color in [
        (eta_s_polys, r"$\eta_s$", COLORS["cold"]),
        (eta_vol_polys, r"$\eta_{vol}$", COLORS["hot"])]:
    support = [calc_EN12900(polys[f], 2, 57) for f in frequencies]
    axes[2].plot(f_fine, np.interp(f_fine, frequencies, support), "-",
                 color=color, label=label)
    axes[2].plot(frequencies, support, "o", ms=6, color=color)
axes[2].set_xlabel("frequency in Hz")
axes[2].set_ylabel("efficiency")
axes[2].set_title("speed dependence at $t_e$ = 2 °C, $t_c$ = 57 °C")
axes[2].grid(alpha=0.3)
axes[2].legend()
fig.tight_layout()
plt.show()

Both efficiencies drop with lower speed. Between the fitted frequencies we interpolate the polynomial coefficients linearly.

Implementation of the map in PolynomialCompressor

The PolynomialCompressor class available in TESPy only works with single polynomials. In this example, we create a subclass of that to allow providing polynomials per frequency:

• Two dictionaries eta_s_polys and eta_vol_polys mapping the frequencies as keys to the respective coefficients.

• Two new equations evaluating the interpolated coefficients at the current frequency

@component_registry
class MultiFrequencyPolynomialCompressor(PolynomialCompressor):

    def get_parameters(self):
        params = super().get_parameters()
        params.update({
            "eta_s_polys": dc_simple(
                dtype="dict",
                description="frequency -> eta_s polynomial coefficients"
            ),
            "eta_vol_polys": dc_simple(
                dtype="dict",
                description="frequency -> eta_vol polynomial coefficients"
            ),
            "eta_s_polys_group": dc_gcp(
                elements=["eta_s_polys", "dissipation_ratio", "frequency"],
                func=self.eta_s_polys_group_func,
                dependents=self.eta_s_polys_group_dependents,
                num_eq_sets=1,
                description="isentropic efficiency from frequency interpolated polynomials"
            ),
            "eta_vol_polys_group": dc_gcp(
                elements=["reference_state", "eta_vol_polys", "frequency"],
                func=self.eta_vol_polys_group_func,
                dependents=self.eta_vol_polys_group_dependents,
                num_eq_sets=1,
                description="displacement equation from frequency interpolated polynomials"
            ),
        })
        return params

    def _interpolate_coefficients(self, polys):
        frequencies = sorted(polys)
        coefficients = np.array([polys[key] for key in frequencies])
        return np.array([
            np.interp(self.frequency.val_SI, frequencies, coefficients[:, i])
            for i in range(coefficients.shape[1])
        ])

    def eta_s_polys_group_func(self):
        i, o = self.inl[0], self.outl[0]
        t_evap = T_dew_p(i.p.val_SI, i.fluid_data)
        t_cond = T_bubble_p(o.p.val_SI, o.fluid_data)
        c = self._interpolate_coefficients(self.eta_s_polys.val)
        eta_s = calc_EN12900(c, t_evap - 273.15, t_cond - 273.15)
        h_out_s = isentropic(
            i.p.val_SI, i.h.val_SI, o.p.val_SI, i.fluid_data, i.mixing_rule,
            T0=i.T.val_SI, T0_out=o.T.val_SI
        )
        return eta_s * (self._calc_h2() - i.h.val_SI) - (h_out_s - i.h.val_SI)

    def eta_s_polys_group_dependents(self):
        return [
            self.inl[0].p, self.inl[0].h, self.outl[0].p, self.outl[0].h,
            self.frequency
        ]

    def eta_vol_polys_group_func(self):
        i, o = self.inl[0], self.outl[0]
        t_evap = T_dew_p(i.p.val_SI, i.fluid_data)
        t_cond = T_bubble_p(o.p.val_SI, o.fluid_data)
        c = self._interpolate_coefficients(self.eta_vol_polys.val)
        eta_vol = calc_EN12900(c, t_evap - 273.15, t_cond - 273.15)
        displacement = (
            self.reference_state.val["swept_volume"] * self.frequency.val_SI
        )
        return i.m.val_SI - eta_vol * displacement / i.calc_vol()

    def eta_vol_polys_group_dependents(self):
        return [
            self.inl[0].m, self.inl[0].p, self.inl[0].h, self.outl[0].p,
            self.frequency
        ]

Validation of the subclass

Before implementing a cycle model, we verify that we can reproduce the original data with an isolated compressor model.

• T_dew on the suction connection sets the evaporation temperature/pressure.

• td_dew sets the 10 K superheat.

• T_bubble on the discharge connection sets the condensation temperature/pressure.

• To validate the cooling capacity, the enthalpy of saturated liquid at discharge temperature is evaluated.

• A Motor with the same efficiency as assumed during the data fitting is included to reproduce the original data.

from tespy.tools.fluid_properties import h_mix_pQ

reference_state = {
    "T_sh": T_SH, "T_sc": T_SC,
    "frequency_poly": 50.0,
    "displacement": DISPLACEMENT,
    "frequency_displacement": FREQ_DISPLACEMENT,
}

nw_validation = Network(iterinfo=False)
nw_validation.units.set_defaults(
    temperature="degC", pressure="bar", power="kW", pressure_difference="kPa"
)
compressor_validation = MultiFrequencyPolynomialCompressor("compressor")
v1 = Connection(Source("suction state"), "out1", compressor_validation, "in1", label="v1")
v2 = Connection(compressor_validation, "out1", Sink("discharge state"), "in1", label="v2")
nw_validation.add_conns(v1, v2)

motor_validation = Motor("motor")
e1_validation = PowerConnection(
    PowerSource("grid"), "power", motor_validation, "power_in", label="ev1"
)
e2_validation = PowerConnection(
    motor_validation, "power_out", compressor_validation, "power", label="ev2"
)
nw_validation.add_conns(e1_validation, e2_validation)
motor_validation.set_attr(eta=ETA_MOTOR)

v1.set_attr(fluid={wf: 1}, td_dew=T_SH)
compressor_validation.set_attr(
    eta_s_polys=eta_s_polys, eta_vol_polys=eta_vol_polys,
    reference_state=reference_state, dissipation_ratio=0.05
)

validation = []
for frequency, data in DATA.items():
    compressor_validation.set_attr(frequency=frequency)
    dev_Q, dev_P = [], []
    for i, t_cond in enumerate(data["t_cond"]):
        for j, t_evap in enumerate(data["t_evap"]):
            if np.isnan(data["Q"][i][j]):
                continue
            v1.set_attr(T_dew=t_evap)
            v2.set_attr(T_bubble=t_cond)
            nw_validation.solve("design")
            nw_validation.assert_convergence()

            h_liquid = h_mix_pQ(v2.p.val_SI, 0, v2.fluid_data)
            Q0 = v1.m.val_SI * (v1.h.val_SI - h_liquid)
            dev_Q.append(Q0 / data["Q"][i][j] - 1)
            dev_P.append(
                e1_validation.E.val_SI / (data["P"][i][j] * 1e3) - 1
            )
    validation.append({
        "frequency in Hz": frequency,
        "capacity deviation mean in %": np.mean(np.abs(dev_Q)) * 100,
        "capacity deviation max in %": np.max(np.abs(dev_Q)) * 100,
        "electrical power deviation mean in %": np.mean(np.abs(dev_P)) * 100,
        "electrical power deviation max in %": np.max(np.abs(dev_P)) * 100,
    })
pd.DataFrame(validation).round(2)

	frequency in Hz	capacity deviation mean in %	capacity deviation max in %	electrical power deviation mean in %	electrical power deviation max in %
0	30	0.01	0.04	0.04	0.12
1	50	0.01	0.03	0.04	0.08
2	90	0.01	0.02	0.04	0.12

Heat pump model

We implement the cycle shown in the flowsheet above in the ModelTemplate class of TESPy.

class HeatPumpModel(ModelTemplate):

    def _parameter_lookup(self):
        return {
            "condenser heat": ["Components", "condenser", "Q"],
            "frequency": ["Components", "compressor", "frequency"],
            "compressor eta_s": ["Components", "compressor", "eta_s"],
            "compressor eta_vol": ["Components", "compressor", "eta_vol"],
            "evaporator pinch": ["Components", "evaporator", "td_pinch"],
            "condenser pinch": ["Components", "condenser", "td_pinch"],
            "source temperature": ["Connections", "11", "T"],
            "source outlet temperature": ["Connections", "12", "T"],
            "sink return temperature": ["Connections", "21", "T"],
            "sink feed temperature": ["Connections", "22", "T"],
            "sink mass flow": ["Connections", "21", "m"],
            "electrical power": ["Connections", "e1", "E"],
            "evaporation temperature": ["Connections", "01", "T_dew"],
            "condensation temperature": ["Connections", "05", "T_bubble"],
            "subcooling before valve": ["Connections", "06", "td_bubble"],
            "suction superheat": ["Connections", "02", "td_dew"],
            "suction volume flow": ["Connections", "02", "v"],
            "COP": {"get": self._calc_cop},
        }

    def _calc_cop(self):
        return (
            abs(self.nw.get_comp("condenser").Q.val)
            / self.nw.get_conn("e1").E.val
        )

    def _create_network(self):
        super()._create_network()
        self.nw.iterinfo = False
        self.nw.units.set_defaults(
            temperature="degC", pressure="bar", pressure_difference="bar",
            power="kW", heat="kW", volumetric_flow="m3/h"
        )

        cc = CycleCloser("cycle closer")
        compressor = MultiFrequencyPolynomialCompressor("compressor")
        condenser = MovingBoundaryHeatExchanger("condenser")
        internal_hx = MovingBoundaryHeatExchanger("internal heat exchanger")
        valve = Valve("expansion valve")
        evaporator = MovingBoundaryHeatExchanger("evaporator")

        src_in = Source("source water inlet")
        src_out = Sink("source water outlet")
        snk_in = Source("sink water inlet")
        snk_out = Sink("sink water outlet")

        c1 = Connection(evaporator, "out2", internal_hx, "in2", label="01")
        c2 = Connection(internal_hx, "out2", compressor, "in1", label="02")
        c3 = Connection(compressor, "out1", cc, "in1", label="03")
        c4 = Connection(cc, "out1", condenser, "in1", label="04")
        c5 = Connection(condenser, "out1", internal_hx, "in1", label="05")
        c6 = Connection(internal_hx, "out1", valve, "in1", label="06")
        c7 = Connection(valve, "out1", evaporator, "in2", label="07")

        c11 = Connection(src_in, "out1", evaporator, "in1", label="11")
        c12 = Connection(evaporator, "out1", src_out, "in1", label="12")
        c21 = Connection(snk_in, "out1", condenser, "in2", label="21")
        c22 = Connection(condenser, "out2", snk_out, "in1", label="22")

        self.nw.add_conns(c1, c2, c3, c4, c5, c6, c7, c11, c12, c21, c22)

        grid = PowerSource("grid")
        motor = Motor("motor")
        e1 = PowerConnection(grid, "power", motor, "power_in", label="e1")
        e2 = PowerConnection(motor, "power_out", compressor, "power", label="e2")
        self.nw.add_conns(e1, e2)

        # invariant specifications
        motor.set_attr(eta=ETA_MOTOR)
        compressor.set_attr(dissipation_ratio=0.05)
        for hx in (condenser, internal_hx, evaporator):
            hx.set_attr(pr1=1, pr2=1)

        evaporator.set_attr(td_pinch=5)
        condenser.set_attr(td_pinch=5)

        c11.set_attr(fluid={"water": 1}, p=2)
        c21.set_attr(fluid={"water": 1}, p=2)
        c1.set_attr(fluid={wf: 1}, td_dew=5)  # evaporator superheat
        c5.set_attr(x=0)  # saturated liquid at condenser outlet

    def set_compressor_maps(self, eta_s_polys, eta_vol_polys, reference_state):
        self.nw.get_comp("compressor").set_attr(
            eta_s=None,
            eta_s_polys=eta_s_polys,
            eta_vol_polys=eta_vol_polys,
            reference_state=reference_state,
        )
        self.set_parameters(**{"condenser heat": None})

    def solve_design(self, **kwargs):
        self.solve_model_design(**kwargs)
        self.nw.assert_convergence()
        self._design_path = self.nw.save(as_dict=True)
        self._stable_solution = self._design_path

    def configure_offdesign(self, ua_scaling):
        for label, config in ua_scaling.items():
            hx = self.nw.get_comp(label)
            design = [] if label == "internal heat exchanger" else ["td_pinch"]
            hx.set_attr(design=design, offdesign=["UA_cecchinato_hc"], **config)
        self.nw.get_conn("02").set_attr(design=["td_dew"])
        self.nw.get_conn("12").set_attr(design=["T"], offdesign=["m"])

    def solve_offdesign(self, **kwargs):
        self.solve_model_offdesign(**kwargs)
        return self._solved

Design simulation

The following boundary conditions are assumed at design:

quantity	value
target heating capacity	300 kW
source water	12 °C -> 7 °C
sink water	30 °C -> 40 °C
evaporator and condenser pinch	5 K
evaporator outlet superheat	5 K
suction superheat after IHX	15 K
compressor design frequency	to match heating capacity

model = HeatPumpModel()

model.set_compressor_maps(eta_s_polys, eta_vol_polys, reference_state)
model.solve_design(
    **{
        "condenser heat": -300,
        "source temperature": 12,
        "source outlet temperature": 7,
        "sink return temperature": 30,
        "sink feed temperature": 40,
        "suction superheat": 15,
        "evaporator pinch": 5,
        "condenser pinch": 5,
        "frequency": "var"  # find frequency to match target heating capacity
    }
)

results_design = model.get_results([
    "evaporation temperature", "condensation temperature",
    "subcooling before valve", "suction superheat",
    "compressor eta_s", "compressor eta_vol",
    "suction volume flow", "electrical power", "frequency", "COP"
])
for key, value in results_design.items():
    print(f"{key + ':':28} {value:10.2f}")

F_DESIGN = results_design["frequency"]

evaporation temperature:           2.00
condensation temperature:         43.24
subcooling before valve:           5.98
suction superheat:                15.00
compressor eta_s:                  0.64
compressor eta_vol:                0.96
suction volume flow:             388.02
electrical power:                 69.76
frequency:                        56.07
COP:                               4.30

The design point can be visualized inside the compressor map, which is interpolated for the respective frequency of about 56 Hz. This design value provides quite a lot of flexibility towards part load operation as well as running above design capacity:

def interpolate_polys(polys, frequency):
    # same linear coefficient interpolation the compressor class applies
    support = sorted(polys)
    coefficients = np.array([polys[f] for f in support])
    return np.array([
        np.interp(frequency, support, coefficients[:, i])
        for i in range(coefficients.shape[1])
    ])


results_map = model.get_results([
    "evaporation temperature", "condensation temperature", "frequency"
])

T_evap, T_cond = np.meshgrid(np.linspace(-10, 10, 60), np.linspace(30, 60, 60))
eta_s_map = calc_EN12900(interpolate_polys(eta_s_polys, F_DESIGN), T_evap, T_cond)

fig, ax = plt.subplots(figsize=(10.5, 5.4))
cs = ax.contourf(T_evap, T_cond, eta_s_map, levels=12, cmap="Blues")
cl = ax.contour(
    T_evap, T_cond, eta_s_map, levels=cs.levels, colors="white",
    linewidths=0.6
)
ax.clabel(cl, fontsize=8, fmt="%.2f")

# catalogue data grid and envelope of the 50 and 90 Hz support maps
data = DATA[50]
grid = [
    (te, tc)
    for i, tc in enumerate(data["t_cond"])
    for j, te in enumerate(data["t_evap"])
    if not np.isnan(data["Q"][i][j])
]
ax.plot(
    results_map["evaporation temperature"],
    results_map["condensation temperature"],
    "*", ms=16, color="#eb6834", markeredgecolor="white",
    label="design point"
)
patch = shade_outside_envelope(ax, T_evap, T_cond, ENVELOPE[50])

ax.set_xlabel("evaporation temperature in °C")
ax.set_ylabel("condensation temperature in °C")
ax.set_title(rf"$\eta_s$ at {F_DESIGN:.1f} Hz (interpolated)")
handles, labels = ax.get_legend_handles_labels()
ax.legend(handles=handles + [patch], loc="lower right")
fig.tight_layout()
plt.show()

The operating point is far up in the map, with an isentropic efficiency which is lower than the optimal efficiency since the lift is quite high for what the compressor is optimized for. A different compressor would be much better in this context, however, these data are I was able to pull out and I could not find another datasheet with a better fitting compressor.

We can also plot the The $\log p$-$h$ diagram of the refrigerant cycle and the $Q$-$T$ diagrams of the heat exchangers.

fig, ax = plt.subplots(figsize=(10.5, 5))
model.plot_logph_diagram_matplotlib("01", ax=ax)
fig.tight_layout()
plt.show()

fig, axes = plt.subplots(1, 3, figsize=(13, 4))
for ax, label in zip(axes, ["condenser", "evaporator", "internal heat exchanger"]):
    model.plot_QT_diagram_matplotlib(label, ax=ax)
    ax.set_title(label)
fig.tight_layout()
plt.show()

Offdesign model

The compressor utilizing the map already is the offdesign model. For the heat exchangers UA scaling is applied following the methodology shown in Cecchinato et al. [5].

$$ f_{UA}=\frac{1 + x} {\left(\frac{\dot m_\text{cold}}{\dot m_\text{cold,design}}\right)^{-n_\text{cold}} + x \cdot \left(\frac{\dot m_\text{hot}}{\dot m_\text{hot,design}}\right)^{-n_\text{hot}}} \qquad 0 = UA_\text{design} \cdot f_{UA} - \sum_j \frac{\dot Q_j}{\Delta T_{log,j}} $$

$x$ is the product of alpha_ratio and area_ratio. The exponents for the mass flow are Reynolds exponents, that are supplied by [5] as well. Quoilin also compiled Reynolds exponents for different configurations in his PhD-thesis [6].

To apply the UA scaling other design boundary conditions are replaced:

• condenser and evaporator pinch

• suction gas superheat at the IHX outlet

• for the heat sink feed and return temperatures are fixed and the mass flow adapts to the available heat production

• on the source side the inlet temperature and design mass flow are fixed and the outlet temperature results from that

• the compressor frequency controls the part load ratio

UA_SCALING = {
    "condenser": {
        "re_exp_hot": 0.4,    # condensing refrigerant
        "re_exp_cold": 0.8,   # sink water side
        "alpha_ratio": 2.0, "area_ratio": 1.0,
    },
    "evaporator": {
        "re_exp_hot": 0.8,    # source water side
        "re_exp_cold": 0.5,   # evaporating refrigerant
        "alpha_ratio": 0.7, "area_ratio": 1.0,
    },
    "internal heat exchanger": {
        "re_exp_hot": 0.8,    # liquid refrigerant
        "re_exp_cold": 0.8,   # suction gas
        "alpha_ratio": 0.3, "area_ratio": 1.0,
    },
}

model.configure_offdesign(UA_SCALING)

# consistency check: off-design at the design point reproduces the design
model.set_parameters(**{"condenser heat": None})
model.solve_offdesign(frequency=F_DESIGN)
check = model.get_results(["condenser heat", "COP", "condenser pinch", "evaporator pinch"])
for key, value in check.items():
    print(f"{key + ':':20} {value:10.2f}")

condenser heat:         -300.00
COP:                       4.30
condenser pinch:           5.00
evaporator pinch:          5.00

Parameter sweep

We sweep the compressor frequency from 30 Hz to 80 Hz around the design across two operating dimensions:

• three sink temperature levels with the design spread of 10 K (20 -> 30 °C, 25 -> 35 °C and 30 -> 40 °C) at the design source temperature of 12 °C

• three source temperatures (8, 12, 16 °C) at the design sink level

# design frequency first, then outward on a 5 Hz grid (grid points
# overlapping the design frequency are dropped)
grid = [f for f in np.arange(30.0, 85.0, 5.0) if abs(f - F_DESIGN) > 2.5]
frequencies = (
    [F_DESIGN]
    + sorted([f for f in grid if f < F_DESIGN], reverse=True)
    + sorted([f for f in grid if f > F_DESIGN])
)
scenarios = (
    # sink temperature levels at design source temperature
    [{"T_source": 12, "T_return": T_return} for T_return in [20, 25, 30]]
    # source temperatures at design sink level
    + [{"T_source": T_source, "T_return": 30} for T_source in [8, 16]]
)
tracked = [
    "condenser heat", "electrical power", "COP",
    "compressor eta_s", "compressor eta_vol",
    "evaporation temperature", "condensation temperature",
    "subcooling before valve", "suction superheat",
    "evaporator pinch", "condenser pinch"
]

records = []
for scenario in scenarios:
    label = (
        f"source {scenario['T_source']} °C, "
        f"sink {scenario['T_return']} -> {scenario['T_return'] + 10} °C"
    )
    for frequency in frequencies:
        solved = model.solve_offdesign(**{
            "frequency": frequency,
            "source temperature": scenario["T_source"],
            "sink return temperature": scenario["T_return"],
            "sink feed temperature": scenario["T_return"] + 10,
        })
        if not solved:
            print(f"skipped: {label} at {frequency} Hz (not converged)")
            continue
        pinches = model.get_results(["evaporator pinch", "condenser pinch"])
        if min(pinches.values()) <= 0:
            print(f"skipped: {label} at {frequency} Hz (temperature cross)")
            continue
        records.append({
            "T_source": scenario["T_source"],
            "T_return": scenario["T_return"],
            "frequency": frequency,
            **{key: model.get_parameter(key) for key in tracked},
        })
    # back to the design point for stable starting values
    model.solve_offdesign(**{
        "frequency": F_DESIGN,
        "source temperature": 12,
        "sink return temperature": 30,
        "sink feed temperature": 40,
    })

results = pd.DataFrame(records).sort_values(["T_source", "T_return", "frequency"])
results["heating capacity"] = results.pop("condenser heat").abs()
# the frame keeps every tracked quantity, the table shows the main ones
overview = [
    "T_source", "T_return", "frequency", "heating capacity",
    "electrical power", "COP", "compressor eta_s"
]
results[overview].round(3).head()

	T_source	T_return	frequency	heating capacity	electrical power	COP	compressor eta_s
38	8	30	30.0	158.525	36.030	4.400	0.628
37	8	30	35.0	181.706	42.201	4.306	0.630
36	8	30	40.0	204.329	48.372	4.224	0.633
35	8	30	45.0	226.481	54.531	4.153	0.636
34	8	30	50.0	248.227	60.669	4.091	0.640

Results

The PQ diagram shows the heat production over the power drawn by the motor at the different heat sink temperature levels.

sink_results = results.loc[results["T_source"] == 12]

fig, ax = plt.subplots(figsize=(10.5, 5.5))

P_max = sink_results["electrical power"].max() * 1.09
Q_max = sink_results["heating capacity"].max() * 1.12
for cop in [4.0, 4.5, 5.0, 5.5, 6.0]:
    ax.plot([0, P_max], [0, P_max * cop], color="#999", lw=0.8, ls="--", zorder=1)
    if P_max * cop < Q_max:
        ax.annotate(
            f"COP {cop}", (P_max, P_max * cop), color="#777", fontsize=8,
            ha="right", va="bottom", xytext=(-4, 2), textcoords="offset points"
        )
    else:
        ax.annotate(
            f"COP {cop}", (Q_max / cop, Q_max), color="#777", fontsize=8,
            ha="center", va="top", xytext=(8, -4), textcoords="offset points"
        )

for T_return, group in sink_results.groupby("T_return"):
    color = COLORS[str(T_return)]
    ax.plot(group["electrical power"], group["heating capacity"], "-o",
            ms=4.5, color=color, zorder=3,
            label=f"sink {T_return} → {T_return + 10} °C")

# one frequency label per end of the curve bundle, at the middle curve
middle = sink_results.loc[sink_results["T_return"] == 25]
for point, offset, ha in [
    (middle.iloc[0], (0, -16), "center"), (middle.iloc[-1], (8, 4), "left")
]:
    ax.annotate(
        f'{point["frequency"]:.0f} Hz',
        (point["electrical power"], point["heating capacity"]),
        color="#555", fontsize=9, ha=ha, xytext=offset,
        textcoords="offset points"
    )

ax.set_xlim(sink_results["electrical power"].min() * 0.72, P_max)
ax.set_ylim(sink_results["heating capacity"].min() * 0.55, Q_max * 1.05)
ax.set_xlabel("electrical power in kW")
ax.set_ylabel("heating capacity in kW")
ax.grid(alpha=0.3)
ax.legend(loc="upper left")
fig.tight_layout()
plt.show()

COP over the relative load:

design_point_mask = (
    (results["T_source"] == 12)
    & (results["T_return"] == 30)
    & (results["frequency"] == F_DESIGN)
)
Q_design_actual = results.loc[design_point_mask, "heating capacity"].iloc[0]

fig, ax = plt.subplots(figsize=(10.5, 4.8))
for T_return, group in sink_results.groupby("T_return"):
    color = COLORS[str(T_return)]
    ax.plot(
        group["heating capacity"] / Q_design_actual * 100, group["COP"],
        "-o", ms=4.5, color=color,
        label=f"sink {T_return} → {T_return + 10} °C"
    )

ax.set_xlabel("relative heating capacity in % of design")
ax.set_ylabel("COP")
ax.grid(alpha=0.3)
ax.legend()
fig.tight_layout()
plt.show()

The isentropic and volumetric efficiencies at the same operating points. The sharp peak at 50 Hz is a result of the compressor maps being set up based on 30 Hz, 50 Hz and 90 Hz only.

fig, axes = plt.subplots(1, 2, figsize=(11, 4.2))
for T_return, group in sink_results.groupby("T_return"):
    color = COLORS[str(T_return)]
    axes[0].plot(group["frequency"], group["compressor eta_s"], "-o", ms=4,
                 color=color, label=f"sink {T_return} → {T_return + 10} °C")
    axes[1].plot(group["frequency"], group["compressor eta_vol"], "-o", ms=4,
                 color=color)
axes[0].set_ylabel(r"isentropic efficiency $\eta_s$")
axes[1].set_ylabel(r"volumetric efficiency $\eta_{vol}$")
for ax in axes:
    ax.set_xlabel("frequency in Hz")
    ax.grid(alpha=0.3)
axes[0].legend()
fig.tight_layout()
plt.show()

All operating points at 30 Hz, 56 Hz (design frequency) and 80 Hz on the compressor performance map for isentropic efficiency:

map_frequencies = [30.0, F_DESIGN, 80.0]
eta_maps_op = {
    f: calc_EN12900(interpolate_polys(eta_s_polys, f), T_evap, T_cond)
    for f in map_frequencies
}
levels = np.linspace(
    min(field.min() for field in eta_maps_op.values()),
    max(field.max() for field in eta_maps_op.values()),
    13
)

fig, axes = plt.subplots(1, 3, figsize=(13.5, 4.2), sharey=True)
for ax, f in zip(axes, map_frequencies):
    cs = ax.contourf(T_evap, T_cond, eta_maps_op[f], levels=levels, cmap="Blues")
    cl = ax.contour(
        T_evap, T_cond, eta_maps_op[f], levels=levels, colors="white",
        linewidths=0.6
    )
    ax.clabel(cl, fontsize=8, fmt="%.2f")

    at_frequency = results.loc[results["frequency"] == f]
    sink_points = at_frequency.loc[at_frequency["T_source"] == 12]
    for T_return, group in sink_points.groupby("T_return"):
        ax.scatter(group["evaporation temperature"],
                   group["condensation temperature"], s=30,
                   color=COLORS[str(T_return)], zorder=3,
                   label=f"sink {T_return} → {T_return + 10} °C")
    source_points = at_frequency.loc[at_frequency["T_source"] != 12]
    ax.scatter(source_points["evaporation temperature"],
               source_points["condensation temperature"], s=30, marker="D",
               color="#777", zorder=3, label="source variation")

    shade_outside_envelope(
        ax, T_evap, T_cond, ENVELOPE[30] if f == 30 else ENVELOPE[50]
    )
    ax.set_xlabel("evaporation temperature in °C")
    ax.set_title(rf"$\eta_s$ at {f:.0f} Hz")

axes[1].plot(results_map["evaporation temperature"],
             results_map["condensation temperature"], "*", ms=15,
             color="white", markeredgecolor="black", label="design point")
axes[1].set_title(rf"$\eta_s$ at {F_DESIGN:.0f} Hz (design)")
axes[0].set_ylabel("condensation temperature in °C")
axes[1].legend(loc="lower right", fontsize=8)
fig.tight_layout()
plt.show()

Finally, the source temperature variation. Here the temperature level changes the achievable heat production range with each frequency.

source_results = results.loc[results["T_return"] == 30]
source_colors = {8: "#2a78d6", 12: "#eb6834", 16: "#1baf7a"}

fig, ax = plt.subplots(figsize=(10.5, 4.8))
for T_source, group in source_results.groupby("T_source"):
    ax.plot(
        group["heating capacity"] / Q_design_actual * 100, group["COP"],
        "-o", ms=4.5, color=source_colors[T_source],
        label=f"source {T_source} °C"
    )

ax.set_xlabel("relative heating capacity in % of design")
ax.set_ylabel("COP")
ax.grid(alpha=0.3)
ax.legend()
fig.tight_layout()
plt.show()