Process Heat Integration

Version 1.0.1 (4.35 KB) by Julian
Calculates the minimum heating and cooling demand of an industrial process, as well as useful excess heat.
Updated 1 Mar 2024

View License

The function HeatIntegration.m can be used to calculate the minimum heating and cooling demand of an industrial process, as well as the amount of useful excess heat the process produces. The model is based on the work of Papoulias and Grossmann ( and I extended it to determine the amount of useful excess heat from the process.
Short Description of Functions
  • HeatIntegration() calculates the minimum heating demand and cooling demand, as well as useful excess heat. Input arguments are the process streams and latent process heats, as well as the minimum temperature difference for heat exchange and a threshold temperature above which excess heat is considered useful.
  • getEnthalpy() Calculates the enthalpy for a liquid or gaseous process stream. For liquid process streams, only constant heat capacity values are supported, while for gaseous process streams, constant heat capacity values and the DIPPR equation 107 are supported. Adding other heat capacity models is easy.
Steam Methane Reforming
This is a simplified version of a steam methane reforming process, assuming methane and liquid water enter the process at 300 K and 10 bar, are heated to 1,000 K, and then react in an isothermal reactor to reach chemical equilibrium. The product mixture is cooled to 700 K to undergo water gas shift reaction (not included).
% Each species is specified by [heat of formation at 298.15 K in kJ/mol,
% [Parameters A, B, C, D, and E for ideal gas heat capacity according to DIPPR equation 107 in kJ/(mol-K)]
% or a constant heat capacity value
properties = struct("Hydrogen", struct("h_0", 0, "Cp_model", "DIPPR107", "Cp_params", [0.027617, 0.00956, 2466, 0.00376, 567.6]), ...
"Methane", struct("h_0", -74.6, "Cp_model", "DIPPR107", "Cp_params", [0.033298, 0.079933, 2086.9, 0.041602, 991.96]), ...
"Carbon_Monoxide", struct("h_0", -110.53, "Cp_model", "DIPPR107", "Cp_params", [0.029108, 0.008773, 3085.1, 0.0084553, 1538.2]), ...
"Water_g", struct("h_0", -241.81, "Cp_model", "DIPPR107", "Cp_params", [0.033363, 0.02679, 2610.5, 0.008896, 1169]), ...
"Water_l", struct("h_0", -285.83, "Cp_model", "const", "Cp_params", 0.07538));
% Inlet stream 1
components_1 = [struct("Name", 'Methane', "x", 1, "Props", properties.Methane)];
stream_1 = struct("n", 1, "T_in", 298.15, "T_out", 1000, "Components", components_1, "State", "gas");
% Inlet stream 2
% Assume 10 bar pressure (boiling point of water is 179.91°C) and 1000 K for the reactor
components_2_l = [struct("Name", 'Water', "x", 1, "Props", properties.Water_l)];
stream_2_l = struct("n", 1, "T_in", 298.15, "T_out", 453.06, "Components", components_2_l, "State", "liquid");
components_2_g = [struct("Name", 'Water', "x", 1, "Props", properties.Water_g)];
stream_2_g = struct("n", 1, "T_in", 453.06, "T_out", 1000, "Components", components_2_g, "State", "gas");
% Outlet stream
% ...the equilibrium conversion is 0.41 under these conditions
X = 0.41; % [-]
% ...calculate stream composition
n_CH4_3 = stream_1.n*(1-X);
n_H2O_3 = stream_2_g.n*(1-X);
n_CO_3 = stream_1.n*X;
n_H2_3 = 3*stream_1.n*X;
n_tot_3 = n_CH4_3 + n_H2O_3 + n_CO_3 + n_H2_3;
components_3 = [struct("Name", 'Methane', "x", n_CH4_3/n_tot_3, "Props", properties.Methane), struct("Name", 'Water', "x", n_H2O_3/n_tot_3, "Props", properties.Water_g), struct("Name", 'Carbon_Monoxide', "x", n_CO_3/n_tot_3, "Props", properties.Carbon_Monoxide), struct("Name", 'Hydrogen', "x", n_H2_3/n_tot_3, "Props", properties.Hydrogen)];
stream_3 = struct("n", n_tot_3, "T_in", 1000, "T_out", 700, "Components", components_3, "State", "gas");
streams = [stream_1, stream_2_l, stream_2_g, stream_3];
% Water evaporation at 179.91°C
Q_vap = getEnthalpy(struct("n", stream_2_g.n, "T", 453.06, "Components", components_2_g, "State", "gas")) - getEnthalpy(struct("n", stream_2_l.n, "T", 453.06, "Components", components_2_l, "State", "liquid"));
% Heat of reaction
Q_R = getEnthalpy(struct("n", stream_3.n, "T", 1000, "Components", components_3, "State", "gas")) - getEnthalpy(struct("n", stream_2_g.n, "T", 1000, "Components", components_2_g, "State", "gas")) - getEnthalpy(struct("n", stream_1.n, "T", 1000, "Components", components_1, "State", "gas"));
latentHeats = [struct("Q", Q_vap, "T", 453.06), struct("Q", -Q_R, "T", 1000)];
dT_min = 10; % [K]
T_use = 373.15; % [K]
[Q_H, Q_C, Q_use] = HeatIntegration(streams, latentHeats, dT_min, T_use);
Where to Find Data
  • Standard enthalpy values for gaseous species can be obtained from, DIPPR 107 parameters for the heat capacities from Aspen Plus or databases.
  • Heat capacity and standard enthalpy values for liquids are reported in the NIST Chemistry WebBook.
Feel free to reach out to me under
The functions are provided 'as is', without warranty of any kind, express or implied.

Cite As

Julian (2024). Process Heat Integration (, MATLAB Central File Exchange. Retrieved .

MATLAB Release Compatibility
Created with R2023b
Compatible with any release
Platform Compatibility
Windows macOS Linux
Tags Add Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!
Version Published Release Notes

Added Image