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

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

Update (files can be found here: http://www.mathworks.com/matlabcentral/fileexchange/53588-build-a-jet-pump-model-in-simscape-s-thermal-liquid-domain-files

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

[Andreas]

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]

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.

Answer 1

[Andreas]

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

Answer 2

[Sohrab Forouzan Mehr]

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

[Andreas]

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.

Answer 3

[miracle Nkwocha]

Hi

How do you create an .ssc file please?! I can't figure it out

[Andreas]

.ssc files are plain text files which can be compiled through ssc_build afterwards.

Have a look here: https://de.mathworks.com/help/physmod/simscape/ref/ssc_build.html?s_tid=srchtitle

Answer 4

[Aram Amouzandeh]

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

[Andreas]

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

?

[Omer Sariyildiz]

Did you find any solution for this problem?

I'm using 2022a and I couldn't find any idea to fix.

[Andreas]

@Omer Sariyildiz @Aram Amouzandeh

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.