Problem with Custom Heat Exchanger (G-TL) model

9 views (last 30 days)
Rohan Kokate on 18 Nov 2020
Commented: Rohan Kokate on 29 Jul 2021
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.
Yifeng Tang on 27 Apr 2021
Hey Rohan,
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.
If you are still interested, I can take a look at the custom block you wrote. Let me know by comment if this is still of interest.

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.
Rohan Kokate on 29 Jul 2021
I figured out the issue. Thanks for your suggestions!

R2020b

Community Treasure Hunt

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

Start Hunting!