Heat Exchanger Model

clc
clear
close all

%Estimate outlet temperature and solve to calculate Q using the 3 equations

%Equation
%Q=m1*cp1*(T1i-T1o)=m2*cp2*(T2o-T2i)=U*A*F*deltaTlm

%deltaTlm=((T1i-T2i)-(T1o-T2o))/log((T1i-T2i)/(T1o-T2o)) for parallel
%deltaTlm=((T1i-T2o)-(T1o-T2i))/log((T1i-T2o)/(T1o-T2i)) for counter flow

%-------------------------------------------------------------------------
%Design requirement
%Engine oil outlet mues be less than or equal to 115 degrees
%Pressure drop on each side must be deltaP<=10kPa

%engine oil
mdot1=0.23;
cp1=2307;%J/kgK
T1i=120+273; %Kelvin
Required_T1o=115+273;
k1=0.135; %W/m.K
miu1=1.027*(10^-2); %N.s/m^2
PR1=175;
rho1=828; %kg/m^3
Rf1=0.176*10^-3; %m^2*k/W
upperrange1=1;
lowerrange1=0.5;

%-------------------------------------------------------------------------

%Engine coolant
mdot2=0.47;
cp2=3650;%J/kgK
T2i=90+273; %Kelvin (100 degrees)
k2=0.442; %W/m.K
miu2=0.08*(10^-2); %N.s/m^2
PR2=6.6;
rho2=1020; %kg/m^3
Rf2=0.353*10^-3;%m^2*k/W
upperrange2=2.4;
lowerrange2=1.2;

%-------------------------------------------------------------------------

%Parameters
pitch_shape=2; %*** Square pitch=1 / Triangular Pitch =2
kw=42.7; %W/m.K
Ds=50.8*10^-3; %Shell inside diameter (m)
Lt=254*10^-3; %Tube length (m)
do=3.175*10^-3; %Tube outside diameter
dr=1.3; %Given Diameter ratio where dr=do/di
di=(1/dr)*do; %Tube inner diameter (m)
fprintf('Tube inner diameter(m)=%',di)
disp(di)
Np=2; %Number of passes (defined)
B=25.4*10^-3; %Baffling space (Assumed)(m)
Nb=Lt/B-1; %Number of Baffles
fprintf('Number of Baffles=%',Nb)
disp(Nb)

%Tube pitch ratio Pr where Pr=Pt/do
%Pt= tube pitch
Pr=1.25; %Given tube pitch ratio
Pt=Pr*do; %Obtain tube pitch
fprintf('Tube pitch(m)=%',Pt)
disp(Pt)

%-------------------------------------------------------------------------

%Tube Clearance
Ct=Pt-do; %Tube clearance (mm)
fprintf('Tube clearance(m)=%',Ct)
disp(Ct)

%tube count calculation constants CTP
%To estimate the number of tubes
Nt=[0 0 0]';
CTP=[0.93 0.9 0.85]';
if pitch_shape==1 %Square pitch layout
    CL=1;
end
if pitch_shape==2 %triangular pitch layout
    CL=0.866;
end
%Constant pitch layout
for i=1:3
    ShadeArea=CL*Pt^2;
    Nt(i)=CTP(i)*((pi*Ds^2/4)/ShadeArea);
end
%disp(Nt)
Nt=round(Nt(Np));
fprintf('Number of tubes needed=%',Nt)
disp(Nt)

%-------------------------------------------------------------------------

%Tube side cross flow area / velocity/ reynolds number for engine oil
Ac1=(pi*(di)^2/4)*(Nt/Np); %Cross-flow area (m^2)
v1=mdot1/(rho1*Ac1); %velocity ms-1
RE1=(rho1*v1*(di))/miu1; %Reynolds number
fprintf('Engine oil cross flow area=%',Ac1)
disp(Ac1)
fprintf('Engine oil velocity=%',v1)
disp(v1)
fprintf('Engine oil Reynolds number=%',RE1)
disp(RE1)

%Varify
if RE1<2300
    disp('Laminar flow')
else
    disp('turbulent flow')
end

if v1<upperrange1
    disp('Velocity within range')
    if v1<lowerrange1
        disp('Velocity lower than range')
        if v1<lowerrange1/2
            dis('Velcity is too low')
        end
    end
    if v1>upperrange1
        disp('Velocity higher than range')
    end
