MATLAB Answers

Build a jet pump model in Simscape's thermal liquid domain.

9 views (last 30 days)
Andreas
Andreas on 13 Oct 2015
Commented: Andreas on 9 Jan 2020
Hello everyone,
we are currently working with the thermal liquid domain in Simscape and try to implement a model of a jet pump.
There is no real workaround by using Standard elements and lookup tables, so we are looking forward to port the "Jet pump" block from the simhydraulics domain: http://de.mathworks.com/help/physmod/hydro/ref/jetpump.html
ssc-file working now, we had a negative term in sqrt before, so initialisation would fail, no matter there is no flow reversal (and thus negative mdot) in a Jet pump. Code:
component jet_pump
% Thermal Liquid Jet Pump
% Defines a thermal liquid domain component with external nodes A, B and C
% located on the left, left and right hand side of the block respectively.
% Also defines associated through variables as well as thermodynamic
% properties and flux scheme. The properties are computed using table
% lookup functions.
% Modification of the simscape two_port.ssc-file
% Goal: Jet pump model
inputs
Anozzle = { 1e-4, 'm^2' }; % ctrl:left
end
nodes
A = foundation.thermal_liquid.thermal_liquid; % A:left
B = foundation.thermal_liquid.thermal_liquid; % B:left
C = foundation.thermal_liquid.thermal_liquid; % C:right
end
parameters
length = { 1e-1, 'm' }; % Characteristic longitudinal length
area = { 1e-2, 'm^2' }; % Pipe cross-sectional area
An_min = { 1e-10, 'm^2' }; % Minimum nozzle area
Ath = { 1e-4, 'm^2' }; % Throat area
Ad = { 1e-3, 'm^2' }; % Diffuser outlet area
Kn = { 0.05, '1' }; % Nozzle hydraulic loss coefficient
Ken = { 0.005, '1' }; % Throat entry hydraulic loss coefficient
Kth = { 0.1, '1' }; % Throat hydraulic loss coefficient
Kdi = { 0.1, '1' }; % Diffuser hydraulic loss coefficient
end
variables
% Component variables
mdot_A = { 0, 'kg/s' }; % Mass flow rate into A
mdot_B = { 0, 'kg/s' }; % Mass flow rate into B
mdot_C = { 0, 'kg/s' }; % Mass flow rate into C
Phi_A = { 0, 'J/s' }; % Thermal flux through A
Phi_B = { 0, 'J/s' }; % Thermal flux through B
Phi_C = { 0, 'J/s' }; % Thermal flux through C
p = { 1.01325,'bar' }; % Pressure
end
variables(Conversion=absolute)
T = { 293.15, 'K' }; % Temperature
end
variables(Access = protected)
Phi_convection_A= { 0, 'J/s' }; % Convective part of thermal flux through A
Phi_convection_B= { 0, 'J/s' }; % Convective part of thermal flux through B
Phi_convection_C= { 0, 'J/s' }; % Convective part of thermal flux through C
% Fluid properties - Default values are for water at p = 1 atmosphere and T = 293.15 Kelvin
%%Internal energy
u = { 84, 'J/g' }; % Internal energy
u_A = { 84, 'J/g' }; % Internal energy at A
u_B = { 84, 'J/g' }; % Internal energy at B
u_C = { 84, 'J/g' }; % Internal energy at C
%%Density
rho = {998.2, 'kg/m^3' }; % Density
rho_A = {998.2, 'kg/m^3' }; % Density at A
rho_B = {998.2, 'kg/m^3' }; % Density at B
rho_C = {998.2, 'kg/m^3' }; % Density at C
%%Specific heat at constant pressure, isobaric thermal expansion, and conductivity
cp = { 4.16, 'J/(g*K)' }; % Specific heat at constant pressure
alpha = { -2.0691e-04, '1/K' }; % Isobaric thermal expansion
k = { 598.5, 'mW/(m*K)' }; % Thermal conductivity
%%Isothermal bulk modulus and viscosity
beta = { 2.1791, 'GPa' }; % Isothermal bulk modulus
nu = { 1, 'mm^2/s' }; % Viscosity
end
branches
mdot_A : A.mdot -> *;
mdot_B : B.mdot -> *;
mdot_C : C.mdot -> *;
Phi_A : A.Phi -> *;
Phi_B : B.Phi -> *;
Phi_C : C.Phi -> *;
end
equations
let
% The nozzle area must be larger than An_min
An = if Anozzle < An_min, An_min else Anozzle end;
% Pressures and temperatures at fluid boundaries
p_A = A.p;
T_A = A.T;
p_B = B.p;
T_B = B.T;
p_C = C.p;
T_C = C.T;
% Thermal conductance for each part of the source
Gth = if k*(area/2)/(length/2) <= A.G_min, A.G_min else k*(area/2)/(length/2) end;
% Thermal energy equation sources and sinks
p_dv_A = A.p * mdot_A * if gt(mdot_A, 0), 1/rho_A - 1/rho else 1/rho - 1/rho end;
p_dv_B = B.p * mdot_B * if gt(mdot_B, 0), 1/rho_B - 1/rho else 1/rho - 1/rho end;
p_dv_C = C.p * mdot_C * if gt(mdot_C, 0), 1/rho_C - 1/rho else 1/rho - 1/rho end;
% Internal energy at the internal temperature and corresponding to boundary pressure
u_out_A = tablelookup(A.T_TLU, A.p_TLU, A.u_TLU, T, p_A, interpolation = linear, extrapolation = nearest);
u_out_B = tablelookup(B.T_TLU, B.p_TLU, A.u_TLU, T, p_B, interpolation = linear, extrapolation = nearest);
u_out_C = tablelookup(C.T_TLU, C.p_TLU, A.u_TLU, T, p_C, interpolation = linear, extrapolation = nearest);
% jet pump specific
b = An / Ath;
c = (1 - b) / b;
Z = rho_A * (((mdot_A / rho_A)^2) / (2 * An^2));
M = mdot_B / mdot_A;
a = Ath / Ad;
in
% Convective part of Thermal Fluxes
Phi_convection_A == mdot_A * if gt(mdot_A, 0), u_A else u_out_A end;
Phi_convection_B == mdot_B * if gt(mdot_B, 0), u_B else u_out_B end;
Phi_convection_C == mdot_C * if gt(mdot_C, 0), u_C else u_out_C end;
% Conservation of mass
% mdot_A == (An/sqrt(1+Kn))*(sqrt((2/rho_A)*(A.p-p)))*rho_A;
mdot_A == (An/sqrt(1+Kn))*(sqrt((2/rho_A)*(abs(A.p-p))))*rho_A*(A.p-p)/abs(A.p-p); % (A.p-p)/abs(A.p-p) for avoiding negative terms, here was the mistake before
% mdot_B == (An*c/sqrt(1+Ken))*(sqrt((2/rho_B)*(B.p-p)))*rho_B;
mdot_B == (An*c/sqrt(1+Ken))*(sqrt((2/rho_B)*(abs(B.p-p))))*rho_B*(B.p-p)/abs(B.p-p); % (A.p-p)/abs(A.p-p) for avoiding negative terms
0 == mdot_A + mdot_B+ mdot_C;
% Conservation of momentum
C.p-p == Z*(b^2)*((2/b)+(2/1-b)*M^2-((1+M)^2)*(1+Kth+Kdi+a^2)); %assuming, that sign has to be minus, thus pressure would be okay
% Conservation of Energy
%%Thermal fluxes
Phi_A == Phi_convection_A + Gth * (A.T - T);
Phi_B == Phi_convection_B + Gth * (B.T - T);
Phi_C == Phi_convection_C + Gth * (C.T - T);
% Phi_A == Phi_convection_A;
% Phi_B == Phi_convection_B;
% Phi_C == Phi_convection_C;
%%Balance
0 == Phi_A + Phi_B + p_dv_A + p_dv_B+Phi_C + p_dv_C; %before: type 0 == a + b + c instead of c == a + b
% 0 == Phi_A + Phi_B + Phi_C;
% Fluid properties table lookup functions
u == tablelookup(A.T_TLU, A.p_TLU, A.u_TLU, T, p, interpolation = linear, extrapolation = nearest);
u_A == tablelookup(A.T_TLU, A.p_TLU, A.u_TLU, T_A, p_A, interpolation = linear, extrapolation = nearest);
u_B == tablelookup(A.T_TLU, A.p_TLU, A.u_TLU, T_B, p_B, interpolation = linear, extrapolation = nearest);
u_C == tablelookup(A.T_TLU, A.p_TLU, A.u_TLU, T_C, p_C, interpolation = linear, extrapolation = nearest);
rho == tablelookup(A.T_TLU, A.p_TLU, A.rho_TLU, T, p, interpolation = linear, extrapolation = nearest);
rho_A == tablelookup(A.T_TLU, A.p_TLU, A.rho_TLU, T_A , p_A, interpolation = linear, extrapolation = nearest);
rho_B == tablelookup(A.T_TLU, A.p_TLU, A.rho_TLU, T_B , p_B, interpolation = linear, extrapolation = nearest);
rho_C == tablelookup(A.T_TLU, A.p_TLU, A.rho_TLU, T_C , p_C, interpolation = linear, extrapolation = nearest);
nu == tablelookup(A.T_TLU, A.p_TLU, A.nu_TLU, T, p, interpolation = linear, extrapolation = nearest);
cp == tablelookup(A.T_TLU, A.p_TLU, A.cp_TLU, T, p, interpolation = linear, extrapolation = nearest);
alpha == tablelookup(A.T_TLU, A.p_TLU, A.alpha_TLU, T, p, interpolation = linear, extrapolation = nearest);
k == tablelookup(A.T_TLU, A.p_TLU, A.k_TLU, T, p, interpolation = linear, extrapolation = nearest);
beta == tablelookup(A.T_TLU, A.p_TLU, A.beta_TLU, T, p, interpolation = linear, extrapolation = nearest);
end
end
end
  2 Comments
