Таблица Б.3 - Параметры бинарного взаимодействия xij

j

i

СН4

C2H6

С3Н8

н-C4H10

и4Н10

N2

CO2

H2S

СН4

0,0

0,036

0,076

0,121

0,129

0,060

0,074

0,089

C2H6

-

0,0

0,0

0,0

0,0

0,106

0,093

0,079

С3Н8

-

-

0,0

0,0

0,0

0,0

0,0

0,0

н-C4Н10

-

-

-

0,0

0,0

0,0

0,0

0,0

и4Н10

-

-

-

-

0,0

0,0

0,0

0,0

N2

-

-

-

-

-

0,0

0,022

0,211

CO2

-

-

-

-

-

-

0,0

0,089

H2S

-

-

-

-

-

-

-

0,0

Таблица Б.4 - Параметры бинарного взаимодействия lij

j

i

СН4

С2Н6

С3Н8

н4Н10

и-C4H10

N2

СО2

H2S

СН4

0,0

-0,074

-0,146

-0,258

-0,222

-0,023

-0,086

0,0

С2Н6

-

0,0

0,0

0,0

0,0

0,0

0,0

0,0

С3Н8

-

-

0,0

0,0

0,0

0,0

0,0

0,0

н-C4H10

-

-

-

0,0

0,0

0,0

0,0

0,0

и4Н10

-

-

-

-

0,0

0,0

0,0

0,0

N2

-

-

-

-

-

0,0

-0,064

0,0

СО2

-

-

-

-

-

-

0,0

-0,062

H2S

-

-

-

-

-

-

-

0,0

ПРИЛОЖЕНИЕ В

(рекомендуемое)

Листинг программы расчета коэффициента сжимаемости природного газа

C **********************************************************

C * *

С * Программа расчета коэффициента сжимаемости природного газа *

С * (основной модуль) *

С * *

C **********************************************************

IMPLICIT REAL*8(A-H,O-Z)

CHARACTER*26 AR(25)

DIMENSION PI(100),TI(100),ZP(100,100)

COMMON/P/P/T/T/RON/RON/YI/YC(25)/Z/Z/NPR/NPR

DATA AR/’ метана (СН4)’,’ этана (С2Н6)’,’ пропана (С3Н8)’,

*’ н-бутана (н-С4Н10)’,’ и-бутана (и-С4Н10)’,’ азота (N2)’,

*’ диоксида углерода (СO2)’,’ сероводорода (H2S)’,

*’ ацетилена (С2Н2)’,’ этилена (С2Н4)’,’ пропилена (С3Н6)’,

*’ н-пентана (н-С5Н12)’,’ и-пентана (и-C5H12)’,

*’ нео-пентана (нео-С5Н12)’,’ н-гексана (н-С6Н14)’,

*’ бензола (С6Н6)’,’ н-гептана (н-С7Н16)’,’ толуола (С7Н8)’,

*’ н-октана (н-С8Н18)’,’ н-нонана (н-С9Н20)’,

*’ н-декана (н-С10Н22)’,’ гелия (Не)’,’ водорода (Н2)’,

*’ моноксида углерода (СО)’,’ кислорода (О2)’/

200 WRITE(*,100)

CALL VAR(NVAR)

IF(NVAR.EQ.5) GO TO 134

WRITE(*,l00)

100 FORMAT(25(/))

WRITE(*,1)

1 FORMAT(’ Введите исходные данные для расчета.’/)

IF(NVAR.LE.2) THEN

WRITE(*,’(A)’)

*’ Плотность при 293.15 К и 101.325 кПа, в кг/куб.м ’

READ(*,*)RON

WRITE(*,53)

53 FORMAT(’ Введите 0, если состав азота и диоксида углерода’,

*’ задан в молярных долях’/

*’ или 1, если состав этих компонентов задан’,

*’ в объемных долях ’)

READ(*,*)NPR

IF(NPR.EQ.0) WRITE(*,3)

3 FORMAT (’ Значение молярной доли, в мол. %’)

IF(NPR.EQ.l) WRITE(*,33)

33 FORMAT(’ Значение объемной доли, в об. %’)

WRITE(*,’(A)’) ’ азота (N2)

READ(*,*)YA

YA = YA/100.

WRITE(*,’(A)’) ’ диоксида углерода (С02) ’

