Example for the equation oriented approach

Example for the equation oriented approach#

For the equation oriented approach we implement the equations and respective derivatives with respect to the variables

mass flow,
pressure and
enthalpy.

The intention is to show the idea behind the approach as it is implemented in the open source software TESPy.

Equation setup#

import numpy as np
from CoolProp.CoolProp import PropsSI
import pandas as pd

The equations implemented are

pressure ratio \(pr=\frac{p_\text{out}}{p_\text{in}}\)
isentropic efficiency \(\eta_\text{s}=\frac{h_\text{2s}-h_1}{h_2 - h_1}\)
energy equation \(\dot E = \dot m \cdot \left(h_2 - h_1\right)\)
specified temperature \(T=T\left(h, p\right)\)
specified saturated state \(Q=Q\left(h, p\right)\)

The functions return a simple value, the derivatives return a dictionary indicating the partial derivatives to each of the relevant variables.

def pr_func(pr, p_1, p_2):
    return pr * p_1 - p_2

def pr_deriv(pr, p_1, p_2):
    return {"p_1": pr, "p_2": -1}

def eta_s_func(eta_s, h_1, p_1, h_2, p_2, fluid):
    h_2s = PropsSI("H", "P", p_2, "S", PropsSI("S", "H", h_1, "P", p_1, fluid), fluid)
    return (h_2 - h_1) * eta_s - (h_2s - h_1)

def eta_s_deriv(eta_s, h_1, p_1, h_2, p_2, fluid):
    d = 1e-2
    return {
        "h_1": (eta_s_func(eta_s, h_1 + d, p_1, h_2, p_2, fluid) - eta_s_func(eta_s, h_1 - d, p_1, h_2, p_2, fluid)) / (2 * d),
        "h_2": eta_s,
        "p_1": (eta_s_func(eta_s, h_1, p_1 + d, h_2, p_2, fluid) - eta_s_func(eta_s, h_1, p_1 - d, h_2, p_2, fluid)) / (2 * d),
        "p_2": (eta_s_func(eta_s, h_1, p_1, h_2, p_2 + d, fluid) - eta_s_func(eta_s, h_1, p_1, h_2, p_2 - d, fluid)) / (2 * d),
    }

def energy_func(E, m, h_1, h_2):
    return E - m * (h_2 - h_1)

def energy_deriv(E, m, h_1, h_2):
    return {
        "m": -(h_2 - h_1),
        "h_2": -m,
        "h_1": m
    }

def temperature_func(T, h, p, fluid):
    return PropsSI("T", "H", h, "P", p, fluid) - T

def temperature_deriv(h, p, fluid):
    d = 1e-2
    return {
        "h": (PropsSI("T", "H", h + d, "P", p, fluid) - PropsSI("T", "H", h - d, "P", p, fluid)) / (2 * d),
        "p": (PropsSI("T", "H", h, "P", p + d, fluid) - PropsSI("T", "H", h, "P", p - d, fluid)) / (2 * d)
    }

def saturation_func(Q, h, p, fluid):
    return PropsSI("Q", "H", h, "P", p, fluid) - Q

def saturation_deriv(h, p, fluid):
    d = 1e-2
    return {
        "h": (PropsSI("Q", "H", h + d, "P", p, fluid) - PropsSI("Q", "H", h - d, "P", p, fluid)) / (2 * d),
        "p": (PropsSI("Q", "H", h, "P", p + d, fluid) - PropsSI("Q", "H", h, "P", p - d, fluid)) / (2 * d)
    }

With these equations employed we can model operation of a heat pump. First, we set some input values:

Name of the working fluid
evaporation temperature \(T_1\)
condensation temperature \(T_3\)
compressor efficiency \(\eta_\text{s}\)
heat supplied by the condenser \(\dot Q\)

On top of that, we have to set some guess values for the variables. These are:

the mass flow,
enthalpy and pressure after evaporation,
enthalpy and pressure after compression,
enthalpy after condensation (pressure is the same as pressure after compression) and
enthalpy and pressure after throttling are also known already (enthalpy is the same as after condensation and pressure is the evaporation pressure).

Note

With the process knowledge we could also calculate a couple more states already. E.g. in TESPy such a presolving step is also included, which generally decreases the problem size. However to showcase the general idea behind the method this step is left out in this context.

fluid = "R290"
t_1 = 10 + 273.15
t_3 = 60 + 273.15
heat = -1000e3
eta_s = 0.8