Andreas
Andreas on 15 Oct 2015
Progress: Forgot, that two_port.ssc is no independent block but only object / foundation for top level blocks. So many of the missing conserving equations can be found in other blocks... I'll try to get this working.

Sign in to comment.

Accepted Answer

Andreas
Andreas on 2 Nov 2015
ssc-file working now, we had a negative term in sqrt before, so initialisation would fail, no matter there is no flow reversal (and thus negative mdot) in a Jet pump

More Answers (3)

Sohrab Forouzan Mehr
Sohrab Forouzan Mehr on 15 Oct 2015
Hi Andreas,
I tried this variant (see below). You won't get the variable exceeds equations error anymore, but running a model won't work either.
Any idea what went wrong?
component three_port
% Thermal Liquid Three Port
% Defines a thermal liquid domain component with external nodes A, B and C
% located on the left, left and right hand side of the block respectively.
% Also defines associated through variables as well as thermodynamic
% properties and flux scheme. The properties are computed using table
% lookup functions.
nodes
A = foundation.thermal_liquid.thermal_liquid; % A:left
B = foundation.thermal_liquid.thermal_liquid; % B:left
C = foundation.thermal_liquid.thermal_liquid; % C:right
end
parameters
% commanded_pressure = { 0, 'Pa' }; % Pressure differential
length = { 1e-1, 'm' }; % Characteristic longitudinal length
area = { 1e-2, 'm^2' }; % Pipe cross-sectional area
end
variables
% Component variables
mdot_A = { 0, 'kg/s' }; % Mass flow rate into A
mdot_B = { 0, 'kg/s' }; % Mass flow rate into B
mdot_C = { 0, 'kg/s' }; % Mass flow rate into C
Phi_A = { 0, 'J/s' }; % Thermal flux through A
Phi_B = { 0, 'J/s' }; % Thermal flux through B
Phi_C = { 0, 'J/s' }; % Thermal flux through C
p = { 1.01325,'bar' }; % Pressure
end
variables(Conversion=absolute)
T = { 293.15, 'K' }; % Temperature
end
variables(Access = protected)
Phi_convection_A= { 0, 'J/s' }; % Convective part of thermal flux through A
Phi_convection_B= { 0, 'J/s' }; % Convective part of thermal flux through B
Phi_convection_C= { 0, 'J/s' }; % Convective part of thermal flux through C
% Fluid properties - Default values are for water at p = 1 atmosphere and T = 293.15 Kelvin
%%Internal energy
u = { 84, 'J/g' }; % Internal energy
u_A = { 84, 'J/g' }; % Internal energy at A
u_B = { 84, 'J/g' }; % Internal energy at B
u_C = { 84, 'J/g' }; % Internal energy at C
%%Density
rho = {998.2, 'kg/m^3' }; % Density
rho_A = {998.2, 'kg/m^3' }; % Density at A
rho_B = {998.2, 'kg/m^3' }; % Density at B
rho_C = {998.2, 'kg/m^3' }; % Density at C
%%Specific heat at constant pressure, isobaric thermal expansion, and conductivity
cp = { 4.16, 'J/(g*K)' }; % Specific heat at constant pressure
alpha = { -2.0691e-04, '1/K' }; % Isobaric thermal expansion
k = { 598.5, 'mW/(m*K)' }; % Thermal conductivity
%%Isothermal bulk modulus and viscosity
beta = { 2.1791, 'GPa' }; % Isothermal bulk modulus
nu = { 1, 'mm^2/s' }; % Viscosity
end
branches
mdot_A : A.mdot -> *;
mdot_B : B.mdot -> *;
mdot_C : C.mdot -> *;
Phi_A : A.Phi -> *;
Phi_B : B.Phi -> *;
Phi_C : C.Phi -> *;
end
equations
let
% Pressures and temperatures at fluid boundaries
p_A = A.p;
T_A = A.T;
p_B = B.p;
T_B = B.T;
p_C = C.p;
T_C = C.T;
% Thermal conductance for each half of the source
Gth = if k*area/(length/2) <= A.G_min, A.G_min else k*area/(length/2) end;
% Thermal energy equation sources and sinks
p_dv = p * mdot_A * if gt(mdot_A, 0), 1/rho_A - 1/rho else 1/rho - 1/rho_B end;
% Internal energy at the internal temperature and corresponding to boundary pressure
u_out_A = tablelookup(A.T_TLU, A.p_TLU, A.u_TLU, T, p_A, interpolation = linear, extrapolation = nearest);
u_out_B = tablelookup(A.T_TLU, A.p_TLU, A.u_TLU, T, p_B, interpolation = linear, extrapolation = nearest);
u_out_C = tablelookup(A.T_TLU, A.p_TLU, A.u_TLU, T, p_C, interpolation = linear, extrapolation = nearest);
% % Thermal energy equation sources and sinks
% p_dv = p * mdot_A * if gt(mdot_A, 0), 1/rho_A - 1/rho else 1/rho - 1/rho_B end;
in
% Convective part of Thermal Fluxes
Phi_convection_A == mdot_A * if gt(mdot_A, 0), u_A else u_out_A end;
Phi_convection_B == mdot_B * if gt(mdot_B, 0), u_B else u_out_B end;
Phi_convection_C == mdot_C * if gt(mdot_C, 0), u_C else u_out_C end;
% Conservation of mass
% mdot_A == -mdot_B;
mdot_C == mdot_A + mdot_B;
% Conservation of momentum
p == (A.p + B.p + C.p)/3;
B.p - A.p == 0;
B.p - C.p == 0;
% Conservation of Energy
%%Thermal fluxes
Phi_A == Phi_convection_A + Gth * (A.T - T);
Phi_B == Phi_convection_B + Gth * (B.T - T);
Phi_C == Phi_convection_C + Gth * (C.T - T);
% Phi_A == Phi_convection_A;
% Phi_B == Phi_convection_B;
% Phi_C == Phi_convection_C;
%%Balance
Phi_C == Phi_A + Phi_B + p_dv;
% Fluid properties table lookup functions
u == tablelookup(A.T_TLU, A.p_TLU, A.u_TLU, T, p, interpolation = linear, extrapolation = nearest);
u_A == tablelookup(A.T_TLU, A.p_TLU, A.u_TLU, T_A, p_A, interpolation = linear, extrapolation = nearest);
u_B == tablelookup(A.T_TLU, A.p_TLU, A.u_TLU, T_B, p_B, interpolation = linear, extrapolation = nearest);
u_C == tablelookup(A.T_TLU, A.p_TLU, A.u_TLU, T_C, p_C, interpolation = linear, extrapolation = nearest);
rho == tablelookup(A.T_TLU, A.p_TLU, A.rho_TLU, T, p, interpolation = linear, extrapolation = nearest);
rho_A == tablelookup(A.T_TLU, A.p_TLU, A.rho_TLU, T_A , p_A, interpolation = linear, extrapolation = nearest);
rho_B == tablelookup(A.T_TLU, A.p_TLU, A.rho_TLU, T_B , p_B, interpolation = linear, extrapolation = nearest);
rho_C == tablelookup(A.T_TLU, A.p_TLU, A.rho_TLU, T_C , p_C, interpolation = linear, extrapolation = nearest);
nu == tablelookup(A.T_TLU, A.p_TLU, A.nu_TLU, T, p, interpolation = linear, extrapolation = nearest);
cp == tablelookup(A.T_TLU, A.p_TLU, A.cp_TLU, T, p, interpolation = linear, extrapolation = nearest);
alpha == tablelookup(A.T_TLU, A.p_TLU, A.alpha_TLU, T, p, interpolation = linear, extrapolation = nearest);
k == tablelookup(A.T_TLU, A.p_TLU, A.k_TLU, T, p, interpolation = linear, extrapolation = nearest);
beta == tablelookup(A.T_TLU, A.p_TLU, A.beta_TLU, T, p, interpolation = linear, extrapolation = nearest);
end
end
end
  1 Comment