READ(*,*)YY

YY = YY/100.

ELSE

WRITE(*,35)

35 FORMAT(’ Введите 0, если состав задан в молярных долях’/

*’ или 1, если состав задан в объемных долях ’)

READ(*,*)NPR

IF(NPR.EQ.0) WRITE(*,3)

IF(NPR.EQ.l) WRITE(*,33)

DO 5 I=1,25

WRITE(*,’(A)’) AR(I)

READ(*,*)YC(I)

5 YC(I) = YC(I)/100.

ENDIF

WRITE(*,’(A)’)

*’ Введите количество точек по давлению: ’

READ(*,*)NP

WRITE(*,’(A)’)

*’ Введите количество точек по температуре: ’

READ(*,*)NT

WRITE(*,’(A)’)

*’ Введите значения давлений в МПа: ’

READ(*,*)(PI(I),I=1,NP)

WRITE(*,’(A)’)

*’ Введите значения температур в К: ’

READ(*,*)(TI(I),I=1,NT)

WRITE(*,’(A)’)

*’ Ввод исходных данных завершен. ’

P=.101325D0

T=293.15D0

ICALC=1

GO TO (10,20,30,40) NVAR

10 CALL NX19(YA,YY)

ZN=Z

GO TO 50

20 CALL GERG2(ICALC,YA,YY)

ZN=Z

GO TO 50

30 CALL AGA8DC(ICALC)

ZN=Z

GO TO 50

40 CALL VNIC(ICALC)

ZN=Z

50 CONTINUE

IF(Z.EQ.0D0) THEN

CALL RANGE(NRANGE)

IF(NRANGE) 134,134,200

ENDIF

ICALC=2

NTS=0

DO 7 I=1,NP

P=PI(I)

D07 J=1,NT

T=TI(J)

IF(NVAR.EQ.l) CALL NX19(YA,YY)

IF(NVAR.EQ.2) CALL GERG2(ICALC,YA,YY)

IF(NVAR.EQ.3) CALL AGA8DC(ICALC)

IF(NVAR.EQ.4) CALL VNIC(ICALC)

IF(Z.NE.0D0) NTS=NTS+1

ZP(I,J)=Z/ZN

7 CONTINUE

IF(NTS.EQ.0) THEN

CALL RANGE(NRANGE)

IF (NRANGE) 134,134,200

ELSE

I=1

9 IС=0

DO 11 J=1,NT

IF(ZP(I,J).EQ.0D0)

IC=IC+1

11 CONTINUE

IF(IC.EQ.NT) THEN

IF(I.NE.NP) THEN

DO 13 J=I,NP-1

PI(J)=PI(J+1)

DO 13 K=1,NT

13 ZP(J,K)=ZP(J+1,K)

ENDIF

NP=NP-1

ELSE

I=I+1

ENDIF

IF(I.LE.NP) GO TO 9

J=l

15 JS=0

DO 17 I=1,NP

IF(ZP(I,J).EQ.0D0) JS=JS+1

17 CONTINUE

IF(JS.EQ.NP) THEN

IF(J.NE.NT) THEN

DO 19 I=J,NT-1

ТI(I)=ТI(I+1)

DO 19 K=1,NP

19 ZP(K,I)=ZP(K,I+1)

ENDIF

NT=NT-1

ELSE

J=J+1

ENDIF

IF(J.LE.NT) GO TO 15

CALL TABL(YA,YY,PI,TI,ZP,NP,NT,NVAR,AR)

ENDIF

GO TO 200

134 STOP

END

SUBROUTINE VAR(NVAR)

WRITE(*,1)

1 FORMAT(//

*10X,’ Расчет коэффициента сжимаемости природного газа’//

*10Х,’ ----------------Метод расчета----------------- ’/

*10Х,’ ’/

*10Х,’ 1. Модифицированный метод NX 19 ’/

*10Х,’ ’/

*10Х,’ 2. Уравнение состояния GERG-91 ’/

*10Х,’ ’/

*10Х,’ 3. Уравнение состояния AGA8-92DC ’/

*10Х,’ ’/

*10Х,’ 4. Уравнение состояния ВНИЦ СМВ ’/

*10Х,’ ’/

*10Х,’---------------------------------------------------’/)

WRITE(*,5)

