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

[Tcout Pcondtubein balance1 xdesuper xcond econd Vair]=fintubecondenser(mdoti,Vair,kW,Qa,Tx,Tdis,Pdis,Vairsolve);
```%{
Input:
mdoti=refrigerant mass flow rate (kg/s)
Vair=condenser air volumetric flow rate (initial guess if solving for Vair) (m3/s)
kW=compressor power input (kW)
(kW+Qa is the amount of heat that needs to be rejected by the condenser. used in energy balance equation)
Tx=outdoor temperature (K)
Tdis=compressor discharge temperature (K)
Pdis=compressor discharge pressure (kPa)
Vairsolve=flag to solve for minimum Vair that satisfies energy balance for a given set of parameters (1 or 0)

Output:
Tcout=condenser outlet temperature (K)
Pcondtubein=condenser outlet pressure (kPa)
balance1=energy balance convergence (kW)
xdesuper=desuperheating portion (%)
xcond=condensation portion (%)
econd= condenser effectiveness (%)
Vair=condenser air volumetric flow rate (m3/s)
%}
function [Tcout Pcondtubein balance1 xdesuper xcond econd Vair]=fintubecondenser(mdoti,Vair,kW,Qa,Tx,Tdis,Pdis,Vairsolve);
p=Tchparm;p.dP_HX=1;%flag if want to include pressure drop
%%
disp('condenser');
numair='Gray and Webb';numref='Gnielinski'; % correlation names
tube.Tx=Tx;p.Q=Qa;Pcondtubein=Pdis;%outside HX conditions
Tsat=tube.Tx;tubecheck=0;[airrho aircp]=refpropm('DC','T',tube.Tx,'P',p.Pamb,p.airname);

% tube parameters
tube.vertical=p.cond.vertical;tube.width=p.cond.width;
tube.N=(p.cond.Ntot-p.cond.tube.Nliq)/p.cond.tube.div+p.cond.tube.Nliq;tube.div=(p.cond.Ntot-p.cond.tube.Nliq)/p.cond.tube.div;
tube.row=p.cond.tube.row;
tube.pitchtransverse=p.cond.tube.pitchtransverse;
tube.length=p.cond.tube.length;tube.lengthtotal=tube.N*tube.length+pi*tube.pitchtransverse/2*(tube.N-1);
tube.Do=p.cond.tube.Do;
tube.t=p.cond.tube.t;tube.Di=tube.Do-2*tube.t;
tube.k=p.cond.tube.k;
tube.elementlength=tube.length;

% fin parameters
fin.t=p.cond.fin.t;tube.Dc=tube.Do+2*fin.t;%wang paper
fin.w=p.cond.fin.w/tube.row;tube.pitchlong=fin.w/2;% fin.w=longitudnal pitch assumption
fin.p=p.cond.fin.p;
fin.k=p.cond.fin.k;

if Vairsolve==1
options_c = optimset('Algorithm','interior-point','Display','off','TolFun', 0.001);
problem=createOptimProblem('fmincon','objective',@mair_solver,'x0',Vair,'lb',.2,'ub',.67,'options',options_c);
bigstart = CustomStartPointSet([.67:-.5:.2]');
gs=MultiStart('Display','off','UseParallel','always');
Vair=run(gs,problem,bigstart);
end

%% condenser loop
balance=mair_solver(Vair);
function balance=mair_solver(Vair);
Pcondtubein=Pdis;Ttubein=Tdis;Tsat=tube.Tx;mdot=mdoti/p.cond.tube.div;tube.elementlength=tube.length;tubelooplength=tube.elementlength;Ttubeout=Ttubein;%resetting variables in loop
i=1;j=1;k=1;xin=.999;Pcondboil=zeros(tube.N*2,1);Pcond=zeros(tube.N*2,1);Qcond=zeros(tube.N*2,1);Qcondboil=zeros(tube.N*2,1);jj=1;minc=0;flage=0;tpcalc=0;tubecheck=0;trcheck=0;trchecksub=0;flagmerge=0;%resetting variables in loop
%trcheck means transition check, tpcalc means two-phase calc
mflux=mdot/(pi/4*(tube.Di)^2);xoillimit=.9;mair=Vair.*airrho;

% condenser loop
while tubelooplength<tube.lengthtotal

% 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 liquid line as two streams are merging so mdot will change
if tubelooplength>tube.length*tube.div && minc==0
mdot=mdot*p.cond.tube.div;mflux=mdot/(pi/4*(tube.Di)^2);minc=1;flagmerge=1;
end

% for calc of local oil concentration
if p.oilcon~=0
if tpcalc==0 && trcheck==0
oilconlocal=0;
elseif xin>0 && xin<xoillimit
oilconlocal=p.oilcon/(1-xin);
else
oilconlocal=p.oilcon/(1-0);
end
else oilconlocal=0;
end

% for calc of saturation temperature for the condensation part ***
if tpcalc==0 %&& trcheck==0
Tsat=refpropm('T','P',Pcondtubein,'Q',1,p.Refname);
else
if xin>0 && xin<xoillimit
Tsat=refpropm('T','P',Pcondtubein,'Q',xin,p.Refname);%due to temperature glide in R410A
T1=refpropm('T','P',Pcondtubein,'Q',xin,p.Refname);T2=refpropm('T','P',Pcondtubein+1,'Q',xin,p.Refname);
equations=[1 T1;1 T2];values=[log(Pcondtubein)*T1;log(Pcondtubein+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(Pcondtubein)-B);
else
Tsat=refpropm('T','P',Pcondtubein,'Q',0,p.Refname);
if oilconlocal~=0
T1=refpropm('T','P',Pcondtubein,'Q',0,p.Refname);T2=refpropm('T','P',Pcondtubein+1,'Q',0,p.Refname);
equations=[1 T1;1 T2];values=[log(Pcondtubein)*T1;log(Pcondtubein+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(Pcondtubein)-B);
end
end
end %***

% check whether to do desuperheating HT or condensation HT
if (round(Ttubein(i)*100)/100>(round((Tsat+.01)*100)/100) && tpcalc==0) || tpcalc==1
% refrigerant properties calculation**
if tpcalc==1 && trchecksub==1 %only true when just came out of condensation phase
[ref.k ref.u]=refpropm('LV','T',round(Ttubein(i)),'Q',0,p.Refmix1,p.Refmix2,p.Refmixcomp);
[ref.rho,ref.cp]=refpropm('DC','T',Ttubein(i),'Q',0,p.Refname);trchecksub=0;
else

[ref.k,ref.u]=refpropm('LV','T',round(Ttubein(i)*10)/10,'P',Pcondtubein,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',0,p.Refmix1,p.Refmix2,p.Refmixcomp);
if ref.k<0
[ref.k ref.u]=refpropm('LV','T',round(Ttubein(i)*100)/100+1,'Q',0,p.Refmix1,p.Refmix2,p.Refmixcomp);end
end
[ref.rho,ref.cp]=refpropm('DC','T',Ttubein(i),'P',Pcondtubein,p.Refname);
if ref.cp<0
[ref.rho ref.cp]=refpropm('DC','T',round(Ttubein(i)),'P',Pcondtubein,p.Refname);
if ref.cp<0
[ref.rho ref.cp]=refpropm('DC','T',round(Ttubein(i)*100)/100,'Q',0,p.Refname);end
end
end

% ref-oil mixture properties calculation
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 desuperheating 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*(Ttubein(i)-tube.Tx);Ttubeout(i)=Ttubein(i)-Q/Cref;

% for accurately identifying the location of transition from superheat to condensation
while (Ttubeout(i))<Tsat  && trcheck==0
tube.elementlength=tube.elementlength/2; %for refinement of discretization to get more accurately the point of end of condensation
%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*(Ttubein(i)-tube.Tx);Ttubeout(i)=Ttubein(i)-Q/Cref;

if p.dP_HX==1
% Pressure drop calculation
Pcondtube=(Pcondtubein*1e3-2*f*tube.elementlength*(mflux)^2/(tube.Di*ref.rho))/1e3;
Tsat=refpropm('T','P',Pcondtubein,'Q',1,p.Refname);
if oilconlocal~=0
T1=refpropm('T','P',Pcondtubein,'Q',1,p.Refname);T2=refpropm('T','P',Pcondtubein+1,'Q',1,p.Refname);
equations=[1 T1;1 T2];values=[log(Pcondtubein)*T1;log(Pcondtubein+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(Pcondtubein)-B);
end
else Pcondtube=Pcondtubein;
end
% for marking end of transition region accurate to within .1K
if (Ttubeout(i)+.01)>Tsat && (Ttubeout(i)-.01)<Tsat
trcheck=1;
Pcondtubein=Pcondtube;
% for calculation of ref quality at beginning of condensation (This is due to element length resolution limit)
if  (Ttubeout(i))<Tsat
ref.hliq=refpropm('H','T',Tsat,'Q',0,p.Refmix1,p.Refmix2,p.Refmixcomp);
ref.hg=refpropm('H','T',Tsat,'Q',1,p.Refmix1,p.Refmix2,p.Refmixcomp);ref.hfg=(ref.hg-ref.hliq);
if tpcalc==0
xin=1-Q/(mdot*ref.hfg);end
if xin>.99
xin=.99;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
Pcondtubein=(Pcondtubein*1e3-2*f*tube.elementlength*(mflux)^2/(tube.Di*ref.rho))/1e3;
end

% for resetting element length after calc of transition from desuperheating to condensation
if tube.elementlength<tube.length && tpcalc~=1 && trcheck==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 subcooling after calc. of point of condensation end and 1 calc. of subcool
end

%for streams operating in same conditions so multiply by
%the number of streams to get total Q for use in energy balance satisfaction
if flagmerge==0
Q=Q*p.cond.tube.div;end

% loop safety
if Pcondtubein<p.Ptp || (Ttubein(i)<tube.Tx && tpcalc==1)
Qcond(k)=Q;
break;
end

%condensation part
elseif tpcalc==0
xdesuper=tubelooplength/tube.lengthtotal;
while xin>0

% loop safety
if Pcondtubein<p.Ptp
break; end

% for stream merging while condensation HT haven't finished
if tubelooplength>tube.length*tube.div && minc==0
mdot=mdot*p.cond.tube.div;mflux=mdot/(pi/4*(tube.Di)^2);minc=1;flagmerge=1;
end

Tsat=refpropm('T','P',Pcondtubein,'Q',xin,p.Refname); %due to temperature glide in R410A
if p.oilcon~=0
if xin<xoillimit
oilconlocal=p.oilcon/(1-xin);
else oilconlocal=0;
end
if oilconlocal~=0
T1=refpropm('T','P',Pcondtubein,'Q',xin,p.Refname);T2=refpropm('T','P',Pcondtubein+1,'Q',xin,p.Refname);
equations=[1 T1;1 T2];values=[log(Pcondtubein)*T1;log(Pcondtubein+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(Pcondtubein)-B);
end
else oilconlocal=0;
end
% loop safety
if Tsat<tube.Tx
break; end

% refrigerant liquid and gas phase properties
[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);

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
ref.Prliq=(ref.uliq*ref.cpliq)/ref.kliq;ref.Prg=(ref.ug*ref.cpg)/ref.kg;

if p.dP_HX==1
beforepdrop=Pcondtubein;pdropcorr='flowtype based';%flowtype based for pdrop
else
pdropcorr='None';
end

% loop to calculate two-phase heat transfer and pressure drop as qflux not known
qfluxc=(Cair*(Tsat-tube.Tx))/(pi/4*tube.Di^2);qfluxcn=qfluxc*100;qc=[1 1 1];flowp='a';flowpp='a';
for qloop=1:30
if qc(2)~=qc(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
qfluxcn=qfluxc;
[dPcondtube twophasecoefflocalcond flowtype]=h2ph(xin,mflux,qfluxc,Pcondtubein,pdropcorr,tube,ref);
UA=(1/(twophasecoefflocalcond*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);
qfluxc=eff*Cmin*(Tsat-tube.Tx)/(pi/4*tube.Di^2);
flowpp=flowp;flowp=flowtype;
if strcmp(flowp,flowpp) && qc(1)~=1
qc(2)=1;qc(3)=1;
elseif qc(1)~=1
qc(2)=2;qc(3)=1;
end
if abs(qfluxc-qfluxcn)/qfluxcn<.01
break;
end
qc(1)=2;
end

if p.oilcon~=0
if ~strcmp('mist',flowtype)
Pcondtubein=Pcondtubein-dPcondtube*(oil.u/ref.uliq)^(.184*oilconlocal); % oil correction from Thome paper
else
Pcondtubein=Pcondtubein-dPcondtube;
end
else
Pcondtubein=Pcondtubein-dPcondtube;
end

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

% for accurately identifying the location of transition from condensation to subcooling
if (xout+.01)<=0
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

%for streams operating in same conditions so multiply by
%the number of streams to get total Q
if flagmerge==0
Q=Q*p.cond.tube.div;end
xin=xout;Pcondboil(j)=Pcondtubein;Qcondboil(j)=Q;j=j+1;Ttubeout(i)=Tsat;Ttubein(i+1)=Ttubeout(i);i=i+1;
% proceed to next element
tubelooplength=tubelooplength+tube.elementlength;
end

% for exit from condensation loop
if xin<=0
if Pcondtubein<p.Ptp || Tsat<tube.Tx
break; end;
tubecheck=1;xcond=(tubelooplength-xdesuper*tube.lengthtotal)/tube.lengthtotal;
tpcalc=1; % to identify have been in condensation loop
trcheck=1;trchecksub=1;%indicates done calc. of transition region from condensation to subcooling

% for resetting element length after calc of transition from condensation to subcooling
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;xcond=(tubelooplength-xdesuper*tube.lengthtotal)/tube.lengthtotal;tpcalc=1;
break;
end

% for estimating evaporator pressure drop for guess of Pevapin
if xin>.5 && p.dP_HX==1
pdropevapsafe=beforepdrop-Pcondtubein;end
if xin<=0.5 && flage==0 && p.dP_HX==1
pdropevap=beforepdrop-Pcondtubein;flage=1;%guess for pdrop for evaporator
if pdropevap~=0
pdropevapsafe=pdropevap;end
if pdropevap==0
pdropevap=pdropevapsafe;end
if minc==1
pdropevap=pdropevap/p.evap.tube.div;end%in evaporator two inlets
end
end
end
% for HT and pressure drop balance
if tpcalc==1 && jj==1
jj=0;Pcond(k:k+length(Pcondboil)-1)=Pcondboil;
Qcond(k:k+length(Qcondboil)-1)=Qcondboil;k=length(Qcond(Qcond>0))+1;
else
Pcond(k)=Pcondtubein;Qcond(k)=Q;k=k+1;
end

% loop safety
if Tsat<tube.Tx || Pcondtubein<p.Ptp
disp('Tsat<Tx || Pcondtubein<0');
break;
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+kW-sum(Qcond)/1e3);
balance=abs(p.Q+kW-sum(Qcond)/1e3);
else balance=p.Q+kW+1-sum(Qcond)/1e3;balance1=balance; %to know from balance value if parameters are not valid i.e. proper condenser phenomenon didn't occur
end
else
balance=(p.Q+kW-sum(Qcond)/1e3);balance1=balance;
end
end
Tcout=Ttubeout(end);
%{
% include Psuc in input if want to use this portion of code
% for calculation of hein
% if balance*1e3<100
if Tcout==Tsat
hcout=refpropm('h','T',Tcout,'Q',0,p.Refmix1,p.Refmix2,p.Refmixcomp);
else
hcout=refpropm('h','T',Tcout,'P',Pcondtubein,p.Refmix1,p.Refmix2,p.Refmixcomp);
end
if p.dP_HX==1
pdropevaptot=pdropevap*p.evap.Ntot/p.evap.tube.div;Pein=Psuc+pdropevaptot;
else Pein=Psuc;pdropevaptot=0;
end
tref=refpropm('t','h',200e3,'s',1e3,p.Refname); %IIR reference
funh=@(T) 4.186.*(.388+.00045.*(1.8.*T+32))./((.97386-6.91473e-4.*T)*1e3/998.5).^.5; %kJ/kg K
hcout=p.oilcon*hoilcout+(1-p.oilcon)*hcout;
hein=hcout; %assume isenthalpic expansion
[Tei xini]=refpropm('TQ','H',hein,'P',Pein,p.Refname);
%}
if balance*1e3<(p.Q+kW)*1e3*.1
disp('converged');end
econd=(p.Q+kW-balance1)./(mair.*aircp.*(Tdis-Tx)/1e3)*100;xdesuper=xdesuper*100;xcond=xcond*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);%Af=2(w*(L+t/2)-pi*r^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;

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 temp and tempwater is dewpoint temp
case 'Fischer Rice'%for 3000<Rew<15000 not for low mass flowtype rate of air
Rew=mairflux*tube.pitchtransverse/air.u;
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
j= 0.14*Reout^-0.328*(tube.pitchtransverse/tube.pitchlong)^-0.502*(fin.p/tube.Dc)^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 some what 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 [dPcondtube twophasecoefflocalcond flowtype]=h2ph(xin,mflux,qfluxc,Pcondtubein,pdropcorr,tube,ref);
Pcritical=4.9012e3;molwt=72.585;Preduced=Pcondtubein/Pcritical;%for R410a
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*(qfluxc/qcrit)^.7);
if xdi<.01
xdi=.01;end
xde=.61*exp(.57-5.8e-3*Weg^.38*Frg^.15*(ref.rhog/ref.rholiq)^-.09*(qfluxc/qcrit)^.27);
if xde<.01
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*(qfluxc/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*(qfluxc/qcrit)^-.27)^.943;

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

% if (strcmp(flowtype,'dryout') || strcmp(flowtype,'mist')) && xde<xdi
%     xde=.99;flowtype='annular';thetastratwavy=0;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
Reliq=4*mflux*(1-xin)*liqfilmthick/((1-e)*ref.uliq);fi=1+(vg/vliq)^.5*((ref.rholiq-ref.rhog)*9.8*liqfilmthick^2/ref.surften)^.25;
if mflux<mfluxstrat
fi=1+(vg/vliq)^.5*((ref.rholiq-ref.rhog)*9.8*liqfilmthick^2/ref.surften)^.25*(mflux/mfluxstrat);
end
switch flowtype
case 'dryout'
Reliq=4*mflux*(1-xdi)*liqfilmthick/((1-e)*ref.uliq);
coeffnucboil=55*Preduced^.12*(-log10(Preduced))^-.55*molwt^-.5*qfluxc^.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;
twophasecoefflocalcond=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;
twophasecoefflocalcond=.0117*ReHomo^.79*ref.Prg^1.06*Y^-1.83*ref.kg/tube.Di;
otherwise
coeffannular=0.003*Reliq^0.74*ref.Prliq^0.5*ref.kliq/liqfilmthick*fi;
coefffilmcond=0.655*(ref.rholiq*(ref.rholiq-ref.rhog)*9.8*ref.hfg*ref.kliq^3/(ref.uliq*tube.Di*qfluxc))^(1/3);
twophasecoefflocalcond=(coefffilmcond*thetastratwavy+(2*pi-thetastratwavy)*coeffannular)/(2*pi);
end

switch pdropcorr
case 'Friedel'
ref.rhohomo=(xin/ref.rhog+(1-xin)/ref.rholiq)^-1;Reliq=mflux*tube.Di/ref.uliq;fliq=.079/(Reliq^.25);Reg=mflux*tube.Di/ref.ug;fg=.079/(Reg^.25);
E=(1-xin)^2+xin^2*ref.rholiq*fg/(ref.rhog*fliq);F=xin^.78*(1-xin)^.224;H=(ref.rholiq/ref.rhog)^.91*(ref.ug/ref.uliq)^.19*(1-(ref.ug/ref.uliq))^.7;
Fr=mflux^2/(ref.rhohomo^2*9.8*tube.Di);Weliq=mflux^2*tube.Di/(ref.surften*ref.rhohomo);twophasefricfact=E+3.24*F*H/(Fr^.045*Weliq^.035);
dPcondtube=2*fliq*tube.elementlength/tube.Di/ref.rholiq*(mflux)^2*twophasefricfact^2/1e3;
case 'flowtype 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 flowtype
case 'slug'
dPcondtube=(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'
dPcondtube=(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'
dPcondtube=2*fstratwavy*tube.elementlength*ref.rhog*vg^2/tube.Di/1e3;
case 'strat'
fstratwavy=thetastrat/(2*pi)*fg+(1-thetastrat/(2*pi))*fannular;
dPcondtube=2*fstratwavy*tube.elementlength*ref.rhog*vg^2/tube.Di/1e3;
if xin<xIA
dPcondtube=(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'
dPcondtube=(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'
dPcondtube=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);
dPcondtube=(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;
dPcondtube=2*fmist*tube.elementlength*mflux^2/(tube.Di*ref.rhohomo)/1e3;
end
end
end```