Heat Pump Model

Here is a simple example of how to do a heat pump cycle calculation for a simple four-component system, with very simple models for each component

from dataclasses import dataclass
import os

import numpy as np
import scipy

import teqp

R = 8.31446261815324
k_B = 1.380649e-23
N_A = R/k_B

@dataclass
class ModelCombo:

    model : object
    aig: object
    name: str

    def get_h(self, T, rhomolar, molefrac):
        Atot10 = self.model.get_Ar10(T, rhomolar, molefrac) + self.aig.get_Aig10(T, rhomolar, molefrac)
        return R*T*(1+self.model.get_Ar01(T, rhomolar, molefrac) + Atot10)

    def get_s(self, T, rhomolar, molefrac):
        Atot10 = self.model.get_Ar10(T, rhomolar, molefrac) + self.aig.get_Aig10(T, rhomolar, molefrac)
        Atot00 = self.model.get_Ar00(T, rhomolar, molefrac) + self.aig.get_Aig00(T, rhomolar, molefrac)
        return R*(Atot10 - Atot00)

    def get_p(self, T, rhomolar, molefrac):
        return rhomolar*self.model.get_R(molefrac)*T*(1 + self.model.get_Ar01(T, rhomolar, molefrac))

def cycle(combo, anc, *, Tevap, Tcond, DELTAT_sh, DELTAT_sc, eta_comp):
    """
    combo(ModelCombo): The joined model with ideal-gas and residual portions
    anc : A set of ancillary functions implementing rhoL(T) and rhoV(T) methods
    Tevap : Saturated vapor temperature in evaporator, in K
    Tcond : Saturated vapor temperature in condenser, in K
    DELTAT_sh : superheat, in K
    DELTAT_sc : subcooling, in K
    eta_comp : compressor efficiency
    """

    model = combo.model
    z = np.array([1.0])

    # VLE densities,
    # w/ guess values from the ancillary
    rhomolar1satL, rhomolar1sat = model.pure_VLE_T(Tevap, anc.rhoL(Tevap), anc.rhoV(Tevap), 10)
    rhomolar3sat, rhomolar3satV = model.pure_VLE_T(Tcond, anc.rhoL(Tcond), anc.rhoV(Tcond), 10)
    p1 = combo.get_p(Tevap, rhomolar1sat, z)
    p2 = combo.get_p(Tcond, rhomolar3sat, z)

    # Evaporator outlet & compressor inlet @ state point 1
    T1 = Tevap + DELTAT_sh
    rhomolar1 = scipy.optimize.newton(lambda rho_: combo.get_p(T1, rho_, z)-p1, rhomolar1sat)
    h1 = combo.get_h(T1, rhomolar1, z)
    s1 = combo.get_s(T1, rhomolar1, z)

    # Solve for isentropic compressor outlet
    res = lambda x: [combo.get_p(x[0], x[1], z)-p2, combo.get_s(x[0], x[1], z)-s1]
    T2s, rho2s = scipy.optimize.fsolve(res, [Tcond, rhomolar1sat])
    h2s = combo.get_h(T2s, rho2s, z)
    h2 = h1 + (h2s-h1)/eta_comp # @ state point 2

    # Condenser outlet and expansion valve inlet @ state point 3
    T3 = Tcond - DELTAT_sc
    rhomolar3 = scipy.optimize.newton(lambda rho_: combo.get_p(T3, rho_, z)-p2, rhomolar3sat)
    h3 = combo.get_h(T3, rhomolar3, z)

    COP = (h1-h3)/(h2-h1)

    return {
        'name': combo.name,
        'COP': COP,
        'pevap / kPa': p1/1e3,
        'pcond / kPa': p2/1e3,
        'rho1 / mol/m^3': rhomolar1,
        'Qvol / kJ/m^3': (h1-h3)*rhomolar1/1e3,
    }

# Build the model (ideal-gas and residual)
FLD = 'R125'
path = teqp.get_datapath()+f'/dev/fluids/{FLD}.json'
assert(os.path.exists(path))
jig = teqp.convert_CoolProp_idealgas(path, 0)
combo = ModelCombo(
    model=teqp.build_multifluid_model([path], ''),
    aig=teqp.IdealHelmholtz([jig]),
    name=FLD
)

# Generic ancillary functions from on-the-fly ancillary construction
anc = teqp.build_ancillaries(combo.model, 360, 6000, 250)

cycle(combo, anc=anc, Tevap=270, Tcond=313, DELTAT_sh=5, DELTAT_sc=5, eta_comp=0.7)

{'name': 'R125',
 'COP': np.float64(3.238624644077594),
 'pevap / kPa': np.float64(606.2416455945954),
 'pcond / kPa': np.float64(2001.2735289000334),
 'rho1 / mol/m^3': np.float64(306.9219167132696),
 'Qvol / kJ/m^3': np.float64(3287.7920883249403)}

Exercise for the reader: plot the points on a P-H diagram, showing the saturated liquid and vapor states