Engineering

Refrigeration Cycle Analysis in Python with CoolProp: A Practical Tutorial

Oktober 26, 2026 / Aldy Dwi Prasetyo

Every refrigeration or heat pump project I have touched starts with the same quick question: what pressures will this thing run at, how much compressor power does it need, and what COP can I realistically expect? You can answer that with a refrigerant chart and a ruler, but the chart reading gets old fast once you start changing temperatures or comparing refrigerants. In this tutorial we will build a small Python script with CoolProp, the free open-source thermophysical property library, that solves a basic vapor-compression cycle in about thirty lines and draws the p-h diagram for you.

The goal is not a full system design tool. It is a sanity-check calculator you can trust, that you understand line by line, and that you can extend later.

What you need

Python 3 and two packages. Install them with pip:

pip install CoolProp matplotlib numpy

Check that CoolProp imports and see which version you have:

python -c "import CoolProp; print(CoolProp.__version__)"

The script below was run with the current CoolProp release from PyPI. The high-level function we use, PropsSI, has been stable for many years, so older 6.x versions will work the same way.

How PropsSI works in one paragraph

PropsSI takes the output you want, two independent input properties with their values, and the fluid name. Everything is in SI units: temperature in K, pressure in Pa, enthalpy in J/kg, entropy in J/kg/K. For example, the saturation pressure of R134a at -10 °C is:

import CoolProp.CoolProp as CP
p_evap = CP.PropsSI("P", "T", 263.15, "Q", 1, "R134a")   # Pa

Q is vapor quality, so Q = 1 means saturated vapor and Q = 0 saturated liquid. That is the only trick you need for the two-phase states. For superheated or subcooled states you simply give temperature and pressure, and for the compressor outlet you give pressure and entropy or pressure and enthalpy.

One habit that saves a lot of debugging: always convert units at the edges of your script (input in °C and bar for humans, SI inside) and never mix kJ and J in the middle.

The cycle we are going to solve

A standard single-stage cycle has four state points:

  1. State 1, compressor inlet: evaporator pressure, slightly superheated vapor.
  2. State 2, compressor outlet: condenser pressure, hot superheated vapor.
  3. State 3, condenser outlet: condenser pressure, slightly subcooled liquid.
  4. State 4, evaporator inlet: after the expansion valve, two-phase mixture at evaporator pressure.

The example inputs I will use are typical of a small air-conditioning or chiller duty. They are example inputs, not measured data from a real unit:

  • Refrigerant: R134a
  • Evaporating temperature: -10 °C
  • Condensing temperature: 40 °C
  • Suction superheat: 5 K
  • Liquid subcooling: 3 K
  • Compressor isentropic efficiency: 0.75
  • Refrigerant mass flow: 0.05 kg/s

Assumptions: no pressure drop in the heat exchangers or lines, adiabatic compressor, and an isenthalpic expansion valve (h4 = h3). These are the usual textbook simplifications and they are fine for a first estimate.

Step 1: pressures from the saturation temperatures

import CoolProp.CoolProp as CP

fluid = "R134a"
Te = -10 + 273.15     # evaporating temperature, K
Tc = 40 + 273.15      # condensing temperature, K
SH, SC = 5.0, 3.0     # superheat and subcooling, K
eta_is = 0.75         # compressor isentropic efficiency
mdot = 0.05           # kg/s

pe = CP.PropsSI("P", "T", Te, "Q", 1, fluid)
pc = CP.PropsSI("P", "T", Tc, "Q", 0, fluid)

For a pure fluid like R134a it does not matter whether you use Q = 0 or Q = 1 here. For a zeotropic blend such as R407C it does matter, because dew and bubble temperatures differ. Using dew point (Q = 1) for the evaporator and bubble point (Q = 0) for the condenser is a sensible convention.

Step 2: compressor inlet and outlet

T1 = Te + SH
h1 = CP.PropsSI("H", "T", T1, "P", pe, fluid)
s1 = CP.PropsSI("S", "T", T1, "P", pe, fluid)

h2s = CP.PropsSI("H", "P", pc, "S", s1, fluid)   # isentropic outlet
h2 = h1 + (h2s - h1) / eta_is                    # real outlet
T2 = CP.PropsSI("T", "P", pc, "H", h2, fluid)    # discharge temperature

The isentropic efficiency definition is the usual one: ideal work divided by actual work. The discharge temperature T2 is worth printing every time, because it is the number that tells you whether oil breakdown or compressor limits will become a problem.

Step 3: condenser outlet, valve and evaporator

T3 = Tc - SC
h3 = CP.PropsSI("H", "T", T3, "P", pc, fluid)
h4 = h3                                          # isenthalpic expansion
x4 = CP.PropsSI("Q", "P", pe, "H", h4, fluid)   # quality entering evaporator

q_evap = h1 - h4      # refrigerating effect, J/kg
w_comp = h2 - h1      # compressor work, J/kg
q_cond = h2 - h3      # heat rejected, J/kg
COP = q_evap / w_comp

Note the subcooled state 3. Giving T and P slightly below saturation works fine. If you ever set exactly zero subcooling, use "Q", 0 instead of temperature and pressure, because a state sitting exactly on the saturation line is ambiguous when specified by T and P.

Step 4: results

Running the script with the inputs above gives these numbers (computed with CoolProp, rounded):

  • Evaporator pressure: 2.01 bar(a); condenser pressure: 10.17 bar(a); pressure ratio 5.07
  • h1 = 396.9 kJ/kg, h2s = 431.6 kJ/kg, h2 = 443.2 kJ/kg, h3 = h4 = 251.9 kJ/kg
  • Discharge temperature: 61.9 °C
  • Quality at evaporator inlet: 0.317
  • Refrigerating effect: 145.0 kJ/kg; compressor work: 46.3 kJ/kg; heat rejected: 191.3 kJ/kg
  • COP = 3.13
  • At 0.05 kg/s: cooling capacity 7.25 kW, compressor power 2.31 kW, condenser duty 9.56 kW