Andreas
Andreas on 19 Oct 2015
We indeed missed some equations, but now the basic model works. The matlab support suggested a graphical solution for the jet pump - we'll see, what works out for us.
Updated code in the top post.

Sign in to comment.


miracle Nkwocha
miracle Nkwocha on 8 Dec 2016
Hi
How do you create an .ssc file please?! I can't figure it out

Aram Amouzandeh
Aram Amouzandeh on 13 Dec 2019
Dear All,
I tried to build the model in MATLAB R2018b and get the following error:
Failed to generate 'TL_addon2_lib'
Caused by:
Error using TL_addon2.elements/jet_pump>equations (line 102)
No function or value named A.G_min.
Updating Model Advisor cache…
Model Advisor cache updated. For new customizations, to update the cache, use the Advisor.Manager.refresh_customizations method.
Multiple compilation errors detected while compiling TL_addon2_test_model.
Multiple compilation errors detected while compiling TL_addon2_test_model.
['TL_addon2_test_model/Thermal Liquid Jet Pump']: Cannot find parameter 'Phi_A_nominal_specify'. Please run ssc_build if you have made changes to Simscape file or if you are upgrading to a new version of Simscape.
In jet_pump.ssc A.G_min is used in a singe equation as argument:
% Thermal conductance for each part of the source
Gth = if k*(area/2)/(length/2) <= A.G_min, A.G_min else k*(area/2)/(length/2) end;
However, the variable is not defined anywhere else. Could the variable be a global MATLAB variable not valid anymore in version 2018?
Thanky for your input!
Cheers,
Aram
  1 Comment
Andreas
Andreas on 9 Jan 2020
This is possible. I build the model in 2015b and 2016a.
I'll have a look on it the next weeks and port it to 2018/2019.
What happens, when you open:
foundation.thermal_liquid.thermal_liquid
?

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!

Translated by