5 FORMAT(/,3X,

*’Введите порядковый номер метода расчета или 5 для выхода в ДОС’,

*)

READ(*,*)NVAR

RETURN

END

SUBROUTINE RANGE(NRANGE)

IMPLICIT REAL*8(A-H,О-Z)

COMMON/Z/Z

WRITE(*,1)

1 FORMAT(//

*’ Выбранная Вами методика при заданных параметрах «не работает»’/

*’ Продолжить работу программы ? 0 - нет, 1 - да ’)

READ(*,*)NRANGE

RETURN

END

SUBROUTINE TABL(YA,YY,PI,TI,ZP,NP,NT,NVAR,AR)

IMPLICIT REAL*8(A-H,О-Z)

CHARACTER*26 AR(25), FNAME

CHARACTER METH(4)*31,A*6,LIN1(5)*9,LIN2(5)*9,LIN3(6)*9,LIN4*9,

*AT(06)*28

CHARACTER*70 F,FZ(11,2)

DIMENSION PI(100),TI(100),ZP(100,100),ZPP(6)

COMMON/RON/RON/YI/YC(25)/NPR/NPR

DATA METH/

*’(модифицированный метод NX19)’,

*’(уравнение состояния GERG-91)’,

*’(уравнение состояния AGA8-92DC)’,

*’(уравнение состояния ВНИЦ СМВ)’/

DATA LIN1/5*’------’/,LIN2/5*’------’/,LIN3/6*’------’/,

*LIN4/’------’/,A/’ - ’/

DATA AT/

*’ T, K’,’ T, K’,’ T, K’,’ T,K’,

*’ T, K’,’ T, K’/

DATA FZ/

*’(3X,F5.2,2X,6(3X,F6.4))’,’(3X,F5.2,5X,A6,5(3X,F6.4))’,

*’(3X,F5.2,2X,2(3X,A6),4(3X,F6.4))’,’(3X,F5.2,2X,3(3X,A6),

*3(3X,F6.4))’,

*’(3X,F5.2,2X,4(3X,A6),2(3X,F6.4))’,’(3X,F5.2,2X,5(3X,A6),

*3X,F6.4)’,

*’(3X,F5.2,2X,5(3X,F6.4),3X,A6)’,’(3X,F5.2,2X,4(3X,F6.4),

*2(3X,A6))’,

*’(3X,F5.2,2X,3(3X,F6.4),3(3X,A6))’,’(3X,F5.2,2X,2(3X,F6.4),

*4(3X,A6))’,

*’(3X,F5.2,5X,F6.4,5(3X,A6))’,’(3X,F9.6,1X,F6.4,5(3X,F6.4))’,

*’(ЗX,F9.6,lX,A6,5(3X,F6.4))’,’(3X,F9.б,lX,A6,3X,A6,4(3X,F6.4))’,

*’(3X,F9.6,1X,A6,2(3X,A6),3(3X,F6.4))’,’(3X,F9.6,1X,A6,3(3X,A6),

*2(3X,F6.4))’,

*’(3X,F9.6,1X,A6,4(3X,A6),3X,F6.4)’,’(3X,F9.6,1X,F6.4,4(3X,F6.4),

*3X,A6)’,

*’(3X,F9.6,1X,F6.4,3(3X,F6.4),2(3X,A6))’,’(3X,F9.6,1X,F6.4),

*2(3X,F6.4),3(3X,A6))’,

*’(3X,F9.6,1X,F6.4,3X,F6.4,4(3X,A6))’,’(3X,F9.6,1X,F6.4,5(3X,A6))’/

22 WRITE(*,44)

