Vapor compression cycle component models

by

 

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)
Qa=evaporator load (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
    hoilcout=200e3+quadl(funh,tref-273.15,Tcout-273.15)*1e3;
    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

Contact us