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;
Best Answer