# guess values to start
h_2 = 400e3
p_2 = 2e5
h_3 = 420e3
p_3 = 20e5
h_4 = 400e3
m = 1
variables = np.array([h_2, p_2, h_3, p_3, h_4, m])
residual = np.ones(6)

Next we set up the Newon-Raphson method by calculating the residuals and the Jacobian matrix and updating the input variables after every iteration. This is done until the norm of the residual vector is sufficiently low.

while np.linalg.norm(residual) > 1e-4:

    ev_outlet_sat = saturation_func(1, variables[0], variables[1], fluid)
    ev_temp = temperature_func(t_1, variables[0], variables[1], fluid)
    cp_eff = eta_s_func(eta_s, variables[0], variables[1], variables[2], variables[3], fluid)
    cd_outlet_sat = saturation_func(0, variables[4], variables[3], fluid)
    cd_temp = temperature_func(t_3, variables[4], variables[3], fluid)
    # cp_power = energy_func(power, variables[5], variables[0], variables[2])
    cd_heat = energy_func(heat, variables[5], variables[2], variables[4])

    residual = np.array([ev_outlet_sat, ev_temp, cp_eff, cd_outlet_sat, cd_temp, cd_heat])
    jacobian = np.zeros((6, 6))

    ev_outlet_sat_j = saturation_deriv(variables[0], variables[1], fluid)
    ev_temp_j = temperature_deriv(variables[0], variables[1], fluid)
    cp_eff_j = eta_s_deriv(eta_s, variables[0], variables[1], variables[2], variables[3], fluid)
    cd_outlet_sat_j = saturation_deriv(variables[4], variables[3], fluid)
    cd_temp_j = temperature_deriv(variables[4], variables[3], fluid)
    cd_heat_j = energy_deriv(heat, variables[5], variables[2], variables[4])

    jacobian[0, 0] = ev_outlet_sat_j["h"]
    jacobian[0, 1] = ev_outlet_sat_j["p"]
    jacobian[1, 0] = ev_temp_j["h"]
    jacobian[1, 1] = ev_temp_j["p"]
    jacobian[2, 0] = cp_eff_j["h_1"]
    jacobian[2, 1] = cp_eff_j["p_1"]
    jacobian[2, 2] = cp_eff_j["h_2"]
    jacobian[2, 3] = cp_eff_j["p_2"]
    jacobian[3, 3] = cd_outlet_sat_j["p"]
    jacobian[3, 4] = cd_outlet_sat_j["h"]
    jacobian[4, 3] = cd_temp_j["p"]
    jacobian[4, 4] = cd_temp_j["h"]
    jacobian[5, 2] = cd_heat_j["h_1"]
    jacobian[5, 4] = cd_heat_j["h_2"]
    jacobian[5, 5] = cd_heat_j["m"]

    variables -= np.linalg.inv(jacobian).dot(residual)

    h_sat = PropsSI("H", "P", variables[1], "Q", 1, fluid)
    if variables[0] > h_sat:
        variables[0] = h_sat * 0.9

    h_sat = PropsSI("H", "P", variables[3], "Q", 0, fluid)
    if variables[4] < h_sat:
        variables[4] = h_sat * 1.2

In post-processing we can check the results:

t_1 = PropsSI("T", "H", variables[0], "P", variables[1], fluid)
t_2 = PropsSI("T", "H", variables[2], "P", variables[3], fluid)
t_3 = PropsSI("T", "H", variables[4], "P", variables[3], fluid)
t_4 = PropsSI("T", "H", variables[4], "P", variables[1], fluid)

p_1 = variables[1]
p_2 = variables[3]
p_3 = p_2
p_4 = p_1

h_1 = variables[0]
h_2 = variables[2]
h_3 = variables[4]
h_4 = h_3

df = pd.DataFrame(index=[1, 2, 3, 4], columns=["m", "T", "h", "p"])
df.loc[1] = [variables[5], t_1, h_1, p_1]
df.loc[2] = [variables[5], t_2, h_2, p_2]
df.loc[3] = [variables[5], t_3, h_3, p_3]
df.loc[4] = [variables[5], t_4, h_4, p_4]
df

	m	T	h	p
1	3.483096	283.15	585672.587117	636601.564181
2	3.483096	343.582866	655236.564663	2116752.919654
3	3.483096	333.15	368135.6923	2116752.919654
4	3.483096	283.15	368135.6923	636601.564181

df.loc[2, "m"] * (df.loc[2, "h"] - df.loc[1, "h"])

242298.03613476685

Now, let’s see, what happens, if we change some parameters. For example, instead of the heat we specify the compressor power and instead of the efficiency we specify the compressor outlet temperature:

power = 200e3
t_2 = 360
residual = np.ones(6)

while np.linalg.norm(residual) > 1e-4:

    ev_outlet_sat = saturation_func(1, variables[0], variables[1], fluid)
    ev_temp = temperature_func(t_1, variables[0], variables[1], fluid)
    cp_temp_out = temperature_func(t_2, variables[2], variables[3], fluid)
    cd_outlet_sat = saturation_func(0, variables[4], variables[3], fluid)
    cd_temp = temperature_func(t_3, variables[4], variables[3], fluid)
    cp_power = energy_func(power, variables[5], variables[0], variables[2])

    residual = np.array([ev_outlet_sat, ev_temp, cp_temp_out, cd_outlet_sat, cd_temp, cp_power])
    jacobian = np.zeros((6, 6))

    ev_outlet_sat_j = saturation_deriv(variables[0], variables[1], fluid)
    ev_temp_j = temperature_deriv(variables[0], variables[1], fluid)
    cp_temp_out_j = temperature_deriv(variables[2], variables[3], fluid)
    cd_outlet_sat_j = saturation_deriv(variables[4], variables[3], fluid)
    cd_temp_j = temperature_deriv(variables[4], variables[3], fluid)
    cp_power_j = energy_deriv(power, variables[5], variables[0], variables[2])

    jacobian[0, 0] = ev_outlet_sat_j["h"]
    jacobian[0, 1] = ev_outlet_sat_j["p"]
    jacobian[1, 0] = ev_temp_j["h"]
    jacobian[1, 1] = ev_temp_j["p"]
    jacobian[2, 2] = cp_temp_out_j["h"]
    jacobian[2, 3] = cp_temp_out_j["p"]
    jacobian[3, 3] = cd_outlet_sat_j["p"]
    jacobian[3, 4] = cd_outlet_sat_j["h"]
    jacobian[4, 3] = cd_temp_j["p"]
    jacobian[4, 4] = cd_temp_j["h"]
    jacobian[5, 0] = cp_power_j["h_1"]
    jacobian[5, 2] = cp_power_j["h_2"]
    jacobian[5, 5] = cp_power_j["m"]

    variables -= np.linalg.inv(jacobian).dot(residual)

    h_sat = PropsSI("H", "P", variables[1], "Q", 1, fluid)
    if variables[0] > h_sat:
        variables[0] = h_sat * 0.9

    h_sat = PropsSI("H", "P", variables[3], "Q", 0, fluid)
    if variables[4] < h_sat:
        variables[4] = h_sat * 1.2

t_1 = PropsSI("T", "H", variables[0], "P", variables[1], fluid)
t_2 = PropsSI("T", "H", variables[2], "P", variables[3], fluid)
t_3 = PropsSI("T", "H", variables[4], "P", variables[3], fluid)
t_4 = PropsSI("T", "H", variables[4], "P", variables[1], fluid)

p_1 = variables[1]
p_2 = variables[3]
p_3 = p_2
p_4 = p_1

h_1 = variables[0]
h_2 = variables[2]
h_3 = variables[4]
h_4 = h_3

df = pd.DataFrame(index=[1, 2, 3, 4], columns=["m", "T", "h", "p"])
df.loc[1] = [variables[5], t_1, h_1, p_1]
df.loc[2] = [variables[5], t_2, h_2, p_2]
df.loc[3] = [variables[5], t_3, h_3, p_3]
df.loc[4] = [variables[5], t_4, h_4, p_4]
df

	m	T	h	p
1	1.820021	283.15	585672.587117	636601.564181
2	1.820021	360.0	695561.432352	2116752.919654
3	1.820021	333.15	368135.6923	2116752.919654
4	1.820021	283.15	368135.6923	636601.564181

df.loc[2, "m"] * (df.loc[2, "h"] - df.loc[1, "h"])

199999.99999999997

(PropsSI("H", "P", df.loc[2, "p"], "S", PropsSI("S", "P", df.loc[1, "p"], "H", df.loc[1, "h"], fluid), fluid) - df.loc[1, "h"]) / (df.loc[2, "h"] - df.loc[1, "h"])

0.5064315847295852

We can see that there is little change to the previous setup of the system. If we are able to abstract guessing starting values, ordering the equations and variables and connect the variables and equations with the Jacobian automatically, we can solve such systems no matter how the user specifications are and how the topology looks like (in principle) without having to extend or tweak or solution algorithm for every application. All this has been implemented in TESPy, as you will learn in the next chapter.

Example for the equation oriented approach

Contents

Example for the equation oriented approach#

Equation setup#