44 FORMAT(//’ Устройство вывода результатов расчета ?,’)

WRITE(*,’(A)’)

*’ 0 - дисплей, 1 - принтер, 2 - файл на диске ’

READ(*,*)NYST

IF(NYST.EQ.0) OPEN(1,FILE=’CON’)

IF(NYST.EQ.l) OPEN(1,FILE=’PRN’)

IF(NYST.EQ.2) WRITE(*,’(A)’) ’ Введите имя файла ’

IF(NYST.EQ.2) READ(*,’(A)’)FNAME

IF(NYST.EQ.2) OPEN(1,FILE=FNAME)

IF(NYST.EQ.0) WRITE(*,100)

100 FORMAT(25(/))

IF(NYST.EQ.l) PAUSE

*’ Включите принтер, вставьте бумагу и нажмите <ВВОД> ’

WRITE(1,88)METH(NVAR)

88 FORMAT(

*13X,’Коэффициент сжимаемости природного газа.’/

*18Х,А31/)

NW=3

IF(NVAR.LE.2) THEN

WRITE(1,1)RON

1 FORMAT(’ Плотность при 293.15 К и 101.325 кПа ’,F6.4,’ кг/куб.м’)

NW=NW+1

IF(YA.NE.0D0.OR.YY.NE.0D0) THEN

IF(NPR.EQ.0) WRITE(1,3)

3 FORMAT(’ Содержание в мол. %’)

IF(NPR.EQ.l) WRITE(1,33)

33 FORMAT(’ Содержание в об.%’)

NW=NW+1

IF(YA.NE.0D0) THEN

WRITE(1,5)AR(6),YA* 100.

5 FORMAT(2(A26,F7.4))

NW=NW+1

ENDIF

IF(YY.NE.0D0) THEN

WRITE(1,5)AR(7),YY*100.

NW=NW+1

ENDIF

ENDIF

ELSE

IF(NPR.EQ.0) WRITE(1,3)

IF(NPR.EQ.l) WRITE(1,33)

NW=NW+1

I=1

9 J=I+1

13 CONTINUE

IF(YC(J).NE.0D0) THEN

WRITE(1,5)AR(I),YC(I)*100.,AR(J),YC(J)*100.

NW=NW+1

DO 11 I=J+1,25

IF(YC(I).NE.0D0.AND.I.NE.25) GO TO 9

IF(YC(I).NE.0D0.AND.I.EQ.25) THEN

WRITE(1,5)AR(I),YC(I)*100.

nw=nw+1

GO TO 99

ENDIF

11 CONTINUE

ELSE

J=J+1

IF(J.LE.25) THEN

GO TO 13

ELSE

WRITE(1,5)AR(I),YC(I)*100.

NW=NW+1

ENDIF

ENDIF

ENDIF

99 CONTINUE

IF(NW.GT.12.AND.NYST.EQ.0) THEN

WRITE(*,7)

7 FORMAT(/)

PAUSE ’ Для продолжения вывода нажмите <ВВОД> ’

WRITE(*,100)

NW=0

ENDIF

DO 15 I=1,NT,6

IF(NW.GT.12.AND.NYST.EQ.0) THEN

WRITE(*,7)

PAUSE ’ Для продолжения вывода нажмите <ВВОД> ’

WRITE(*,100)

NW=0

ENDIF

IF(NW.GT.46.AND.NYST.NE.O) THEN

WRITE(1,7)

WRITE(*,7)

IF(NYST.EQ.l)

PAUSE

*’ Для продолжения вывода вставьте бумагу и нажмите <ВВОД> ’

NW=0

ENDIF

IF(I+5.LE.NT) THEN

NL=6

ELSE

NL=NT-I+1

ENDIF

WRITE(1,7)

IF(NL.GT.1) WRITE(1,17)LIN2(1),(LIN1(K),K=1,NL-1)

IF(NL.EQ.l) WRITE(1,17)LIN2(1)

17 FORMAT(’ ------’,6A9)

WRITE(1,19)AT(NL)

19 FORMAT(’ ------’,A28)

IF(NL.GT.1)WRITE(1,21)LIN4,(LIN2(K),K=1,NL-1)

IF(NL.EQ.l) WRITE(1,21)LIN4

21 FORMAT(’ p, МПа ’,6А9)

WRITE(1,23)(TI(K),K=I,I+NL-1)

23 FORMAT(10X,6(:,’|’,F6.2))

WRITE(1,17)(LIN3(K),K=1,NL)

NW=NW+6

DO 25 J=1,NP

JP=1

IF(PI(J).EQ.0.101325D0) JP=2

NL1=0

NLN=0

DO 27K=I,I+NL-1

NL1=NL1+1

IF(ZP(J,K).EQ.0D0) THEN

ZPP(NL1)=A

NLN=NLN+1

ELSE

ZPP(NL1)=ZP(J,K)

ENDIF

27 CONTINUE

IF(NLN.EQ.NL) GO TO 133

IF(NLN.EQ.0) THEN

F=FZ(1,JP)

ELSE

IF(ZP(J,I).EQ.0D0) F=FZ(NLN+1,JP)

IF(ZP(J,I+NL-1).EQ.0D0) F=FZ(NLN+12-NL,JP)

ENDIF

IF(NLI.EQ.1)WRITE(1,F)PI(J),ZPP(1)

IF(NL1.EQ.2)WRITE(1,F)PI(J),ZPP(1),ZPP(2)

IF(NL1.EQ.3)WRITE(1,F)PI(J),ZPP(1),ZPP(2),ZPP(3)

IF(NL1.EQ.4)WRITE(1,F)PI(J),ZPP(1),ZPP(2),ZPP(3),ZPP(4)

IF(NL1.EQ.5)

*WRITE(1,F)PI(J),ZPP(1),ZPP(2),ZPP(3),ZPP(4),ZPP(5)

IF(NL1.EQ.6)

*WRITE(1,F)PI(J),ZPP(1),ZPP(2),ZPP(3),ZPP(4),ZPP(5),ZPP(6)

NW=NW+1

133 CONTINUE

IF(NW.EQ.20.AND.NYST.EQ.0) THEN

IF(J.EQ.NP.AND.I+NL-1.EQ.NT) GO TO 29

WRITE(*,7)

PAUSE ’ Для продолжения вывода нажмите <ВВОД> ’

WRITE(*,100)

NW=0

WRITE(1,7)

IF(NL.GT.1)WRITE(1,17)LIN2(1),(LIN1(K),K=1,NL-1)

IF(NL.EQ.l) WRITE(1,17)LIN2(1)

WRITE(1,19)AT(NL)

IF(NL.GT.1)WRITE(1,21)LIN4,(LIN2(K),K=1,NL-1)

IF(NL.EQ.l) WRITE(1,21)LIN4

WRITE(1,23)(TI(K),K=I,I+NL-1)

WRITE(1,17)(LIN3(K),K=1,NL)

NW=NW+6

ENDIF

IF(NW.EQ.54.AND.NYST.NE.0) THEN

IF(J.EQ.NP.AND.I+NL-1.EQ.NT) GO TO 29

WRITE(1,7)

WRITE(*,7)

IF(NYST.EQ.l) PAUSE

*’ Для продолжения вывода вставьте бумагу и нажмите <ВВОД> ’

NW=0

IF(NL.GT.1) WRITE(1,17)LIN2(1),(LIN1(K),K=1,NL-1)

IF(NL.EQ.l) WRITE(1,17)LIN2(1)

WRITE(1,19)AT(NL)

IF(NL.GT.1) WRITE(1,21)LIN4,(LIN2(K),K=1,NL-1)

IF(NL.EQ.l) WRITE(1,21)LIN4

WRITE(1,23)(TI(K),K=I,I+NL-1)

WRITE(1,17)(LIN3(K),K=1,NL)

NW=NW+6

ENDIF

25 CONTINUE

15 CONTINUE

29 CLOSE(l)

WRITE(*,7)

PAUSE ’ Вывод завершен, для продолжения работы нажмите <ВВОД> ’

WRITE(*,66)

66 FORMAT(/’ Назначить другое устройство вывода ?’,

*’, 0 - нет, 1 - да ’)

READ(*,*)NBOLB

IF(NBOLB.EQ.l) GO TO 22

RETURN

END

C **********************************************************

С * *

С * Подпрограмма расчета коэффициента сжимаемости природного *

С * газа по модифицированному методу NX19. *

C **********************************************************

SUBROUTINE NX19(YA,YY)

IMPLICIT REAL*8(A-H,O-Z)

COMMON/NCONT/NCONT/YA/Y(2)/RON/RON

Y(1)=YA

Y(2)=YY

CALL PTCONT

IF(NCONT.EQ.l) GO TO 134

CALL EA

CALL PHASEA

134 RETURN

END

SUBROUTINE PTCONT

IMPLICIT REAL*8(A-H,O-Z)

COMMON/NCONT/NCONT/Z/Z/P/P/T/T/YA/Y(2)/RON/RON

NCONT=0

IF(RON.LT.0.66D0.OR.RON.GT.1D0) NCONT=1

IF(Y(1).GT.0.2D0.OR.Y(2).GT.0.15D0) NCONT=l

IF(P.LE.0.D0.OR.T.LE.0.D0) NCONT=1

IF(T.LT.250.D0.OR.T.GT.340.D0) NCONT=1

IF(P.GT.12.D0) NCONT=1

IF(NCONT.EQ.1) Z=0D0

RETURN

END

SUBROUTINE EA

IMPLICIT REAL*8(A-H,O-Z)

COMMON/T/T/YA/Y(2)/RON/RON/P/P/PT/PA,TA/BI/B1,B2/T0/T0

PCM=2.9585*(1.608D0-0.05994*RON+Y(2)-.392*Y(1))

TCM=88.25*(0.9915D0+1.759*RON-Y(2)-1.681*Y(1))

PA=0.6714*P/PCM+0.0147

TA=0.71892*T/TCM+0.0007

DTA=TA-1.09D0

F=0D0

IF(PA.GE.0D0.AND.PA.LT.2D0.AND.DTA.GE.0D0.AND.DTA.LT.0.3D0)

F=75D-5*PA**2.3/DEXP(20.*DTA)+

*11D-4*DTA**0.5*(PA*(2.17D0-PA+1.4*DTA**0.5))**2

IF(PA.GE.0D0.AND.PA.LT.1.3D0.AND.DTA.GE.-0.25D0.AND.DTA.LT.0D0)

*F=75D-5*PA**2.3*(2D0-DEXP(20.*DTA))+

*1.317*PA*(1.69D0-PA**2)*DTA**4

IF(PA.GE.1.3D0.AND.PA.LT.2D0.AND.DTA.GE.-0.21D0.AND.DTA.LT.0D0)

*F=75D-5*PA**2.3*(2D0-DEXP(20.*DTA))+

*0.455*(1.3D0-PA)*(1.69*2.D0**1.25-PA**2)*(DTA*(0.03249D0+

*18.028*DTA**2)+DTA**2*(2.0167D0+DTA**2*(42.844D0+200.*DTA**2)))

T1=:TA**5/(TA**2*(6.60756*TA-4.42646D0)+3.22706D0)

T0=(TA**2*(1.77218D0-0.8879*TA)+0.305131D0)*T1/TA**4

B1=2.*T1/3.-TO**2

B0=T0*(T1-T0**2)+0.1*T1*PA*(F-1D0)

B2=(B0+(B0**2+B1**3)**0.5)**(1D0/3D0)

RETURN

END

SUBROUTINE PHASEA

IMPLICIT REAL*8(A-H,O-Z)

COMMON/Z/Z/PT/PA,TA/BI/B1,B2/T0/T0

Z=(1D0+0.00132/TA**3.25)**2*0.1*PA/(B1/B2-B2+T0)

RETURN

END

C *************************************************************

С * *

С * Подпрограмма расчета коэффициента сжимаемости природного *

С * газа по модифицированному уравнению состояния GERG-91. *

С * *

C *************************************************************

$NOTRUNCATE

SUBROUTINE GERG2(ICALC,YA,YY)

IMPLICIT REAL*8(A-H,O-Z)

COMMON/T/T1/P/PRESS/RON/RON/Z/Z

COMMON/XBLOK/X1,X2,X3,X11,X12,X13,X22,X23,X33

COMMON/MBLOK/GM2,GM3,FA,FB,TO,R

DATABMO/.0838137D0/,BM1/-.00851644D0/,WD0/134.2153D0/,

*WD1/1067.943D0/

Z=-1D0

IF(ICALC.EQ.2) GO TO 3

X2=YA

X3=YY

IF(RON.LT.0.66D0.OR.RON.GT.1D0)Z=0D0

IF(X2.LT.0D0.OR.X2.GT.0.2D0)Z=0D0

IF(X3.LT.0D0.OR.X3.GT.0.15D0) Z=0D0

IF(Z.EQ.0D0) GO TO 133

X1=1D0-X2-X3

Х11=Х1*Х1

X12=X1*X2

X13=X1*X3

X22=X2*X2

X23=X2*X3

X33=X3*X3

Z=1D0-(.0741*RON-.006D0-.063*YA-.0575*YY)**2

BMNG=24.05525*Z*RON

Y1=1D0-YA-YY

BMY=(BMNG-28.0135*YA-44.01*YY)/Y1

С Расчет теплоты сгорания эквивалентного углеводорода (Н)