# kringproces.py - het koelkringproces met propaan (R290), berekend met de toestandsvergelijking van Peng-Robinson
import numpy as np

R, M = 8.314462618, 44.097e-3                    # J/(mol K), kg/mol
Tk, pk, om = 369.89, 42.512e5, 0.1521            # kritisch punt (K, Pa) en acentrische factor van propaan
kap = 0.37464 + 1.54226 * om - 0.26992 * om**2
a0, b = 0.45724 * R**2 * Tk**2 / pk, 0.07780 * R * Tk / pk
C = (-4.224, 0.3063, -1.586e-4, 3.215e-8)        # cp van het ideale gas: C0 + C1 T + C2 T² + C3 T³ (J/(mol K))
w2 = np.sqrt(2)


def a(T):
    return a0 * (1 + kap * (1 - np.sqrt(T / Tk)))**2


def dadT(T):
    return -a0 * kap * (1 + kap * (1 - np.sqrt(T / Tk))) / np.sqrt(T * Tk)


def Z(T, p, fase):
    """Compressibiliteitsfactor: de kleinste wortel is de vloeistof, de grootste de damp."""
    A, B = a(T) * p / (R * T)**2, b * p / (R * T)
    w = np.roots([1, B - 1, A - 3 * B**2 - 2 * B, -(A * B - B**2 - B**3)])
    w = sorted(x.real for x in w if abs(x.imag) < 1e-9 and x.real > B)
    return (w[0] if fase == "vloeistof" else w[-1]), A, B


def L(z, B):                                     # een logaritme die in alle formules terugkomt
    return np.log((z + (1 + w2) * B) / (z + (1 - w2) * B))


def psat(T):
    """Dampspanning: de druk waarbij vloeistof en damp even 'vluchtig' zijn (gelijke fugaciteit)."""
    p = pk * np.exp(5.373 * (1 + om) * (1 - Tk / T))          # beginschatting (Wilson)
    for _ in range(60):
        zl, A, B = Z(T, p, "vloeistof")
        zv, _, _ = Z(T, p, "damp")
        lnf = [z - 1 - np.log(z - B) - A / (2 * w2 * B) * L(z, B) for z in (zl, zv)]
        p = p * np.exp(lnf[0] - lnf[1])
    return p


def hs_ruw(T, p, fase):
    """Enthalpie (J/kg) en entropie (J/(kg K)): ideaal gas vanaf 273,15 K en 1 bar, plus de afwijking van Peng-Robinson."""
    z, A, B = Z(T, p, fase)
    f = lambda t: C[0] * t + C[1] * t**2 / 2 + C[2] * t**3 / 3 + C[3] * t**4 / 4
    g = lambda t: C[0] * np.log(t) + C[1] * t + C[2] * t**2 / 2 + C[3] * t**3 / 3
    h = f(T) - f(273.15) + R * T * (z - 1) + (T * dadT(T) - a(T)) / (2 * w2 * b) * L(z, B)
    s = g(T) - g(273.15) - R * np.log(p / 1e5) + R * np.log(z - B) + dadT(T) / (2 * w2 * b) * L(z, B)
    return h / M, s / M


# referentie: verzadigde vloeistof bij 0 °C heeft h = 200 kJ/kg en s = 1 kJ/(kg K)
_h0, _s0 = hs_ruw(273.15, psat(273.15), "vloeistof")


def hs(T, p, fase):
    h, s = hs_ruw(T, p, fase)
    return h - _h0 + 200e3, s - _s0 + 1000.0


def bisectie(f, lo, hi):
    for _ in range(60):
        mid = (lo + hi) / 2
        lo, hi = (mid, hi) if f(lo) * f(mid) > 0 else (lo, mid)
    return (lo + hi) / 2


def cyclus(te, tc, oververhit=5.0, onderkoel=3.0, eta_is=0.70):
    """Verdampen bij te en condenseren bij tc (°C). Geeft drukken (Pa), enthalpieën (J/kg) en de COP van het kringproces."""
    Te, Tc = te + 273.15, tc + 273.15
    pe, pc = psat(Te), psat(Tc)
    h1, s1 = hs(Te + oververhit, pe, "damp")                            # 1: oververhitte damp naar de compressor
    T2s = bisectie(lambda T: hs(T, pc, "damp")[1] - s1, Tc, Tc + 150)   # isentropisch samenpersen
    h2 = h1 + (hs(T2s, pc, "damp")[0] - h1) / eta_is                   # 2: echte compressor
    T2 = bisectie(lambda T: hs(T, pc, "damp")[0] - h2, Tc, Tc + 200)
    h3 = hs(Tc - onderkoel, pc, "vloeistof")[0]                         # 3: onderkoelde vloeistof na de condensor
    h4 = h3                                                             # 4: smoren in het expansieventiel
    rho1 = pe * M / (Z(Te + oververhit, pe, "damp")[0] * R * (Te + oververhit))
    return {"pe": pe, "pc": pc, "h1": h1, "h2": h2, "h3": h3, "h4": h4, "T2": T2 - 273.15, "rho1": rho1,
            "qe": h1 - h4, "w": h2 - h1, "qc": h2 - h3, "cop": (h2 - h3) / (h2 - h1)}


if __name__ == "__main__":
    for t in (-42.1, 0.0, 40.0):
        print(f"propaan bij {t:6.1f} °C: dampspanning {psat(t + 273.15) / 1e5:6.2f} bar")
    c = cyclus(-4.0, 40.0)                       # buitenlucht 2 °C, water 35 °C (A2/W35)
    print(f"verdampen {c['pe'] / 1e5:.2f} bar, condenseren {c['pc'] / 1e5:.2f} bar, uit de compressor {c['T2']:.1f} °C")
    print(f"h1 = {c['h1'] / 1e3:.1f}, h2 = {c['h2'] / 1e3:.1f}, h3 = h4 = {c['h3'] / 1e3:.1f} kJ/kg")
    print(f"per kg: verdamper {c['qe'] / 1e3:.1f} kJ, compressor {c['w'] / 1e3:.1f} kJ, condensor {c['qc'] / 1e3:.1f} kJ")
    print(f"COP van het kringproces: {c['cop']:.2f}")