Газ природный. Методы расчета физических свойств. Определение коэффициента сжимаемости. Межгосударственный стандарт

-

-

0,0

0,0

0,0

0,0

N₂

-

-

-

-

-

0,0

-0,064

0,0

СО₂

-

-

-

-

-

-

0,0

-0,062

H₂S

-

-

-

-

-

-

-

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

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

H=47.479*BMY+128.64D0

RETURN

3 Т=Т1

ТС=Т1-Т0

P=PRESS

IF(PRESS.LE.0D0.OR.PRESS.GT.12D0)Z=0D0

IF(T1.LT.250D0.OR.T1.GT.340D0)Z=0D0

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

CALL B11BER(T,H,B11)

CALL BBER(T,B11,B,Z)

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

CALL CBER(T,H,C,Z)

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

CALL ITER2(P,T,B,C,Z)

133 RETURN

END

SUBROUTINE B11BER(T,H,B11)

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

COMMON/BBLOK/BR11H0(3),BR11H1(3),BR11H2(3),BR22(3),BR23(3),BR33(3)

T2=T*T

B11=BR11H0(1)+BR11H0(2)*T+BR11H0(3)*T2+

*(BR11H1(1)+BR11H1(2)*T+BR11H1(3)*T2)*H+

*(BR11H2(1)+BR11H2(2)*T+BR11H2(3)*T2)*H*H

END

SUBROUTINE BBER(T,B11,BEFF,Z)

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

COMMON/BBLOK/BR11H0(3),BR11H1(3),BR11H2(3),BR22(3),BR23(3),BR33(3)

COMMON/ZETA/Z12,Z13,Y12,Y13,Y123

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

T2=T*T

B22=BR22(1)+BR22(2)*T+BR22(3)*T2

B23=BR23(1)+BR23(2)*T+BR23(3)*T2

B33=BR33(1)+BR33(2)*T+BR33(3)*T2

BA13=B11*B33

IF(BA13.LT.0D0) THEN

Z=0D0

RETURN

ENDIF

ZZZ=Z12+(320D0-T)**2*1.875D-5

BEFF=X11*B11+X12*ZZZ*(B11+B22)+2.*X13*Z13*DSQRT(BA13)+

*X22*B22+2.*X23*B23+X33*B33

END

SUBROUTINE CBER(T,H,CEFF,Z)

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

COMMON/CBLOK/CR111H0(3),CR111H1(3),CR111H2(3),CR222(3),CR223(3),

*CR233(3),CR333(3)

COMMON/ZETA/Z12,Z13,Y12,Y13,Y123

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

T2=T*T

C111=CR111 H0(1)+CR111H0(2)*T+CR111H0(3)*T2+

*(CR111H1(1)+CR111H1(2)*T+CR111H1(3)*T2)*H+

*(CR111H2(1)+CR111H2(2)*T+CR111H2(3)*T2)H*H

C222=CR222(1)+CR222(2)*T+CR222(3)*T2

C223=CR223(1)+CR223(2)*T+CR223(3)*T2

C233=CR233(1)+CR233(2)*T+CR233(3)*T2

C333=CR333(1)+CR333(2)*T+CR333(3)*T2

CA112=C111*C111*C222

CA113=C111*C111*C333

CA122=C111*C222*C222

CA123=C111*C222*C333

CA133=C111^<C333*C333

IF(CA112.LT.0D0.OR.CA113.LT.0D0.OR.CA122.LT.0DO0.OR.

*CA123.LT.0D0.OR.CA133.LT.0D0)THEN

Z=0D0

RETURN

ENDIF

D3REP=1D0/3D0

CEFF=X1*X11*C111+3D0*X11*X2*(CA112)**D3REP*(Y12+(T-270D0)*.0013D0)

*+3.*X11*X3*(CA113)**D3REP*Y13+

*3.*X1*X22*(CA122)**D3REP*(Y12+(T-270D0)*.0013D0)+

*6.*X1*X2*X3*(CA123)**D3REP*Y123+3.*X1*X33*(CA133)**D3REP*Y13+

*X22*X2*C222+3.*X22*X3*C223+3.*X2*X33*C233+X3*X33*C333

