matlab结合python的CoolProp库来进行热泵热循环仿真

news2024/11/28 3:23:58

前言

需要安装python，不同matlab版本需要下载对用的python版本！！！！，切记！！！！否则程序无法运行，下图是展示了matlab和python之间的版本对应
在这里插入图片描述
安装完python后，打开cmd(mac电脑终端) 输入 pip install CoolProp

实现步骤

多状态点和过程：包括压缩、冷凝、节流、蒸发等过程。
详细的能量平衡和效率计算：包括各个设备的效率（如压缩机、泵、涡轮等），以及整体循环的性能。
详细的状态点计算：使用 CoolProp 计算精确的热力学性质。
火用分析：计算火用损失和火用效率。
动态仿真：通过迭代方法模拟系统的动态行为

核心代码

% Parameters
T0 = 298.15; % Reference temperature (K)
P1 = 100; % Pressure at state point 1 (kPa)
T1 = 30 + 273.15; % Temperature at state point 1 (K)
P2 = 3000; % Pressure at state point 2 (kPa)
P3 = P2; % Pressure at state point 3 (kPa)
P4 = P1; % Pressure at state point 4 (kPa)
eff_pump = 0.85; % Efficiency of the pump
eff_turbine = 0.85; % Efficiency of the turbine
m_dot = 1; % Mass flow rate (kg/s)

% State 1: Liquid water at low pressure and temperature
h1 = py.CoolProp.CoolProp.PropsSI('H', 'P', P1*1000, 'T', T1, 'Water') / 1000; % Enthalpy (kJ/kg)
s1 = py.CoolProp.CoolProp.PropsSI('S', 'P', P1*1000, 'T', T1, 'Water') / 1000; % Entropy (kJ/kg.K)

% State 2: Liquid water at high pressure (after pump)
s2 = s1; % Isentropic process
h2s = py.CoolProp.CoolProp.PropsSI('H', 'P', P2*1000, 'S', s2*1000, 'Water') / 1000; % Isentropic enthalpy (kJ/kg)
h2 = h1 + (h2s - h1) / eff_pump; % Actual enthalpy after pump (kJ/kg)
T2 = py.CoolProp.CoolProp.PropsSI('T', 'P', P2*1000, 'H', h2*1000, 'Water') - 273.15; % Temperature (°C)

% State 3: High pressure steam (after boiler)
T3 = 450 + 273.15; % Temperature at state point 3 (K)
h3 = py.CoolProp.CoolProp.PropsSI('H', 'P', P2*1000, 'T', T3, 'Water') / 1000; % Enthalpy (kJ/kg)
s3 = py.CoolProp.CoolProp.PropsSI('S', 'P', P2*1000, 'T', T3, 'Water') / 1000; % Entropy (kJ/kg.K)

% State 4: Low pressure steam (after turbine)
s4 = s3; % Isentropic process
h4s = py.CoolProp.CoolProp.PropsSI('H', 'P', P1*1000, 'S', s4*1000, 'Water') / 1000; % Isentropic enthalpy (kJ/kg)
h4 = h3 - eff_turbine * (h3 - h4s); % Actual enthalpy after turbine (kJ/kg)
T4 = py.CoolProp.CoolProp.PropsSI('T', 'P', P1*1000, 'H', h4*1000, 'Water') - 273.15; % Temperature (°C)

% Energy balance and exergy analysis
Q_in = h3 - h2; % Heat input in the boiler (kJ/kg)
Q_out = h4 - h1; % Heat rejected in the condenser (kJ/kg)
W_net = (h3 - h4) - (h2 - h1); % Net work output (kJ/kg)
COP = Q_in / W_net; % Coefficient of Performance

% Display results
fprintf('State Points:\n');
fprintf('1: T = %.2f K, P = %.2f kPa, h = %.2f kJ/kg, s = %.2f kJ/kg.K\n', T1, P1, h1, s1);
fprintf('2: T = %.2f °C, P = %.2f kPa, h = %.2f kJ/kg, s = %.2f kJ/kg.K\n', T2, P2, h2, s2);
fprintf('3: T = %.2f °C, P = %.2f kPa, h = %.2f kJ/kg, s = %.2f kJ/kg.K\n', T3-273.15, P2, h3, s3);
fprintf('4: T = %.2f °C, P = %.2f kPa, h = %.2f kJ/kg, s = %.2f kJ/kg.K\n', T4, P1, h4, s4);