A quick check you should always do: the energy balance closes, since 7.25 + 2.31 is about 9.56 kW. Another one: the Carnot COP between -10 °C and 40 °C is 263.15 / 50 = 5.26, so a real cycle at 3.13 sits at roughly 60 percent of Carnot. That is a believable ratio for this kind of estimate. If your script ever returns a COP above Carnot, there is a unit or sign error somewhere.

Step 5: draw the p-h diagram

Plotting the cycle on the pressure-enthalpy diagram makes it much easier to spot mistakes, like a compressor inlet that accidentally landed inside the dome. The saturation dome comes straight from CoolProp:

import numpy as np
import matplotlib.pyplot as plt

Tcrit = CP.PropsSI("Tcrit", fluid)
Ts = np.linspace(233.15, Tcrit - 0.02, 400)
hl = [CP.PropsSI("H", "T", t, "Q", 0, fluid) / 1e3 for t in Ts]
hv = [CP.PropsSI("H", "T", t, "Q", 1, fluid) / 1e3 for t in Ts]
ps = [CP.PropsSI("P", "T", t, "Q", 0, fluid) / 1e5 for t in Ts]

fig, ax = plt.subplots(figsize=(8, 5.5))
ax.plot(hl, ps, "k-")
ax.plot(hv, ps, "k-")
H = np.array([h1, h2, h3, h4, h1]) / 1e3
P = np.array([pe, pc, pc, pe, pe]) / 1e5
ax.plot(H, P, "r-o")
ax.set_yscale("log")
ax.set_xlabel("Specific enthalpy h [kJ/kg]")
ax.set_ylabel("Pressure p [bar]")
ax.grid(True, which="both", alpha=0.3)
fig.savefig("ph_cycle.png", dpi=150)

The log scale on pressure is the convention you see on manufacturer charts, and it makes the low-pressure side readable. Stop the dome slightly below the critical temperature, because properties change very steeply right at the critical point.

R134a vapor-compression refrigeration cycle on a log pressure-enthalpy diagram, evaporating at -10 C and condensing at 40 C

R134a cycle from this tutorial on the p-h diagram (example inputs: Te = -10 °C, Tc = 40 °C, 5 K superheat, 3 K subcooling, isentropic efficiency 0.75). Plot generated with CoolProp and Matplotlib. Image: Bugo Center (own plot generated with CoolProp and Matplotlib; script in cycle.py).

Step 6: compare refrigerants by changing one string

This is where the script pays off. Wrap the calculation in a function and loop over fluid names. With the same temperatures, superheat, subcooling and efficiency, I got:

  • R134a: 2.01 / 10.17 bar, refrigerating effect 145.0 kJ/kg, COP 3.13, discharge 61.9 °C, 0.050 kg/s for 7.25 kW
  • R1234yf: 2.22 / 10.18 bar, refrigerating effect 110.6 kJ/kg, COP 2.98, discharge 49.7 °C, 0.066 kg/s for 7.25 kW
  • R290 (propane): 3.45 / 13.69 bar, refrigerating effect 273.5 kJ/kg, COP 3.09, discharge 61.1 °C, 0.027 kg/s for 7.25 kW

A few engineering observations come out of this immediately. R1234yf runs at almost the same pressures as R134a, which is why it was chosen as a near drop-in in car air conditioning, but it needs about 30 percent more mass flow for the same duty and gives a slightly lower COP under the same assumptions. Propane needs roughly half the mass flow because of its large latent heat, and runs at higher pressures. Keep in mind that this comparison holds the compressor efficiency fixed, which a real compressor will not do across fluids, and it says nothing about flammability, charge limits or volumetric capacity. Add CP.PropsSI("D", "T", T1, "P", pe, fluid) for suction density if you want to compare swept volume as well, which is often the more important number when selecting a compressor.

Common mistakes I see

  • Wrong fluid names. CoolProp uses names like "R134a", "R1234yf", "Propane" or "R290". If a name fails, list the available fluids with CP.get_global_param_string("fluids_list").
  • Gauge versus absolute pressure. CoolProp works in absolute pressure. Service gauges read gauge pressure, so add about 1.013 bar before comparing.
  • Reference state confusion. Absolute enthalpy values depend on the reference state, so your h values may not match a chart from another source. Differences in enthalpy, and therefore COP and duties, will match.
  • Specifying a two-phase state with T and P. Inside the dome, temperature and pressure are not independent. Use quality, or pressure with enthalpy.

Where to take it next

Once the basic cycle works, the natural extensions are pressure drops in the suction and discharge lines, a suction-line heat exchanger, a volumetric efficiency model for the compressor, or a sweep of condensing temperature to see how COP falls on a hot day. If you want a refresher on the heat exchanger side, see the earlier Heat Transfer post on this blog. All of it is just more PropsSI calls, and the script stays readable. For mixtures, CoolProp also supports predefined blends such as "R407C.mix" and "R410A", but check the documentation for how each one is modeled before trusting the results near the critical point.

If you enjoy working through the physics behind the flows inside heat exchangers and pipes, my book Memahami Turbulensi: Konsep, Persamaan, dan Praktik CFD dengan Perangkat Lunak Gratis on Google Play Books (written in Indonesian) goes deeper into turbulence and heat transfer, with hands-on chapters using free tools including OpenFOAM, Cantera and CoolProp. You can read more about it in this overview of the book. And if this kind of practical, open-source walkthrough is useful to you, follow the Bugo Center blog so you catch the next tutorial.

References

Related posts