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 Pevapplatein Twaterout balance1 xevap eevap Vwater]=plateevaporator(mdot,Vwater,Qa,Tchwin,Tevapin,Pei,xini,Vwatersolve)
```%{
Input:
mdot=refrigerant mass flow rate (kg/s)
Vwater=evaporator water volumetric flow rate (initial guess if solving for Vwater) (m3/s)
(Qa is the amount of heat that needs to be gained by the evaporator. used in energy balance equation)
Tchwin=chilled water inlet temperature (K)
Tevapin=evporator inlet temperature (K)
Pei=evaporator inlet pressure (kPa)
xini=evaporator inlet quality of refrigerant (0-1)
Vwatersolve=flag to solve for minimum Vwater that satisfies energy balance for a given set of parameters (1 or 0)

Output:
Tevapout=evaporator outlet temperature (K)
Pevapplatein=evaporator outlet pressure (kPa)
balance1=energy balance convergence (kW)
xevap=evaporation portion (%)
eevap=evaporator effectiveness (%)
Vwater=evaporator water volumetric flow rate (m3/s)
%}
function [Tevapout Pevapplatein Twaterout balance1 xevap eevap Vwater]=plateevaporator(mdot,Vwater,Qa,Tchwin,Tevapin,Pei,xini,Vwatersolve)
p=Tchparm;p.dP_HX=1;%flag if want to include pressure drop
disp('plate_evaporator');

p.Q=Qa;%outside HX conditions
platecheck=0;[rhowater cpwater]=refpropm('DC','T',Tchwin,'P',p.Pamb,p.watername);

% chiller geometry
chiller.plates = p.chiller.plates/2;
chiller.plate_length = p.chiller.plate_length;
chiller.plate_width = p.chiller.plate_width;
chiller.plate_thick=p.chiller.plate_thick;
chiller.plate.k=p.chiller.plate.k;
chiller.developed_length= p.chiller.developed_length;
chiller.protracted_length = p.chiller.protracted_length;
chiller.channel_thickness=p.chiller.channel_thickness;
chiller.corrugation_pitch=p.chiller.corrugation_pitch;
chiller.chevron_angle=p.chiller.chevron_angle;
chiller.dia_water =p.chiller.dia_water;
chiller.dia_ref = p.chiller.dia_ref;

chiller.areaenlarge=chiller.developed_length/chiller.protracted_length;% enlargement factor
chiller.channelarea=chiller.channel_thickness*chiller.plate_width;
chiller.parameterwet=2*(chiller.channel_thickness+chiller.areaenlarge*chiller.plate_width);
chiller.diaeq=4*chiller.channelarea/chiller.parameterwet;
elementdiv=10;chiller.elementlength=chiller.plate_length/elementdiv; %number of starting elements

% pressure loss at port on ref side
oilconlocal=p.oilcon/(1-xini);
T1=refpropm('T','P',Pei,'Q',abs(xini),p.Refname);T2=refpropm('T','P',Pei+1,'Q',abs(xini),p.Refname);
equations=[1 T1;1 T2];values=[log(Pei)*T1;log(Pei+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(Pei)-B);
refrho=refpropm('D','P',Pei,'Q',0,p.Refname);rhog=refpropm('D','P',Pei,'Q',1,p.Refname);
oilrho=(0.97386-6.91473e-4*(Tsat-273))*1e3; %kg/m3
rholiq=(oilconlocal/oilrho+(1-oilconlocal)/refrho)^-1;
rhom=(xini/rhog+(1-xini)/rholiq)^-1;
dpport= 1.4*(mdot/(pi/4*chiller.dia_ref^2))^2/(2*rhom); % empirical correlation of kakac and liu
Pevapplatein=Pei-dpport/1e3;
mdot=mdot/chiller.plates;mflux=mdot/(chiller.channelarea);

if Vwatersolve==1
options_e = optimset('Algorithm','interior-point','Display','off','TolFun', 0.001);
problem=createOptimProblem('fmincon','objective',@mwater_solver,'x0',Vwater,'lb',.12e-3,'ub',.28e-3,'options',options_e);
bigstart = CustomStartPointSet([.28:-.1:.12]'*1e-3);
gs=MultiStart('Display','off','UseParallel','never');
[Vwater,balance]=run(gs,problem,bigstart);
end

%% evaporator loop
balance=mwater_solver(Vwater);
function balance=mwater_solver(Vwater)
elementdiv=10;chiller.elementlength=chiller.plate_length/elementdiv; %number of starting elements
platelooplength=chiller.elementlength;Pevapplatein=Pei;Tplateinref=Tevapin;Tplateinwater=Tchwin;xin=xini;%resetting variables
Pevapboil=zeros(chiller.plates*2,1);Pevap=zeros(chiller.plates*2,1);Qevap=zeros(chiller.plates*2,1);Qevapboil=zeros(chiller.plates*2,1);Tplateoutref=[];Tplateoutwater=[];i=1;j=1;k=1;jj=1;tpcalc=0;trcheck=0;platecheck=0;trchecksub=0;%resetting variables
xinlimit=1-p.oilcon;xoillimit=0.9;mwater=Vwater.*rhowater;mwater=mwater/(chiller.plates);mwaterflux=mwater/chiller.channelarea;
% evaporator loop
while platelooplength<chiller.plate_length

% loop safety
if Pevapplatein<p.Ptp || Tplateinref(i)>Tplateinwater(i)
break;end

% for not exceeding total length and calc. what remains
if (platelooplength+chiller.elementlength)>chiller.plate_length
chiller.elementlength=(chiller.plate_length-platelooplength); % element length changed. need to calc hwater and fin heat transfer
end

% for calc of local oil concentration
if p.oilcon~=0
if tpcalc==0 || xin<1
%oilconlocal=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',Pevapplatein,'Q',abs(xin),p.Refname);
if oilconlocal~=0
T1=refpropm('T','P',Pevapplatein,'Q',abs(xin),p.Refname);T2=refpropm('T','P',Pevapplatein+1,'Q',abs(xin),p.Refname);
equations=[1 T1;1 T2];values=[log(Pevapplatein)*T1;log(Pevapplatein+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(Pevapplatein)-B);
end
else
if xin<xoillimit
Tsat=refpropm('T','P',Pevapplatein,'Q',abs(xin),p.Refname);%this is because of temperature glide in r410a
if oilconlocal~=0
T1=refpropm('T','P',Pevapplatein,'Q',abs(xin),p.Refname);T2=refpropm('T','P',Pevapplatein+1,'Q',abs(xin),p.Refname);
equations=[1 T1;1 T2];values=[log(Pevapplatein)*T1;log(Pevapplatein+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(Pevapplatein)-B);
end
else
Tsat=refpropm('T','P',Pevapplatein,'Q',1,p.Refname);
end
end %***

% for checking whether to do heating if ref. is in liquid phase after expansion
if round(Tplateinref(i))<round(Tsat) && 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(Tplateinref(i)),'Q',1,p.Refmix1,p.Refmix2,p.Refmixcomp);
[ref.rho,ref.cp]=refpropm('DC','T',round(Tplateinref(i)),'Q',1,p.Refname);trchecksub=0;
% Refrigerant Properties **
else
[ref.k,ref.u]=refpropm('LV','T',round(Tplateinref(i)*100)/100,'P',Pevapplatein,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(Tplateinref(i)),'Q',1,p.Refmix1,p.Refmix2,p.Refmixcomp);
if ref.k<0
[ref.k ref.u]=refpropm('LV','T',round(Tplateinref(i)*100)/100+1,'Q',1,p.Refmix1,p.Refmix2,p.Refmixcomp);end
end

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

if p.oilcon~=0
oil=oil_properties(Tplateinref(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; %**

[water.rho,water.k,water.cp,water.u,water.h]=refpropm('DLCVH','T',Tplateinwater(i),'P',p.Pamb,p.watername);
water.uw=refpropm('V','T',Tplateinref(i),'P',p.Pamb,p.watername);water.Pr=(water.u*water.cp)/water.k;
ref.uw=refpropm('V','T',Tplateinwater(i),'P',Pevapplatein,p.watername);
% hwater
[coeffconvwater ~]=hf1ph(mwaterflux,water,chiller);Cwater=mwater*water.cp;

% href
[coeffconvref f]=hf1ph(mflux,ref,chiller);Cref=mdot*ref.cp;

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

% for subcooling/superheating part
%total UA=A*[1/href+t/k+1/hw]^-1;
UA=(1/(coeffconvref)+chiller.plate_thick/(chiller.plate.k)+1/(coeffconvwater))^-1*chiller.parameterwet*chiller.elementlength;
NTU=UA/Cmin;Cr=Cmin/Cmax;eff=(1-exp(-NTU*(1-Cr)))/(1-Cr*exp(-NTU*(1-Cr)));%counterflow
Q=eff*Cmin*(Tplateinwater(i)-Tplateinref(i));Tplateoutref(i)=Tplateinref(i)+Q/Cref;Tplateoutwater(i)=Tplateinwater(i)-Q/Cwater;

% for accurately identifying the location of transition from subcool to evaporation
while (Tplateoutref(i))>Tsat  && trcheck==0 && tpcalc==0
chiller.elementlength=chiller.elementlength/2; %for refinement of discretization to get more accurately the point of start of evaporation

%energy balance
%total UA=A*[1/href+t/k+1/hw]^-1;
UA=(1/(coeffconvref)+chiller.plate_thick/(chiller.plate.k)+1/(coeffconvwater))^-1*chiller.parameterwet*chiller.elementlength;
NTU=UA/Cmin;Cr=Cmin/Cmax;eff=(1-exp(-NTU*(1-Cr)))/(1-Cr*exp(-NTU*(1-Cr)));%counterflow
Q=eff*Cmin*(Tplateinwater(i)-Tplateinref(i));Tplateoutref(i)=Tplateinref(i)+Q/Cref;Tplateoutwater(i)=Tplateinwater(i)-Q/Cwater;

if p.dP_HX==1
% Pressure drop calculation
Pevapplate=(Pevapplatein*1e3-2*f*chiller.elementlength*mflux^2/(ref.rho*chiller.diaeq)-ref.rho*9.8*chiller.elementlength)/1e3;
Tsat=refpropm('T','P',Pevapplatein,'Q',1,p.Refname);
if oilconlocal~=0
T1=refpropm('T','P',Pevapplatein,'Q',1,p.Refname);T2=refpropm('T','P',Pevapplatein+1,'Q',1,p.Refname);
equations=[1 T1;1 T2];values=[log(Pevapplatein)*T1;log(Pevapplatein+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(Pevapplatein)-B);
end
else Pevapplate=Pevapplatein;
end
% for marking end of transition region accurate to within .1K
if (Tplateoutref(i)+.1)>Tsat && (Tplateoutref(i)-.1)<Tsat
trcheck=1;
Pevapplatein=Pevapplate;
% for calculation of ref quality at beginning of evaporation (This is due to element length resolution limit)
if  (Tplateoutref(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

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

% proceed to next element
platelooplength=platelooplength+chiller.elementlength;Tplateinref(i+1)=Tplateoutref(i);Tplateinwater(i+1)=Tplateoutwater(i);i=i+1;

% for resetting element length after calc of transition from subcooling to evaporation
if chiller.elementlength<chiller.plate_length/elementdiv && tpcalc~=1
chiller.elementlength=chiller.plate_length/elementdiv;
else
chiller.elementlength=chiller.elementlength;% for calc. of point of evaporation and for calc. of superheating
end

% loop safety
if Pevapplatein<p.Ptp || (Tplateinref(i)>Tplateinwater(i) && tpcalc==1)
break;end

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

% refrigerant liquid and gas phase properties
Tsat=refpropm('T','P',Pevapplatein,'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',Pevapplatein,'Q',abs(xin),p.Refname);T2=refpropm('T','P',Pevapplatein+1,'Q',abs(xin),p.Refname);
equations=[1 T1;1 T2];values=[log(Pevapplatein)*T1;log(Pevapplatein+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(Pevapplatein)-B);
end
else oilconlocal=0;
end

% loop safety
if Tsat>Tplateinwater(i)
Tsat=refpropm('T','P',Pevapplatein,'Q',abs(xin),p.Refname); 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
ref.rhom=(xin/ref.rhog+(1-xin)/ref.rholiq)^-1;

[water.rho,water.k,water.cp,water.u,water.h]=refpropm('DLCVH','T',Tplateinwater(i),'P',p.Pamb,p.watername);
water.uw=refpropm('V','T',Tplateinref(i),'P',p.Pamb,p.watername);water.Pr=(water.u*water.cp)/water.k;
% hwater
[coeffconvwater ~]=hf1ph(mwaterflux,water,chiller);Cwater=mwater*water.cp;

% loop to calculate two-phase heat transfer and pressure drop as qflux not known
qfluxe=(Cwater*(Tplateinwater(i)-Tsat))/chiller.channelarea;
for qloop=1:30
qfluxen=qfluxe;
[coeffconvref2ph f]=hf2ph(xin,mflux,qfluxe,Pevapplatein,ref,chiller);
UA=(1/(coeffconvref2ph)+chiller.plate_thick/(chiller.plate.k)+1/(coeffconvwater))^-1*chiller.parameterwet*chiller.elementlength;
Cmin=Cwater;NTU=UA/Cmin;eff=1-exp(-1*NTU);
qfluxe=eff*Cmin*(Tplateinwater(i)-Tsat)/chiller.channelarea;
if abs(qfluxe-qfluxen)/qfluxen<.01
break;
end
end
Q=eff*Cmin*(Tplateinwater(i)-Tsat);xout=xin+Q/(mdot*ref.hfg);Tplateoutwater(i)=Tplateinwater(i)-Q/Cwater;

if p.dP_HX==1
% Pressure drop calculation
mfluxeq=mflux*(1-xin+xin*(ref.rholiq/ref.rhog)^.5);dpacc=mfluxeq^2*xout/(ref.rholiq-ref.rhog)-mfluxeq^2*xin/(ref.rholiq-ref.rhog);
Pevapplatein=(Pevapplatein*1e3-2*f*chiller.elementlength*mflux^2/(ref.rhom*chiller.diaeq)-ref.rhom*9.8*chiller.elementlength-dpacc)/1e3;
end
% for accurately identifying the location of transition from evaporation to subcooling
if (xout-.01)>=xinlimit
xout=xin;
chiller.elementlength=chiller.elementlength/2;
else
xin=xout;Pevapboil(j)=Pevapplatein;Qevapboil(j)=Q;j=j+1;Tplateoutref(i)=Tsat;Tplateinref(i+1)=Tplateoutref(i);Tplateinwater(i+1)=Tplateoutwater(i);i=i+1;
% proceed to next element
platelooplength=platelooplength+chiller.elementlength;
if chiller.elementlength<chiller.plate_length/elementdiv
chiller.elementlength=chiller.plate_length/elementdiv;
end
end
tpcalc=1;
% for exit from evaporation loop
if xin>=xinlimit
if Pevapplatein<p.Ptp || Tplateinref(i)>Tplateinwater(i)
break; end;
platecheck=1;xevap=platelooplength/chiller.plate_length;
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 chiller.elementlength<chiller.plate_length/elementdiv
chiller.elementlength=chiller.plate_length/elementdiv;
end
elseif platelooplength>=chiller.plate_length
platecheck=0;xevap=platelooplength/chiller.plate_length;
break;
end
end
end

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)=Pevapplatein;Qevap(k)=Q;k=k+1;
end
end

% energy balance convergence equation
if Vwatersolve==1
% use abs for balance so that fmincon tries to find Vwater that makes balance zero if possible
if platecheck==1
balance1=(p.Q-sum(Qevap)*chiller.plates/1e3);
balance=abs(p.Q-sum(Qevap)*chiller.plates/1e3);
else balance=p.Q+1-sum(Qevap)*chiller.plates/1e3;balance1=balance;
end
else
balance=(p.Q-sum(Qevap)*chiller.plates/1e3);balance1=balance;
end
end

rhog=refpropm('D','T',Tplateoutref(end),'Q',1,p.Refname);
dpport= 1.4*(mdot*chiller.plates/(pi/4*chiller.dia_ref^2))^2/(2*rhog); % empirical correlation of kakac and liu
Pevapplatein=Pevapplatein-dpport/1e3;

Tevapout=Tplateoutref(end);Twaterout=Tplateoutwater(end);pdropevapn=Pei-Pevapplatein;
if balance*1e3<p.Q*1e3*.1
disp('converged');end
eevap=(p.Q-balance1)./(Vwater*rhowater.*cpwater.*(Tchwin-Tevapin)/1e3)*100;xevap=xevap*100;
end

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

%% single-phase heat transfer and pressure drop
function [coeffconv f]=hf1ph(mflux, fluid, chiller);
beta = 90 - chiller.chevron_angle;
phi=chiller.areaenlarge;
Re=mflux*chiller.diaeq/fluid.u;
if Re<10^3 % Wanniarachchi correlation
m = 0.646+0.0011*beta;
Nut = 12.6*beta^(-1.142)*phi^(1-m)*Re^m;
Nul = 3.65*beta^(-0.455)*phi^0.661*Re^0.339;
Nu = (Nul^3+Nut^3)^(1/3)*fluid.Pr^(1/3)*(fluid.u/fluid.uw)^0.17;
coeffconv = Nu*fluid.k/chiller.diaeq; % [W/m2K]

p = 0.00423*beta+0.0000223*beta^2;
ft = 46.6*beta^-1.08*phi^(1+p)*Re^(-p);
fl=1774*beta^-1.026*phi^2*Re^(-1);
f = (fl^3+ft^3)^(1/3);

else % Muley and Manglik correlation
Nu = (0.2668-0.006967*beta+7.244e-5*beta^2)*(20.78-50.94*phi+41.16*phi^2-10.51*phi^3)*...
Re^(0.728+0.0543*sin(pi*beta/45+3.7))*fluid.Pr^(1/3)*(fluid.u/fluid.uw)^0.14;
coeffconv = Nu*fluid.k/chiller.diaeq; % [W/m2K]

f = (2.917-0.1277*beta+2.016e-3*beta^2)*(5.474-19.02*phi+18.93*phi^2-5.341*phi^3)...
*Re^-(0.2+0.0577*sin(pi*beta/45+2.1));
end
end

%% two-phase refrigerant side heat transfer and pressure drop
function [coeffconv2ph f]=hf2ph(xin,mflux,qfluxe,Pevapplatein,ref,chiller);
Pcritical=4.9019e3;molwt=72.585;Preduced=Pevapplatein/Pcritical;%for R410a
Reliq=mflux*chiller.diaeq/ref.uliq;
coeffconvliq= 0.023*Reliq^(4/5)*ref.Prliq^0.4*ref.kliq/(chiller.diaeq);
coeffconvpool= 55*Preduced^0.12*(-1*log10(Preduced))^-0.55*molwt^-0.5*qfluxe^0.67; % Cooper equation
Xtt=((1-xin)/xin)^.9*(ref.rhog/ref.rholiq)^.5*(ref.uliq/ref.ug)^.1;
E=1+24000*(qfluxe/(mflux*ref.hfg))^1.16+1.37*(1/Xtt)^.86;
S=(1+1.15e-6*E^2*Reliq^1.17)^-1;
coeffconv2ph=E*coeffconvliq+S*coeffconvpool;

mfluxeq=mflux*((1-xin)+xin*(ref.rholiq/ref.rhog)^.5);
Reeq=mfluxeq*chiller.diaeq/ref.uliq;
f=23820*Reeq^(-1.12);
end```