END

С Подпрограмма, реализующая схему Кардано для определения

С фактора сжимаемости из уравнения состояния

SUBROUTINE ITER2(P,T,Bm,Cm,Z)

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

B1=1D3*P/2.7715/T

B0=Bl*Bm

C0=Bl**2*Cm

A1=1D0+B0

A0=1D0+1.5*(B0+C0)

A01=A0**2-A1**3

IF(A01.LE.0D0) THEN

Z=0D0

RETURN

ENDIF

A=A0-A01**0.5

A2=DABS(A)**(1D0/3D0)

IF(A-LT.0D0) A2=-A2

Z=(1D0+A2+A1/A2)/3.

END

BLOCK DATA BDGRG2

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

COMMON/BBLOK/BR11H0(3),BR11H1(3),BR11H2(3),BR22(3),BR23(3),

*BR33(3)/CBLOK/CR111H0(3),CR111H1(3),CR111H2(3),CR222(3),

*CR223(3),CR233(3),CR333(3)

COMMON/ZETA/Z12,Z13,Y12,Y13,Y123

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

DATA BR11H0/-.425468D0,.2865D-2,-.462073D-5/,

* BR11H1/.877118D-3,-.556281D-5,.881514D-8/,

* BR11H2/-.824747D-6,.431436D-8,-.608319D-11/,

* BR22/-.1446D0,.74091D-3,-.91195D-6/,

* BR23/-.339693D0,.161176D-2,-.204429D-5/,

* BR33/-.86834D0,.40376D-2,-.51657D-5/

DATA CR111H0/-.302488D0,.195861D-2,-.316302D-5/,

* CR111 H1/.646422D-3,-.422876D-5,.688157D-8/,

* CR111H2/-.332805D-6,.22316D-8,-.367713D-11/,

* CR222/.78498D-2,-.39895D-4,.61187D-7/,

* CR223/.552066D-2,-.168609D-4,.157169D-7/,

* CR233/.358783D-2,.806674D-5,-.325798D-7/,

* CR333/.20513D-2,.34888D-4,-.83703D-7/

DATA Z12/.72D0/,Z13/-.865D0/,Y12/.92D0/,Y13/.92D0/,Y123/1.1D0/

DATA GM2/28.0135D0/,GM3/44.01D0/,

* FA/22.414097D0/,FB/22.710811D0/,

* TO/273.15D0/,R/.0831451D0/

END 46

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

С * *

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

С * газа по уравнению состояния AGA8-92DC. *

C * *

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

SUBROUTINE AGA8DC(ICALC)

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

REAL*8 KI,KIJ,KD

COMMON/RM/RM/Y1/Y(19)/NC1/NC/NI1/NI(19)/EFI/EI(19),KI(19),

*GI(19),QI(19),FI(19)

*/INTER1/EIJ(19,19),UIJ(19,19),KIJ(19,19),GIJ(19,19)

*/EFD/ED(19),KD(19),GD(19),QD(19),FD(19)/Z/Z

RM=8.31448D0

IF(ICALC.NE.l) GO TO 3

CALL COMPO1

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

CALL PARIN1

DO 75 I=1,NC

EI(I)=ED(NI(I))

KI(I)=KD(NI(I))

GI(I)=GD(NI(I))

QI(I)=QD(NI(I))

FI(I)=FD(NI(I))

DO 123 J=1,NC

IF(I.GE.J) GO TO 123

EIJ(I,J)=EIJ(NI(I),NI(J))

UIJ(I,J)=UIJ(NI(I),NI(J))

KIJ(I,J)=KIJ(NI(I),NI(J))

GIJ(I,J)=GIJ(NI(I),NI(J))

123 CONTINUE

75 CONTINUE

CALL PARMI1

3 CALL PHASE1

133 RETURN

END

SUBROUTINE COMPO1

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

DIMENSION ZNI(25),YI(25)

COMMON/YI/Y(19)/YI/YC(25)/NC1/NC/NT1/NI(19)/NPR/NPR

DATA ZNI/.9981D0,.992D0,.9834D0,.9682D0,.971D0,.9997D0,.9947D0,