fprintf('\nCycle Performance:\n');
fprintf('Net Work Output (kJ/kg): %.2f\n', W_net);
fprintf('Heat Input (kJ/kg): %.2f\n', Q_in);
fprintf('Heat Rejected (kJ/kg): %.2f\n', Q_out);
fprintf('Coefficient of Performance (COP): %.2f\n', COP);

% Plotting T-s Diagram
figure;
hold on;
plot([s1 s2], [T1-273.15 T2], 'bo-'); % Process 1-2
plot([s2 s3], [T2 T3-273.15], 'ro-'); % Process 2-3
plot([s3 s4], [T3-273.15 T4], 'go-'); % Process 3-4
plot([s4 s1], [T4 T1-273.15], 'mo-'); % Process 4-1
xlabel('Entropy (kJ/kg.K)');
ylabel('Temperature (°C)');
title('T-s Diagram for Rankine Cycle');
legend('Process 1-2', 'Process 2-3', 'Process 3-4', 'Process 4-1');
grid on;
hold off;

% Total cycle performance
Q_in_total = Q_in * m_dot; % Total heat input (kW)
W_net_total = W_net * m_dot; % Total net work output (kW)
COP_real = Q_in_total / W_net_total; % Real COP

fprintf('\nTotal Cycle Performance:\n');
fprintf('Total Heat Input (kW): %.2f\n', Q_in_total);
fprintf('Total Net Work Output (kW): %.2f\n', W_net_total);
fprintf('Real Coefficient of Performance (COP): %.2f\n', COP_real);

% Exergy analysis
exergy_in = m_dot * (h3 - h2 - T0 * (s3 - s2)); % Exergy input (kW)
exergy_out = m_dot * (h4 - h1 - T0 * (s4 - s1)); % Exergy output (kW)
exergy_efficiency = exergy_out / exergy_in; % Exergy efficiency

fprintf('\nExergy Analysis:\n');
fprintf('Exergy Input (kW): %.2f\n', exergy_in);
fprintf('Exergy Output (kW): %.2f\n', exergy_out);
fprintf('Exergy Efficiency: %.2f\n', exergy_efficiency);

% Detailed component analysis
% Pump work
W_pump = h2 - h1; % Pump work (kJ/kg)
W_pump_total = W_pump * m_dot; % Total pump work (kW)
fprintf('\nPump Work:\n');
fprintf('Pump Work (kJ/kg): %.2f\n', W_pump);
fprintf('Total Pump Work (kW): %.2f\n', W_pump_total);

% Turbine work
W_turbine = h3 - h4; % Turbine work (kJ/kg)
W_turbine_total = W_turbine * m_dot; % Total turbine work (kW)
fprintf('\nTurbine Work:\n');
fprintf('Turbine Work (kJ/kg): %.2f\n', W_turbine);
fprintf('Total Turbine Work (kW): %.2f\n', W_turbine_total);

% Boiler heat input
Q_boiler = h3 - h2; % Boiler heat input (kJ/kg)
Q_boiler_total = Q_boiler * m_dot; % Total boiler heat input (kW)
fprintf('\nBoiler Heat Input:\n');
fprintf('Boiler Heat Input (kJ/kg): %.2f\n', Q_boiler);
fprintf('Total Boiler Heat Input (kW): %.2f\n', Q_boiler_total);
rejection (kW)
fprintf('\nCondenser Heat Rejection:\n');
fprintf('Condenser Heat Rejection (kJ/kg): %.2f\n', Q_condenser);
fprintf('Total Condenser Heat Rejection (kW): %.2f\n', Q_condenser_total);

% Efficiency calculations
thermal_efficiency = W_net_total / Q_in_total; % Thermal efficiency
fprintf('\nEfficiency Calculations:\n');
fprintf('Thermal Efficiency: %.2f%%\n', thermal_efficiency * 100);