Vapor compression cycle component models

by

 

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

Contact us