*.99D0,.993D0,.994D0,985D0,.945D0,.953D0,1D0,.919D0,

*.936D0,.876D0,.892D0,3*1D0,1.0005D0,1.0006D0,.9996D0,.9993D0/

DO 100 I=1,25

100 YI(I)=YC(I)

YI(13)=YI(13)+YI(14)

YI(14)=0D0

IF(NPR.EQ.0D0) GO TO 5

YI(17)=YI(17)+YI(19)+YI(20)+YI(21)

YI(19)=0D0

YI(20)=0D0

YI(21)=0D0

SUM=0D0

DO 7 I=1,25

7 SUM=SUM+YI(I)/ZNI(I)

DO 9 I=1,25

9 YI(I)=YI(I)/ZNI(I)/SUM

5 YI(2)=YI(2)+YI(9)+YI(10)

YI(9)=0D0

YI(10)=O0D0

YI(3)=YI(3)+YI(11)

YI(11)=0D0

YI(15)=YI(15)+YI(16)

YI(16)=0D0

YI(17)=YI(17)+YI(18)

YI(18)=0D0

NC=0

ИС=0

YSUM=0D0

DO 11 1=1,25

IF((I.GE.9.AND.I.LE.11).OR.I.EQ.14.0R.I.EQ.16.0R.I.EQ.18)

*ИС=ИС+1

IF(YI(I).EQ.0D0) GO TO 11

NC=NC+1

NI(NC)=I-ИC

Y(NC)=YI(I)

YSUM=YSUM+Y(NC)

11 CONTINUE

CALL MOLDOl(YI)

DO 13 I=1,NC

13 Y(I)=Y(I)/YSUM

RETURN

END

SUBROUTINE MOLDO1(YI)

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

DIMENSION YI(25)

COMMON/Z/Z

Z=-1D0

YS=0D0

DO 1 I=9,25

1 YS=YS+YI(I)

1F(YI(1).LT.0.65D0.OR.YI(2).GT.0.15D0.OR.YI(3).GT.0.035D0.OR.

*YI(4).GT.0.015D0.OR.YI(5).GT.0.015D0.OR.YS.GT.0.01D0) Z=0D0

IF(YI(6).GT.0.2D0.OR.YI(7).GT.0.15D0.OR.Y1(8).GT.5D-5) Z=O0D0

RETURN

END

SUBROUTINE PARIN1

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

REAL*8 KIJ

COMMON/INTER1/EIJ(19,19),UIJ(19,19),KIJ(19,19),GIJ(19,19)

DO 1 I=1,19

DO 1 J=1,19

EIJ(I,J)=1D0

UIJ(I,J)=1D0

KIJ(I,J)=1D0

1 GIJ(I,J)=1D0

EIJ(1,6)=0.97164D0

UIJ(1,6)=0.886106D0

KIJ(1,6)=1.00363D0

EIJ(1,7)=0.960644D0

UIJ(1,7)=0.963827D0

KIJ(1,7)=0.995933D0

GIJ(1,7)=0.807653D0

EIJ(1,3)=0.99605D0

UIJ(1,3)=1.02396D0

EU(1,17)=1.17052D0

UIJ(1,17)=1.15639D0

KIJ(1,17)=1.02326D0

GIJ(1,17)=1.95731D0

EIJ(1,18)=0.990126D0

EIJ(1,5)=1.01953D0

EIJ(1,4)=0.995474D0

UIJ(1,4)=1.02128D0

EIJ(1,10)=1.00235D0

EIJ(1,9)=1.00305D0

EIJ(1,11)=1.01293D0

EIJ(1,12)=0.999758D0

EIJ(1,13)=0.988563D0

EIJ(6,7)=1.02274D0

UIJ(6,7)=0.835058D0

KIJ(6,7)=0.982361D0

GIJ(6,7)=0.982746D0

EIJ(2,6)=0.97012D0

UIJ(2,6)=0.816431D0

KIJ(2,6)=1.00796D0

EIJ(3,6)=0.945939D0

UIJ(3,6)=0.915502D0

EIJ(6,17)=1.08632D0

UIJ(6,17)=0.408838D0