end

%Friction factor
if RE1<2300
    fricfac=16/RE1;
else
    fricfac=(1.58*log(RE1)-3.28)^2;
end
fprintf('Friction Factor for engine oil=%',fricfac)
disp(fricfac);

%Nusselt number
if RE1>2300
    Nu1=(fricfac/2)*((RE1-1000)*PR1)/(1+12.7*((fricfac/2)^0.5)*(PR1^(2/3)-1));
else
    Nu1=1.86*((di*RE1*PR1)/Lt)^(1/3);
end
fprintf('Nusselt number=%',Nu1)
disp(Nu1);

%Convection heat tranfer coefficient
h1=(Nu1*k1)/di;
fprintf('Convection heat transfer coefficient h (W/m^2*k)=%',h1)
disp(h1);

%-------------------------------------------------------------------------

%Shell side free flow area / velocity
Ac2=(Ds*Ct*B)/Pt;
fprintf('free flow area(m^2)=%',Ac2);
disp(Ac2)
v2=mdot2/(rho2*Ac2);
fprintf('Velocity of Coolant(ms-1)=%',v2);
disp(v2)

%Varify
if v2<upperrange2
    disp('Velocity within range')
    if v2<lowerrange2
        disp('Velocity lower than range')
        if v2<lowerrange2/2
            dis('Velcity is too low')
        end
    end
    if v2>upperrange2
        disp('Velocity higher than range')
    end
end

%Equivalent diameter
 if pitch_shape==1 %Square pitch layout
     De=4*(Pt^2-pi*(do^2)/4)/pi*do;
 end
 if pitch_shape==2 %triangular pitch layout
     De=4*(((Pt^2*sqrt(3)/4)-(pi*do^2/8))/(pi*do/2));
 end
 fprintf('Equivalent Diameter (m)=%',De);
 disp(De)

%Reynolds number for engine coolant
RE2=(rho2*v2*De)/miu2;
fprintf('Reynolds number for engine coolant=%',RE2);
disp(RE2)

%Friction factor 2
fricfac2= exp(0.576-0.19*log(RE2));
fprintf('Friction Factor2 for coolant=%',fricfac2)
disp(fricfac2);


%Nusselt number
if RE2>2000
    if RE2<1*10^6
    Nu2=0.36*(RE2^0.55)*(PR2^(1/3));
    end
end
fprintf('Nusselt number2=%',Nu2)
disp(Nu2);

%Heat transfer coefficient
h2=(Nu2*k2)/De;
fprintf('Heat transfer coefficient for engine coolant h (W/m^2*k)=%',h2);
disp(h2)

%-------------------------------------------------------------------------

%Total heat transfer areas for both fluids
Ai=pi*di*Lt*Nt;
Ao=pi*do*Lt*Nt;
fprintf('Total heat transfer inner area (m^2) =%',Ai);
disp(Ai)
fprintf('Total heat transfer outer area (m^2)=%',Ao);
disp(Ao)

%The overall heat transfer coefficient UAo
a=1/(h1*Ai);
b=Rf1/Ai;
c=log(do/di)/(2*pi*kw*Lt*Nt);
d=Rf2/Ao;
e=1/(h2*Ao);

UAo=1/(a+b+c+d+e);
fprintf('Overall heat transfer coefficient (W/K)=%',UAo);
disp(UAo)

%-------------------------------------------------------------------------

%Heat capacities for both fluids
C1=mdot1*cp1;
C2=mdot2*cp2;
Cmin=C1;
Cmax=C2;
fprintf('Heat capacity for Engine oil(C minimum)=%',C1);
disp(C1)
fprintf('Heat capacity for Engine coolant(C minimum)=%',C2);
disp(C2)
%Heat capacity ratio
Cr=Cmin/Cmax;
fprintf('Heat capacity ratio=%',Cr);
disp(Cr)
%Number of transfer units
NTU=UAo/Cmin;
fprintf('Number of transfer units(NTU)=%',NTU);
disp(NTU)
%Effectiveness for a shell and tube heat exchanger
NTU1=NTU/Np;
Ehx=2*(1+Cr+(1+Cr^2)^0.5*((1+exp(-NTU1*(1+Cr^2)^0.5))/(1-exp(-NTU1*(1+Cr^2)^0.5))))^-1;
fprintf('effectiveness Ehx=%',Ehx);
disp(Ehx)

