clc
double deltapGSaez(i)
%read in all the values by experiment from excel-file
in = xlsread('C:\Dokumente und Einstellungen\ben\Eigene Dateien\VT\Diplomarbeit\models\holdup-pressure_drop\pressure-drop\draft_pressure_drop_10ppi_1.xls','10ppi','A6:R21');
%coefficient
m1=3;
m2=3;
f1=20;
f2=0.9;
sigmaL=in(:,1);
A=in(:,2);
B=in(:,3);
ds=in(:,4);
rohG=in(:,5);
rohL=in(:,6);
muG=in(:,7);
muL=in(:,8);
uG=in(:,9);
uL=in(:,10);
dw=in(:,11);
eps=in(:,12);
ReL=in(:,13);
ReG=in(:,14);
ssv=in(:,15);
pexp=in(:,16);
pemp=in(:,17);
hL=in(:,18);
g=9.81;
%calculate the hydrodynamics for each velocity combination
for i=1:1:16
m3(i)=-3.6;%???
%equivalent diameter calculation    
de(i)=1.5.*(1-eps(i)).*(dw(i)+ds(i))./eps(i);
%Galileo number calculation
GaG(i)=rohG(i).^2.*g.*dw(i).^3.*eps(i).^3./(muG(i).^2.*(1-eps(i)).^3);
GaL(i)=rohL(i).^2.*g.*dw(i).^3.*eps(i).^3./(muL(i).^2.*(1-eps(i)).^3);
%saturation
SL(i)=hL(i)./eps(i);
SG(i)=1-hL(i)./eps(i);
%Eoetoes Number
eoe(i)=rohL(i).*g.*dw(i).^2.*eps(i)./(sigmaL(i).*(1-eps(i)).^2);
%saturation based on static holdup
SL0(i)=1./(eps(i).*(f1+f2.*eoe(i)));
%permeability
kL(i)=((SL(i)-SL0(i))./(1-SL0(i))).^m1;
kG(i)=SG(i).^m2;
%gas holdup
hG(i)=eps(i)-hL(i);

%Saez&Carbonell 1986 AIchE Journal March 32, no.3 "The hydrodynamics of trickling flow in packed beds; Part I: Conduit models"
%from equation 4 
%-deltap
dimlessdeltapGSaez(i)=(A(i).*ReG(i)./GaG(i)+B(i).*ReG(i).^2./GaG(i))./(SG(i).^m2);
%inserted in equation 2
deltapGSaez(i)=-hG(i).*g.*rohG(i)+dimlessdeltapGSaez(i).*rohG(i).*g.*hG(i);
%durch hG teilen auf der rechten Seite?? evtl. zu viel im Vergleich zu Holub und Edouard??
%not taken

%modified Saez& Carbonell 1986 Journal March 32, no.3 "The hydrodynamics of trickling flow in packed beds; Part I: Conduit models"
%from equation 4 
% with normal momentum equation
deltapGSaezmod(i)=-g.*rohG(i)+dimlessdeltapGSaez(i).*rohG(i).*g;
% leads to the same as per Edouard 2008

%Edouard 2008
%-deltap
%liquid drag force from equation 5
FLEdouard(i)=hL(i).*(A(i).*ReL(i)./GaL(i)+B(i).*ReL(i).^2./GaL(i)).*rohL(i).*g./kL(i);
%gas drag force
FGEdouard(i)=hG(i).*(A(i).*ReG(i)./GaG(i)+B(i).*ReG(i).^2./GaG(i)).*rohG(i).*g./kG(i);
%liquid pressure drop delta p/H
%from equation 7
deltapLEdouard(i)=-rohL(i).*g+FLEdouard(i)./hL(i);
%gas pressure drop delta p/H
deltapGEdouard(i)=-rohG(i).*g+FGEdouard(i)./hG(i);

%Al-Dahhan, Iliuta 1998 and Holub 1993
%from equation 2 Al-Dahhan basing on Holub slit model equation 16, p. 310
deltapGHolub(i)=(eps(i)./(eps(i)-hL(i))).^3.*(A(i).*ReG(i)./GaG(i)+B(i).*ReG(i).^2./GaG(i)).*rohG(i).*g-rohG(i).*g;
%extended slit model by Al-Dahhan
% expected to be not applicable as at atmospheric pressure: fv and fs are 0.
%Al-Dahhan 1998 p.796 extended slit model even though not foreseen for atmospheric
%dimensionless body force
%interaction Reynoldsnumber equation 5, p.796
Rei(i)=uL(i).*de(i)./(muL(i).*(1-eps(i)));
%interaction parameter exponents originally 0.05 and -0.05 
%(uncertain - to be modified)
fv(i)=-2.3.*ReG(i).^0.05.*ReL(i).^-0.05;
bodyforceGAlDahhan(i)=(eps(i)./(eps(i)-hL(i))).^3.*(A(i).*(ReG(i)-fv(i).*hG(i).*Rei(i))./GaG(i)+B(i).*(ReG(i)-fv(i).*hG(i).*Rei(i)).^2./GaG(i));
deltapGAlDahhan(i)=bodyforceGAlDahhan(i).*g.*rohG(i)-rohG(i).*g;
%hardly any change by modifying the fv exponents. Model is far away from
%experimental results, - not to be used

%Nemec&Levec model of gas pressure drop
%relative permeability
kGNemec(i)=kG(i).*m3(i);
%pressure drop
deltapGNemec(i)=(A(i).*ReG(i)./GaG(i)+B(i).*ReG(i).^2./GaG(i)).*rohG(i).*g./kGNemec(i);

%modified Nemec&Levec model of gas pressure drop
%pressure drop
deltapGNemecModified(i)=-rohG(i).*g+deltapGNemec(i)

end
%xlswrite('C:\Dokumente und Einstellungen\ben\Eigene Dateien\VT\Diplomarbeit\models\holdup-pressure_drop\pressure-drop\draft_pressure_drop_10ppi_1.xls','10ppi','S6:S21',peltapGSaez(i));