KIJ(6,17)=1.03227D0

EIJ(6,18)=1.00571D0

EIJ(5,6)=0.946914D0

EIJ(4,6)=0.973384D0

UIJ(4,6)=0.993556D0

EIJ(6,10)=0.95934D0

EIJ(6,9)=0.94552D0

EIJ(6,11)=0.93788D0

EIJ(6,12)=0.935977D0

EIJ(6,13)=0.933269D0

EIJ(2,7)=0.925053D0

UIJ(2,7)=0.96987D0

KIJ(2,7)=1.00851D0

GIJ(2,7)=0.370296D0

EIJ(3,7)=0.960237D0

EIJ(7,17)=1.28179D0

EIJ(7,18)=1.5D0

UIJ(7,18)=0.9D0

EIJ(5,7)=0.906849D0

EIJ(4,7)=0.897362D0

EIJ(7,10)=0.726255D0

EIJ(7,9)=0.859764D0

EIJ(7,11)=0.766923D0

EIJ(7,12)=0.782718D0

EIJ(7,13)=0.805823D0

EIJ(2,3)=1.03502D0

UIJ(2,3)=1.0805D0

KIJ(2,3)=1.00046D0

EIJ(2,17)=1.16446D0

UIJ(2,17)=1.61666D0

KIJ(2,17)=1.02034D0

UIJ(2,5)=1.25D0

EIJ(2,4)=1.01306D0

UIJ(2,4)=1.25D0

UIJ(2,10)=1.25D0

EIJ(2,9)=1.00532D0

UIJ(2,9)=1.25D0

EIJ(3,17)=1.034787D0

EIJ(3,4)=1.0049D0

EIJ(17,18)=1.1D0

EIJ(5,17)=1.3D0

EIJ(4,17)=1.3D0

RETURN

END

SUBROUTINE PARMI1

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

REAL*8 KI,KIJ,KM

INTEGER GN,QN,FN

DIMENSION EIJM(19,19),GIJM(19,19)

COMMON/Y1/Y(19)/NC1/NC/EFI/EI(19),KI(19),GI(19),QI(19),FI(19)

*/INTER1/EIJ(19,19),UIJ(19,19),KIJ(19,19),GIJ(19,19)

*/KM/KM/COEF1/B1(13),C1(53)/AN/AN(53)

*/GQFN/GN(53),QN(53),FN(53)/UN/UN(53)

DO 1 I=1,NC

EIJM(I,I)=EI(I)

GIJM(I,I)=GI(I)

DO 1 J=1,NC

IF(I.GE.J) GO TO 1

EIJM(I,J)=EIJ(I,J)*(EI(I)*EI(J))**.5

GIJM(I,J)=GIJ(I,J)*(GI(I)+GI(J))/2.

1 CONTINUE

KM=0D0

UM=0D0

KM=0D0

UM=0D0

GM=0D0

QM=0D0

FM=0D0

DO 3 I=1,NC

KM=KM+Y(I)*KI(I)**2.5

UM=UM+Y(I)*EI(I)**2.5

GM=GM+Y(I)*GI(I)

QM=QM+Y(I)*QI(I)

3 FM=FM+Y(I)**2*FI(I)

KM=KM*KM

UM=UM*UM

DO 5 I=1,NC-1

DO 5 J=I+1,NC

UM=UM+2.*Y(I)*Y(J)*(UIJ(I,J)**5-1D0)*(EI(I)*EI(J))**2.5

GM=GM+2.*Y(I)*Y(J)*(GIJ(I,J)-1D0)*(GI(I)+GI(J))

5 KM=KM+2.*Y(I)*Y(J)*(KIJ(I,J)**5-1D0)*(KI(I)*KI(J))**2.5

KM=KM**.6

UM=UM**.2

DO 7 N=1,13

B1(N)=0D0

DO 9 I=1,NC

9 В1(N)=B1(N)+Y(I)*Y(I)(GIJM(I,I)+ 1D0-GN(N))**GN(N)*

*(QI(I)*QI(I)+1D0-QN(N))**QN(N)*(FI(I)+1D0-FN(N))*FN(N)*