# Problem with Custom Heat Exchanger (G-TL) model

25 views (last 30 days)
Rohan Kokate on 18 Nov 2020
Commented: Rohan Kokate on 21 Apr 2022
Hi,
I am trying to write a code for the heat exchanger (G-TL) component. Earlier, I did a similar thing for the heat exchanger (TL-TL) and was successful. But it seems, with gas components this becomes really challenging. I used the source code for the pipe (G) and two_port_steady (G) components to write my code. I have the basic structure but when I try to run the model, I am getting the following error:
• Cannot solve for one or more variables, including dynamic variable derivatives: Time derivative of 'Simscape_Component.T1' (Temperature of the liquid volume 1)
• Cannot solve for one or more variables, including dynamic variable derivatives: Time derivative of 'Simscape_Component.T2' (Temperature of the liquid volume 2)
• Solver encountered an error while simulating model 'test_model_customHX_G_1' at time 0 and cannot continue. Please check the model for errors.
• Solver could not solve the system of algebraic equations because a singular iteration matrix was encountered. Consider providing more accurate initial conditions. If the problem persists, check the model structure and values of parameters.
Attaching the model and the code below: component (Propagation = blocks) HX_gas
nodes
A1 = foundation.gas.gas; % A1:top
B1 = foundation.gas.gas; % 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 coefficient
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' };
volume = HX_section*HX_length; % Pipe volume
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
% h_A1 = {420, 'kJ/kg' }; % Specific enthalpy at port A
% h_B1 = {420, 'kJ/kg' }; % Specific enthalpy at port A
Q = {1000, 'J/s'}; % Total heat trasnfer rate (absorbed/disipated)
p_A1 = {0.1, 'MPa' }; % Pressure at port A
p_B1 = {0.1, 'MPa' }; % Pressure at port B
end
variables
% Default initial conditions
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
variables (Access=protected, ExternalAccess=none)
rho_in_A1 = {1.2, 'kg/m^3'}; % Density for inflow at port A
rho_in_B1 = {1.2, 'kg/m^3'}; % Density for inflow at port B
ht_out = {420, 'kJ/kg' }; % Specific total enthalpy for outflow
end
% Parameter groups
annotations
UILayout = [UIGroup('physmod:simscape:library:tabs:Geometry', ...
HX_length, HX_section,A)
UIGroup('physmod:simscape:library:tabs:FrictionAndHeatTransfer', ...
U,K1,K2,dynamic_compressibility)]
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
% Domain parameters for look up table only
T2_TLU = A2.T_TLU;
p2_TLU = A2.p_TLU;
cp2_TLU = A2.cp_TLU;
rho2_TLU = A2.rho_TLU;
u2_TLU = A2.u_TLU;
alpha2_TLU = A2.alpha_TLU;
beta2_TLU = A2.beta_TLU;
p_A2 = A2.p;
p_B2 = B2.p;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Upwind energy scheme
k2_cv = A2.k_cv;
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_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
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_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_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_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
pv_A2 = mdot_A2/mdot_abs_A2 * p_A2/rho_A2;
pv_B2 = mdot_B2/mdot_abs_B2 * p_B2/rho_B2;
% GAS SIDE
% Pressure and temperature for inflow
p_in_A1 = A1.p;
p_in_B1 = B1.p;
T_in_A1 = A1.T;
T_in_B1 = B1.T;
% Domain parameters
gas_spec = A1.gas_spec;
R = A1.R;
Z = A1.Z;
T1_ref = A1.T_ref;
h1_ref = A1.h_ref;
cp1_ref = A1.cp_ref;
cv1_ref = A1.cv_ref;
T1_TLU1 = A1.T_TLU1;
h1_TLU1 = A1.h_TLU1;
a1_TLU1 = A1.a_TLU1;
T1_TLU2 = A1.T_TLU2;
p1_TLU2 = A1.p_TLU2;
log_T1_TLU2 = A1.log_T_TLU2;
log_p1_TLU2 = A1.log_p_TLU2;
log_rho1_TLU2 = A1.log_rho_TLU2;
h1_TLU2 = A1.h_TLU2;
a1_TLU2 = A1.a_TLU2;
Mach1_rev = A1.Mach_rev
T1_unit = A1.T_unit;
p1_unit = A1.p_unit;
rho1_unit = A1.rho_unit;
log_ZR = A1.log_ZR;
% Log temperature
log_T_in_A1 = log(T_in_A1/T1_unit);
log_T_in_B1 = log(T_in_B1/T1_unit);
% Log pressure
log_p_in_A1 = log(p_in_A1/p1_unit);
log_p_in_B1 = log(p_in_B1/p1_unit);
% Log density
log_rho_A1 = log(rho_A1 /rho1_unit);
log_rho_B1 = log(rho_B1 /rho1_unit);
log_rho_in_A1 = log(rho_in_A1/rho1_unit);
log_rho_in_B1 = log(rho_in_B1/rho1_unit);
% Specific enthalpy for inflow
[h_in_A1, h_in_B1] = ...
if gas_spec == 1, ...
h1_ref + cp1_ref*(T_in_A1 - T1_ref); ...
h1_ref + cp1_ref*(T_in_B1 - T1_ref) ...
elseif gas_spec == 2, ...
tablelookup(T1_TLU1, h1_TLU1, T_in_A1, interpolation=linear, extrapolation=linear); ...
tablelookup(T1_TLU1, h1_TLU1, T_in_B1, interpolation=linear, extrapolation=linear) ...
else ...
tablelookup(T1_TLU2, p1_TLU2, h1_TLU2, T_in_A1, p_in_A1, interpolation=linear, extrapolation=linear); ...
tablelookup(T1_TLU2, p1_TLU2, h1_TLU2, T_in_B1, p_in_B1, interpolation=linear, extrapolation=linear) ...
end;
% Specific total enthalpy for inflow
% Take absolute value of density to eliminate nonphysical solutions
ht_in_A1 = h_in_A1 + (mdot_A1/HX_section)^2 / rho_in_A1 / abs(rho_in_A1) / 2;
ht_in_B1 = h_in_B1 + (mdot_B1/HX_section)^2 / rho_in_B1 / abs(rho_in_B1) / 2;
% Specific total enthalpy at ports A and B
% Take absolute value of density to eliminate nonphysical solutions
ht_A1 = h_A1 + (mdot_A1/HX_section)^2 / rho_A1 / abs(rho_A1) / 2;
ht_B1 = h_B1 + (mdot_B1/HX_section)^2 / rho_B1 / abs(rho_B1) / 2;
% Speed of sound
[a_in_A1, a_in_B1] = ...
if gas_spec == 1, ...
sqrt(cp1_ref / cv1_ref * Z * R * T_in_A1); ...
sqrt(cp1_ref / cv1_ref * Z * R * T_in_B1) ...
elseif gas_spec == 2, ...
tablelookup(T1_TLU1, a1_TLU1, T_in_A1, interpolation=linear, extrapolation=linear); ...
tablelookup(T1_TLU1, a1_TLU1, T_in_B1, interpolation=linear, extrapolation=linear) ...
else ...
tablelookup(T1_TLU2, p1_TLU2, a1_TLU2, T_in_A1, p_in_A1, interpolation=linear, extrapolation=linear); ...
tablelookup(T1_TLU2, p1_TLU2, a1_TLU2, T_in_B1, p_in_B1, interpolation=linear, extrapolation=linear) ...
end;
% Smoothing threshold for flow reversal
mdot_rev_A1 = Mach1_rev * rho_in_A1 * a_in_A1 * HX_section;
mdot_rev_B1 = Mach1_rev * rho_in_B1 * a_in_B1 * HX_section;
% Smooth absolute value of mass flow rate
mdot_abs_A1 = sqrt(mdot_A1^2 + mdot_rev_A1^2);
mdot_abs_B1 = sqrt(mdot_B1^2 + mdot_rev_B1^2);
% Smooth step functions for energy flow rate during flow reversal
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;
% Average mass flow rate from port A to port B
mdot1 = (mdot_A1 - mdot_B1)/2;
mdot1_abs = sqrt(mdot1^2 + mdot_rev_A1^2/2 + mdot_rev_B1^2/2);
step_plus1 = (1 + mdot1/mdot1_abs)/2;
step_minus1 = (1 - mdot1/mdot1_abs)/2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Liquid properties table lookup upstream
cp2 = tablelookup(T2_TLU, p2_TLU, cp2_TLU, T2, p2, interpolation=linear, extrapolation=linear);
beta2 = tablelookup(T2_TLU, p2_TLU, beta2_TLU, T2, p2, 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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
log_T_A1 = simscape.function.logProtected(T_A1/A1.T_unit, 1);
log_T_B1 = simscape.function.logProtected(T_B1/A1.T_unit, 1);
log_T1 = simscape.function.logProtected(T1 /A1.T_unit, 1);
% Log pressure
log_p_A1 = simscape.function.logProtected(p_A1/A1.p_unit, 1);
log_p_B1 = simscape.function.logProtected(p_B1/A1.p_unit, 1);
log_p1 = simscape.function.logProtected(p1/A1.p_unit, 1);
% Specific heat at constant pressure (for perfect and semiperfect gas)
cp1 = ...
if A1.gas_spec == foundation.enum.gas_spec.perfect_gas, ...
A1.cp_ref ...
elseif A1.gas_spec == foundation.enum.gas_spec.semiperfect_gas, ...
tablelookup(A1.T_TLU1, A1.cp_TLU1, T1, interpolation = linear, extrapolation = linear) ...
else ... % A.gas_spec == foundation.enum.gas_spec.real_gas
{1, 'kJ/(kg*K)'} ...
end;
% Derivative of density and density * internal energy
% with respect to pressure and temperature (for real gas)
% Use log-space to improve accuracy
[drho_dp1, drho_dT1, drhou_dp1, drhou_dT1] = ...
if A1.gas_spec == foundation.enum.gas_spec.real_gas, ...
exp(tablelookup(A1.log_T_TLU2, A1.log_p_TLU2, A1.log_drho_dp_TLU2, log_T1, log_p1, interpolation = linear, extrapolation = linear)) * A1.drho_dp_unit; ...
-exp(tablelookup(A1.log_T_TLU2, A1.log_p_TLU2, A1.log_drho_dT_TLU2, log_T1, log_p1, interpolation = linear, extrapolation = linear)) * A1.drho_dT_unit; ...
tablelookup(A1.T_TLU2, A1.p_TLU2, A1.drhou_dp_TLU2, T1, p1, interpolation = linear, extrapolation = linear); ...
tablelookup(A1.T_TLU2, A1.p_TLU2, A1.drhou_dT_TLU2, T1, p1, interpolation = linear, extrapolation = linear) ...
else ...
A1.drho_dp_unit; ...
A1.drho_dT_unit; ...
1; ...
{1, 'kJ/(m^3*K)'} ...
end;
% Partial derivatives of mass and internal energy
% with respect to pressure and temperature at constant volume
[dMdp1, dMdT1, dUdp1, dUdT1] = ...
if A1.gas_spec ~= foundation.enum.gas_spec.real_gas, ...
volume * rho1 / p1; ...
-volume * rho1 / T1; ...
volume * (h1 / (A1.Z * A1.R * T1) - 1); ...
volume * rho1 * (cp1 - h1/T1) ...
else ...
volume * drho_dp1; ...
volume * drho_dT1; ...
volume * drhou_dp1; ...
volume * drhou_dT1 ...
end;
% Thermal equation of state
% Use log-space to improve accuracy
[rho_A1, rho_B1] = ...
if A1.gas_spec ~= foundation.enum.gas_spec.real_gas, ...
exp(log_p_A1 - A1.log_ZR - log_T_A1) * A1.rho_unit; ...
exp(log_p_B1 - A1.log_ZR - log_T_B1) * A1.rho_unit ...
else ...
exp(tablelookup(A1.log_T_TLU2, A1.log_p_TLU2, A1.log_rho_TLU2, log_T_A1, log_p_A1, interpolation = linear, extrapolation = linear)) * A1.rho_unit; ...
exp(tablelookup(A1.log_T_TLU2, A1.log_p_TLU2, A1.log_rho_TLU2, log_T_B1, log_p_B1, interpolation = linear, extrapolation = linear)) * A1.rho_unit ...
end;
% Caloric equation of state
[h_A1, h_B1] = ...
if A1.gas_spec == foundation.enum.gas_spec.perfect_gas, ...
A1.h_ref + A1.cp_ref*(T_A1 - A1.T_ref); ...
A1.h_ref + A1.cp_ref*(T_B1 - A1.T_ref) ...
elseif A1.gas_spec == foundation.enum.gas_spec.semiperfect_gas, ...
tablelookup(A1.T_TLU1, A1.h_TLU1, T_A1, interpolation = linear, extrapolation = linear); ...
tablelookup(A1.T_TLU1, A1.h_TLU1, T_B1, interpolation = linear, extrapolation = linear) ...
else ... % A.gas_spec == foundation.enum.gas_spec.real_gas
tablelookup(A1.T_TLU2, A1.p_TLU2, A1.h_TLU2, T_A1, p_A1, interpolation = linear, extrapolation = linear); ...
tablelookup(A1.T_TLU2, A1.p_TLU2, A1.h_TLU2, T_B1, p_B1, interpolation = linear, extrapolation = linear) ...
end;
% Speed of sound
T_A1_limited = if ge(T_A1, A1.T_min), T_A1 else A1.T_min end;
T_B1_limited = if ge(T_B1, A1.T_min), T_B1 else A1.T_min end;
[a_A1, a_B1] = ...
if A1.gas_spec == foundation.enum.gas_spec.perfect_gas, ...
sqrt(A1.cp_ref / A1.cv_ref * A1.Z * A1.R * T_A1_limited); ...
sqrt(A1.cp_ref / A1.cv_ref * A1.Z * A1.R * T_B1_limited) ...
elseif A1.gas_spec == foundation.enum.gas_spec.semiperfect_gas, ...
tablelookup(A1.T_TLU1, A1.a_TLU1, T_A1_limited, interpolation = linear, extrapolation = linear); ...
tablelookup(A1.T_TLU1, A1.a_TLU1, T_B1_limited, interpolation = linear, extrapolation = linear) ...
else ... % A.gas_spec == foundation.enum.gas_spec.real_gas
tablelookup(A1.T_TLU2, A1.p_TLU2, A1.a_TLU2, T_A1_limited, p_A1, interpolation = linear, extrapolation = linear); ...
tablelookup(A1.T_TLU2, A1.p_TLU2, A1.a_TLU2, T_B1_limited, p_B1, interpolation = linear, extrapolation = linear) ...
end;
volume = HX_length * HX_section;
h1 = ...
if A1.gas_spec == foundation.enum.gas_spec.perfect_gas, ...
A1.h_ref + A1.cp_ref*(T1 - A1.T_ref) ...
elseif A1.gas_spec == foundation.enum.gas_spec.semiperfect_gas, ...
tablelookup(A1.T_TLU1, A1.h_TLU1, T1, interpolation = linear, extrapolation = linear) ...
else ... % A.gas_spec == foundation.enum.gas_spec.real_gas
tablelookup(A1.T_TLU2, A1.p_TLU2, A1.h_TLU2, T1, p1, interpolation = linear, extrapolation = linear) ...
end; % Specific enthalpy of gas 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;
%cp based on TLU again
% Calculation 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;
% Calculation of the thermal capacity ratio(Cr = C_min/C_max)
Cr = min(C1,C2)/max(C1,C2)
% Calculation of Number of Transfer Units (NTU)
NTU = U*(A/2)/min(C1,C2)
epsilon = (1-exp(-NTU*(1+Cr)))/(1+Cr); % parallel flow heat exchanger
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Upwind energy scheme
in
% Density table lookup
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_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_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_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
rho2 == tablelookup(T2_TLU, p2_TLU, rho2_TLU, T2, p2, interpolation=linear, extrapolation=linear);
u2 == tablelookup(T2_TLU, p2_TLU, u2_TLU, T2, p2, interpolation=linear, extrapolation=linear);
% Upwind energy flow rate
% Normalized by mass flow rate to improve scaling
Phi_A1/mdot_abs_A1 == step_plus_A1*ht_in_A1 - step_minus_A1*ht_out;
Phi_B1/mdot_abs_B1 == step_plus_B1*ht_in_B1 - step_minus_B1*ht_out;
% Thermal equation of state
% Use log-space to improve accuracy
if (gas_spec == 1) || (gas_spec == 2)
log_rho_A1 == log_p_A1 - log_ZR - log_T_A1;
log_rho_B1 == log_p_B1 - log_ZR - log_T_B1;
log_rho_in_A1 == log_p_in_A1 - log_ZR - log_T_in_A1;
log_rho_in_B1 == log_p_in_B1 - log_ZR - log_T_in_B1;
else
log_rho_A1 == tablelookup(log_T1_TLU2, log_p1_TLU2, log_rho1_TLU2, log_T_A1, log_p_A1, interpolation=linear, extrapolation=linear);
log_rho_B1 == tablelookup(log_T1_TLU2, log_p1_TLU2, log_rho1_TLU2, log_T_B1, log_p_B1, interpolation=linear, extrapolation=linear);
log_rho_in_A1 == tablelookup(log_T1_TLU2, log_p1_TLU2, log_rho1_TLU2, log_T_in_A1, log_p_in_A1, interpolation=linear, extrapolation=linear);
log_rho_in_B1 == tablelookup(log_T1_TLU2, log_p1_TLU2, log_rho1_TLU2, log_T_in_B1, log_p_in_B1, interpolation=linear, extrapolation=linear);
end
% Caloric equation of state
if gas_spec == 1
h_A1 == h1_ref + cp1_ref*(T_A1 - T1_ref);
h_B1 == h1_ref + cp1_ref*(T_B1 - T1_ref);
elseif gas_spec == 2
h_A1 == tablelookup(T1_TLU1, h1_TLU1, T_A1, interpolation=linear, extrapolation=linear);
h_B1 == tablelookup(T1_TLU1, h1_TLU1, T_B1, interpolation=linear, extrapolation=linear);
else
h_A1 == tablelookup(T1_TLU2, p1_TLU2, h1_TLU2, T_A1, p_A1, interpolation=linear, extrapolation=linear);
h_B1 == tablelookup(T1_TLU2, p1_TLU2, h1_TLU2, T_B1, p_B1, interpolation=linear, extrapolation=linear);
end
% Specific total enthalpy at ports A or B is equal to the outflow value
% depending on flow direction
step_plus1*ht_B1 + step_minus1*ht_A1 == ht_out;
% Thermal equation of state
% Use log-space to improve accuracy
rho1 == ...
if A1.gas_spec ~= foundation.enum.gas_spec.real_gas, ...
exp(log_p1 - A1.log_ZR - log_T1) * A1.rho_unit ...
else ...
exp(tablelookup(A1.log_T_TLU2, A1.log_p_TLU2, A1.log_rho_TLU2, log_T1, log_p1, interpolation = linear, extrapolation = linear)) * A1.rho_unit ...
end; % Density of gas volume
if dynamic_compressibility == simscape.enum.onoff.off
% Mass balance equation
mdot_A1 + mdot_B1 == der(p1)*dMdp1 + der(T1)*dMdT1;
mdot_A2 + mdot_B2 == 0;
% Energy balance equation for each fluid (Q_1 = - Q_2 = epsilon * Q_max)
Q == epsilon * min(C1,C2) * (T_A1 - T_A2);
Q + Phi_A1 + Phi_B1 == der(p1)*dUdp1 + der(T1)*dUdT1;
T2.der * cv_2 * rho2* volume == -Q + Phi_A2 + Phi_B2;
else % dynamic_compressibility == simscape.enum.onoff.on
% Mass balance equation
mdot_A1 + mdot_B1 == der(p1)*dMdp1 + der(T1)*dMdT1;
mdot_A2 + mdot_B2 == (p2.der/ beta2 - T2.der * alpha2) * rho2 * volume;
% Energy balance equation for each fluid
Q == epsilon * min(C1,C2) * (T_A1 - T_A2);
Q + Phi_A1 + Phi_B1 == der(p1)*dUdp1 + der(T1)*dUdT1;
-Q + Phi_A2 + Phi_B2 == p2.der * dUdp2 + T2.der * dUdT2;
end
% p_A1 - p1 == K1 * mdot_A1;
% p_B1 - p1 == K1 * mdot_B1;
p_A2 - p2 == K2 * mdot_A2;
p_B2 - p2 == K2 * mdot_B2;
% % Adiabatic process for each half of pipe
h_A1 + (mdot_A1/HX_section/rho_A1)^2/2 == h1 + (mdot_A1/HX_section/rho1)^2/2;
h_B1 + (mdot_B1/HX_section/rho_B1)^2/2 == h1 + (mdot_B1/HX_section/rho1)^2/2;
end
end
end
I don't understand what am I doing wrong. Any help would be greatly appreciated.
Regards,
Rohan Kokate.
Rohan Kokate on 21 Apr 2022
Hi Maria, try to look at the the Pipe(TL) and Pipe(G) library components. They have the source codes available. You will need two differents sets of conservation equations (one for TL and one for G). The equation for Q [Q == epsilon * min(C1,C2) * (T_A1 - T_A2)] will be the connection between the two sets of conservation equations. What I have tried to do in the code above is have the conservation equations together in one section which causes the problem. Hope that helps.

Yifeng Tang on 14 Jun 2021
There is a G-TL heat exchanger block available in Simscape Fluids, starting in R2019a. The source code is not visible since it's an add-on library block, so the education value isn't there. If you have access to Simscape Fluids, it may solve your modeling problem.
##### 2 CommentsShow 1 older commentHide 1 older comment
Maria Coutinho on 21 Apr 2022
Can you share the solution? I'm having the same problem! Thanks

### Categories

Find more on Sources in Help Center and File Exchange

R2020b

### Community Treasure Hunt

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

Start Hunting!