Gas Diffusion
clc
clear
KB=1.380658*10^(-23);
M=100;
dx=1;
dy=1;
dz=1;
VV=dx*dy*dz;
AVGADO=6.022*10^23;
P=1.01*10^(5);
RA=0.7;
RB=0.3;
PA=RA*P;
PB=RB*P;
% B=CH4(3.758) + A=AIR(3.711) OR CO2 + AIR
% pV = nRT ;
% T=31+273=304K
% p=1.4*1.01*10^5
% V= 1m^3
% n= 1.4*1.01*10^5*1/8.314*304=55.95moles
% 1 mole contains 6.022*10^23 molecules
% So result 3.369*10^25molecules
% 2) in 1.8m^3 there are in moles 55.95*1.8=100.8moles
% and so 25.2 moles of O2 and 75.6 moles of N2
% mass of O2 25.2*32=806g
% mass of N2 75.6*28=2117g
a=29;
b=44;
for i=1:M
MAAAA(i)=a;
MBBBB(i)=b;
MAA(i)=806;
MBB(i)=2117;
MAAA(i)=MAA(i)/MAAAA(i);
MBBB(i)=MBB(i)/MBBBB(i);
NA(i)=MAAA(i)/RA;
NB(i)=MBBB(i)/RB;
end
for i=1:M
PP(i)=PA;
TT(i)=273+20;
V(i)=VV;
R(i)=8.314;
NMOLES(i)=(PP(i)*V(i))/(R(i)*TT(i));
TOT_NMOLES(i)=NMOLES(i)*AVGADO;
end
SEKMAA=3.711;
SEKMAB=3.758;
SKMA=0.5*(SEKMAA+SEKMAB);
for i=1:M;
T(i)=273+0.00001*randn(1,1);
P(i)=1.01*10^5;
MA(i)=28.97;
NA(i)=1;
VAVG(i)=(8*KB*T(i)/(pi*MA(i)))^(0.5);
ZA(i)=0.25*(TOT_NMOLES(i)/V(i))*VAVG(i);
LAMDA(i)=(((2)^(0.5))*pi*(TOT_NMOLES(i)/V(i))*(SKMA)^2)^(-1);
DAB(i)=(2/3)*((((KB)^3)*T(i)/((pi^3)*MA(i)))^0.5)*(T(i)/(((SKMA)^2)*P(i)));
end
scrsz = get(0,'ScreenSize');
figure('Position',[10 10 1500 700])
subplot 231
plot(DAB,'r')
Xlabel('Number of Samples')
ylabel('DAB (m/s)')
grid on
subplot 232
plot(TOT_NMOLES)
Xlabel('Number of Samples')
ylabel('TOTAL NUMBER OF MOLECULES (Moles)')
grid on
subplot 233
plot(NMOLES)
Xlabel('Number of Samples')
ylabel('NUMBER OF MOLECULES (Moles)')
grid on
subplot 234
plot(T)
Xlabel('Number of Samples')
ylabel('T (k)')
grid on
subplot 235
plot(VAVG)
Xlabel('Number of Samples')
ylabel('VAVG (m/s)')
grid on
pause(10)
close all
clear
KB=1.380658*10^(-23);
M=100;
dx=1;
dy=1;
dz=1;
VV=dx*dy*dz;
AVGADO=6.022*10^23;
P=1.01*10^(5);
RA=0.7;
RB=0.3;
PA=RA*P;
PB=RB*P;
% B=CH4(3.758) + A=AIR(3.711) OR CO2 + AIR
% pV = nRT ;
% T=31+273=304K
% p=1.4*1.01*10^5
% V= 1m^3
% n= 1.4*1.01*10^5*1/8.314*304=55.95moles
% 1 mole contains 6.022*10^23 molecules
% So result 3.369*10^25molecules
% 2) in 1.8m^3 there are in moles 55.95*1.8=100.8moles
% and so 25.2 moles of O2 and 75.6 moles of N2
% mass of O2 25.2*32=806g
% mass of N2 75.6*28=2117g
a=29;
b=44;
for i=1:M
MAAAA(i)=a;
MBBBB(i)=b;
MAA(i)=806;
MBB(i)=2117;
MAAA(i)=MAA(i)/MAAAA(i);
MBBB(i)=MBB(i)/MBBBB(i);
NA(i)=MAAA(i)/RA;
NB(i)=MBBB(i)/RB;
end
for i=1:M
PP(i)=PA;
TT(i)=273+20;
V(i)=VV;
R(i)=8.314;
NMOLES(i)=(PP(i)*V(i))/(R(i)*TT(i));
TOT_NMOLES(i)=NMOLES(i)*AVGADO;
end
SEKMAA=3.711;
SEKMAB=3.758;
SKMA=0.5*(SEKMAA+SEKMAB);
for i=1:M;
T(i)=273+0.00001*randn(1,1);
P(i)=1.01*10^5;
MA(i)=28.97;
NA(i)=1;
VAVG(i)=(8*KB*T(i)/(pi*MA(i)))^(0.5);
ZA(i)=0.25*(TOT_NMOLES(i)/V(i))*VAVG(i);
LAMDA(i)=(((2)^(0.5))*pi*(TOT_NMOLES(i)/V(i))*(SKMA)^2)^(-1);
DAB(i)=(2/3)*((((KB)^3)*T(i)/((pi^3)*MA(i)))^0.5)*(T(i)/(((SKMA)^2)*P(i)));
end
scrsz = get(0,'ScreenSize');
figure('Position',[10 10 1500 700])
subplot 231
plot(DAB,'r')
Xlabel('Number of Samples')
ylabel('DAB (m/s)')
grid on
subplot 232
plot(TOT_NMOLES)
Xlabel('Number of Samples')
ylabel('TOTAL NUMBER OF MOLECULES (Moles)')
grid on
subplot 233
plot(NMOLES)
Xlabel('Number of Samples')
ylabel('NUMBER OF MOLECULES (Moles)')
grid on
subplot 234
plot(T)
Xlabel('Number of Samples')
ylabel('T (k)')
grid on
subplot 235
plot(VAVG)
Xlabel('Number of Samples')
ylabel('VAVG (m/s)')
grid on
pause(10)
close all
Gas Diffusion in Porous Media
Diffusion of gas molecules.
The following code shows different types of gases diffusing
clc
clear
M=100;
t0=0;
t1=30;
T0=273+t0;
T1=300+t1;
DT=(T1-T0)/M;
KB=1.380658e12-23;
R=8.314;
T(1)=T0;
ekCO2=78.6;
ekAIR=195.2;
SQCO2=3.711;
SQAIR=3.941;
SQ12=0.5*(SQCO2+SQAIR);
MCO2=44.01;
MAIR=29;
MW=1/MCO2+1/MAIR;
a=(ekCO2*ekAIR)^-0.5;
for i=1:M;
P(i)=1;
T(i+1)=T(i)+DT;
TSCO2(i)=a*T(i+1);
end
A=1.06036;
B=0.1561;
C=0.193;
D=0.47635;
E=1.03587;
F=1.52996;
G=1.76474;
H=3.89411;
for i=1:M;
MDCO2(i)=A/(TSCO2(i)^B)+C/exp(D*TSCO2(i))+E/exp(F*TSCO2(i))+G/exp(H*TSCO2(i));
end
for i=1:M;
DCO2(i)=(1.86e-3*(T(i)^1.5)*(MW^0.5))/((P(i)*MDCO2(i)*(SQ12)^2));
end
DL=1;
TOR=1.8;
for i=1:M;
FYI(i)=0.45;
DeCO2(i)=(DCO2(i)*DL*FYI(i))/TOR;
end
plot(DeCO2,T(1:M),'^-')
title('Mass Diffusion in Pouros Meduim vs Tempreture')
ylabel('Tempreture in Kelvin')
hold on
%CH4 SECTION
ekCH4=78.6;
ekAIR=195.2;
SQCH4=3.711;
SQAIR=3.941;
SQ12=0.5*(SQCH4+SQAIR);
MCH4=16;
MAIR=29;
MW=1/MCH4+1/MAIR;
a=(ekCO2*ekAIR)^-0.5;
for i=1:M;
P(i)=1;
T(i+1)=T(i)+DT;
TSCH4(i)=a*T(i+1);
end
A=1.06036;
B=0.1561;
C=0.193;
D=0.47635;
E=1.03587;
F=1.52996;
G=1.76474;
H=3.89411;
for i=1:M;
MDCH4(i)=A/(TSCH4(i)^B)+C/exp(D*TSCH4(i))+E/exp(F*TSCH4(i))+G/exp(H*TSCH4(i));
end
for i=1:M;
DCH4(i)=(1.86e-3*(T(i)^1.5)*(MW^0.5))/((P(i)*MDCH4(i)*(SQ12)^2));
end
DL=1;
TOR=1.8;
for i=1:M;
FYI(i)=0.45;
DeCH4(i)=(DCH4(i)*DL*FYI(i))/TOR;
end
plot(DeCH4,T(1:M),'+-r')
%N2 SECTION
ekN2=71.4;
ekAIR=195.2;
SQN2=3.798;
SQAIR=3.941;
SQ12=0.5*(SQN2+SQAIR);
MN2=28;
MAIR=29;
MW=1/MN2+1/MAIR;
a=(ekN2*ekAIR)^-0.5;
for i=1:M;
P(i)=1;
T(i+1)=T(i)+DT;
TSN2(i)=a*T(i+1);
end
A=1.06036;
B=0.1561;
C=0.193;
D=0.47635;
E=1.03587;
F=1.52996;
G=1.76474;
H=3.89411;
for i=1:M;
MDN2(i)=A/(TSN2(i)^B)+C/exp(D*TSN2(i))+E/exp(F*TSN2(i))+G/exp(H*TSN2(i));
end
for i=1:M;
DN2(i)=(1.86e-3*(T(i)^1.5)*(MW^0.5))/((P(i)*MDN2(i)*(SQ12)^2));
end
DL=1;
TOR=1.8;
for i=1:M;
FYI(i)=0.45;
DeN2(i)=(DN2(i)*DL*FYI(i))/TOR;
end
plot(DeN2,T(1:M),'+-g')
%O2 SECTION
ekO2=106.7;
ekAIR=195.2;
SQO2=3.467;
SQAIR=3.941;
SQ12=0.5*(SQO2+SQAIR);
MO2=32;
MAIR=29;
MW=1/MO2+1/MAIR;
a=(ekO2*ekAIR)^-0.5;
for i=1:M;
P(i)=1;
T(i+1)=T(i)+DT;
TSO2(i)=a*T(i+1);
end
A=1.06036;
B=0.1561;
C=0.193;
D=0.47635;
E=1.03587;
F=1.52996;
G=1.76474;
H=3.89411;
for i=1:M;
MDO2(i)=A/(TSO2(i)^B)+C/exp(D*TSO2(i))+E/exp(F*TSO2(i))+G/exp(H*TSO2(i));
end
for i=1:M;
DO2(i)=(1.86e-3*(T(i)^1.5)*(MW^0.5))/((P(i)*MDO2(i)*(SQ12)^2));
end
DL=1;
TOR=1.8;
for i=1:M;
FYI(i)=0.45;
DeO2(i)=(DO2(i)*DL*FYI(i))/TOR;
end
plot(DeO2,T(1:M),'--r')
%Ar SECTION
ekAr=93.3;
ekAIR=195.2;
SQAr=3.542;
SQAIR=3.941;
SQ12=0.5*(SQAr+SQAIR);
MAr=39.948;
MAIR=29;
MW=1/MAr+1/MAIR;
a=(ekAr*ekAIR)^-0.5;
for i=1:M;
P(i)=1;
T(i+1)=T(i)+DT;
TSAr(i)=a*T(i+1);
end
A=1.06036;
B=0.1561;
C=0.193;
D=0.47635;
E=1.03587;
F=1.52996;
G=1.76474;
H=3.89411;
for i=1:M;
MDAr(i)=A/(TSAr(i)^B)+C/exp(D*TSAr(i))+E/exp(F*TSAr(i))+G/exp(H*TSAr(i));
end
for i=1:M;
DAr(i)=(1.86e-3*(T(i)^1.5)*(MW^0.5))/((P(i)*MDAr(i)*(SQ12)^2));
end
DL=1;
TOR=1.8;
for i=1:M;
FYI(i)=0.45;
DeAr(i)=(DAr(i)*DL*FYI(i))/TOR;
end
plot(DeAr,T(1:M),'*-r')
%Ne SECTION
ekNe=32.8;
ekAIR=195.2;
SQNe=2.820;
SQAIR=3.941;
SQ12=0.5*(SQCH4+SQAIR);
MNe=20.1797;
MAIR=29;
MW=1/MNe+1/MAIR;
a=(ekNe*ekAIR)^-0.5;
for i=1:M;
P(i)=1;
T(i+1)=T(i)+DT;
TSNe(i)=a*T(i+1);
end
A=1.06036;
B=0.1561;
C=0.193;
D=0.47635;
E=1.03587;
F=1.52996;
G=1.76474;
H=3.89411;
for i=1:M;
MDNe(i)=A/(TSNe(i)^B)+C/exp(D*TSNe(i))+E/exp(F*TSNe(i))+G/exp(H*TSNe(i));
end
for i=1:M;
DNe(i)=(1.86e-3*(T(i)^1.5)*(MW^0.5))/((P(i)*MDNe(i)*(SQ12)^2));
end
DL=1;
TOR=1.8;
for i=1:M;
FYI(i)=0.45;
DeNe(i)=(DNe(i)*DL*FYI(i))/TOR;
end
plot(DeNe,T(1:M),'^-k')
%He SECTION
ekHe=10.2;
ekAIR=195.2;
SQHe=2.551;
SQAIR=3.941;
SQ12=0.5*(SQHe+SQAIR);
MHe=4.002602;
MAIR=29;
MW=1/MHe+1/MAIR;
a=(ekHe*ekAIR)^-0.5;
for i=1:M;
P(i)=1;
T(i+1)=T(i)+DT;
TSHe(i)=a*T(i+1);
end
A=1.06036;
B=0.1561;
C=0.193;
D=0.47635;
E=1.03587;
F=1.52996;
G=1.76474;
H=3.89411;
for i=1:M;
MDHe(i)=A/(TSHe(i)^B)+C/exp(D*TSHe(i))+E/exp(F*TSHe(i))+G/exp(H*TSHe(i));
end
for i=1:M;
DHe(i)=(1.86e-3*(T(i)^1.5)*(MW^0.5))/((P(i)*MDHe(i)*(SQ12)^2));
end
DL=1;
TOR=1.8;
for i=1:M;
FYI(i)=0.45;
DeHe(i)=(DHe(i)*DL*FYI(i))/TOR;
end
plot(DeHe,T(1:M),'^-r')
%H2 SECTION
ekH2=59.7;
ekAIR=195.2;
SQH2=2.827;
SQAIR=3.941;
SQ12=0.5*(SQH2+SQAIR);
MH2=1.00794;
MAIR=29;
MW=1/MH2+1/MAIR;
a=(ekH2*ekAIR)^-0.5;
for i=1:M;
P(i)=1;
T(i+1)=T(i)+DT;
TSH2(i)=a*T(i+1);
end
A=1.06036;
B=0.1561;
C=0.193;
D=0.47635;
E=1.03587;
F=1.52996;
G=1.76474;
H=3.89411;
for i=1:M;
MDH2(i)=A/(TSH2(i)^B)+C/exp(D*TSH2(i))+E/exp(F*TSH2(i))+G/exp(H*TSH2(i));
end
for i=1:M;
DH2(i)=(1.86e-3*(T(i)^1.5)*(MW^0.5))/((P(i)*MDH2(i)*(SQ12)^2));
end
DL=1;
TOR=1.8;
for i=1:M;
FYI(i)=0.45;
DeH2(i)=(DH2(i)*DL*FYI(i))/TOR;
end
plot(DeH2,T(1:M),'-b')
%Xe SECTION
ekXe=229;
ekAIR=195.2;
SQXe=4.055;
SQAIR=3.941;
SQ12=0.5*(SQXe+SQAIR);
MXe=131.293;
MAIR=29;
MW=1/MXe+1/MAIR;
a=(ekXe*ekAIR)^-0.5;
for i=1:M;
P(i)=1;
T(i+1)=T(i)+DT;
TSXe(i)=a*T(i+1);
end
A=1.06036;
B=0.1561;
C=0.193;
D=0.47635;
E=1.03587;
F=1.52996;
G=1.76474;
H=3.89411;
for i=1:M;
MDXe(i)=A/(TSXe(i)^B)+C/exp(D*TSXe(i))+E/exp(F*TSXe(i))+G/exp(H*TSXe(i));
end
for i=1:M;
DXe(i)=(1.86e-3*(T(i)^1.5)*(MW^0.5))/((P(i)*MDXe(i)*(SQ12)^2));
end
DL=1;
TOR=1.8;
for i=1:M;
FYI(i)=0.45;
DeXe(i)=(DXe(i)*DL*FYI(i))/TOR;
end
plot(DeXe,T(1:M),'.k')
%H2O SECTION
ekH2O=809.1;
ekAIR=195.2;
SQH2O=2.641;
SQAIR=3.941;
SQ12=0.5*(SQH2O+SQAIR);
MH2O=18.01528;
MAIR=29;
MW=1/MH2O+1/MAIR;
a=(ekH2O*ekAIR)^-0.5;
for i=1:M;
P(i)=1;
T(i+1)=T(i)+DT;
TSH2O(i)=a*T(i+1);
end
A=1.06036;
B=0.1561;
C=0.193;
D=0.47635;
E=1.03587;
F=1.52996;
G=1.76474;
H=3.89411;
for i=1:M;
MDH2O(i)=A/(TSH2O(i)^B)+C/exp(D*TSH2O(i))+E/exp(F*TSH2O(i))+G/exp(H*TSH2O(i));
end
for i=1:M;
DH2O(i)=(1.86e-3*(T(i)^1.5)*(MW^0.5))/((P(i)*MDH2O(i)*(SQ12)^2));
end
DL=1;
TOR=1.8;
for i=1:M;
FYI(i)=0.45;
DeH2O(i)=(DH2O(i)*DL*FYI(i))/TOR;
end
plot(DeH2O,T(1:M),'--')
%N2O SECTION
ekN2O=232.4;
ekAIR=195.2;
SQN2O=3.828;
SQAIR=3.941;
SQ12=0.5*(SQN2O+SQAIR);
MN2O=44.0128;
MAIR=29;
MW=1/MN2O+1/MAIR;
a=(ekN2O*ekAIR)^-0.5;
for i=1:M;
P(i)=1;
T(i+1)=T(i)+DT;
TSN2O(i)=a*T(i+1);
end
A=1.06036;
B=0.1561;
C=0.193;
D=0.47635;
E=1.03587;
F=1.52996;
G=1.76474;
H=3.89411;
for i=1:M;
MDN2O(i)=A/(TSN2O(i)^B)+C/exp(D*TSN2O(i))+E/exp(F*TSN2O(i))+G/exp(H*TSN2O(i));
end
for i=1:M;
DN2O(i)=(1.86e-3*(T(i)^1.5)*(MW^0.5))/((P(i)*MDN2O(i)*(SQ12)^2));
end
DL=1;
TOR=1.8;
for i=1:M;
FYI(i)=0.45;
DeN2O(i)=(DN2O(i)*DL*FYI(i))/TOR;
end
plot(DeN2O,T(1:M),'c')
%CO SECTION
ekCO=91.7;
ekAIR=195.2;
SQCO=3.690;
SQAIR=3.941;
SQ12=0.5*(SQCO+SQAIR);
MCO=28.0101;
MAIR=29;
MW=1/MCO+1/MAIR;
a=(ekCO*ekAIR)^-0.5;
for i=1:M;
P(i)=1;
T(i+1)=T(i)+DT;
TSCO(i)=a*T(i+1);
end
A=1.06036;
B=0.1561;
C=0.193;
D=0.47635;
E=1.03587;
F=1.52996;
G=1.76474;
H=3.89411;
for i=1:M;
MDCO(i)=A/(TSCO(i)^B)+C/exp(D*TSCO(i))+E/exp(F*TSCO(i))+G/exp(H*TSCO(i));
end
for i=1:M;
DCO(i)=(1.86e-3*(T(i)^1.5)*(MW^0.5))/((P(i)*MDCO(i)*(SQ12)^2));
end
DL=1;
TOR=1.8;
for i=1:M;
FYI(i)=0.45;
DeCO(i)=(DCO(i)*DL*FYI(i))/TOR;
end
plot(DeCO,T(1:M),'-+')
%SO2 SECTION
ekSO2=335.2;
ekAIR=195.2;
SQSO2=4.112;
SQAIR=3.941;
SQ12=0.5*(SQSO2+SQAIR);
MSO2=64.0638;
MAIR=29;
MW=1/MSO2+1/MAIR;
a=(ekSO2*ekAIR)^-0.5;
for i=1:M;
P(i)=1;
T(i+1)=T(i)+DT;
TSSO2(i)=a*T(i+1);
end
A=1.06036;
B=0.1561;
C=0.193;
D=0.47635;
E=1.03587;
F=1.52996;
G=1.76474;
H=3.89411;
for i=1:M;
MDSO2(i)=A/(TSSO2(i)^B)+C/exp(D*TSSO2(i))+E/exp(F*TSSO2(i))+G/exp(H*TSSO2(i));
end
for i=1:M;
DSO2(i)=(1.86e-3*(T(i)^1.5)*(MW^0.5))/((P(i)*MDSO2(i)*(SQ12)^2));
end
DL=1;
TOR=1.8;
for i=1:M;
FYI(i)=0.45;
DeSO2(i)=(DSO2(i)*DL*FYI(i))/TOR;
end
plot(DeSO2,T(1:M),'+g')
%O3 SECTION
% ekO3=78.6;
% ekAIR=195.2;
% SQO3=3.711;
% SQAIR=3.941;
% SQ12=0.5*(SQO3+SQAIR);
% MO3=47.9982;
% MAIR=29;
% MW=1/MO3+1/MAIR;
%
% a=(ekO3*ekAIR)^-0.5;
%
% for i=1:M;
% P(i)=1;
% T(i+1)=T(i)+DT;
% TSO3(i)=a*T(i+1);
% end
%
% A=1.06036;
% B=0.1561;
% C=0.193;
% D=0.47635;
% E=1.03587;
% F=1.52996;
% G=1.76474;
% H=3.89411;
% for i=1:M;
% MDO3(i)=A/(TSO3(i)^B)+C/exp(D*TSO3(i))+E/exp(F*TSO3(i))+G/exp(H*TSO3(i));
% end
%
%
% for i=1:M;
% DO3(i)=(1.86e-3*(T(i)^1.5)*(MW^0.5))/((P(i)*MDO3(i)*(SQ12)^2));
% end
%
% DL=1;
% TOR=1.8;
%
% for i=1:M;
% FYI(i)=0.45;
% DeO3(i)=(DO3(i)*DL*FYI(i))/TOR;
% end
%
% plot(DeO3,T(1:M),'+-k')
% ylabel('Tempreture in Kelvin')
%NO2 SECTION
% ekNO2=78.6;
% ekAIR=195.2;
% SQNO2=3.711;
% SQAIR=3.941;
% SQ12=0.5*(SQNO2+SQAIR);
% MNO2=46.0055;
% MAIR=29;
% MW=1/MNO2+1/MAIR;
%
% a=(ekNO2*ekAIR)^-0.5;
%
% for i=1:M;
% P(i)=1;
% T(i+1)=T(i)+DT;
% TSNO2(i)=a*T(i+1);
% end
%
% A=1.06036;
% B=0.1561;
% C=0.193;
% D=0.47635;
% E=1.03587;
% F=1.52996;
% G=1.76474;
% H=3.89411;
% for i=1:M;
% MDNO2(i)=A/(TSNO2(i)^B)+C/exp(D*TSNO2(i))+E/exp(F*TSNO2(i))+G/exp(H*TSNO2(i));
% end
%
%
% for i=1:M;
% DNO2(i)=(1.86e-3*(T(i)^1.5)*(MW^0.5))/((P(i)*MDNO2(i)*(SQ12)^2));
% end
%
% DL=1;
% TOR=1.8;
%
% for i=1:M;
% FYI(i)=0.45;
% DeNO2(i)=(DNO2(i)*DL*FYI(i))/TOR;
% end
%
% plot(DeNO2,T(1:M),'+-b')
grid on
legend('Mass Diffusion in Pouros Meduim CO_2','Mass Diffusion in Pouros Meduim CH_4','Mass Diffusion in Pouros Meduim N_2','Mass Diffusion in Pouros Meduim O_2','Mass Diffusion in Pouros Meduim Ar','Mass Diffusion in Pouros Meduim Ne','Mass Diffusion in Pouros Meduim He','Mass Diffusion in Pouros Meduim H2','Mass Diffusion in Pouros Meduim Xe','Mass Diffusion in Pouros Meduim H2O','Mass Diffusion in Pouros Meduim N2O','Mass Diffusion in Pouros Meduim CO','Mass Diffusion in Pouros Meduim SO2');%,'Mass Diffusion in Pouros Meduim O3')%,'Mass Diffusion in Pouros Meduim NO2')
xlabel('Mass Diffusion in Pouros Meduim cm^2/sec')
ylabel('Tempreture in Kelvin')
clear
M=100;
t0=0;
t1=30;
T0=273+t0;
T1=300+t1;
DT=(T1-T0)/M;
KB=1.380658e12-23;
R=8.314;
T(1)=T0;
ekCO2=78.6;
ekAIR=195.2;
SQCO2=3.711;
SQAIR=3.941;
SQ12=0.5*(SQCO2+SQAIR);
MCO2=44.01;
MAIR=29;
MW=1/MCO2+1/MAIR;
a=(ekCO2*ekAIR)^-0.5;
for i=1:M;
P(i)=1;
T(i+1)=T(i)+DT;
TSCO2(i)=a*T(i+1);
end
A=1.06036;
B=0.1561;
C=0.193;
D=0.47635;
E=1.03587;
F=1.52996;
G=1.76474;
H=3.89411;
for i=1:M;
MDCO2(i)=A/(TSCO2(i)^B)+C/exp(D*TSCO2(i))+E/exp(F*TSCO2(i))+G/exp(H*TSCO2(i));
end
for i=1:M;
DCO2(i)=(1.86e-3*(T(i)^1.5)*(MW^0.5))/((P(i)*MDCO2(i)*(SQ12)^2));
end
DL=1;
TOR=1.8;
for i=1:M;
FYI(i)=0.45;
DeCO2(i)=(DCO2(i)*DL*FYI(i))/TOR;
end
plot(DeCO2,T(1:M),'^-')
title('Mass Diffusion in Pouros Meduim vs Tempreture')
ylabel('Tempreture in Kelvin')
hold on
%CH4 SECTION
ekCH4=78.6;
ekAIR=195.2;
SQCH4=3.711;
SQAIR=3.941;
SQ12=0.5*(SQCH4+SQAIR);
MCH4=16;
MAIR=29;
MW=1/MCH4+1/MAIR;
a=(ekCO2*ekAIR)^-0.5;
for i=1:M;
P(i)=1;
T(i+1)=T(i)+DT;
TSCH4(i)=a*T(i+1);
end
A=1.06036;
B=0.1561;
C=0.193;
D=0.47635;
E=1.03587;
F=1.52996;
G=1.76474;
H=3.89411;
for i=1:M;
MDCH4(i)=A/(TSCH4(i)^B)+C/exp(D*TSCH4(i))+E/exp(F*TSCH4(i))+G/exp(H*TSCH4(i));
end
for i=1:M;
DCH4(i)=(1.86e-3*(T(i)^1.5)*(MW^0.5))/((P(i)*MDCH4(i)*(SQ12)^2));
end
DL=1;
TOR=1.8;
for i=1:M;
FYI(i)=0.45;
DeCH4(i)=(DCH4(i)*DL*FYI(i))/TOR;
end
plot(DeCH4,T(1:M),'+-r')
%N2 SECTION
ekN2=71.4;
ekAIR=195.2;
SQN2=3.798;
SQAIR=3.941;
SQ12=0.5*(SQN2+SQAIR);
MN2=28;
MAIR=29;
MW=1/MN2+1/MAIR;
a=(ekN2*ekAIR)^-0.5;
for i=1:M;
P(i)=1;
T(i+1)=T(i)+DT;
TSN2(i)=a*T(i+1);
end
A=1.06036;
B=0.1561;
C=0.193;
D=0.47635;
E=1.03587;
F=1.52996;
G=1.76474;
H=3.89411;
for i=1:M;
MDN2(i)=A/(TSN2(i)^B)+C/exp(D*TSN2(i))+E/exp(F*TSN2(i))+G/exp(H*TSN2(i));
end
for i=1:M;
DN2(i)=(1.86e-3*(T(i)^1.5)*(MW^0.5))/((P(i)*MDN2(i)*(SQ12)^2));
end
DL=1;
TOR=1.8;
for i=1:M;
FYI(i)=0.45;
DeN2(i)=(DN2(i)*DL*FYI(i))/TOR;
end
plot(DeN2,T(1:M),'+-g')
%O2 SECTION
ekO2=106.7;
ekAIR=195.2;
SQO2=3.467;
SQAIR=3.941;
SQ12=0.5*(SQO2+SQAIR);
MO2=32;
MAIR=29;
MW=1/MO2+1/MAIR;
a=(ekO2*ekAIR)^-0.5;
for i=1:M;
P(i)=1;
T(i+1)=T(i)+DT;
TSO2(i)=a*T(i+1);
end
A=1.06036;
B=0.1561;
C=0.193;
D=0.47635;
E=1.03587;
F=1.52996;
G=1.76474;
H=3.89411;
for i=1:M;
MDO2(i)=A/(TSO2(i)^B)+C/exp(D*TSO2(i))+E/exp(F*TSO2(i))+G/exp(H*TSO2(i));
end
for i=1:M;
DO2(i)=(1.86e-3*(T(i)^1.5)*(MW^0.5))/((P(i)*MDO2(i)*(SQ12)^2));
end
DL=1;
TOR=1.8;
for i=1:M;
FYI(i)=0.45;
DeO2(i)=(DO2(i)*DL*FYI(i))/TOR;
end
plot(DeO2,T(1:M),'--r')
%Ar SECTION
ekAr=93.3;
ekAIR=195.2;
SQAr=3.542;
SQAIR=3.941;
SQ12=0.5*(SQAr+SQAIR);
MAr=39.948;
MAIR=29;
MW=1/MAr+1/MAIR;
a=(ekAr*ekAIR)^-0.5;
for i=1:M;
P(i)=1;
T(i+1)=T(i)+DT;
TSAr(i)=a*T(i+1);
end
A=1.06036;
B=0.1561;
C=0.193;
D=0.47635;
E=1.03587;
F=1.52996;
G=1.76474;
H=3.89411;
for i=1:M;
MDAr(i)=A/(TSAr(i)^B)+C/exp(D*TSAr(i))+E/exp(F*TSAr(i))+G/exp(H*TSAr(i));
end
for i=1:M;
DAr(i)=(1.86e-3*(T(i)^1.5)*(MW^0.5))/((P(i)*MDAr(i)*(SQ12)^2));
end
DL=1;
TOR=1.8;
for i=1:M;
FYI(i)=0.45;
DeAr(i)=(DAr(i)*DL*FYI(i))/TOR;
end
plot(DeAr,T(1:M),'*-r')
%Ne SECTION
ekNe=32.8;
ekAIR=195.2;
SQNe=2.820;
SQAIR=3.941;
SQ12=0.5*(SQCH4+SQAIR);
MNe=20.1797;
MAIR=29;
MW=1/MNe+1/MAIR;
a=(ekNe*ekAIR)^-0.5;
for i=1:M;
P(i)=1;
T(i+1)=T(i)+DT;
TSNe(i)=a*T(i+1);
end
A=1.06036;
B=0.1561;
C=0.193;
D=0.47635;
E=1.03587;
F=1.52996;
G=1.76474;
H=3.89411;
for i=1:M;
MDNe(i)=A/(TSNe(i)^B)+C/exp(D*TSNe(i))+E/exp(F*TSNe(i))+G/exp(H*TSNe(i));
end
for i=1:M;
DNe(i)=(1.86e-3*(T(i)^1.5)*(MW^0.5))/((P(i)*MDNe(i)*(SQ12)^2));
end
DL=1;
TOR=1.8;
for i=1:M;
FYI(i)=0.45;
DeNe(i)=(DNe(i)*DL*FYI(i))/TOR;
end
plot(DeNe,T(1:M),'^-k')
%He SECTION
ekHe=10.2;
ekAIR=195.2;
SQHe=2.551;
SQAIR=3.941;
SQ12=0.5*(SQHe+SQAIR);
MHe=4.002602;
MAIR=29;
MW=1/MHe+1/MAIR;
a=(ekHe*ekAIR)^-0.5;
for i=1:M;
P(i)=1;
T(i+1)=T(i)+DT;
TSHe(i)=a*T(i+1);
end
A=1.06036;
B=0.1561;
C=0.193;
D=0.47635;
E=1.03587;
F=1.52996;
G=1.76474;
H=3.89411;
for i=1:M;
MDHe(i)=A/(TSHe(i)^B)+C/exp(D*TSHe(i))+E/exp(F*TSHe(i))+G/exp(H*TSHe(i));
end
for i=1:M;
DHe(i)=(1.86e-3*(T(i)^1.5)*(MW^0.5))/((P(i)*MDHe(i)*(SQ12)^2));
end
DL=1;
TOR=1.8;
for i=1:M;
FYI(i)=0.45;
DeHe(i)=(DHe(i)*DL*FYI(i))/TOR;
end
plot(DeHe,T(1:M),'^-r')
%H2 SECTION
ekH2=59.7;
ekAIR=195.2;
SQH2=2.827;
SQAIR=3.941;
SQ12=0.5*(SQH2+SQAIR);
MH2=1.00794;
MAIR=29;
MW=1/MH2+1/MAIR;
a=(ekH2*ekAIR)^-0.5;
for i=1:M;
P(i)=1;
T(i+1)=T(i)+DT;
TSH2(i)=a*T(i+1);
end
A=1.06036;
B=0.1561;
C=0.193;
D=0.47635;
E=1.03587;
F=1.52996;
G=1.76474;
H=3.89411;
for i=1:M;
MDH2(i)=A/(TSH2(i)^B)+C/exp(D*TSH2(i))+E/exp(F*TSH2(i))+G/exp(H*TSH2(i));
end
for i=1:M;
DH2(i)=(1.86e-3*(T(i)^1.5)*(MW^0.5))/((P(i)*MDH2(i)*(SQ12)^2));
end
DL=1;
TOR=1.8;
for i=1:M;
FYI(i)=0.45;
DeH2(i)=(DH2(i)*DL*FYI(i))/TOR;
end
plot(DeH2,T(1:M),'-b')
%Xe SECTION
ekXe=229;
ekAIR=195.2;
SQXe=4.055;
SQAIR=3.941;
SQ12=0.5*(SQXe+SQAIR);
MXe=131.293;
MAIR=29;
MW=1/MXe+1/MAIR;
a=(ekXe*ekAIR)^-0.5;
for i=1:M;
P(i)=1;
T(i+1)=T(i)+DT;
TSXe(i)=a*T(i+1);
end
A=1.06036;
B=0.1561;
C=0.193;
D=0.47635;
E=1.03587;
F=1.52996;
G=1.76474;
H=3.89411;
for i=1:M;
MDXe(i)=A/(TSXe(i)^B)+C/exp(D*TSXe(i))+E/exp(F*TSXe(i))+G/exp(H*TSXe(i));
end
for i=1:M;
DXe(i)=(1.86e-3*(T(i)^1.5)*(MW^0.5))/((P(i)*MDXe(i)*(SQ12)^2));
end
DL=1;
TOR=1.8;
for i=1:M;
FYI(i)=0.45;
DeXe(i)=(DXe(i)*DL*FYI(i))/TOR;
end
plot(DeXe,T(1:M),'.k')
%H2O SECTION
ekH2O=809.1;
ekAIR=195.2;
SQH2O=2.641;
SQAIR=3.941;
SQ12=0.5*(SQH2O+SQAIR);
MH2O=18.01528;
MAIR=29;
MW=1/MH2O+1/MAIR;
a=(ekH2O*ekAIR)^-0.5;
for i=1:M;
P(i)=1;
T(i+1)=T(i)+DT;
TSH2O(i)=a*T(i+1);
end
A=1.06036;
B=0.1561;
C=0.193;
D=0.47635;
E=1.03587;
F=1.52996;
G=1.76474;
H=3.89411;
for i=1:M;
MDH2O(i)=A/(TSH2O(i)^B)+C/exp(D*TSH2O(i))+E/exp(F*TSH2O(i))+G/exp(H*TSH2O(i));
end
for i=1:M;
DH2O(i)=(1.86e-3*(T(i)^1.5)*(MW^0.5))/((P(i)*MDH2O(i)*(SQ12)^2));
end
DL=1;
TOR=1.8;
for i=1:M;
FYI(i)=0.45;
DeH2O(i)=(DH2O(i)*DL*FYI(i))/TOR;
end
plot(DeH2O,T(1:M),'--')
%N2O SECTION
ekN2O=232.4;
ekAIR=195.2;
SQN2O=3.828;
SQAIR=3.941;
SQ12=0.5*(SQN2O+SQAIR);
MN2O=44.0128;
MAIR=29;
MW=1/MN2O+1/MAIR;
a=(ekN2O*ekAIR)^-0.5;
for i=1:M;
P(i)=1;
T(i+1)=T(i)+DT;
TSN2O(i)=a*T(i+1);
end
A=1.06036;
B=0.1561;
C=0.193;
D=0.47635;
E=1.03587;
F=1.52996;
G=1.76474;
H=3.89411;
for i=1:M;
MDN2O(i)=A/(TSN2O(i)^B)+C/exp(D*TSN2O(i))+E/exp(F*TSN2O(i))+G/exp(H*TSN2O(i));
end
for i=1:M;
DN2O(i)=(1.86e-3*(T(i)^1.5)*(MW^0.5))/((P(i)*MDN2O(i)*(SQ12)^2));
end
DL=1;
TOR=1.8;
for i=1:M;
FYI(i)=0.45;
DeN2O(i)=(DN2O(i)*DL*FYI(i))/TOR;
end
plot(DeN2O,T(1:M),'c')
%CO SECTION
ekCO=91.7;
ekAIR=195.2;
SQCO=3.690;
SQAIR=3.941;
SQ12=0.5*(SQCO+SQAIR);
MCO=28.0101;
MAIR=29;
MW=1/MCO+1/MAIR;
a=(ekCO*ekAIR)^-0.5;
for i=1:M;
P(i)=1;
T(i+1)=T(i)+DT;
TSCO(i)=a*T(i+1);
end
A=1.06036;
B=0.1561;
C=0.193;
D=0.47635;
E=1.03587;
F=1.52996;
G=1.76474;
H=3.89411;
for i=1:M;
MDCO(i)=A/(TSCO(i)^B)+C/exp(D*TSCO(i))+E/exp(F*TSCO(i))+G/exp(H*TSCO(i));
end
for i=1:M;
DCO(i)=(1.86e-3*(T(i)^1.5)*(MW^0.5))/((P(i)*MDCO(i)*(SQ12)^2));
end
DL=1;
TOR=1.8;
for i=1:M;
FYI(i)=0.45;
DeCO(i)=(DCO(i)*DL*FYI(i))/TOR;
end
plot(DeCO,T(1:M),'-+')
%SO2 SECTION
ekSO2=335.2;
ekAIR=195.2;
SQSO2=4.112;
SQAIR=3.941;
SQ12=0.5*(SQSO2+SQAIR);
MSO2=64.0638;
MAIR=29;
MW=1/MSO2+1/MAIR;
a=(ekSO2*ekAIR)^-0.5;
for i=1:M;
P(i)=1;
T(i+1)=T(i)+DT;
TSSO2(i)=a*T(i+1);
end
A=1.06036;
B=0.1561;
C=0.193;
D=0.47635;
E=1.03587;
F=1.52996;
G=1.76474;
H=3.89411;
for i=1:M;
MDSO2(i)=A/(TSSO2(i)^B)+C/exp(D*TSSO2(i))+E/exp(F*TSSO2(i))+G/exp(H*TSSO2(i));
end
for i=1:M;
DSO2(i)=(1.86e-3*(T(i)^1.5)*(MW^0.5))/((P(i)*MDSO2(i)*(SQ12)^2));
end
DL=1;
TOR=1.8;
for i=1:M;
FYI(i)=0.45;
DeSO2(i)=(DSO2(i)*DL*FYI(i))/TOR;
end
plot(DeSO2,T(1:M),'+g')
%O3 SECTION
% ekO3=78.6;
% ekAIR=195.2;
% SQO3=3.711;
% SQAIR=3.941;
% SQ12=0.5*(SQO3+SQAIR);
% MO3=47.9982;
% MAIR=29;
% MW=1/MO3+1/MAIR;
%
% a=(ekO3*ekAIR)^-0.5;
%
% for i=1:M;
% P(i)=1;
% T(i+1)=T(i)+DT;
% TSO3(i)=a*T(i+1);
% end
%
% A=1.06036;
% B=0.1561;
% C=0.193;
% D=0.47635;
% E=1.03587;
% F=1.52996;
% G=1.76474;
% H=3.89411;
% for i=1:M;
% MDO3(i)=A/(TSO3(i)^B)+C/exp(D*TSO3(i))+E/exp(F*TSO3(i))+G/exp(H*TSO3(i));
% end
%
%
% for i=1:M;
% DO3(i)=(1.86e-3*(T(i)^1.5)*(MW^0.5))/((P(i)*MDO3(i)*(SQ12)^2));
% end
%
% DL=1;
% TOR=1.8;
%
% for i=1:M;
% FYI(i)=0.45;
% DeO3(i)=(DO3(i)*DL*FYI(i))/TOR;
% end
%
% plot(DeO3,T(1:M),'+-k')
% ylabel('Tempreture in Kelvin')
%NO2 SECTION
% ekNO2=78.6;
% ekAIR=195.2;
% SQNO2=3.711;
% SQAIR=3.941;
% SQ12=0.5*(SQNO2+SQAIR);
% MNO2=46.0055;
% MAIR=29;
% MW=1/MNO2+1/MAIR;
%
% a=(ekNO2*ekAIR)^-0.5;
%
% for i=1:M;
% P(i)=1;
% T(i+1)=T(i)+DT;
% TSNO2(i)=a*T(i+1);
% end
%
% A=1.06036;
% B=0.1561;
% C=0.193;
% D=0.47635;
% E=1.03587;
% F=1.52996;
% G=1.76474;
% H=3.89411;
% for i=1:M;
% MDNO2(i)=A/(TSNO2(i)^B)+C/exp(D*TSNO2(i))+E/exp(F*TSNO2(i))+G/exp(H*TSNO2(i));
% end
%
%
% for i=1:M;
% DNO2(i)=(1.86e-3*(T(i)^1.5)*(MW^0.5))/((P(i)*MDNO2(i)*(SQ12)^2));
% end
%
% DL=1;
% TOR=1.8;
%
% for i=1:M;
% FYI(i)=0.45;
% DeNO2(i)=(DNO2(i)*DL*FYI(i))/TOR;
% end
%
% plot(DeNO2,T(1:M),'+-b')
grid on
legend('Mass Diffusion in Pouros Meduim CO_2','Mass Diffusion in Pouros Meduim CH_4','Mass Diffusion in Pouros Meduim N_2','Mass Diffusion in Pouros Meduim O_2','Mass Diffusion in Pouros Meduim Ar','Mass Diffusion in Pouros Meduim Ne','Mass Diffusion in Pouros Meduim He','Mass Diffusion in Pouros Meduim H2','Mass Diffusion in Pouros Meduim Xe','Mass Diffusion in Pouros Meduim H2O','Mass Diffusion in Pouros Meduim N2O','Mass Diffusion in Pouros Meduim CO','Mass Diffusion in Pouros Meduim SO2');%,'Mass Diffusion in Pouros Meduim O3')%,'Mass Diffusion in Pouros Meduim NO2')
xlabel('Mass Diffusion in Pouros Meduim cm^2/sec')
ylabel('Tempreture in Kelvin')
Example Methane Diffusion
clc
clear
M=100;
t0=0;
t1=30;
T0=273+t0;
T1=300+t1;
aAIR=78.6;
aCH4=148.6;
a=aAIR*aCH4;
DT=(T1-T0)/M;
T(1)=T0;
for i=1:M;
P(i)=1+0.01*randn(1,1);
T(i+1)=T(i)+DT;
TS(i)=a*T(i+1);
end
A=1.06036;
B=0.1561;
C=0.193;
D=0.47635;
E=1.03587;
F=1.52996;
G=1.76474;
H=3.89411;
for i=1:M;
MD(i)=A/(TS(i)^B)+C/exp(D*TS(i))+E/exp(F*TS(i))+G/exp(H*TS(i));
end
SQ1=3.758;
SQ2=3.941;
SQ12=0.5*(SQ1+SQ2);
M1=16;
M2=29;
MW=2/(1/M1+1/M2);
KB=1.380658e-23;
R=8.314;
for i=1:M;
D(i)=(3/16)*((4*pi*KB*T(i))/MW)^0.5/((P(i)/(R*T(i)))*MD(i)*(SQ12)^2);
end
plot(D,T(1:M),'.-');
hold on
DL=1;
TOR=1.8;
for i=1:M;
FYI(i)=0.45;%+0.1*randn(1,1);
De(i)=(D(i)*DL*FYI(i))/TOR;
end
plot(De,T(1:M),'+-')
legend('Mass Diffusion','Mass Diffusion in Pouros Meduim')
title('Mass Diffusion vs Tempreture')
grid on
xlabel('CH_4 Mass Diffusion in Air (m^2/s)')
ylabel('Tempreture in Kelvin')
clear
M=100;
t0=0;
t1=30;
T0=273+t0;
T1=300+t1;
aAIR=78.6;
aCH4=148.6;
a=aAIR*aCH4;
DT=(T1-T0)/M;
T(1)=T0;
for i=1:M;
P(i)=1+0.01*randn(1,1);
T(i+1)=T(i)+DT;
TS(i)=a*T(i+1);
end
A=1.06036;
B=0.1561;
C=0.193;
D=0.47635;
E=1.03587;
F=1.52996;
G=1.76474;
H=3.89411;
for i=1:M;
MD(i)=A/(TS(i)^B)+C/exp(D*TS(i))+E/exp(F*TS(i))+G/exp(H*TS(i));
end
SQ1=3.758;
SQ2=3.941;
SQ12=0.5*(SQ1+SQ2);
M1=16;
M2=29;
MW=2/(1/M1+1/M2);
KB=1.380658e-23;
R=8.314;
for i=1:M;
D(i)=(3/16)*((4*pi*KB*T(i))/MW)^0.5/((P(i)/(R*T(i)))*MD(i)*(SQ12)^2);
end
plot(D,T(1:M),'.-');
hold on
DL=1;
TOR=1.8;
for i=1:M;
FYI(i)=0.45;%+0.1*randn(1,1);
De(i)=(D(i)*DL*FYI(i))/TOR;
end
plot(De,T(1:M),'+-')
legend('Mass Diffusion','Mass Diffusion in Pouros Meduim')
title('Mass Diffusion vs Tempreture')
grid on
xlabel('CH_4 Mass Diffusion in Air (m^2/s)')
ylabel('Tempreture in Kelvin')
Concentration of Gas in Soil With Relation to Height
clc
clear
Constructivity=1;
Porosity=0.45;
Tortisity=2;
DiffCo2=12.8E-6;
EDiffCO2=(Porosity*DiffCo2*Constructivity)/Tortisity;
AlphaCo2=0.4e-6;
M=100;
L=0.6;
dz=L/M;
Z(1)=0;
for i=1:M;
Z(i+1)=Z(i)+dz;
NmCO2(i)=(0.5*AlphaCo2/EDiffCO2)*((0.6)^2-Z(i)^2)
end
plot(NmCO2,Z(1:100),'-')
grid on
title('CO2 Moles VS Height')
xlabel('XCO2 g/m^3')
ylabel('Z Height in Meters')
legend('XCO2 %')
clear
Constructivity=1;
Porosity=0.45;
Tortisity=2;
DiffCo2=12.8E-6;
EDiffCO2=(Porosity*DiffCo2*Constructivity)/Tortisity;
AlphaCo2=0.4e-6;
M=100;
L=0.6;
dz=L/M;
Z(1)=0;
for i=1:M;
Z(i+1)=Z(i)+dz;
NmCO2(i)=(0.5*AlphaCo2/EDiffCO2)*((0.6)^2-Z(i)^2)
end
plot(NmCO2,Z(1:100),'-')
grid on
title('CO2 Moles VS Height')
xlabel('XCO2 g/m^3')
ylabel('Z Height in Meters')
legend('XCO2 %')
Diffusion for Different Gases
clc
clear
%
% i SPECIES XI MWI SEKMAI EPSLONI/KB
% 1 N2 0.78 28 71.4
% 2 O2 0.2095 16 106.7
% 3 Ar 0.0093 39.95 93.3
% 4 Ne 0.000018 20.18
% 5 He 0.000005 4.00 10.22
% 6 H2 0.0000005 2 59.7
% 7 Xe 0.000000009 131.29
% 8 H2O 0.000789181 20.16 809.1
% 9 CO2 0.00035 44.01 195.2
% 10 CH4 0.0000017 16 148.6
% 11 N2O 0.0000003 44.01
% 12 CO 0.000035 28 91.7
% 13 SO2 0.00000014 64.06 335.4
% 14 O3 0.00000012 48.00
% 15 NO2 0.00000005 28.01
% 16 AIR
%note that 14 and 15 have been approximated by me.
NAME(1,1:2)='N2';
NAME(2,1:2)='O2';
NAME(3,1:2)='Ar';
NAME(4,1:2)='Ne';
NAME(5,1:2)='He';
NAME(6,1:2)='H2';
NAME(7,1:2)='Xe';
NAME(8,1:3)='H2O';
NAME(9,1:3)='CO2';
NAME(10,1:3)='CH4';
NAME(11,1:3)='N2O';
NAME(12,1:2)='CO';
NAME(13,1:3)='SO2';
NAME(14,1:2)='O3';
NAME(15,1:3)='NO2';
X(1)=0.78;
X(2)=0.2095;
X(3)=0.0093;
X(4)=0.000018;
X(5)=0.000005;
X(6)=0.0000005;
X(7)=0.000000009;
X(8)=0.000789181;
X(9)=0.00035;
X(10)=0.0000017;
X(11)=0.0000003;
X(12)=0.000035;
X(13)=0.00000014;
X(14)=0.00000012;
X(15)=0.00000005;
%X(16)=
MW(1)=28;
MW(2)=16;
MW(3)=39.95;
MW(4)=20.18;
MW(5)=4;
MW(6)=2;
MW(7)=131.29;
MW(8)=20.16;
MW(9)=44.01;
MW(10)=16;
MW(11)=44.01;
MW(12)=28;
MW(13)=64.06;
MW(14)=48;
MW(15)=28.01;
SEKMA(1)=3.798;
SEKMA(2)=3.467;
SEKMA(3)=3.542;
SEKMA(4)=0.282;
SEKMA(5)=2.551;
SEKMA(6)=59.17;
SEKMA(7)=0.405;
SEKMA(8)=2.641;
SEKMA(9)=3.941;
SEKMA(10)=3.758;
SEKMA(11)=0.3828;
SEKMA(12)=3.69;
SEKMA(13)=4.112;
SEKMA(14)=0.3884;
SEKMA(15)=0.3686;
EP_KB(1)=71.4;
EP_KB(2)=106.7;
EP_KB(3)=93.3;
EP_KB(4)=32.8;
EP_KB(5)=10.22;
EP_KB(6)=59.7;
EP_KB(7)=229;
EP_KB(8)=809.1;
EP_KB(9)=195.2;
EP_KB(10)=148.6;
EP_KB(11)=232.4;
EP_KB(12)=91.7;
EP_KB(13)=335.4;
EP_KB(14)=106.7;
EP_KB(15)=162;
A=1.06036;
B=0.15610;
C=0.19300;
D=0.47635;
E=1.03587;
F=1.52996;
G=1.76474;
H=3.89411;
KB=1.380658*10^(-23);
for i=1:15;
for j=1:15;
if (i==j)
EP_BOTH(i,j)=0;
else
EP_BOTH(i,j)=(EP_KB(i)*EP_KB(j))/KB;
end
end
end
for i=1:15;
for j=1:15;
T(i,j)=600;
P(i,j)=101.325;
T_STAR(i,j)=KB*T(i,j)/EP_BOTH(i,j);
end
end
for i=1:15;
for j=1:15;
if (i==j)
SEKMA_BOTH(i,j)=0;
else
SEKMA_BOTH(i,j)=(SEKMA(i)+SEKMA(j))/2;
end
end
end
for i=1:15;
for j=1:15;
if (i==j)
MW_BOTH(i,j)=0;
else
MW_BOTH(i,j)=2/(1/MW(i)+1/MW(j));
end
end
end
for i=1:15;
for j=1:15;
if (i==j)
MU_D(i,j)=0;
else
MU_D(i,j)=A/((T_STAR(i,j))^(B))+B/exp(D*T_STAR(i,j))+E/exp(F*T_STAR(i,j))+G/exp(H*T_STAR(i,j));
end
end
end
for i=1:15;
for j=1:15;
if (i==j)
D(i,j)=0;
else
DD(i,j)=(0.0266*T(i,j)^1.5)/(P(i,j)*((MW_BOTH(i,j))^(0.5))*(SEKMA_BOTH(i,j))^(2)*MU_D(i,j));
end
end
end
for i=1:15;
for j=1:15;
if (i==j)
L(i,j)=0;
elseif (i<=14)&(j<=14)
L(i,j)=(X(j)*(MW(j)*X(j)+MW(i)*X(i)))/(DD(i,j)*MW(i))+(X(j+1)*(MW(j)*X(j)))/(MW(i)*DD(i,j+1));
elseif (i==15)&(j==15)
L(i,j)=X(j)*(MW(j)*X(j)+MW(i)*X(i))+X(j-1)*((MW(j)*X(j)))/(MW(i)*DD(i,j));
end
end
end
for i=1:15;
for j=1:15;
F(i,j)=0.01*randn(1,1);
end
end
for i=1:15;
for j=1:15;
delta_F(i,j)=F(i,j)-F(i,i);
end
end
for i=1:15;
z(i)=X(i)*MW(i);
end
MW_MIX=sum(z);
for i=1:15;
for j=1:15;
D(i,j)=((X(i)*MW_MIX)/MW(j))*delta_F(i,j);
end
end
DD
surf(D)
for i=1:15
for j=1:15;
if (D(i,j)==max(max(D)))
imax=i
jmax=j
NAME(imax,:)
NAME(jmax,:)
else
display('no');
end
end
end
clear
%
% i SPECIES XI MWI SEKMAI EPSLONI/KB
% 1 N2 0.78 28 71.4
% 2 O2 0.2095 16 106.7
% 3 Ar 0.0093 39.95 93.3
% 4 Ne 0.000018 20.18
% 5 He 0.000005 4.00 10.22
% 6 H2 0.0000005 2 59.7
% 7 Xe 0.000000009 131.29
% 8 H2O 0.000789181 20.16 809.1
% 9 CO2 0.00035 44.01 195.2
% 10 CH4 0.0000017 16 148.6
% 11 N2O 0.0000003 44.01
% 12 CO 0.000035 28 91.7
% 13 SO2 0.00000014 64.06 335.4
% 14 O3 0.00000012 48.00
% 15 NO2 0.00000005 28.01
% 16 AIR
%note that 14 and 15 have been approximated by me.
NAME(1,1:2)='N2';
NAME(2,1:2)='O2';
NAME(3,1:2)='Ar';
NAME(4,1:2)='Ne';
NAME(5,1:2)='He';
NAME(6,1:2)='H2';
NAME(7,1:2)='Xe';
NAME(8,1:3)='H2O';
NAME(9,1:3)='CO2';
NAME(10,1:3)='CH4';
NAME(11,1:3)='N2O';
NAME(12,1:2)='CO';
NAME(13,1:3)='SO2';
NAME(14,1:2)='O3';
NAME(15,1:3)='NO2';
X(1)=0.78;
X(2)=0.2095;
X(3)=0.0093;
X(4)=0.000018;
X(5)=0.000005;
X(6)=0.0000005;
X(7)=0.000000009;
X(8)=0.000789181;
X(9)=0.00035;
X(10)=0.0000017;
X(11)=0.0000003;
X(12)=0.000035;
X(13)=0.00000014;
X(14)=0.00000012;
X(15)=0.00000005;
%X(16)=
MW(1)=28;
MW(2)=16;
MW(3)=39.95;
MW(4)=20.18;
MW(5)=4;
MW(6)=2;
MW(7)=131.29;
MW(8)=20.16;
MW(9)=44.01;
MW(10)=16;
MW(11)=44.01;
MW(12)=28;
MW(13)=64.06;
MW(14)=48;
MW(15)=28.01;
SEKMA(1)=3.798;
SEKMA(2)=3.467;
SEKMA(3)=3.542;
SEKMA(4)=0.282;
SEKMA(5)=2.551;
SEKMA(6)=59.17;
SEKMA(7)=0.405;
SEKMA(8)=2.641;
SEKMA(9)=3.941;
SEKMA(10)=3.758;
SEKMA(11)=0.3828;
SEKMA(12)=3.69;
SEKMA(13)=4.112;
SEKMA(14)=0.3884;
SEKMA(15)=0.3686;
EP_KB(1)=71.4;
EP_KB(2)=106.7;
EP_KB(3)=93.3;
EP_KB(4)=32.8;
EP_KB(5)=10.22;
EP_KB(6)=59.7;
EP_KB(7)=229;
EP_KB(8)=809.1;
EP_KB(9)=195.2;
EP_KB(10)=148.6;
EP_KB(11)=232.4;
EP_KB(12)=91.7;
EP_KB(13)=335.4;
EP_KB(14)=106.7;
EP_KB(15)=162;
A=1.06036;
B=0.15610;
C=0.19300;
D=0.47635;
E=1.03587;
F=1.52996;
G=1.76474;
H=3.89411;
KB=1.380658*10^(-23);
for i=1:15;
for j=1:15;
if (i==j)
EP_BOTH(i,j)=0;
else
EP_BOTH(i,j)=(EP_KB(i)*EP_KB(j))/KB;
end
end
end
for i=1:15;
for j=1:15;
T(i,j)=600;
P(i,j)=101.325;
T_STAR(i,j)=KB*T(i,j)/EP_BOTH(i,j);
end
end
for i=1:15;
for j=1:15;
if (i==j)
SEKMA_BOTH(i,j)=0;
else
SEKMA_BOTH(i,j)=(SEKMA(i)+SEKMA(j))/2;
end
end
end
for i=1:15;
for j=1:15;
if (i==j)
MW_BOTH(i,j)=0;
else
MW_BOTH(i,j)=2/(1/MW(i)+1/MW(j));
end
end
end
for i=1:15;
for j=1:15;
if (i==j)
MU_D(i,j)=0;
else
MU_D(i,j)=A/((T_STAR(i,j))^(B))+B/exp(D*T_STAR(i,j))+E/exp(F*T_STAR(i,j))+G/exp(H*T_STAR(i,j));
end
end
end
for i=1:15;
for j=1:15;
if (i==j)
D(i,j)=0;
else
DD(i,j)=(0.0266*T(i,j)^1.5)/(P(i,j)*((MW_BOTH(i,j))^(0.5))*(SEKMA_BOTH(i,j))^(2)*MU_D(i,j));
end
end
end
for i=1:15;
for j=1:15;
if (i==j)
L(i,j)=0;
elseif (i<=14)&(j<=14)
L(i,j)=(X(j)*(MW(j)*X(j)+MW(i)*X(i)))/(DD(i,j)*MW(i))+(X(j+1)*(MW(j)*X(j)))/(MW(i)*DD(i,j+1));
elseif (i==15)&(j==15)
L(i,j)=X(j)*(MW(j)*X(j)+MW(i)*X(i))+X(j-1)*((MW(j)*X(j)))/(MW(i)*DD(i,j));
end
end
end
for i=1:15;
for j=1:15;
F(i,j)=0.01*randn(1,1);
end
end
for i=1:15;
for j=1:15;
delta_F(i,j)=F(i,j)-F(i,i);
end
end
for i=1:15;
z(i)=X(i)*MW(i);
end
MW_MIX=sum(z);
for i=1:15;
for j=1:15;
D(i,j)=((X(i)*MW_MIX)/MW(j))*delta_F(i,j);
end
end
DD
surf(D)
for i=1:15
for j=1:15;
if (D(i,j)==max(max(D)))
imax=i
jmax=j
NAME(imax,:)
NAME(jmax,:)
else
display('no');
end
end
end
Unless otherwise noted, all content on this site is @Copyright by Ahmed Al Makky 2012-2013 - http://cfd2012.com