Problem with Custom Heat Exchanger (G-TL) model

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
              % instead of temperature 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.

Thanks in advance for your time and efforts!

Regards,

Rohan Kokate.

[Yifeng Tang]

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.

[Maria Coutinho]

Can you share the solution? I'm having the same problme! Thanks in advance

[Rohan Kokate]

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.

Answer 1

[Yifeng Tang]

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]

I figured out the issue. Thanks for your suggestions!

[Maria Coutinho]

Can you share the solution? I'm having the same problem! Thanks