Energy balance issue with Simscape Heat Exchanger model

Hello,

I am trying to write a code of simple heat exchanger system based on the E-NTU method in simscape environment. I have a general code which does the heat transfer based on other questions (https://www.mathworks.com/matlabcentral/answers/index/?s_tid=gn_mlc_an) here. However, when I try to validate my model against the Simscape component of Heat Exchanger (TL-TL) it doesn't work. I am not even getting proper mass balance which I think should be the basic first step. I have worked long enough on this but still don't see the issue with the code. I am attaching my code and model for reference.

Any help would be greatly appreciated.

Thanks in advance!

component (Propagation = blocks) HX
% Heat Exchanger module
% 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   
          U             = {200, 'W/(m^2)/K'}; %Heat Transfer co-efficient
          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' }
      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;

% pv_A1 = 0;
% pv_A2 = 0;
% pv_B1 = 0;
% pv_B2 = 0;
              
              %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
              
              % 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;
              
               % Specific heat at constant volume
              
              cv_1 = cp1 - beta1 * alpha1^2 * T1 / rho1;
              cv_2 = cp2 - beta2 * alpha2^2 * T2 / rho2;
              
              % 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 = U*A/min(C1,C2)
          
%            Calcul of "efficiency" cross flow
%            if flow_type == flow.direct
           
               
              epsilon = if Cr == 1,  NTU/(1+NTU) else (1-exp(-NTU*(1-Cr)))/(1-Cr*exp(-NTU*(1-Cr))) end;
              
%            else
%               
%                 epsilon == 1 - exp(-1/Cr*(1-exp(-Cr*NTU))); 
%            end;  
              
          
              %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
              % Upwind energy scheme
          in
              % 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);
              
              if dynamic_compressibility == simscape.enum.onoff.off
                  
              
              % Mass balance equation
              mdot_A1 + mdot_B1 == 0;
              mdot_A2 + mdot_B2 == 0;
               
              % Energy balance equation for each fluid
              Q  == epsilon * min(C1,C2) * (T1 - T2);
              
              
              
              -Q + Phi_A1 + Phi_B1 == T1.der * cv_1 * rho1* volume;
               Q + Phi_A2 + Phi_B2 == T2.der * cv_2 * rho2* volume;
                
      
              
              else % dynamic_compressibility == simscape.enum.onoff.on
              
                  % Mass balance equation
              mdot_A1 + mdot_B1 == (p1.der/ beta1 - T1.der * alpha1) * rho1 * volume*0.00000001;
              mdot_A2 + mdot_B2 == (p2.der/ beta2 - T2.der * alpha2) * rho2 * volume*0.00000001;
              
              % 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;
               
              end
              % 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

Answer 1

[Rohan Kokate]

UPDATE: I resolved the issue with the energy balance. The heat transfer rate should be calculated as:

Q == epsilon * min(C1,C2) * (T_A1 - T_A2);

And also, I was using counter flow heat exchanger effectiveness formula instead of parallel flow heat exchanger effectiveness. That caused the unphyscial temperature profile of the fluids at outlet. The system and code are both working as expected now. I validated the code and system with the library component of Heat Exchanger (TL-TL).

Hope this helps anyone trying to do similar thing!

[Niveditha Rajendrakumar]

Hi Rohan,

I came across your solutions, and I figured you might have ideas as to where I'm going wrong with my heat exchanger block. You can find my problem statement in the link below.

https://nl.mathworks.com/matlabcentral/answers/1652275-why-does-the-simscape-heat-exchanger-block-not-reduce-in-effectiveness-as-the-mass-flow-rate-increas?s_tid=prof_contriblnk

I hope the explanation is clear and hope you can help!

Thanks!

Niveditha