Energy balance issue with Simscape Heat Exchanger model

10 views (last 30 days)
Hello,
I am trying to write a code of simple heat exchanger system based on the E-NTU method in simscape environment. I have a general code which does the heat transfer based on other questions (https://www.mathworks.com/matlabcentral/answers/index/?s_tid=gn_mlc_an) here. However, when I try to validate my model against the Simscape component of Heat Exchanger (TL-TL) it doesn't work. I am not even getting proper mass balance which I think should be the basic first step. I have worked long enough on this but still don't see the issue with the code. I am attaching my code and model for reference.
Any help would be greatly appreciated.
Thanks in advance!
component (Propagation = blocks) HX
% Heat Exchanger module
% Based on the epsilon-NTU method
nodes
A1 = foundation.thermal_liquid.thermal_liquid; % A1:top
B1 = foundation.thermal_liquid.thermal_liquid; % B1:top
B2 = foundation.thermal_liquid.thermal_liquid; % B2:bottom
A2 = foundation.thermal_liquid.thermal_liquid; % A2:bottom
end
outputs
heat = { 0.0, 'J/s' }; % H:bottom
end
parameters
U = {200, 'W/(m^2)/K'}; %Heat Transfer co-efficient
A = {10, 'm^2' }; % Total surface area
HX_length = {5, 'm' }; % Pipe length
HX_section = {0.01, 'm^2'}; % Cross-sectional area
K1 = {0.01, 'bar*s/kg'}; % Coefficient of pressure drop fluid 1
K2 = {0.01, 'bar*s/kg'}; % Coefficient of pressure drop fluid 2
dynamic_compressibility = simscape.enum.onoff.on; % Fluid dynamic compressibility
% 1 - on
% 0 - off
end
parameters(Access=protected)
mdot_min = { 0.0001, 'kg/s' }
end
variables (Access=protected)
% Through variables
mdot_A1 = { 0.01, 'kg/s' }; % Mass flow branch 1
mdot_B1 = { 0.01, 'kg/s' }; % Mass flow branch 2
mdot_A2 = { 0.01, 'kg/s' }; % Mass flow branch 1
mdot_B2 = { 0.01, 'kg/s' }; % Mass flow branch 2
Phi_A1 = {0, 'J/s'}; % Heat flow branch 1 port A
Phi_A2 = {0, 'J/s'}; % Heat flow branch 2 port A
Phi_B1 = {0, 'J/s'}; % Heat flow branch 1 port B
Phi_B2 = {0, 'J/s'}; % Heat flow branch 2 port B
rho_A1 = {1000, 'kg/m^3'}; % Density at port A1
rho_A2 = {1000, 'kg/m^3'}; % Density at port A2
rho_B1 = {1000, 'kg/m^3'}; % Density at port B1
rho_B2 = {1000, 'kg/m^3'}; % Density at port B2
T_A1 = {293.15, 'K'}; % Temperature at port A1
T_B1 = {293.15, 'K'}; % Temperature at port B1
T_A2 = {293.15, 'K'}; % Temperature at port A2
T_B2 = {293.15, 'K'}; % Temperature at port B2
rho1 = {1000, 'kg/m^3'}; % Density of fluid 1
rho2 = {1000, 'kg/m^3'}; % Density of fluid 2
u1 = {85, 'kJ/kg'};
u2 = {85, 'kJ/kg'};
% Internal nodes
u_A1 = {85, 'kJ/kg' }; % Specific internal energy at port A1
u_A2 = {85, 'kJ/kg' }; % Specific internal energy at port A2
u_B1 = {85, 'kJ/kg' }; % Specific internal energy at port B1
u_B2 = {85, 'kJ/kg' }; % Specific internal energy at port B2
Q = {1000, 'J/s'};
end
variables
% Internal nodes
p1 = {value = {1, 'bar'}, priority = priority.high}; % Pressure of the liquid volume 1
T1 = {value = {293.15, 'K'}, priority = priority.high}; % Temperature of the liquid volume 1
p2 = {value = {1, 'bar'}, priority = priority.high}; % Pressure of the liquid volume 2
T2 = {value = {293.15, 'K'}, priority = priority.high}; % Temperature of the liquid volume 2
end
branches
mdot_A1 : A1.mdot -> *;
mdot_B1 : B1.mdot -> *;
mdot_A2 : A2.mdot -> *;
mdot_B2 : B2.mdot -> *;
Phi_A1 : A1.Phi -> *;
Phi_A2 : A2.Phi -> *;
Phi_B1 : B1.Phi -> *;
Phi_B2 : B2.Phi -> *;
end
equations
let
% Across variables
p_A1 = A1.p;
p_A2 = A2.p;
p_B1 = B1.p;
p_B2 = B2.p;
% Domain parameters for look up table only
T1_TLU = A1.T_TLU;
p1_TLU = A1.p_TLU;
T2_TLU = A2.T_TLU;
p2_TLU = A2.p_TLU;
cp1_TLU = A1.cp_TLU;
cp2_TLU = A2.cp_TLU;
rho1_TLU = A1.rho_TLU;
rho2_TLU = A2.rho_TLU;
u1_TLU = A1.u_TLU;
u2_TLU = A2.u_TLU;
alpha1_TLU = A1.alpha_TLU;
alpha2_TLU = A2.alpha_TLU;
beta1_TLU = A1.beta_TLU;
beta2_TLU = A2.beta_TLU;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Upwind energy scheme
k1_cv = A1.k_cv;
k2_cv = A2.k_cv;
max_aspect_ratio_A1 = A1.max_aspect_ratio;
max_aspect_ratio_B1 = B1.max_aspect_ratio;
max_aspect_ratio_A2 = A2.max_aspect_ratio;
max_aspect_ratio_B2 = B2.max_aspect_ratio;
% Max length for conduction based on the max component aspect ratio (length/diameter)
max_conduction_length_A1 = max_aspect_ratio_A1 * sqrt(4*HX_section/pi);
max_conduction_length_B1 = max_aspect_ratio_B1 * sqrt(4*HX_section/pi);
max_conduction_length_A2 = max_aspect_ratio_A2 * sqrt(4*HX_section/pi);
max_conduction_length_B2 = max_aspect_ratio_B2 * sqrt(4*HX_section/pi);
% Thermal conductance coefficient
% Approximate conduction using specific internal energy differential
% instead of temperature differential
G_A1 = ...
if le(HX_length/2, max_conduction_length_A1), ...
k1_cv*HX_section/(HX_length/2) ...
else ...
k1_cv*HX_section/max_conduction_length_A1 ...
end;
G_B1 = ...
if le(HX_length/2, max_conduction_length_B1), ...
k1_cv*HX_section/(HX_length/2) ...
else ...
k1_cv*HX_section/max_conduction_length_B1 ...
end;
G_A2 = ...
if le(HX_length/2, max_conduction_length_A2), ...
k2_cv*HX_section/(HX_length/2) ...
else ...
k2_cv*HX_section/max_conduction_length_A2 ...
end;
G_B2 = ...
if le(HX_length/2, max_conduction_length_B2), ...
k2_cv*HX_section/(HX_length/2) ...
else ...
k2_cv*HX_section/max_conduction_length_B2 ...
end;
% Smooth absolute value of mass flow rate for energy flow rate calculations
% There can be non-zero energy flow even at zero mass flow due to conduction
mdot_abs_A1 = sqrt(mdot_A1^2 + 4*G_A1^2);
mdot_abs_B1 = sqrt(mdot_B1^2 + 4*G_B1^2);
mdot_abs_A2 = sqrt(mdot_A2^2 + 4*G_A2^2);
mdot_abs_B2 = sqrt(mdot_B2^2 + 4*G_B2^2);
% Smoothed step functions for energy flow rate during flow reversal
% Smoothing is based on conduction heat flow rate which dominates near zero flow
% and is negligible otherwise
step_plus_A1 = (1 + mdot_A1/mdot_abs_A1)/2;
step_plus_B1 = (1 + mdot_B1/mdot_abs_B1)/2;
step_minus_A1 = (1 - mdot_A1/mdot_abs_A1)/2;
step_minus_B1 = (1 - mdot_B1/mdot_abs_B1)/2;
step_plus_A2 = (1 + mdot_A2/mdot_abs_A2)/2;
step_plus_B2 = (1 + mdot_B2/mdot_abs_B2)/2;
step_minus_A2 = (1 - mdot_A2/mdot_abs_A2)/2;
step_minus_B2 = (1 - mdot_B2/mdot_abs_B2)/2;
% Upstream specific internal energy for inflow and outflow
u_in_A1 = tablelookup(T1_TLU, p1_TLU, u1_TLU, A1.T, p_A1, interpolation=linear, extrapolation=linear);
u_out_A1 = u1;
u_in_B1 = tablelookup(T1_TLU, p1_TLU, u1_TLU, B1.T, p_B1, interpolation=linear, extrapolation=linear);
u_out_B1 = u1;
u_in_A2 = tablelookup(T2_TLU, p2_TLU, u2_TLU, A2.T, p_A2, interpolation=linear, extrapolation=linear);
u_out_A2 = u2;
u_in_B2 = tablelookup(T2_TLU, p2_TLU, u2_TLU, B2.T, p_B2, interpolation=linear, extrapolation=linear);
u_out_B2 = u2;
% Flow work (modif rho_B1 = rho1 et rho_B2 = Rho2)
pv_A1 = mdot_A1/mdot_abs_A1 * p_A1/rho_A1;
pv_B1 = mdot_B1/mdot_abs_B1 * p_B1/rho_B1;
pv_A2 = mdot_A2/mdot_abs_A2 * p_A2/rho_A2;
pv_B2 = mdot_B2/mdot_abs_B2 * p_B2/rho_B2;
% pv_A1 = 0;
% pv_A2 = 0;
% pv_B1 = 0;
% pv_B2 = 0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Liquid properties table lookup upstream
cp1 = tablelookup(T1_TLU, p1_TLU, cp1_TLU, T1, p1, interpolation=linear, extrapolation=linear);
cp2 = tablelookup(T2_TLU, p2_TLU, cp2_TLU, T2, p2, interpolation=linear, extrapolation=linear);
beta1 = tablelookup(T1_TLU, p1_TLU, beta1_TLU, T1, p1, interpolation=linear, extrapolation=linear);
beta2 = tablelookup(T2_TLU, p2_TLU, beta2_TLU, T2, p2, interpolation=linear, extrapolation=linear);
alpha1 = tablelookup(T1_TLU, p1_TLU, alpha1_TLU, T1, p1, interpolation=linear, extrapolation=linear);
alpha2 = tablelookup(T2_TLU, p2_TLU, alpha2_TLU, T2, p2, interpolation=linear, extrapolation=linear);
% Partial derivatives of internal energy
% with respect to pressure and temperature at constant volume
volume = HX_length * HX_section;
h1 = u1 + p1/rho1;
dUdT1 = (cp1 - h1*alpha1) * rho1 * volume;
dUdp1 = (rho1*h1/beta1 - T1*alpha1) * volume;
h2 = u2 + p2/rho2;
dUdT2 = (cp2 - h2*alpha2) * rho2 * volume;
dUdp2 = (rho2*h2/beta2 - T2*alpha2) * volume;
% Specific heat at constant volume
cv_1 = cp1 - beta1 * alpha1^2 * T1 / rho1;
cv_2 = cp2 - beta2 * alpha2^2 * T2 / rho2;
% Calcul of heat capacity flow (important to avoid 0 flow)
C1 = (abs(mdot_A1) + mdot_min) / abs(abs(mdot_A1) + mdot_min) * (abs(mdot_A1) + mdot_min) * cp1;
C2 = (abs(mdot_A2) + mdot_min) / abs(abs(mdot_A2) + mdot_min) * (abs(mdot_A2) + mdot_min) * cp2;
% Calcul of the ratio of minimum /maximum
Cr = min(C1,C2)/max(C1,C2)
% Calcul of Number of Transfer Unit
NTU = U*A/min(C1,C2)
% Calcul of "efficiency" cross flow
% if flow_type == flow.direct
epsilon = if Cr == 1, NTU/(1+NTU) else (1-exp(-NTU*(1-Cr)))/(1-Cr*exp(-NTU*(1-Cr))) end;
% else
%
% epsilon == 1 - exp(-1/Cr*(1-exp(-Cr*NTU)));
% end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Upwind energy scheme
in
% Density table lookup
rho_A1 == tablelookup(T1_TLU, p1_TLU, rho1_TLU, T_A1, p_A1, interpolation=linear, extrapolation=linear);
rho_B1 == tablelookup(T1_TLU, p1_TLU, rho1_TLU, T_B1, p_B1, interpolation=linear, extrapolation=linear);
rho_A2 == tablelookup(T2_TLU, p2_TLU, rho2_TLU, T_A2, p_A2, interpolation=linear, extrapolation=linear);
rho_B2 == tablelookup(T2_TLU, p2_TLU, rho2_TLU, T_B2, p_B2, interpolation=linear, extrapolation=linear);
% Upwinded energy flow rate
% (internal energy advection + thermal conduction + flow work)
% Normalized by mass flow rate to improve scaling
Phi_A1/mdot_abs_A1 == step_plus_A1*u_in_A1 - step_minus_A1*u_out_A1 + pv_A1;
Phi_B1/mdot_abs_B1 == step_plus_B1*u_in_B1 - step_minus_B1*u_out_B1 + pv_B1;
Phi_A2/mdot_abs_A2 == step_plus_A2*u_in_A2 - step_minus_A2*u_out_A2 + pv_A2;
Phi_B2/mdot_abs_B2 == step_plus_B2*u_in_B2 - step_minus_B2*u_out_B2 + pv_B2;
% Upwind specific internal energy
u_A1 == step_plus_A1*u_in_A1 + step_minus_A1*u_out_A1;
u_B1 == step_plus_B1*u_in_B1 + step_minus_B1*u_out_B1;
u_A2 == step_plus_A2*u_in_A2 + step_minus_A2*u_out_A2;
u_B2 == step_plus_B2*u_in_B2 + step_minus_B2*u_out_B2;
% Solve for temperature corresponding to upwind specific internal energy
u_A1 == tablelookup(T1_TLU, p1_TLU, u1_TLU, T_A1, p_A1, interpolation=linear, extrapolation=linear);
u_B1 == tablelookup(T1_TLU, p1_TLU, u1_TLU, T_B1, p_B1, interpolation=linear, extrapolation=linear);
u_A2 == tablelookup(T2_TLU, p2_TLU, u2_TLU, T_A2, p_A2, interpolation=linear, extrapolation=linear);
u_B2 == tablelookup(T2_TLU, p2_TLU, u2_TLU, T_B2, p_B2, interpolation=linear, extrapolation=linear);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% table look up for solving
rho1 == tablelookup(T1_TLU, p1_TLU, rho1_TLU, T1, p1, interpolation=linear, extrapolation=linear);
rho2 == tablelookup(T2_TLU, p2_TLU, rho2_TLU, T2, p2, interpolation=linear, extrapolation=linear);
u1 == tablelookup(T1_TLU, p1_TLU, u1_TLU, T1, p1, interpolation=linear, extrapolation=linear);
u2 == tablelookup(T2_TLU, p2_TLU, u2_TLU, T2, p2, interpolation=linear, extrapolation=linear);
if dynamic_compressibility == simscape.enum.onoff.off
% Mass balance equation
mdot_A1 + mdot_B1 == 0;
mdot_A2 + mdot_B2 == 0;
% Energy balance equation for each fluid
Q == epsilon * min(C1,C2) * (T1 - T2);
-Q + Phi_A1 + Phi_B1 == T1.der * cv_1 * rho1* volume;
Q + Phi_A2 + Phi_B2 == T2.der * cv_2 * rho2* volume;
else % dynamic_compressibility == simscape.enum.onoff.on
% Mass balance equation
mdot_A1 + mdot_B1 == (p1.der/ beta1 - T1.der * alpha1) * rho1 * volume*0.00000001;
mdot_A2 + mdot_B2 == (p2.der/ beta2 - T2.der * alpha2) * rho2 * volume*0.00000001;
% Energy balance equation for each fluid
% Real heat exchanged (Q = epsilon * Q_max)
%Q == epsilon * min(C1,C2) * (T_B1 - T_B2);
Q == epsilon * min(C1,C2) * (T1 - T2);
-Q + Phi_A1 + Phi_B1 == p1.der * dUdp1 + T1.der * dUdT1;
Q + Phi_A2 + Phi_B2 == p2.der * dUdp2 + T2.der * dUdT2;
end
% Pressure loss equation for each fluid
p_A1 - p1 == K1 * mdot_A1;
p_B1 - p1 == K1 * mdot_B1;
p_A2 - p2 == K2 * mdot_A2;
p_B2 - p2 == K2 * mdot_B2;
% Sortie
heat == Q;
end
end

Accepted Answer

Rohan Kokate
Rohan Kokate on 5 Nov 2020
UPDATE: I resolved the issue with the energy balance. The heat transfer rate should be calculated as:
Q == epsilon * min(C1,C2) * (T_A1 - T_A2);
And also, I was using counter flow heat exchanger effectiveness formula instead of parallel flow heat exchanger effectiveness. That caused the unphyscial temperature profile of the fluids at outlet. The system and code are both working as expected now. I validated the code and system with the library component of Heat Exchanger (TL-TL).
Hope this helps anyone trying to do similar thing!

More Answers (0)

Community Treasure Hunt

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

Start Hunting!

Translated by