Код программы
program termodinamika;
var Pkr1,Pkr2,Pkr3,Tkr1,Tkr2,Tkr3,P,T,y,k,q,akr1,akr2,akr3,a,a1,a2,a3,b,b1,b2,b3,R,Qq,m1,m2,m3,w1,w2,w3,D,F,g1,g2,g3:real;
mol1,mol2,mol3,v,dv,let,let1,Z,ph,bj,aj,bg,ag:real;
S,fg1,fg2,fg3,ff1,ff2,ff3,vf,vg:real;
y1,y2,y3,alfa,beta,C,osn,pow, Vmol, Lmol, Ps1,Ps2,Ps3,K1,K2,K3,Xd1,Xd2,Xd3,Yd1,Yd2,Yd3:real;
begin
R:=8.31;
P:=200*100000;
T:=200;
mol1:=0.4; Vmol:=0.01;
mol2:=0.2; Lmol:=0.99;
mol3:=0.4;
Pkr1:=45.4*100000; {CH4}
Pkr2:=72.8*100000; {C02}
Pkr3:=36.0*100000; {C4H10}
Tkr1:=190.6;
Tkr2:=304.2;
Tkr3:=408.1;
w1:=0.008;
w2:=0.225;
w3:=0.176;
m1:=0.37+(1.5*w1)-(0.269*w1*w1);
m2:=0.37+(1.5*w2)-(0.269*w2*w2);
m3:=0.37+(1.5*w3)-(0.269*w3*w3);
akr1:=(0.45*R*R*Tkr1*Tkr1)/Pkr1;
akr2:=(0.45*R*R*Tkr2*Tkr2)/Pkr2;
akr3:=(0.45*R*R*Tkr3*Tkr3)/Pkr3;
b1:=0.077*R*Tkr1/Pkr1;
b2:=0.077*R*Tkr2/Pkr2;
b3:=0.077*R*Tkr3/Pkr3;
g1:=(1+(m1-(T/Tkr1)))*(1+(m1-(T/Tkr1)));
g2:=(1+(m2-(T/Tkr2)))*(1+(m2-(T/Tkr2)));
g3:=(1+(m3-(T/Tkr3)))*(1+(m3-(T/Tkr3)));
a1:=akr1*g1;
a2:=akr2*g2;
a3:=akr3*g3;
Ps1:=(exp(5.37*(1+w1)*(1-(Tkr1/T)))*Pkr1);
Ps2:=(exp(5.37*(1+w2)*(1-(Tkr2/T)))*Pkr2);
Ps3:=(exp(5.37*(1+w3)*(1-(Tkr3/T)))*Pkr3);
K1:=Ps1/P;
K2:=Ps2/P;
K3:=Ps3/P;
repeat
BEGIN
Xd1:=mol1/(Vmol*(K1-1)+1);
Xd2:=mol2/(Vmol*(K2-1)+1);
Xd3:=mol3/(Vmol*(K3-1)+1);
Yd1:=(mol1*K1)/(Lmol*(K1-1)+1);
Yd2:=(mol2*K2)/(Lmol*(K2-1)+1);
Yd3:=(mol3*K3)/(Lmol*(K3-1)+1);
{module for Z} {fluid}
bj:=(Xd1*b1)+(Xd2*b2)+(Xd3*b3);
aj:=0.2*Xd1*Xd1*a1+0.2*Xd2*Xd2*a2+0.2*Xd3*Xd3*a3+2*(0.2*Xd1*Xd2*sqrt(a1*a2)+0.2*Xd1*Xd3*sqrt(a1*a3)+0.2*Xd2*Xd3*sqrt(a2*a3));
{par}
bg:=(Yd1*b1)+(Yd2*b2)+(Yd3*b3);
ag:=0.2*Yd1*Yd1*a1+0.2*Yd2*Yd2*a2+0.2*Yd3*Yd3*a3+2*(0.2*Yd1*Yd2*sqrt(a1*a2)+0.2*Yd1*Yd3*sqrt(a1*a3)+0.2*Yd2*Yd3*sqrt(a2*a3));
begin {fluid}
b:=bj; a:=aj;
F:=(b*P)/(R*T);
D:=((a*F*F)/(b*b*P))-(3*F*F)-(2*F);
k:=D-((1/3)*(F-1)*(F-1));
q:=(2/27)*((F-1)*(F-1)*(F-1))-((4/3)*F*D)-(2*F*F*F)-(F*F)+((1/3)*D);
Qq:=(k/3)*(k/3)*(k/3)+(q/2)*(q/2);
if (Qq<0) and (k<0) then
begin
C:=(-q)/(2*sqrt((-k)*k*k/27));
alfa:=arctan(sqrt((1-C*C)/(C*C)));
y1:=2*sqrt(-k/3)*cos(alfa/3);
y2:=(-2)*sqrt(-k/3)*cos((alfa/3)+(pi/3));
y3:=(-2)*sqrt(-k/3)*cos((alfa/3)+(pi/3));
if (y1<y2) and (y1<y3) then y:=y1;
if (y2<y1) and (y2<y1) then y:=y2;
if (y3<y2) and (y3<y1) then y:=y3;
end;
if (Qq>=0) and (k>0) then
begin
beta:=arctan((2/q)*sqrt(k*k*k/27));
osn:=sin(beta/2)/cos(beta/2);
if osn<0 then Pow:=(-1)*exp((1/3)*Ln(abs(osn)))
else
if osn>0 then Pow:=exp((1/3)*Ln(abs(osn)))
else pow:=0;
osn:=round(osn);
alfa:=arctan(pow);
y:=2*sqrt(k/3)*(cos(2*alfa)/sin(2*alfa));
end;
if (Qq>=0) and (k<0) then
begin
S:=(2/q)*sqrt((-k)*k*k/27);
beta:=arctan(sqrt(S*S/(1-S*S)));
osn:=sin(beta/2)/cos(beta/2);
if osn<0 then Pow:=(-1)*exp((1/3)*Ln(abs(osn)));
if osn>0 then Pow:=exp((1/3)*Ln(abs(osn)))
else pow:=0;
alfa:=arctan(pow);
y:=(-2)*sqrt((-k)/3)*(1/sin(2*alfa));
end;
V:=((R*T*Y1)/p)+(b/(T*T)); Vf:=V;
Z:=Y+((b*P)/(R*T*3));
end;
dv:=1;
Let:=0;
while dv<1000000 do
begin
Let1:=(((R*T)/V)-P)*dv;
Let:=Let+Let1;
dv:=dv+1;
end;
Let1:=(Z-1)-Ln(Z)+(1/(R*T))*Let;
ff1:=Let1*P*mol1;
ff2:=Let1*P*mol2;
ff3:=Let1*P*mol3;
begin {par}
a:=ag; b:=bg;
F:=(b*P)/(R*T);
D:=((a*F*F)/(b*b*P))-(3*F*F)-(2*F);
k:=D-((1/3)*(F-1)*(F-1));
q:=(2/27)*((F-1)*(F-1)*(F-1))-((4/3)*F*D)-(2*F*F*F)-(F*F)+((1/3)*D);
Qq:=(k/3)*(k/3)*(k/3)+(q/2)*(q/2);
if (Qq<0) and (k<0) then
begin
C:=(-q)/(2*sqrt((-k)*k*k/27));
alfa:=arctan(sqrt((1-C*C)/(C*C)));
y1:=2*sqrt(-k/3)*cos(alfa/3);
y2:=(-2)*sqrt(-k/3)*cos((alfa/3)+(pi/3));
y3:=(-2)*sqrt(-k/3)*cos((alfa/3)+(pi/3));
if (y1>y2) and (y1>y3) then y:=y1;
if (y2>y1) and (y2>y1) then y:=y2;
if (y3>y2) and (y3>y1) then y:=y3;
end;
if (Qq>=0) and (k>0) then
begin
beta:=arctan((2/q)*sqrt(k*k*k/27));
osn:=sin(beta/2)/cos(beta/2);
if osn<0 then Pow:=(-1)*exp((1/3)*Ln(abs(osn)))
else
if osn>0 then Pow:=exp((1/3)*Ln(abs(osn)))
else pow:=0;
osn:=round(osn);
alfa:=arctan(pow);
y:=2*sqrt(k/3)*(cos(2*alfa)/sin(2*alfa));
end;
if (Qq>=0) and (k<0) then
begin
S:=(2/q)*sqrt((-k)*k*k/27);
beta:=arctan(sqrt(S*S/(1-S*S)));
osn:=sin(beta/2)/cos(beta/2);
if osn<0 then Pow:=(-1)*exp((1/3)*Ln(abs(osn)));
if osn>0 then Pow:=exp((1/3)*Ln(abs(osn)))
else pow:=0;
alfa:=arctan(pow);
y:=(-2)*sqrt((-k)/3)*(1/sin(2*alfa));
end;
V:=((R*T*Y1)/p)+(b/(T*T)); Vg:=V;
Z:=Y+((b*P)/(R*T*3));
dv:=1;
Let:=0;
while dv<1000000 do
begin
Let1:=(((R*T)/V)-P)*dv;
Let:=Let+Let1;
dv:=dv+1;
end;
Let1:=(Z-1)-Ln(Z)+(1/(R*T))*Let;
fg1:=Let1*P*mol1;
fg2:=Let1*P*mol2;
fg3:=Let1*P*mol3;
K1:=ff1/fg1;
K2:=ff2/fg2;
K3:=ff3/fg3;
end;
END;
until (1-K1)*(1-K1)+(1-K2)*(1-K2)+(1-K3)*(1-K3)<10E-12;
writeln ('X1= ',Xd1:1:6);
writeln ('X2= ',Xd2:2:6);
writeln ('X3= ',Xd3:2:6);
writeln ('Y1= ',Yd1:2:6);
writeln ('Y2= ',Yd2:2:6);
writeln ('Y3= ',Yd3:2:6);
writeln ('V Luquid ',Vf);
writeln ('V Vapour ',Vg);
readln;
end.
Дата добавления: 2015-10-20 | Просмотры: 541 | Нарушение авторских прав
|