%Effectiveness maybe expressed as
%        q          C1*(T1i-T1o)       C2*(T2o-T2i)
%Ehx= ------- = ------------------ = ----------------
%      qmax        Cmin*(T1i-T2i)      Cmin(T1i-T2i)

%Comparison with the outlet temperatures
T1o=T1i-Ehx*(Cmin/C1)*(T1i-T2i);
T2o=T2i+Ehx*(Cmin/C2)*(T1i-T2i);
fprintf('Calculated outlet temperature of engine oil (Kelvin)=%',T1o);
disp(T1o)
fprintf('Calculated outlet temperature of engine coolant (Kelvin)=%',T2o);
disp(T2o)
fprintf('Calculated outlet temperature of engine oil (degrees C)=%',T1o);
disp(T1o-273)
fprintf('Calculated outlet temperature of engine coolant (degrees C)=%',T2o);
disp(T2o-273)
RndT1o=round(T1o);

if RndT1o>Required_T1o-2
    if RndT1o<Required_T1o+2
        disp('Outlet temperature Requirement satisfied')
    end
else
    disp('Outlet temperature Requirement not satisfied')
end

%Heat transfer rate
q=Ehx*Cmin*(T1i-T2i);
fprintf('Overall heat transfer (W)=%',q);
disp(q)

%Pressure drops
deltaP1=4*(fricfac*Lt/di+1)*Np*0.5*rho1*v1^2;
deltaP2=fricfac2*(Ds/De)*(Nb+1)*0.5*rho2*v2^2;
fprintf('Pressure drop of engine oil delta P1 (Pa)=%',deltaP1);
disp(deltaP1)
fprintf('Pressure drop of engine coolant delta P2 (Pa)=%',deltaP2);
disp(deltaP2)

%Surface density
beta=(Ao+Ai)/((pi*Ds*2/4)*Lt);
fprintf('Surface density=%',beta);
disp(beta)

clc
%-------------------------------------------------------------------------
disp('Default Input')
fprintf('Engine oil inlet temperature (Kelvin)=%',T1i);
disp(T1i)
fprintf('Engine coolant inlet temperature=%',T2i);
disp(T2i)
fprintf('Engine Oil Mass flow rate=%',mdot1);
disp(mdot1)
fprintf('Engine coolant mass flow rate=%',mdot2);
disp(mdot2)
fprintf('Fouling factor of engine oil=%',Rf1);
disp(Rf1)
fprintf('Fouling factor of engine oil=%',Rf2);
disp(Rf2)
disp('---------------------------------------------')
disp('Requirements')
fprintf('Engine oil requirement outlet temperature=%',Required_T1o);
disp(T1o)
disp('---------------------------------------------')
disp('Desgin specification')
fprintf('Number of Passes=%',Np);
disp(Np)
fprintf('Shell inner diameter(m)=%',Ds);
disp(Ds)
fprintf('Tube Outer diameter(m)=%',do);
disp(do)
fprintf('Tube Inner diameter(m)=%',di);
disp(di)
fprintf('Tube length(m)=%',Lt);
disp(Lt)
fprintf('Number of tubes=%',Nt);
disp(Nt)
fprintf('Tube Clearance(m)=%',Ct);
disp(Ct)
fprintf('Baffle Spacing(m)=%',B);
disp(B)
fprintf('Number of Baffles=%',Nb);
disp(Nb)
fprintf('Engine oil outlet temperature(Kelvin)=%',T1o);
disp(T1o)
fprintf('Engine coolant outlet temperature(Kelvin)=%',T2o);
disp(T2o)
fprintf('Heat transfer rate(W)=%',q);
disp(q)
fprintf('Surface area density (m^2/m^3)=%',beta);
disp(beta)
fprintf('Pressure drop of engine oil delta P1 (Pa)=%',deltaP1);
disp(deltaP1)
fprintf('Pressure drop of engine coolant delta P2 (Pa)=%',deltaP2);
disp(deltaP2)
fprintf('Over heat transfer coefficient UA =%',UAo);
disp(UAo)