Code covered by the BSD License

# Vapor compression cycle component models

### Muhammad Tauha Ali (view profile)

13 Sep 2011 (Updated )

Modeling of vapor compression cycle components

[Tevapout Pevaptubein balance1 xevap eevap Vair]=fintubeevaporator(mdot,Vair,Qa,Tz,Tevapin,Pei,xini,Vairsolve)
```%{
Input:
mdot=refrigerant mass flow rate (kg/s)
Vair=evaporator air volumetric flow rate (initial guess if solving for Vair) (m3/s)
(Qa is the amount of heat that needs to be gained by the evaporator. used in energy balance equation)
Tz=zone temperature (K)
Tevapin=evporator inlet temperature (K)
Pei=evaporator inlet pressure (kPa)
xini=evaporator inlet quality of refrigerant (0-1)
Vairsolve=flag to solve for minimum Vair that satisfies energy balance for a given set of parameters (1 or 0)

Output:
Tevapout=evaporator outlet temperature (K)
Pevaptubein=evaporator outlet pressure (kPa)
balance1=energy balance convergence (kW)
xevap=evaporation portion (%)
eevap=evaporator effectiveness (%)
Vair=evaporator air volumetric flow rate (m3/s)
%}
function [Tevapout Pevaptubein balance1 xevap eevap Vair]=fintubeevaporator(mdot,Vair,Qa,Tz,Tevapin,Pei,xini,Vairsolve)
p=Tchparm;p.dP_HX=1;%flag if want to include pressure drop
disp('evaporator');

numair='Gray and Webb';numref='Gnielinski'; %correlation names
tube.Tx=Tz;p.Q=Qa;Pevaptubein=Pei;%outside HX conditions
tubecheck=0;[airrho aircp]=refpropm('DC','T',tube.Tx,'P',p.Pamb,p.airname);

tube.vertical=p.evap.vertical;tube.width=p.evap.width;
%tube parameters
tube.N=p.evap.Ntot/p.evap.tube.div;
tube.row=p.evap.tube.row;
tube.pitchtransverse=p.evap.tube.pitchtransverse;
tube.length=p.evap.tube.length;tube.lengthtotal=tube.N*tube.length+pi*tube.pitchtransverse/2*(tube.N-1);
tube.Do=p.evap.tube.Do;
tube.t=p.evap.tube.t;tube.Di=tube.Do-2*tube.t;
tube.k=p.evap.tube.k;
tube.elementlength=tube.length;%1/fin.p;
%fin parameters
fin.t=p.evap.fin.t;tube.Dc=tube.Do+2*fin.t;%wang paper collar diameter
fin.w=p.evap.fin.w/tube.row;tube.pitchlong=fin.w;% fin.w=longitudnal pitch assumption
fin.p=p.evap.fin.p;
fin.k=p.evap.fin.k;

mdot=mdot/p.evap.tube.div;mflux=mdot/(pi/4*tube.Di^2);

if Vairsolve==1
options_e = optimset('Algorithm','interior-point','Display','off','TolFun', 0.001);
problem=createOptimProblem('fmincon','objective',@mair_solver,'x0',Vair,'lb',.13,'ub',.18,'options',options_e);
bigstart = CustomStartPointSet([.18:-.03:.13]');
gs=MultiStart('Display','off','UseParallel','never');
[Vair,balance]=run(gs,problem,bigstart);
end

%% evaporator loop
balance=mair_solver(Vair);
function balance=mair_solver(Vair)
tube.elementlength=tube.length;tubelooplength=tube.elementlength;Pevaptubein=Pei;Ttubein=Tevapin;xin=xini;%resetting variables
Pevapboil=zeros(tube.N*2,1);Pevap=zeros(tube.N*2,1);Qevap=zeros(tube.N*2,1);Qevapboil=zeros(tube.N*2,1);Ttubeout=[];i=1;j=1;k=1;jj=1;tpcalc=0;trcheck=0;tubecheck=0;trchecksub=0;drycheck=0;%resetting variables
xinlimit=1-p.oilcon;xoillimit=1-p.oilcon;mair=Vair.*airrho;

% evaporator loop
while tubelooplength<tube.lengthtotal

% loop safety
if Pevaptubein<p.Ptp || Ttubein(i)>tube.Tx
break;end

% hair and fin heat transfer
[etasurf At Asurftube coeffconvout Cair]=airHT(numair, mair, tube, fin,p.airname,p.Pamb);

% for not exceeding total length and calc. what remains
if (tubelooplength+tube.elementlength)>tube.lengthtotal
tube.elementlength=tube.lengthtotal-tubelooplength; % element length changed. need to calc hair and fin heat transfer
[etasurf At Asurftube coeffconvout Cair]=airHT(numair, mair, tube, fin,p.airname,p.Pamb);
end

% for calc of local oil concentration
if p.oilcon~=0
if tpcalc==0 || xin<1

if xin<xoillimit
oilconlocal=p.oilcon/(1-xin);
else
oilconlocal=0;
end
elseif xin<xoillimit
oilconlocal=p.oilcon/(1-xin);
else
oilconlocal=0;
end
else oilconlocal=0;
end

if tpcalc==0 || xin<xoillimit
Tsat=refpropm('T','P',Pevaptubein,'Q',abs(xin),p.Refname);
if oilconlocal~=0
T1=refpropm('T','P',Pevaptubein,'Q',abs(xin),p.Refname);T2=refpropm('T','P',Pevaptubein+1,'Q',abs(xin),p.Refname);
equations=[1 T1;1 T2];values=[log(Pevaptubein)*T1;log(Pevaptubein+1)*T2];constants=equations\values;
A=constants(1)+182.5*oilconlocal-724.2*oilconlocal^2+3868*oilconlocal^3-5268.9*oilconlocal^4;
B=constants(2)-.722*oilconlocal+2.391*oilconlocal^2-13.779*oilconlocal^3+17.066*oilconlocal^4;
Tsat=A/(log(Pevaptubein)-B);
end
else
if xin<xoillimit
Tsat=refpropm('T','P',Pevaptubein,'Q',abs(xin),p.Refname);%this is because of temperature glide in r410a
if oilconlocal~=0
T1=refpropm('T','P',Pevaptubein,'Q',abs(xin),p.Refname);T2=refpropm('T','P',Pevaptubein+1,'Q',abs(xin),p.Refname);
equations=[1 T1;1 T2];values=[log(Pevaptubein)*T1;log(Pevaptubein+1)*T2];constants=equations\values;
A=constants(1)+182.5*oilconlocal-724.2*oilconlocal^2+3868*oilconlocal^3-5268.9*oilconlocal^4;
B=constants(2)-.722*oilconlocal+2.391*oilconlocal^2-13.779*oilconlocal^3+17.066*oilconlocal^4;
Tsat=A/(log(Pevaptubein)-B);
end
else
Tsat=refpropm('T','P',Pevaptubein,'Q',1,p.Refname);
end
end %***

% for checking whether to do heating if ref. is in liquid phase after expansion
if round(Ttubein(i)*10)/10<round(Tsat*10)/10 && tpcalc==0 || tpcalc==1
if tpcalc==1 && trchecksub==1 %only true when just came out of evaporation phase
[ref.k ref.u]=refpropm('LV','T',round(Ttubein(i)),'Q',1,p.Refmix1,p.Refmix2,p.Refmixcomp);
[ref.rho,ref.cp]=refpropm('DC','T',round(Ttubein(i)),'Q',1,p.Refname);trchecksub=0;
% Refrigerant Properties **
else
[ref.k,ref.u]=refpropm('LV','T',round(Ttubein(i)*100)/100,'P',Pevaptubein,p.Refmix1,p.Refmix2,p.Refmixcomp);
%due to problem with mixture properties near 2-phase after evaporation
if ref.k<0
[ref.k ref.u]=refpropm('LV','T',round(Ttubein(i)),'Q',1,p.Refmix1,p.Refmix2,p.Refmixcomp);
if ref.k<0
[ref.k ref.u]=refpropm('LV','T',round(Ttubein(i)*100)/100+1,'Q',1,p.Refmix1,p.Refmix2,p.Refmixcomp);end
end

[ref.rho,ref.cp]=refpropm('DC','T',round(Ttubein(i)*100)/100,'P',Pevaptubein,p.Refname);
if ref.cp<0
[ref.rho ref.cp]=refpropm('DC','T',round(Ttubein(i)),'P',Pevaptubein,p.Refname);
if ref.cp<0
[ref.rho ref.cp]=refpropm('DC','T',round(Ttubein(i)*100)/100+1,'Q',1,p.Refname);end
end
end

if p.oilcon~=0
oil=oil_properties(Ttubein(i));
ref.rho=(oilconlocal/oil.rho+(1-oilconlocal)/ref.rho)^-1;
ref.cp=(1-oilconlocal)*ref.cp+oilconlocal*oil.cp;
ref.k=ref.k*(1-oilconlocal)+oil.k*oilconlocal-.72*(oil.k-ref.k)*(1-oilconlocal)*oilconlocal;
ref.u=ref.u^(1-oilconlocal)*oil.u^oilconlocal;
end
ref.Pr=(ref.u*ref.cp)/ref.k; %**

% hrefi
Rein=mflux*tube.Di/ref.u;[f coeffconvin]=href(numref, Rein, ref, tube);Cref=mdot*ref.cp;

if Cair<Cref
Cmin=Cair;Cmax=Cref;
else Cmin=Cref;Cmax=Cair;
end

% for subcooling/superheating part
%total UA=[1/hwAw+ln(Do/Di)/(2pi*k*l)+1/etasurfhaAt]^-1;
UA=(1/(coeffconvin*Asurftube)+log(tube.Dc/tube.Di)/(2*pi*tube.k*tube.elementlength)+1/(etasurf*coeffconvout*At))^-1;
NTU=UA/Cmin;Cr=Cmin/Cmax;eff=1-exp(1/Cr*NTU^.22*(exp(-1*Cr*NTU^.78)-1));%exact for cr=1 but a good approx for 0<cr<1
Q=eff*Cmin*(Tz-Ttubein(i));Ttubeout(i)=Ttubein(i)+Q/Cref;

% for accurately identifying the location of transition from subcool to evaporation
while (Ttubeout(i))>Tsat  && trcheck==0
tube.elementlength=tube.elementlength/2; %for refinement of discretization to get more accurately the point of start of evaporation
%hair and fin HT
[etasurf At Asurftube coeffconvout Cair]=airHT(numair, mair, tube, fin,p.airname,p.Pamb);% element length changed. need to calc hair and fin heat transfer
%hrefi
[f coeffconvin]=href(numref, Rein, ref, tube);
if Cair<Cref
Cmin=Cair;Cmax=Cref;
else Cmin=Cref;Cmax=Cair;
end

%energy balance
UA=(1/(coeffconvin*Asurftube)+log(tube.Dc/tube.Di)/(2*pi*tube.k*tube.elementlength)+1/(etasurf*coeffconvout*At))^-1;
NTU=UA/Cmin;Cr=Cmin/Cmax;eff=1-exp(1/Cr*NTU^.22*(exp(-1*Cr*NTU^.78)-1));%exact for cr=1 but a good approx for 0<cr<1
Q=eff*Cmin*(Tz-Ttubein(i));Ttubeout(i)=Ttubein(i)+Q/Cref;

if p.dP_HX==1
% Pressure drop calculation
Pevaptube=(Pevaptubein*1e3-2*f*tube.elementlength*(mflux)^2/(tube.Di*ref.rho))/1e3;
Tsat=refpropm('T','P',Pevaptubein,'Q',1,p.Refname);
if oilconlocal~=0
T1=refpropm('T','P',Pevaptubein,'Q',1,p.Refname);T2=refpropm('T','P',Pevaptubein+1,'Q',1,p.Refname);
equations=[1 T1;1 T2];values=[log(Pevaptubein)*T1;log(Pevaptubein+1)*T2];constants=equations\values;
A=constants(1)+182.5*oilconlocal-724.2*oilconlocal^2+3868*oilconlocal^3-5268.9*oilconlocal^4;
B=constants(2)-.722*oilconlocal+2.391*oilconlocal^2-13.779*oilconlocal^3+17.066*oilconlocal^4;
Tsat=A/(log(Pevaptubein)-B);
end
else Pevaptube=Pevaptubein;
end
% for marking end of transition region accurate to within .01K
if (Ttubeout(i)+.01)>Tsat && (Ttubeout(i)-.01)<Tsat
trcheck=1;
Pevaptubein=Pevaptube;
% for calculation of ref quality at beginning of evaporation (This is due to element length resolution limit)
if  (Ttubeout(i))>Tsat
ref.hliq=refpropm('H','T',Tsat,'Q',0,p.Refname);
ref.hg=refpropm('H','T',Tsat,'Q',1,p.Refname);ref.hfg=(ref.hg-ref.hliq);
if tpcalc==0
xin=Q/(mdot*ref.hfg);end
end
end
end

% proceed to next element
tubelooplength=tubelooplength+tube.elementlength;Ttubein(i+1)=Ttubeout(i);i=i+1;

% for avoiding double calculation of pressure drop after transition region
if (trcheck==0  || tpcalc==1) && p.dP_HX==1
% Pressure drop calculation
Pevaptubein=(Pevaptubein*1e3-2*f*tube.elementlength/tube.Di/ref.rho*(mflux)^2)/1e3;
end

% for resetting element length after calc of transition from subcooling to evaporation
if tube.elementlength<tube.length && tpcalc~=1
tube.elementlength=tube.length-tube.elementlength;
[etasurf At Asurftube coeffconvout Cair]=airHT(numair, mair, tube, fin,p.airname,p.Pamb);% element length changed. need to calc hair and fin heat transfer
else
tube.elementlength=tube.length;% for calc. of point of evaporation and for calc. of superheating
end

% loop safety
if Pevaptubein<p.Ptp || (Ttubein(i)>tube.Tx && tpcalc==1)
break;end

% for evaporation
elseif tpcalc==0
while xin<xinlimit
% loop safety
if Pevaptubein<p.Ptp
break
end

% refrigerant liquid and gas phase properties
Tsat=refpropm('T','P',Pevaptubein,'Q',abs(xin),p.Refname);

if p.oilcon~=0
if xin<xoillimit
oilconlocal=p.oilcon/(1-xin);
else oilconlocal=0;
end
if oilconlocal~=0
T1=refpropm('T','P',Pevaptubein,'Q',abs(xin),p.Refname);T2=refpropm('T','P',Pevaptubein+1,'Q',abs(xin),p.Refname);
equations=[1 T1;1 T2];values=[log(Pevaptubein)*T1;log(Pevaptubein+1)*T2];constants=equations\values;
A=constants(1)+182.5*oilconlocal-724.2*oilconlocal^2+3868*oilconlocal^3-5268.9*oilconlocal^4;
B=constants(2)-.722*oilconlocal+2.391*oilconlocal^2-13.779*oilconlocal^3+17.066*oilconlocal^4;
Tsat=A/(log(Pevaptubein)-B);
end
else oilconlocal=0;
end

% loop safety
if Tsat>Tz
break; end

[ref.rholiq,ref.kliq,ref.cpliq,ref.uliq,ref.hliq,ref.surften]=refpropm('DLCVHI','T',Tsat,'Q',0,p.Refmix1,p.Refmix2,p.Refmixcomp);
[ref.rhog,ref.kg,ref.cpg,ref.ug,ref.hg]=refpropm('DLCVH','T',Tsat,'Q',1,p.Refmix1,p.Refmix2,p.Refmixcomp);ref.hfg=(ref.hg-ref.hliq);
ref.Prliq=(ref.uliq*ref.cpliq)/ref.kliq;ref.Prg=(ref.ug*ref.cpg)/ref.kg;

if p.oilcon~=0
oilconlocalliq=p.oilcon;
oil=oil_properties(Tsat);
ref.rholiq=(oilconlocalliq/oil.rho+(1-oilconlocalliq)/ref.rholiq)^-1;
ref.cpliq=(1-oilconlocalliq)*ref.cpliq+oilconlocalliq*oil.cp;
ref.kliq=ref.kliq*(1-oilconlocalliq)+oil.k*oilconlocalliq-.72*(oil.k-ref.kliq)*(1-oilconlocalliq)*oilconlocalliq;
ref.uliq=ref.uliq^(1-oilconlocalliq)*oil.u^oilconlocalliq;
ref.surften=ref.surften+(oil.surften-ref.surften)*oilconlocalliq^.51;
end

if p.dP_HX==1
numrefevap='Wojtan';pdropcorr='Flow based';
else
numrefevap='Wojtan';pdropcorr='None';
end

% loop to calculate two-phase heat transfer and pressure drop as qflux not known
qfluxe=(Cair*(Tz-Tsat))/(pi/4*tube.Di^2);qfluxen=qfluxe*100;qe=[1 1 1];flowp='a';flowpp='a';effpp=[];effp=[];
for qloop=1:30
if qe(2)~=qe(3)
tube.elementlength=tube.elementlength/2;
[etasurf At Asurftube coeffconvout Cair]=airHT(numair, mair, tube, fin,p.airname,p.Pamb);% element length changed. need to calc hair and fin heat transfer
end
qfluxen=qfluxe;
[dPevaptube twophasecoefflocalevap flow]=h2ph(xin,mflux,qfluxe,Pevaptubein,numrefevap,pdropcorr,tube,ref);
UA=(1/(twophasecoefflocalevap*Asurftube)+log(tube.Do/tube.Di)/(2*pi*tube.k*tube.elementlength)+1/(etasurf*coeffconvout*At))^-1;
Cmin=Cair;NTU=UA/Cmin;eff=1-exp(-1*NTU);
qfluxe=eff*Cmin*(Tz-Tsat)/(pi/4*tube.Di^2);
flowpp=flowp;flowp=flow;effpp=effp;effp=eff;
if strcmp(flowp,flowpp) && qe(1)~=1
qe(2)=1;qe(3)=1;
elseif qe(1)~=1
qe(2)=2;qe(3)=1;
end
if (strcmp('dryout',flowpp) || strcmp('mist',flowpp)) && ~(strcmp(flow,'dryout') || strcmp(flow,'mist'))
qfluxe=effpp*Cmin*(Tz-Tsat)/(pi/4*tube.Di^2);break;
end
if abs(qfluxe-qfluxen)/qfluxen<.01
break;
end
qe(1)=2;
end
if p.oilcon~=0
if ~(strcmp('mist',flow) || strcmp('dryout',flow))
Pevaptubein=Pevaptubein-dPevaptube*(oil.u/ref.uliq)^(.184*oilconlocal); % oil correction from Thome paper
else
Pevaptubein=Pevaptubein-dPevaptube;
end
else
Pevaptubein=Pevaptubein-dPevaptube;
end

Q=eff*Cmin*(Tz-Tsat);xout=xin+Q/(mdot*ref.hfg);

% for accurately identifying the location of transition from evapensation to subcooling
if (xout-.01)>=xinlimit
xout=xin;
tube.elementlength=tube.elementlength/2;
[etasurf At Asurftube coeffconvout Cair]=airHT(numair, mair, tube, fin,p.airname,p.Pamb);% element length changed. need to calc hair and fin heat transfer
else
if (strcmp('mist',flow) || strcmp('dryout',flow)) && drycheck==0
xoillimit=xin-.0001;
xout=xin;drycheck=1;
else
xin=xout;Pevapboil(j)=Pevaptubein;Qevapboil(j)=Q;j=j+1;Ttubeout(i)=Tsat;Ttubein(i+1)=Ttubeout(i);i=i+1;
% proceed to next element
tubelooplength=tubelooplength+tube.elementlength;
if tube.elementlength<tube.length
tube.elementlength=tube.length-tube.elementlength;
[etasurf At Asurftube coeffconvout Cair]=airHT(numair, mair, tube, fin,p.airname,p.Pamb);% element length changed. need to calc hair and fin heat transfer
end
end
end
tpcalc=1;
% for exit from evaporation loop
if xin>=xinlimit
if Pevaptubein<p.Ptp || Tsat>tube.Tx
break; end;
tubecheck=1;xevap=tubelooplength/tube.lengthtotal;
tpcalc=1; % to identify have been in evapensation loop
trcheck=1;trchecksub=1;%indicates done calc. of transition region from evaporation to superheating

% for resetting element length after calc of transition from evaporation to superheating
if tube.elementlength<tube.length
tube.elementlength=tube.length-tube.elementlength;
[etasurf At Asurftube coeffconvout Cair]=airHT(numair, mair, tube, fin,p.airname,p.Pamb);% element length changed. need to calc hair and fin heat transfer
end
elseif tubelooplength>=tube.lengthtotal
tubecheck=0;xevap=tubelooplength/tube.lengthtotal;tpcalc=1;
break;
end
end
end

% for HT and pressure drop balance
if tpcalc==1 && jj==1
jj=0;Pevap(k:k+length(Pevapboil)-1)=Pevapboil;
Qevap(k:k+length(Qevapboil)-1)=Qevapboil;k=length(Qevap(Qevap>0))+1;
else
Pevap(k)=Pevaptubein;Qevap(k)=Q;k=k+1;
end
end

% energy balance convergence equation
if Vairsolve==1
% use abs for balance so that fmincon tries to find Vair that makes balance zero if possible
if tubecheck==1
balance1=p.Q-sum(Qevap)*p.evap.tube.div/1e3;
balance=abs(p.Q-sum(Qevap)*p.evap.tube.div/1e3);
else balance=p.Q+1-sum(Qevap)*p.evap.tube.div/1e3;balance1=balance;
end
else
balance=p.Q-sum(Qevap)*p.evap.tube.div/1e3;balance1=balance;
end
end

Tevapout=Ttubeout(end);p.Te=Tevapout;
if balance*1e3<p.Q*1e3*.1
disp('converged');end
eevap=(p.Q-balance1)./(mair.*aircp.*(Tz-Tevapin)/1e3)*100;xevap=xevap*100;
end

%% oil properties equation
function oil=oil_properties(Ttubein)
oil.rho=(0.97386-6.91473e-4*(Ttubein-273))*1e3; %kg/m3
oil.cp=(4.186.*(.388+.00045.*(1.8.*(Ttubein-273)+32))./(oil.rho/998.5).^.5)*1e3;%J/kgK
oil.k=.1172/(oil.rho/998.5)*(1-.0054*(Ttubein-273)); %W/mK
oil.u=(1062.075*exp(-(Ttubein-273)/32.29)+4.90664)*1e-6*oil.rho;%Pa.s
oil.surften=(29-.4*(Ttubein-273))*1e-3;
end

%% air side heat transfer calculation
function [etasurf At Asurftube coeffconvout Cair]=airHT(numair, mair, tube, fin,airname,Pamb);
%area calculation
Aminairside=tube.pitchtransverse*tube.elementlength-tube.elementlength/fin.p*fin.t*tube.pitchtransverse-tube.Dc*tube.elementlength;
Af=2*(tube.pitchtransverse*(tube.pitchlong+fin.t/2)-pi*(tube.Dc/2)^2);
Ao=pi*tube.Dc*tube.elementlength;
At=Af*tube.elementlength/fin.p+Ao;Asurftube=pi*tube.Di*tube.elementlength;
%fin UA(schimdt method)
Xl=tube.pitchlong/2;%for inline and 1 row coil
%Xl=sqrt((tube.pitchlong)^2+(tube.pitchtransverse/2)^2)/2; for staggered layout ref(heat transfer enhancement of HX by Sadik Kakac page 160)
Xm=tube.pitchtransverse/2;
reqr0=1.28*Xm/(tube.Dc/2)*sqrt(Xl/Xm-0.2);%for one-row or inline coil layout ref(heat transfer enhancement of HX by Sadik Kakac page 145)
phi=(reqr0-1)*(1+0.35*log(reqr0));

%element mass calculation
mairel=mair*(tube.elementlength*tube.pitchtransverse)/(tube.vertical*tube.width);
mairflux=mair/Aminairside;

% Air Properties
[air.rho,air.k,air.cp,air.u]=refpropm('DLCV','T',tube.Tx,'P',Pamb,airname);air.Pr=(air.u*air.cp)/air.k;
Cair=mairel*air.cp;Reout=mairflux*tube.Dc/air.u;

%hair
switch numair
case 'Wang'
Dh=4*Aminairside*fin.w/Ao;
if tube.row==1
P1=1.9-.23*log(Reout);P2=-.236+.126*log(Reout);
j=.108*Reout^-.29*(tube.pitchtransverse/tube.pitchlong)^P1*(fin.p/tube.Dc)^-1.084*(fin.p/Dh)^-.786*(fin.p/tube.pitchtransverse)^P2;
coeffconvout=j*Reout*air.Pr^(1/3)*air.k/tube.Dc;
else
P6=-5.735+1.21*log(Reout/tube.row);
P5=-.083+0.058*tube.row/log(Reout);P4=-1.224-.076*(tube.pitchlong/Dh)^1.42/log(Reout);
P3=-.361-.042*tube.row/log(Reout)+.158*log(tube.row*(fin.p/tube.Dc)^.41);
j=.086*Reout^P3*tube.row^P4*(fin.p/tube.Dc)^P5*(fin.p/Dh)^P6*(fin.p/tube.pitchtransverse)^-.93;
coeffconvout=j*Reout*air.Pr^(1/3)*air.k/tube.Dc;
end
case 'domanski'%need wetbulb temp for using this
%correlation of young and briggs 1962 refernced in domanski
hc=0.134*air.k/tube.Do*Reout^0.681*air.Pr^0.333*(fin.p/fin.w)^0.2*(fin.p/fin.t)^0.1134;
%correction for latent heat by domanski
coeffconvout=hc*(1+(hfgwater*(wair-wwater))/(air.cp*(tempair-tempwater)));%wair can be found out by the dry bulb and wetbulb tem and tempwater is dewpoint temp
case 'Fischer Rice'
Rew=mairflux*tube.pitchtransverse/air.u;%for 3000<Rew<15000 not for low mass flow rate of air
j=0.0014+0.2618*(1/(1-(fin.w*(2*Xl))/(2*pi*tube.Do/2)))^-0.15*(mairflux*tube.Do/air.u)^-0.4;
%Co=1.0, 1.45 or 1.75 depending on whether the fins are smooth,wavy or louvered
Co=1;
coeffconvout=Co*mairflux*air.cp*air.Pr^(2/3)*j*((1-1280*tube.row*Reout^-1.2)/(1-5120*Rew^-1.2));
%3000<Rew<15000
case 'Gray and Webb'%accounting for tube geometry in detail (don't have original paper)
j= 0.14*Reout^-0.328*(tube.pitchtransverse/fin.w)^-0.502*(fin.p/tube.Do)^0.0312;%1row so tube.pitchlong=fin.w
j1=.991*(2.24*Reout^-.092*(tube.row/4)^-.031)^.607*(4-tube.row)*j;
coeffconvout=j1*mairflux*air.cp/air.Pr^(2/3);
%500 < ReD < 24700 ; 1.97 < Pt/D < 2.55; 1.7 < Pl/D < 2.58 ; 0.08 < Fp/D < 0.64
case 'lee'
%comparison of air side ht evaporators***.pdf%experimental samples
%are not close to our case
coeffconvout=0.162*Reout^0.61*air.Pr^(1/3)*air.k/tube.Do;
case 'Hakan Karatas'%experimental samples somewhat close to our case
Reout=mairflux*tube.pitchtransverse/(air.u);
j= 0.138*Reout^-0.281*(At/Ao)^-0.407;
coeffconvout=j*Reout*air.Pr^(1/3)*air.k/tube.Do;
%300<Re<1000 and 1<At/Ao<6
case 'Jack R. Barbosa Jr.'%for low fin density
j= 0.5685*Reout^-0.4446*(At/Ao)^-0.3824;
coeffconvout=j*mairflux*air.cp/air.Pr^(2/3);
%320<Reair<1200, 2.6<At/Ao<5.8 and 2<N/2<5
end

% fin HT calculation
%eta fin
m=sqrt(2*coeffconvout/(fin.k*fin.t));
etafin=tanh(phi*m*(tube.Dc/2))/(phi*m*(tube.Dc/2));
etasurf=1-tube.elementlength/fin.p*Af/At*(1-etafin);
end

%% single-phase refrigerant side heat transfer and pressure drop
function [f coeffconvin]=href(numref, Rein, ref, tube);
f=(1.58*log(Rein)-3.28)^(-2);
switch numref
case 'Gnielinski'
coeffconvin=(f/2*(Rein-1000)*ref.Pr)/(1+12.7*(f/2)^0.5*(ref.Pr^(2/3)-1))*ref.k/(tube.Di);
case 'Petukhov'
coeffconvin=(f/2*Rein*ref.Pr)/(1.07+12.7*(f/2)^0.5*(ref.Pr^(2/3)-1))*ref.k/(tube.Doi);
case 'Dittus & Boelter'
coeffconvin= 0.023*Rein^(4/5)*ref.Pr^0.4*ref.k/(tube.Doi);%0.7<=Pr<=160;Re>=10000;L/D>=10 (for small to moderate temperature difference)
end
end

%% two-phase refrigerant side heat transfer and pressure drop
function [dPevaptube twophasecoefflocalevap flow]=h2ph(xin,mflux,qfluxe,Pevaptubein,numrefevap,pdropcorr,tube,ref);
Pcritical=4.9012e3;molwt=72.585;Preduced=Pevaptubein/Pcritical;%for R410a
switch numrefevap
%   case'Gungor Winterton'
%    Bo=qfluxe/(mflux*ref.hfg);
%    f=(0.790*log(Reliq)-1.64)^(-2);
%    coeffliq=(f/8*(Reliq-1000)*ref.Prliq)/(1+12.7*(f/8)^0.5*(ref.Prliq^(2/3)-1))*ref.kliq/(tube.Di);
%    coefftwophaselocalevap=1+3000*Bo^0.86+1.12*(xin/(1-xin))^0.75*(ref.rholiq/ref.rhog)^0.41*coeffliq;
case 'Wojtan'
e=xin/ref.rhog*((1+0.12*(1-xin))*(xin/ref.rhog+(1-xin)/ref.rholiq)+(1.18*(1-xin)*(9.8*ref.surften*(ref.rholiq-ref.rhog))^0.25)/(mflux*ref.rholiq^0.5))^-1;
Al=(1-e)*pi/4*(tube.Di)^2;Ag=e*pi/4*(tube.Di)^2;Ald=Al/tube.Di^2;Agd=Ag/tube.Di^2;
vliq=mflux/ref.rholiq*(1-xin)/(1-e);vg=mflux/ref.rhog*xin/e;
xIA=((.2914*(ref.rhog/ref.rholiq)^-(1/1.75)*(ref.uliq/ref.ug)^-(1/7))+1)^-1;
eIA=xIA/ref.rhog*((1+0.12*(1-xIA))*(xIA/ref.rhog+(1-xIA)/ref.rholiq)+(1.18*(1-xIA)*(9.8*ref.surften*(ref.rholiq-ref.rhog))^0.25)/(mflux*ref.rholiq^0.5))^-1;
if eIA<e
eIA=e;end
AldIA=(1-eIA).*pi./4;AgdIA=eIA.*pi./4;
thetastrat=2*pi-2*(pi*(1-e)+(3*pi/2)^(1/3)*(1-2*(1-e)+(1-e)^(1/3)-e^(1/3))-1/200*(1-e)*e*(1-2*(1-e))*(1+4*((1-e)^2+e^2)));
thetastratIA=2*pi-2*(pi*(1-eIA)+(3*pi/2)^(1/3)*(1-2*(1-eIA)+(1-eIA)^(1/3)-eIA^(1/3))-1/200*(1-eIA)*eIA*(1-2*(1-eIA))*(1+4*((1-eIA)^2+eIA^2)));
hLD=0.5*(1-cos((2*pi-thetastrat)/2));inctheta=0;%inclination of tube  hliq=hLD*tube.Di;
hLDIA=0.5.*(1-cos((2.*pi-thetastratIA)./2));inctheta=0;%inclination of tube  hliq=hLD.*tube.Di;
if xin<xIA
xin1=xIA;
mfluxstrat=((226.3^2*AldIA*AgdIA^2*ref.rhog*(ref.rholiq-ref.rhog)*ref.uliq*9.8*cos(inctheta))/(xin1^2*(1-xin1)*pi^3))^(1/3);
else
mfluxstrat=((226.3^2*Ald*Agd^2*ref.rhog*(ref.rholiq-ref.rhog)*ref.uliq*9.8*cos(inctheta))/(xin^2*(1-xin)*pi^3))^(1/3);
end
qcrit=0.131*ref.rhog^0.5*ref.hfg*(9.8*(ref.rholiq-ref.rhog)*ref.surften)^0.25;
WeFrl=9.8*tube.Di^2*ref.rholiq/ref.surften;
mfluxwavyIA=(((16*AgdIA^3*9.8*tube.Di*ref.rholiq*ref.rhog)/(xIA^2*pi^2*(1-(2*hLDIA-1)^2)^0.5))*(pi^2/(25*hLDIA^2)*(WeFrl)^-1+1/cos(inctheta)))^0.5+50;
mfluxwavy=(((16*Agd^3*9.8*tube.Di*ref.rholiq*ref.rhog)/(xin^2*pi^2*(1-(2*hLD-1)^2)^0.5))*(pi^2/(25*hLD^2)*(WeFrl)^-1+1/cos(inctheta)))^0.5+50;
Weg=mflux^2*tube.Di/(ref.rhog*ref.surften);Frg=mflux^2/(9.8*tube.Di*ref.rhog^2);
xdi=.58*exp(.52-.235*Weg^.17*Frg^.37*(ref.rhog/ref.rholiq)^.25*(qfluxe/qcrit)^.7);
if xdi<.01
xdi=.01;end
xde=.61*exp(.57-5.8e-3*Weg^.38*Frg^.15*(ref.rhog/ref.rholiq)^-.09*(qfluxe/qcrit)^.27);
if xde<.01 ||xde>1
xde=0.99;end
mfluxdryout=(1/.235*(log(.58/xin)+.52)*(tube.Di/(ref.rhog*ref.surften))^-.17*(9.8*tube.Di*ref.rhog*(ref.rholiq-ref.rhog))^.37*(ref.rhog/ref.rholiq)^-.25*(qfluxe/qcrit)^-.7)^.926;
mfluxmist=(1/.0058*(log(.61/xin)+.57)*(tube.Di/(ref.rhog*ref.surften))^-.38*(9.8*tube.Di*ref.rhog*(ref.rholiq-ref.rhog))^.15*(ref.rhog/ref.rholiq)^.09*(qfluxe/qcrit)^-.27)^.943;

if mflux>mfluxwavyIA
thetastratwavy=0;flow='slug';
elseif mflux<mfluxwavyIA && xin<xIA && mflux>mfluxstrat
thetastratwavy=thetastrat*xin/xIA*((mfluxwavy-mflux)/(mfluxwavy-mfluxstrat))^0.61;flow='slug-strat-wavy';
elseif xin>=xIA && xin<1
thetastratwavy=thetastrat*((mfluxwavy-mflux)/(mfluxwavy-mfluxstrat))^0.61;flow='strat-wavy';% x due to multiple use in diff flows in pressure drop
elseif mflux>mfluxstrat
thetastratwavy=thetastrat;flow='strat';
end
if mflux>mfluxwavy && xin<xIA
thetastratwavy=0;flow='intermittent';
elseif mflux>mfluxwavy && xin>xIA
thetastratwavy=0;flow='annular';
end
if mfluxstrat>mfluxdryout
mfluxdryout=mfluxstrat;end
if mfluxwavy>mfluxdryout
mfluxdryout=mfluxwavy;
end
if mfluxdryout>mfluxmist
mfluxmist=mfluxdryout;
end
if mflux>mfluxdryout || xin>xdi
flow='dryout';
end
if mflux>mfluxmist || xin>xde
flow='mist';
end

liqfilmthick=tube.Di/2-sqrt((tube.Di/2)^2-2*Al/(2*pi-thetastratwavy));
if liqfilmthick>tube.Di/2
liqfilmthick=tube.Di/2;
end

switch flow
case 'dryout'
Reliq=4*mflux*(1-xdi)*liqfilmthick/((1-e)*ref.uliq);
coeffnucboil=55*Preduced^.12*(-log10(Preduced))^-.55*molwt^-.5*qfluxe^.67;
coeffconvboil=0.0133*Reliq^0.69*ref.Prliq^0.4*ref.kliq/liqfilmthick;coeffwet=((.8*coeffnucboil)^3+coeffconvboil^3)^(1/3);
coeffvap=0.023*(mflux*xdi*tube.Di/(e*ref.ug))^.8*ref.Prg^.4*ref.kg/tube.Di;twophasecoeffxdi=(coeffvap*thetastratwavy+(2*pi-thetastratwavy)*coeffwet)/(2*pi);
ReHomo=mflux*tube.Di/ref.ug*(xde+ref.rhog/ref.rholiq*(1-xde));Y=1-.01*((ref.rholiq/ref.rhog-1)*(1-xde))^.4;mistcoeffxde=.0117*ReHomo^.79*ref.Prg^1.06*Y^-1.83*ref.kg/tube.Di;
twophasecoefflocalevap=twophasecoeffxdi-(xin-xdi)/(xde-xdi)*(twophasecoeffxdi-mistcoeffxde);
case 'mist'
ReHomo=mflux*tube.Di/ref.ug*(xin+ref.rhog/ref.rholiq*(1-xin));Y=1-.01*((ref.rholiq/ref.rhog-1)*(1-xin))^.4;
twophasecoefflocalevap=.0117*ReHomo^.79*ref.Prg^1.06*Y^-1.83*ref.kg/tube.Di;
otherwise
Reliq=4*mflux*(1-xin)*liqfilmthick/((1-e)*ref.uliq);
coeffnucboil=55*Preduced^.12*(-log10(Preduced))^-.55*molwt^-.5*qfluxe^.67;
coeffconvboil=0.0133*Reliq^0.69*ref.Prliq^0.4*ref.kliq/liqfilmthick;
coeffwet=((.8*coeffnucboil)^3+coeffconvboil^3)^(1/3);
coeffvap=0.023*(mflux*xin*tube.Di/(e*ref.ug))^.8*ref.Prg^.4*ref.kg/tube.Di;
twophasecoefflocalevap=(coeffvap*thetastratwavy+(2*pi-thetastratwavy)*coeffwet)/(2*pi);
end
end
switch pdropcorr
case 'Flow based'
Rel=mflux*tube.Di/ref.uliq;fliq=(1.58*log(Rel)-3.28)^(-2);%for liquid pressure drop evaluation
Weliq=ref.rholiq*vliq^2*tube.Di/ref.surften;liqfilmthickannular=pi*tube.Di*(1-e)/(4*pi);fannular=.67*(liqfilmthickannular/tube.Di)^1.2*((ref.rholiq-ref.rhog)*9.8*liqfilmthickannular^2/ref.surften)^-.4*(ref.ug/ref.uliq)^.08*Weliq^-.034;%for annular pressure drop evaluation
Reg=mflux*tube.Di/ref.ug*xin/e;fg=.079/(Reg^.25);
fstratwavy=thetastratwavy/(2*pi)*fg+(1-thetastratwavy/(2*pi))*fannular;
switch flow
case 'slug'
dPevaptube=(2*fliq*tube.elementlength*(mflux)^2/(tube.Di*ref.rholiq)*(1-e/eIA)^.25+2*fannular*tube.elementlength*ref.rhog*vg^2/tube.Di*(e/eIA)^.25)/1e3;
case 'slug-strat-wavy'
dPevaptube=(2*fliq*tube.elementlength/(tube.Di*ref.rholiq)*(mflux)^2*(1-e/eIA)^.25+2*fstratwavy*tube.elementlength*ref.rhog*vg^2/tube.Di*(e/eIA)^.25)/1e3;
case 'strat-wavy'
dPevaptube=2*fstratwavy*tube.elementlength*ref.rhog*vg^2/tube.Di/1e3;
case 'strat'
fstratwavy=thetastrat/(2*pi)*fg+(1-thetastrat/(2*pi))*fannular;
dPevaptube=2*fstratwavy*tube.elementlength*ref.rhog*vg^2/tube.Di/1e3;
if xin<xIA
dPevaptube=(2*fliq*tube.elementlength/tube.Di/ref.rholiq*(mflux)^2*(1-e/eIA)^.25+2*fstratwavy*tube.elementlength*ref.rhog*vg^2/tube.Di*(e/eIA)^.25)/1e3;end
case 'intermittent'
dPevaptube=(2*fliq*tube.elementlength*(mflux)^2/(tube.Di*ref.rholiq)*(1-e/eIA)^.25+2*fannular*tube.elementlength*ref.rhog*vg^2/tube.Di*(e/eIA)^.25)/1e3;%intermittent and slug are treated in same manner
case 'annular'
dPevaptube=2*fannular*tube.elementlength*ref.rhog*vg^2/tube.Di/1e3;
case 'dryout'
edi=xdi/ref.rhog*((1+0.12*(1-xdi))*(xdi/ref.rhog+(1-xdi)/ref.rholiq)+(1.18*(1-xdi)*(9.8*ref.surften*(ref.rholiq-ref.rhog))^0.25)/(mflux*ref.rholiq^0.5))^-1;
vliq=mflux/ref.rholiq*(1-xdi)/(1-edi);vg=mflux*xdi/(ref.rhog*edi);
Weliq=ref.rholiq*vliq^2*tube.Di/ref.surften;
fannular=.67*(liqfilmthickannular/tube.Di)^1.2*((ref.rholiq-ref.rhog)*9.8*liqfilmthickannular^2/ref.surften)^-.4*(ref.ug/ref.uliq)^.08*Weliq^-.034;
dpdi=2*fannular*tube.elementlength*ref.rhog*vg^2/tube.Di;
eh=1/(1+(1-xde)/xde*(ref.rhog/ref.rholiq));ref.rhohomo=ref.rholiq*(1-eh)+ref.rhog*eh;
Rede=mflux*tube.Di/(xde*ref.ug+(1-xde)*ref.uliq);fde=.079/(Rede^.25);dpde=2*fde*tube.elementlength*mflux^2/(tube.Di*ref.rhohomo);
dPevaptube=(dpdi-(xin-xdi)/(xde-xdi)*(dpdi-dpde))/1e3;
case 'mist'
Remist=mflux*tube.Di/(xin*ref.ug+(1-xin)*ref.uliq);fmist=.079/(Remist^.25);eh=1/(1+((1-xin)/xin*(ref.rhog/ref.rholiq)));ref.rhohomo=ref.rholiq*(1-eh)+ref.rhog*eh;
dPevaptube=2*fmist*tube.elementlength*mflux^2/(tube.Di*ref.rhohomo)/1e3;
end
end
end```