Vapor compression cycle component models

by

 

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=evaporator load (kW)
(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

Contact us