3. Блок схема

Программные обозначения:

Z_i – moli (мольная доля i-ой компоненты в смеси);

P – P (пластовое давление);

T – T (температура в пласте);

P_кр – Pkri (критическое давление i-ой компоненты);

Т_кр – Тkri (критическая температура i-ой компоненты);

– wi (ацентрический фактор i-ой компоненты);

x_i – Xdi (молярная доля жидкой фазы в смеси i-ой компоненты);

y_i – Ydi (молярная доля газовой фазы в смеси i-ой компоненты);

НАЧАЛО

Вводим критические параметры, ацентрические факторы, P, T.

Вычисляем а_i и b_i используя формулы (4)-(5).

Находим P_si по эмпирическим формулам (8) для каждого отдельного компонента.

Рассчитываем пробные значения коэффициентов распределения K_i по формуле (7).

Рассчитываем y_i , x_i используя (9),(10)

Вычисляем a и b по правилам смешения (13), (14) для y_i

Вычисляем a и b по правилам смешения (13), (14) для x_i

Вычисляем Q, k используя (15),(15’).

Вычисляем Q, k используя (15),(15’).

Q<0

k<0

Q<0

да нет да нет

Используем (16),(17)

Используем (16),(17)

k<0

да нет да нет

Используем (18)

Используем (19)

Используем (18)

Используем (19)

Выбираем максимальное y

Выбираем минимальный y

Считаем Z

Считаем Z

Вычисляем φ по формуле (21)

Вычисляем φ по формуле (20)

Пересчитываем K_i по формуле (23)

нет

если

да

Выводим на экран x_i, y_i, V^L, V^V

КОНЕЦ

3. Код программы

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.