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

11 Ansichten (letzte 30 Tage)
Andreas
Andreas am 13 Okt. 2015
Bearbeitet: Andreas am 27 Jun. 2022
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 Kommentare
Andreas
Andreas am 14 Okt. 2015
Bearbeitet: Andreas am 28 Okt. 2015
Update:
I tried to modify the two-port-block at first, but get the "number of variables exceeds number of equations" error when testing. The block should deliver the same results as a normal 3-way-crossing of thermal liquid paths.
Edit: removed old ssc code
Andreas
Andreas am 15 Okt. 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.

Melden Sie sich an, um zu kommentieren.

Akzeptierte Antwort

Andreas
Andreas am 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

Weitere Antworten (3)

Sohrab Forouzan Mehr
Sohrab Forouzan Mehr am 15 Okt. 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 Kommentar
Andreas
Andreas am 19 Okt. 2015
Bearbeitet: Andreas am 28 Okt. 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.

Melden Sie sich an, um zu kommentieren.


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

Aram Amouzandeh
Aram Amouzandeh am 13 Dez. 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
  3 Kommentare
Omer Sariyildiz
Omer Sariyildiz am 27 Jun. 2022
Did you find any solution for this problem?
I'm using 2022a and I couldn't find any idea to fix.
Andreas
Andreas am 27 Jun. 2022
Bearbeitet: Andreas am 27 Jun. 2022
Somehow the structure of the thermal liquid domain has changed. A simple solution would be, to build the jetpump model based on lookup tables often delivered from manufacturers like BAELZ and so on:
Works as well as the older ssc-variant.

Melden Sie sich an, um zu kommentieren.

Produkte

Community Treasure Hunt

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

Start Hunting!

Translated by