Simple heat exchanger model problem on simscape

Hello,

I am trying to do a simple model of heat exchanger following the epsilon-NTU method. I tried different approach but none of them work. My last attempt is the following (with the help of other replies https://fr.mathworks.com/matlabcentral/answers/302891-problem-with-custom-simscape-heat-exchanger-model).

My model can simulate (in a wrong way) under special condition. If I change the initial value, nothing works anymore and even if I just change the priority level (without changing the associated value), nothing works anymore.

I looked and double check everything but I cannot find what is wrong with it...

Any help would be greatly appreciated !

PS: i attached the simulation file

component (Propagation = blocks) HX_model
% Heat Exchanger mod
% Based on the epsilon-NTU method
      nodes
          A1 = foundation.thermal_liquid.thermal_liquid; % A1:top
          B1 = foundation.thermal_liquid.thermal_liquid; % B1:top
          B2 = foundation.thermal_liquid.thermal_liquid; % B2:bottom
          A2 = foundation.thermal_liquid.thermal_liquid; % A2:bottom
      end
      
      parameters   
          UA            = {20000, 'J/s/K'};
          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
      end
      
      parameters(Access=protected)
          mdot_min = { 0.0001, 'kg/s' }
      end
      
      variables (Access=protected)
          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
          
          Q          = {0, 'J/s'}; % Heat flow between branches
          rho1       = {1000, 'kg/m^3'};
          rho2       = {1000, 'kg/m^3'};
          u1         = {85,  'kJ/kg'};
          u2         = {85,  'kJ/kg'};
          
      end
      
      variables
          % Internal nodes
          p1          = {1, 'bar'  };    % Pressure of the liquid volume 1
          T1          = {293.15, 'K'    }; % Temperature of the liquid volume 1
          p2          = {1, 'bar'  };    % Pressure of the liquid volume 2
          T2          = {293.15, 'K'    }; % Temperature of the liquid volume 2
      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
              p_A1 = A1.p;
              p_A2 = A2.p;
              T_A1 = A1.T;
              T_A2 = A2.T;
              p_B1 = B1.p;
              p_B2 = B2.p;
              T_B1 = B1.T;
              T_B2 = B2.T;
              
              % Domain parameters for look up table only
              T1_TLU         = A1.T_TLU;
              p1_TLU         = A1.p_TLU;
              T2_TLU         = A2.T_TLU;
              p2_TLU         = A2.p_TLU;
              
              cp1_TLU        = A1.cp_TLU;
              cp2_TLU        = A2.cp_TLU;
              rho1_TLU       = A1.rho_TLU;
              rho2_TLU       = A2.rho_TLU;
              u1_TLU         = A1.u_TLU;
              u2_TLU         = A2.u_TLU;
              alpha1_TLU     = A1.alpha_TLU;
              alpha2_TLU     = A2.alpha_TLU;
              beta1_TLU      = A1.beta_TLU;
              beta2_TLU      = A2.beta_TLU;
              
              % Liquid properties table lookup upstream
              cp1    = tablelookup(T1_TLU, p1_TLU, cp1_TLU,  T1, p1, interpolation=linear, extrapolation=linear);
              cp2    = tablelookup(T2_TLU, p2_TLU, cp2_TLU,  T2, p2, interpolation=linear, extrapolation=linear);
              
              beta1  = tablelookup(T1_TLU, p1_TLU, beta1_TLU,  T1, p1, interpolation=linear, extrapolation=linear);
              beta2  = tablelookup(T2_TLU, p2_TLU, beta2_TLU,  T2, p2, interpolation=linear, extrapolation=linear);
              
              alpha1 = tablelookup(T1_TLU, p1_TLU, alpha1_TLU, T1, p1, 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
              volume = HX_length * HX_section;
              
              h1 = u1 + p1/rho1;
              dUdT1 = (cp1 - h1*alpha1) * rho1 * volume;
              dUdp1 = (rho1*h1/beta1 - T1*alpha1) * volume;
              
              h2 = u2 + p2/rho2;
              dUdT2 = (cp2 - h2*alpha2) * rho2 * volume;
              dUdp2 = (rho2*h2/beta2 - T2*alpha2) * volume;
              
              % Calcul of heat capacity 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;
              
              % Calcul of the ratio of minimum /maximum
              Cr = min(C1,C2)/max(C1,C2)
              
              % Calcul of Number of Transfer Unit
              NTU = UA/min(C1,C2)
              
              % Calcul of "efficiency"
              epsilon = if Cr == 1,  NTU/(1+NTU) else (1-exp(-NTU*(1-Cr)))/(1-Cr*exp(-NTU*(1-Cr))) end;
              
          in
              
              % table look up for solving 
              rho1   == tablelookup(T1_TLU, p1_TLU, rho1_TLU, T1, p1, interpolation=linear, extrapolation=linear);
              rho2   == tablelookup(T2_TLU, p2_TLU, rho2_TLU, T2, p2, interpolation=linear, extrapolation=linear);
              
              u1     == tablelookup(T1_TLU, p1_TLU, u1_TLU,   T1, p1, interpolation=linear, extrapolation=linear);
              u2     == tablelookup(T2_TLU, p2_TLU, u2_TLU,   T2, p2, interpolation=linear, extrapolation=linear);
              
             % Mass balance equation
             mdot_A1 + mdot_B1 == (p1.der/ beta1 - T1.der * alpha1) * rho1 * volume;
             mdot_A2 + mdot_B2 == (p2.der/ beta2 - T2.der * alpha2) * rho2 * volume;
             
              % Real heat exchanged (Q = epsilon Q_max)
             Q  == epsilon * min(C1,C2) * (T1 - T2);
              
             % Heat transfer equation
             -Q  == C1 * (T_A1 - T_B1);
             Q   == C2 * (T_A2 - T_B2);
             
             % Temperature fluid
             T1 == (T_A1 + T_B1) / 2;
             T2 == (T_A2 + T_B2) / 2;
              
             % Energy balance equation for each fluid
             Q + Phi_A1 + Phi_A2 == p1.der * dUdp1 + T1.der * dUdT1;
             Q + Phi_B1 + Phi_B2 == p2.der * dUdp2 + T2.der * dUdT2;
             
             % Pressure loss equation for each fluid
             p_A1 - p_B1 == K1 * mdot_A1;
             (p_B1+p_A1)/2 == p1;

             p_A2 - p_B2 == K2 * mdot_A2;
             (p_B2+p_A2)/2 == p2;

          end
      end
      end

[Sylvain Mathonniere]

So,

I managed to fix my issues. Among other little details, the main problem was that I did not implement the "upwind energy scheme" ( https://fr.mathworks.com/help/physmod/hydro/ug/energy-flow-calculations.html ) that allow simscape to know in which direction the fluid is flowing and using the correct value of internal energy for the energy flow. In the end, the code is fairly complicated...

I wrote the code below for people like me who are struggling making a "simple" model work. If there is anything to learn is that there is no such thing as a "simple model" with thermal fluid, since you need to take into account many things as soon as heat exchange is involved.

I ran into another issue while testing my code however. My simulation (attached) runs only for specific value of total simulation time. For example, If I run my simulation for 1060s then I have no issue. If I run it for 1062s, simscape will never manage to find the first point, it will keep decreasing the size step without result. This is really strange to me as I do not see why the total simulation time can influence the resolution of the first point ?

If anybody has any information about it, I would be grateful.

component (Propagation = blocks) HX_model
% Heat Exchanger mod
% Based on the epsilon-NTU method
      nodes
          A1 = foundation.thermal_liquid.thermal_liquid; % A1:top
          B1 = foundation.thermal_liquid.thermal_liquid; % 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   
          UA            = {2000, 'J/s/K'};
          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
      end
      
      parameters(Access=protected)
          mdot_min = { 0.0001, 'kg/s' }
      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
          
          Q          = {1000, 'J/s'};
          
      end
      
      variables
          % Internal nodes
          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

      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
              p_A1 = A1.p;
              p_A2 = A2.p;
              p_B1 = B1.p;
              p_B2 = B2.p;
              
              % Domain parameters for look up table only
              T1_TLU         = A1.T_TLU;
              p1_TLU         = A1.p_TLU;
              T2_TLU         = A2.T_TLU;
              p2_TLU         = A2.p_TLU;
              
              cp1_TLU        = A1.cp_TLU;
              cp2_TLU        = A2.cp_TLU;
              rho1_TLU       = A1.rho_TLU;
              rho2_TLU       = A2.rho_TLU;
              u1_TLU         = A1.u_TLU;
              u2_TLU         = A2.u_TLU;
              alpha1_TLU     = A1.alpha_TLU;
              alpha2_TLU     = A2.alpha_TLU;
              beta1_TLU      = A1.beta_TLU;
              beta2_TLU      = A2.beta_TLU;
              
              %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
              % Upwind energy scheme
              k1_cv             = A1.k_cv;
              k2_cv             = A2.k_cv;
              
              max_aspect_ratio_A1 = A1.max_aspect_ratio;
              max_aspect_ratio_B1 = B1.max_aspect_ratio;
              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_A1 = max_aspect_ratio_A1 * sqrt(4*HX_section/pi);
              max_conduction_length_B1 = max_aspect_ratio_B1 * sqrt(4*HX_section/pi);
              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_A1 = ...
                    if le(HX_length/2, max_conduction_length_A1), ...
                        k1_cv*HX_section/(HX_length/2) ...
                    else ...
                        k1_cv*HX_section/max_conduction_length_A1 ...
                    end;

              G_B1 = ...
                    if le(HX_length/2, max_conduction_length_B1), ...
                        k1_cv*HX_section/(HX_length/2) ...
                    else ...
                        k1_cv*HX_section/max_conduction_length_B1 ...
                    end;
                
              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_A1 = sqrt(mdot_A1^2 + 4*G_A1^2);
              mdot_abs_B1 = sqrt(mdot_B1^2 + 4*G_B1^2);
              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_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;
              
              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_A1  = tablelookup(T1_TLU, p1_TLU, u1_TLU, A1.T, p_A1, interpolation=linear, extrapolation=linear);
              u_out_A1 = u1;
              u_in_B1  = tablelookup(T1_TLU, p1_TLU, u1_TLU, B1.T, p_B1, interpolation=linear, extrapolation=linear);
              u_out_B1 = u1;
              
              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 (modif rho_B1 = rho1  et rho_B2 = Rho2)
              pv_A1 = mdot_A1/mdot_abs_A1 * p_A1/rho_A1;
              pv_B1 = mdot_B1/mdot_abs_B1 * p_B1/rho_B1;
              pv_A2 = mdot_A2/mdot_abs_A2 * p_A2/rho_A2;
              pv_B2 = mdot_B2/mdot_abs_B2 * p_B2/rho_B2;              
              
              %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
              
              % Liquid properties table lookup upstream
              cp1    = tablelookup(T1_TLU, p1_TLU, cp1_TLU,  T1, p1, interpolation=linear, extrapolation=linear);
              cp2    = tablelookup(T2_TLU, p2_TLU, cp2_TLU,  T2, p2, interpolation=linear, extrapolation=linear);
              
              beta1  = tablelookup(T1_TLU, p1_TLU, beta1_TLU,  T1, p1, interpolation=linear, extrapolation=linear);
              beta2  = tablelookup(T2_TLU, p2_TLU, beta2_TLU,  T2, p2, interpolation=linear, extrapolation=linear);
              
              alpha1 = tablelookup(T1_TLU, p1_TLU, alpha1_TLU, T1, p1, 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
              volume = HX_length * HX_section;
              
              h1 = u1 + p1/rho1;
              dUdT1 = (cp1 - h1*alpha1) * rho1 * volume;
              dUdp1 = (rho1*h1/beta1 - T1*alpha1) * volume;
              
              h2 = u2 + p2/rho2;
              dUdT2 = (cp2 - h2*alpha2) * rho2 * volume;
              dUdp2 = (rho2*h2/beta2 - T2*alpha2) * volume;
              
              % Calcul 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;
              
              % Calcul of the ratio of minimum /maximum
              Cr = min(C1,C2)/max(C1,C2)
              
              % Calcul of Number of Transfer Unit
              NTU = UA/min(C1,C2)
              
              % Calcul of "efficiency"
              epsilon = if Cr == 1,  NTU/(1+NTU) else (1-exp(-NTU*(1-Cr)))/(1-Cr*exp(-NTU*(1-Cr))) end;
              
              
          in
              %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
              % Upwind energy scheme
              
              % Density table lookup
              rho_A1 == tablelookup(T1_TLU, p1_TLU, rho1_TLU, T_A1, p_A1, interpolation=linear, extrapolation=linear);
              rho_B1 == tablelookup(T1_TLU, p1_TLU, rho1_TLU, T_B1, p_B1, interpolation=linear, extrapolation=linear);
              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_A1/mdot_abs_A1 == step_plus_A1*u_in_A1 - step_minus_A1*u_out_A1 + pv_A1;
              Phi_B1/mdot_abs_B1 == step_plus_B1*u_in_B1 - step_minus_B1*u_out_B1 + pv_B1;
              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_A1 == step_plus_A1*u_in_A1 + step_minus_A1*u_out_A1;
              u_B1 == step_plus_B1*u_in_B1 + step_minus_B1*u_out_B1;
              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_A1 == tablelookup(T1_TLU, p1_TLU, u1_TLU, T_A1, p_A1, interpolation=linear, extrapolation=linear);
              u_B1 == tablelookup(T1_TLU, p1_TLU, u1_TLU, T_B1, p_B1, interpolation=linear, extrapolation=linear);
              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 
              rho1   == tablelookup(T1_TLU, p1_TLU, rho1_TLU, T1, p1, interpolation=linear, extrapolation=linear);
              rho2   == tablelookup(T2_TLU, p2_TLU, rho2_TLU, T2, p2, interpolation=linear, extrapolation=linear);
              
              u1     == tablelookup(T1_TLU, p1_TLU, u1_TLU,   T1, p1, interpolation=linear, extrapolation=linear);
              u2     == tablelookup(T2_TLU, p2_TLU, u2_TLU,   T2, p2, interpolation=linear, extrapolation=linear);
              
              % Mass balance equation
              mdot_A1 + mdot_B1 == (p1.der/ beta1 - T1.der * alpha1) * rho1 * volume;
              mdot_A2 + mdot_B2 == (p2.der/ beta2 - T2.der * alpha2) * rho2 * volume;
              
              % Energy balance equation for each fluid
              
              % Real heat exchanged (Q = epsilon * Q_max)
              %Q  == epsilon * min(C1,C2) * (T_B1 - T_B2);
              Q  == epsilon * min(C1,C2) * (T1 - T2);
              
              -Q + Phi_A1 + Phi_B1 == p1.der * dUdp1 + T1.der * dUdT1;
               Q + Phi_A2 + Phi_B2 == p2.der * dUdp2 + T2.der * dUdT2;
             
              % Pressure loss equation for each fluid
              p_A1 - p1 == K1 * mdot_A1;
              p_B1 - p1 == K1 * mdot_B1;
             
              p_A2 - p2 == K2 * mdot_A2;
              p_B2 - p2 == K2 * mdot_B2;
             
              % Sortie
              heat == Q;

          end
      end
  end

[Sylvain Mathonniere]

Ok, I manage to simulate better using the stiff ODE solver : ode15s. Automatically, simscape was selecting ode23t which did not perform well on my problem. Hope it helps other people as well.

[Rohan Kokate]

I was trying to write a similar code and that was really helpful. Thanks!

However, I tried to validate your model with the HeatExchanger(TL) component and the results are not the same. I looked at the equations and everything looks perfect to me. I wonder what the reason might be? Also, the simple heat exhcanger energy balance doesn't seem to work for this model.

Q = mhCph(T1h-T2h) = mcCpc(T2c-T1c)

Can you please correct me if I'm missing something?

Thanks in advance!