Откройте актуальную версию документа прямо сейчас
Если вы являетесь пользователем интернет-версии системы ГАРАНТ, вы можете открыть этот документ прямо сейчас или запросить по Горячей линии в системе.
Вход
Купить систему ГАРАНТ
Получить демо-доступ
Узнать стоимость
Информационный банк
Подобрать комплект
Семинары
- Главная
- Межгосударственный стандарт ГОСТ 30319.3-96 "Газ природный. Методы расчета физических свойств. Определение физических свойств по уравнению состояния" (введен в действие постановлением Государственного комитета РФ по стандартизации, метрологии и сертификации от 30 декабря 1996 г. N 723)
- Приложение А (рекомендуемое). Листинг программы расчета физических свойств природного газа
Приложение А (рекомендуемое). Листинг программы расчета физических свойств природного газа
Приложение А
(рекомендуемое)
Листинг программы расчета физических свойств природного газа
Расчет физических свойств природного газа по уравнению состояния (1) и по уравнению для вязкости (15) реализован на ПЭВМ, совместимых с IBM PC/AT/XT, на языке программирования ФОРТРАН-77.
С *******************************************************************
С * *
С * Программа расчета физических свойств (плотности, показате- *
С * ля адиабаты, скорости звука и вязкости) природного газа по *
С * уравнению состояния ВНИЦ СМВ. *
С * *
С *******************************************************************
IMPLICIT REAL*8(A-H,O-Z)
CHARACTER*26 AR
DIMENSION PI(100),TI(100),ROP(100,100),PAP(100,100),
*WP(100,100),ETAP(100,100)
COMMON/P/P/T/T/RON/RON/YI/YC(25)/NPR/NPR/Z/Z/TS/RO,PA,W
*/ETA/ETA/AR/AR(25)
200 WRITE(*,300)
300 FORMAT(18(/))
WRITE(*,400)
400 FORMAT(
*' Расчет физических свойств природного газа'/
*' по уравнению состояния'/////)
WRITE(*,1)
1 FORMAT(' Введите исходные данные для расчета.'/)
WRITE(*,35)
35 FORMAT(' Введите 0, если состав задан в молярных долях'/
*' или 1, если состав задан в объемных долях '\)
READ(*,*)NPR
IF(NPR.EQ.l) THEN
WRITE(*,'(A\)')
*' Плотность при 293.15 К и 101.325 кПа, в кг/куб.м '
READ(*,*)RON
WRITE(*,33)
33 FORMAT(' Значение объемной доли, в об.%')
ELSE
RON=0D0
WRITE(*,3)
3 FORMAT(' Значение молярной доли, в мол.%')
ENDIF
DO 5 I=1,25
WRITE(*,'(A\)') AR(I)
READ(*,*)YC(I)
5 YC(I)=YC(I)/100.
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
CALL EOSVNIC(ICALC)
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)
DO 7 J=1,NT
T=TI(J)
CALL EOSVNIC(ICALC)
IF(Z.NE.0D0) THEN
NTS=NTS+1
ROP(I,J)=RO
PAP(I,J)=PA
WP(I,J)=W
ETAP(I,J)=ETA
ELSE
ROP(I,J)=0D0
PAP(I,J)=0D0
WP(I,J)=0D0
ETAP(I,J)=0D0
ENDIF
7 CONTINUE
500 WRITE(*,100)
100 FORMAT(25(/))
IF(NTS.EQ.0) THEN
CALL RANGE(NRANGE)
IF (NRANGE) 134,134,200
ELSE
I=1
9 IS=0
DO 11 J=1,NT
IF(ROP(I,J).EQ.0D0) IS=IS+1
11 CONTINUE
IF(IS.EQ.NT) THEN
IF(I.NE.NP) THEN
DO 13 J=I,NP-1
PI(J)=PI(J+1)
DO 13 K=1,NT
ROP(J,K)=ROP(J+1,K)
PAP(J,K)=PAP(J+1,K)
WP(J,K)=WP(J+1,K)
13 ETAP(J,K)=ETAP(J+1 ,K)
ENDIF
NP=NP-1
ELSE
I=I+1
ENDIF
lF(I.LE.NP) GO TO 9
J=l
15 JS=0
DO 17 I=1,NP
IF(ROP(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
TI(I)=TI(I+1)
DO 19 K=1,NP
ROP(K,I)=ROP(K,I+1)
PAP(K,I)=PAP(K,I+1)
WP(K,I)=WP(K,I+1)
19 ETAP(K,I)=ETAP(K,I+1)
ENDIF
NT=NT-1
ELSE
J=J+1
ENDIF
IF(J.LE.NT) GO TO 15
CALL PROP(NPROP)
IF(NPROP.EQ.5) GO TO 134
IF(NPROP.EQ.l) CALL TABL(PI,TI,ROP,NP,NT,NPROP)
IF(NPROP.EQ.2) CALL TABL(PI,TI,PAP,NP,NT,NPROP)
IF(NPROP.EQ.3) CALL TABL(PI,TI,WP,NP,NT,NPROP)
IF(NPROP.EQ.4) CALL TABL(PI,TI,ETAP,NP,NT,NPROP)
WRITE(*,'(A\)')
*' Продолжить вывод рассчитанных свойств ? 0 - нет, 1 - да '
READ(*,*)NCONT
IF(NCONT.EQ.l) GO TO 500
ENDIF
134 STOP
END
SUBROUTINE PROP(NPROP)
WRITE(*,1)
1 FORMATS(//
*10X,'-----Рассчитаны следующие физические свойства-------------'/
*10Х,' '/
*10Х,' 1. Плотность '/
*10Х,' '/
*10Х,' 2. Показатель адиабаты '/
*10Х,' '/
*10Х,' 3. Скорость звука '/
*10Х,' '/
*10Х,' 4. Коэффициент динамической вязкости '/
*10Х,' '/
*10Х,'----------------------------------------------------------'/)
WRITE(*,5)
5 FORMAT(/,3X,
*'Введите порядковый номер свойства для вывода результатов расче',
*'та'/
*' или 5 для выхода в ДОС '\)
READ(*,*)NPROP
RETURN
END
SUBROUTINE RANGE(NRANGE)
IMPLICIT REAL*8(A-H,O-Z)
COMMON/Z/Z
WRITE(*,1)
1 FORMAT(//
*' Метод расчета при заданных параметрах "не работает"'/
*' Продолжить работу программы ? 0 - нет, 1 - да '\)
READ(*,*)NRANGE
RETURN
END
SUBROUTINE TABL(PI,TI,ZP,NP,NT,NPROP)
IMPLICIT REAL*8(A-H,O-Z)
CHARACTER*26 AR,FNAME
CHARACTER PROP(4)*58A*6,LIN1(5)*9,LIN2(5)*9,LIN3(6)*9,LIN4*9,
*AT(6)*28,RAZM(4)*39
CHARACTER*70 F,FZ(11,2),FW(11,2)
DIMENSION PI(100),TI(100),ZP(100,100),ZPP(6)
COMMON/YI/YC(25)/NPR/NPR/AR/AR(25)
DATA PROP/
*' Плотность природного газа.',
*' Показатель адиабаты природного газа.',
*' Скорость звука природного газа.',
*' Коэффициент динамической вязкости природного газа.'/
DATA RAZM/
*' (в кг/куб.м)',' ',
*' (в м/с)',
*' (в мкПа*с)'/
DATA LIN1/5*'----------'/,LIN2/5*'----------'/,LIN3/6*'----------'/,
*LIN4/'-----------'/,А/' - '/
DATA AT/
*' Т, К',' Т, К',' Т, К',' Т, К',
*' Т, К',' Т, К'/
DATA FZ/
*'(3X,F5.2,2X,6(3X,F6.2))','(3X,F5.2,5X,A6,5(3X,F6.2))',
*'(3X,F5.2,2X,2(3X,A6),4(3X,F6.2))','(3X,F5.2,2X,3(3X,A6),
*3(3X,F6.2))',
*'(3X,F5.2,2X,4(3X,A6),2(3X,F6.2))','(3X,F5.2,2X,5(3X,A6),
*3X,F6.2)',
*'(3X,F5.2,2X,5(3X,F6.2),3X,A6)','(3X,F5.2,2X,4(3X,F6.2),
*2(3X,A6))',
*'(3X,F5.2,2X,3(3X,F6.2),3(3X,A6))','(3X,F5.2,2X,2(3X,F6.2),
*4(3X,A6))',
*'(3X,F5.2,5X,F6.2,5(3X,A6))','(3X,F9.6,1X,F6.2,5(3X,F6.2))',
*'(3X,F9.6,1X,A6,5(3X,F6.2))','(3X,F9.6,1X,A6,3X,A6,4(3X,F6.2))',
*'(3X,F9.6,1X,A6,2(3X,A6),3(3X,F6.2))','(3X,F9.6,1X,A6,3(3X,A6),
*2(3X,F6.2))',
*'(3X,F9.6,1X,A6,4(3X,A6),3X,F6.2)','(3X,F9.6,1X,F6.2,4(3X,F6.2),
*3X,A6)',
*'(3X,F9.6,1X,F6.2,3(3X,F6.2),2(3X,A6))','(3X,F9.6,1X,F6.2,
*2(3X,F6.2),3(3X,A6))',
*'(3X,F9.6,1X,F6.2,3X,F6.2,4(3X,A6))','(3X,F9.6,1X,F6.2,5(3X,A6))'/
DATA FW/
'(3X,F5.2,2X,6(4X,F5.1))','(3X,F5.2,5X,A6,5(4X,F5.1))',
*'(3X,F5.2,2X,2(3X,A6),4(4X,F5.1))','(3X,F5.2,2X,3(3X,A6),
*3(4X,F5.1))',
*'(3X,F5.2,2X,4(3X,A6),2(4X,F5.1))','(3X,F5.2,2X,5(3X,A6),
*4X,F5.1)',
*'(3X,F5.2,2X,5(4X,F5.1),3X,A6)','(3X,F5.2,2X,4(4X,F5.1),
*2(3X,A6))',
*'(3X,F5.2,2X,3(4X,F5.1),3(3X,A6))','(3X,F5.2,2X,2(4X,F5.1),
*4(3X,A6))',
*'(3X,F5.2,6X,F5.1,5(3X,A6))','(3X,F9.6,2X,F5.1,5(4X,F5.1))',
*'(3X,F9.6,1X,A6,5(4X,F5.1))','(3X,F9.6,1X,A6,3X,A6,4(4X,F5.1))',
*'(3X,F9.6,1X,A6,2(3X,A6),3(4X,F5.1))','(3X,F9.6,1X,A6,3(3X,A6),
*2(4X,F5.1))',
*'(3X,F9.6,1X,A6,4(3X,A6),4X,F5.1)','(3X,F9.6,2X,F5.1,4(4X,F5.1),
*3X,A6)',
*'(3X,F9.6,2X,F5.1,3(4X,F5.1),2(3X,A6))','(3X,F9.6,2X,F5.1,
*2(4X,F5.1),3(3X,A6))',
*'(3X,F9.6,2X,F5.1,4X,F5.1,4(3X,A6))','(3X,F9.6,2X,F5.1,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.1) PAUSE
*' Включите принтер, вставьте бумагу и нажмите <ВВОД> '
WRITE(1,88)PROP(NPROP),RAZM(NPROP)
88 FORMAT(A58/A39/)
NW=3
IF(NPR.EQ.0) WRITE(l,3)
3 FORMAT(' Содержание в мол.%')
IF(NPR.EQ.l) WRITE(1,33)
33 FORMAT(' Содержание в об.%')
NW=NW+1
I=1
9 J=I+l
13 CONTINUE
IF(YC(J).NE.0D0) THEN
WRITE(l,5)AR(I),YC(I)*100.,AR(J),YC(J)*100.
5 FORMAT(2(A26,F7.4))
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(l,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(l,5)AR(I),YC(I)*100.
NW=NW+1
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.0) 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(' |',А28)
IF(NL.GT.l) 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 27 K=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
IF(NPROP.NE.3) F=FZ(1,JP)
IF(NPROP.EQ3) F=FW(1,JP)
ELSE
IF(ZP(J,I).EQ.0D0.AND.NPROP.NE.3) F=FZ(NLN+1,JP)
IF(ZP(J,I+NL-l).EQ.0D0.AND.NPROP.NE.3)F=FZ(NLN+12-NL,JP)
IF(ZP(J,I).EQ.0D0.AND.NPROP.EQ.3)F=FW(NLN+l,JP)
IF(ZP(J,I+NL-1).EQ.ODO.AND.NPROP.EQ.3)F=FW(NLN+12-NL,JP)
ENDIF
IF(NLl.EQ.l) 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(l),(LINl(K),K=l,NL-l)
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(l,21)LIN4
WRITE(1,23)(TI(K),K=I,I+NL-1)
WRITE(1,17)(LIN3(К),К=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
SUBROUTINE EOSVNIC(ICALC)
IMPLICIT REAL*8(A-H,O-Z)
REAL*8 LIJ(8,8)
DIMENSION VC(8),TC(8),PII(8),DU(8,8)
COMMON/PARCD/VCD(8),TCD(8),PIID(8)/ABIJ/AIJ(10,8),BIJ(10,8)
*/B/B(10,8)/RM/RM/Y/Y(8)/BM/BM(8)/NI/NI(8)/NC/NC/RON/RON/PIM/PIM
СОММОN/СРСI/СРС1(20,5),СРС2(20,3)/IDGFD/TOID(8),МСОD(8),МСРD(8)
*/IDGF/CPC(20,8),TOI(8),MCO(8),MCP(8)
COMMON/P/P/Т/Т/Z/Z/TS/RO,PA,W/ETA/ETA
RM=8.31451D0
IF(ICALC.NE.1) GO TO 1
CALL COMPON
IF(Z.EQ.0D0)GO TO 133
DO 11111 J=l,8
DO 11111 I=1,20
IF(J.LE.5) CPC(I,J)=CPC1(I,J)
IF(J.GT.5) CPC(I,J)=CPC2(I,J-5)
11111 CONTINUE
CALL DDIJ(DIJ,LIJ)
DO 75 1=1,NC
TC(I)=TCD(NI(I))
VC(I)=BM(I)/VCD(NI(I))
PII(I)=PIID(NI(I))
MCO(1)=MCOD(NI(I))
MCP(I)=MCPD(NI(I))
TOI(I)=TOID(NI(I))
MP=MCO(I)+MCP(I)+1
DO 23 J=1,MP
23 CPC(J,I)=CPC(J,NI(I))
DO 123 J=1,NC
IF(I.GE.J) GO TO 123
DIJ(I,J)=DIJ(NI(I),NI(J))
LIJ(I,J)=LIJ(NI(I),NI(J))
123 CONTINUE
75 CONTINUE
CALL PARMIX(DIJ,LIJ,TC,VC,P11)
DO 27 1=1,10
DO 27 J=l,8
27 B(I,J)=AIJ(I,J)+BIJ(I,J)*PIM
IF(RON.NE.0D0) THEN
CALL PHASE
RON=0D0
GO TO 133
ENDIF
1 CALL PHASE
133 RETURN
END
SUBROUTINE COMPON
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION BMI(25),ROI(8),GI(8),YI(25)
COMMON/Y/Y(8)/BMM/BMM/BM/BM(8)/YI/YC(25)/NI/NI(8)/NC/NC/RON/RON
DATA BMI/16.043D0,30.07D0,44.097D0,2*58.123D0,28.0135D0,
*44.01D0,34.082D0,26.038D0,28.054D0,42.081D0,3*72.15D0,
*86.177D0,78.114D0,100.204D0,92.141D0,114.231D0,128.259D0,
*142.286D0,4.0026D0,2.0159D0,28.01D0,31.9988D0/
DATA ROI/0.6682D0,1.2601D0,1.8641D0,2.4956D0,2.488D0,
*1.1649D0,1.8393D0,1.4311 D0/
DO 100 I=1,25
100 YI(I)=YC(I)
IF(RON.NE.0D0) GO TO 333
BMM=0D0
DO 3333 I=1,25
3333 BMM=BMM+YI(I)*BMI(I)
333 YS=0D0
DO 55 I=9,25
55 YS=YS+YI(I)
YS1=0D0
DO 67 I=12,21
67 YS1=YS1+YI(I)
YS2=0D0
DO 69 I=22,25
69 YS2=YS2+YI(I)
YI(2)=YI(2)+YI(9)+YI(10)
YI(3)=YI(3)+YI(11)
YI(4)=YI(4)+YS1
YS3=YI(4)+YI(5)
IF(RON.NE.0D0.AND.YI(5).LT.0.01D0.AND.YS3.LT.0.03D0)THEN
YI(4)=YS3
YI(5)=0D0
ENDIF
IF(RON.EQ.0D0.AND.YI(5).LT.0.01D0.AND.YS3.LE.0.03D0)THEN
YI(4)=YS3
YI(5)=0D0
ENDIF
YI(6)=YI(6)+YS2
IF(RON.EQ.0D0) GO TO 555
ROM=0D0
DO 7 I=1,8
7 ROM=ROM+YI(I)*ROI(I)
DO 9 I=1,8
9 GI(I)=YI(I)*ROI(I)/ROM
SUM=0D0
DO 11 I=1,8
11 SUM=SUM+GI(I)/BMI(I)
SUM=1./SUM
DO 13 I=1,8
13 YI(I)=GI(I)*SUM/BMI(I)
555 NC=0
YSUM=0D0
DO 155 I=1,8
IF(YI(I).EQ.0D0) GO TO 155
NC=NC+1
NI(NC)=I
Y(NС)=YI(I)
YSUM=YSUM+Y(NC)
BM(NC)=BMI(I)
155 CONTINUE
CALL MOLDOL(YI,YS)
DO 551 I=1,NC
551 Y(I)=Y(I)/YSUM
RETURN
END
SUBROUTINE MOLDOL(YI,YS)
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION YI(25)
COMMON/Z/Z
Z=-1D0
IF(YI(l).LT.0.5D0.OR.YI(2).GT.0.2D0.OR.YI(3).GT.0.05D0.OR.
*YI(4).GT.0.03D0.OR.YI(5).GT.0.03D0.OR.YS.GT.0.01D0)Z=0D0
IF(YI(6).GT.0.3D0.OR.YI(7).GT.0.3D0.OR.YI(8).GT.0.3D0)Z=0D0
RETURN
END
SUBROUTINE DDIJ(DIJ,LIJ)
IMPLICIT REAL*8(A-H,0-Z)
REAL*8 LIJ(8,8)
DIMENSION DIJ(8,8)
DO 1 I=1,8
DO 1 J=l,8
LIJ(I,J)=0.D0
1 DIJ(I,J)=0.D0
DIJ(l,2)=0.036D0
DIJ(l,3)=0.076D0
DIJ(l,4)=0.121D0
DIJ(l,5)=0.129D0
DIJ(l,6)=0.06D0
DIJ(l,7)=0.074D0
DIJ(2,6)=0.106D0
DIJ(2,7)=0.093D0
DIJ(6,7)=0.022D0
DIJ(l,8)=0.089D0
DIJ(2,8)=0.079D0
DIJ(6,8)=0.211D0
DIJ(7,8)=0.089D0
LIJ(l,2)=-0.074D0
LIJ(l,3)=-0.146D0
LIJ(l,4)=-0.258D0
LIJ(l,5)=-0.222D0
LIJ(l,6)=-0.023D0
LIJ(l,7)=-0.086D0
LIJ(6,7)=-0.064D0
LIJ(7,8)=-0.062D0
RETURN
END
SUBROUTINE PARMIX(DIJ,LIJ,TC,VC,PII)
IMPLICIT REAL*8(A-H,O-Z)
REAL*8 LIJ(8,8)
DIMENSION Y(8),DIJ(8,8),VCIJ(8,8),TCIJ(8,8),V13(8),TC(8),VC(8),
*PII(8),PIIJ(8,8)
COMMON/PARCM/TCM,VCM/Y/Y/NC/NC/PCM/PCM/PIM/PIM
DO 1 I=1,NC
1 V13(I)=VC(I)**(1.D0/3.D0)
DO 3 I=1,NC
VCIJ(I,I)=VC(I)
PIIJ(I,I)=PII(I)
TCIJ(I,I)=TC(I)
DO 3 J=1,NC
IF(I.GE.J) GO TO 3
VCIJ(I,J)=(l.D0-LIJ(I,J))*((V13(I)+V13(J))/2.)**3
PIIJ(I,J)=(VC(I)*PII(I)+VC(J)*PII(J))/(VC(I)+VC(J))
TCIJ(I,J)=(l.D0-DIJ(I,J))*(TC(I)*TC(J))**0.5
VCIJ(J,I)=VCIJ(I,J)
PIIJ(J,I)=PIIJ(I,J)
TCIJ(J,I)=TCIJ(I,J)
3 CONTINUE
VCM=0.D0
PIM=0.D0
TCM=0.D0
DO 5 I=1,NC
DO 5 J=1,NC
VCM=VCM+Y(I)*Y(J)*VCIJ(I,J)
PIM=PIM+Y(I)*Y(J)*VCIJ(I,J)*PIIJ(I,J)
5 TCM=TCM+Y(I)*Y(J)*VCIJ(I,J)*TCIJ(I,J)**2
PIM=PIM/VCM
ТСМ=(ТСМ/VСМ)**0.5
PCM=8.31451D-3*(0.28707D0-0.05559*PIM)*TCM/VCM
RETURN
END
SUBROUTINE PHASE
IMPLICIT REAL*8(A-H,O-Z)
COMMON/Z/Z/RM/RM/T/T/P/P/PCM/PCM/RON/RON/BMM/BMM
*/АI/АО,А1,А2,А3
IF(T.LT.240D0.OR.T.GT.480D0.OR.P.LE.0D0.OR.P.GT.12D0)THEN
Z=0D0
GO TO 134
ENDIF
PR=P/PCM
RO=9D3*P/(RM*T*(1.1*PR+0.7D0))
CALL FUN(RO)
CALL OMTAU(RO,T)
IF(Z.EQ.0D0) GO TO 134
Z=l.D0+AO
IF(RON.NE.0D0) THEN
BMM=1D-3*Z*RON*RM*T/P
GO TO 134
ENDIF
NPRIZ=2
CALL COMPL(RO,T,NPRIZ)
CALL TP(RO)
CALL ETAS(RO)
134 RETURN
END
С Подпрограмма, реализующая итерационный процесс определения
С плотности из уравнения состояния (метод Ньютона)
SUBROUTINE FUN(X)
IMPLICIT REAL*8(A-H,O-Z)
COMMON/P/P/RM/RM/Т/Т/AI/AO,A1,A2,A3
ITER=1
1 CONTINUE
NPRIZ=0
IF(ITER.NE.l) NPRIZ=1
CALL COMPL(X,T,NPRIZ)
Z=1.D0+AO
FX=1.D6*(P-(1.D-3*RM*T*Z*X))
F=1.D3*RM*T*(1.D0+A1)
DR=FX/F
X=X+DR
IF(ITER.GT.10) GO TO 4
ITER=ITER+1
IF(DABS(DR/X).GT.l.D-6) GO TO 1
4 CALL COMPL(X,T,NPRIZ)
RETURN
END
SUBROUTINE OMTAU(RO,T)
IMPLICIT REAL*8(A-H,O-Z)
COMMON/PARCM/TCM,VCM/Z/Z
Z=-1D0
TR=T/TCM
ROR=RO*VCM
IF(TR.LT.1.05D0) Z=0D0
IF(ROR.LT.0.D0.OR.ROR.GT.3.D0)Z=0D0
RETURN
END
С Подпрограмма определения безразмерных комплексов АО,А1,А2 и A3
SUBROUTINE COMPL(RO,T,NPRIZ)
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION B(10,8),BK(10)
COMMON/PARCM/TCM,VCM/B/B/AI/AO,A1,A2,A3
IF(NPRIZ.NE.0) GO TO 7
TR=T/TCM
DO 1 I=1,10
BK(I)=0
DO 1 J=l,8
1 BK(I)=BK(I)+B(I,J)/TR**(J-1)
7 ROR=RO*VCM
AO=0.D0
A1=0.D0
IF(NPRIZ.EQ.1) GO TO 5
A2=0.D0
A3=0.D0
5 DO 33 I=1,10
D=BK(I)*ROR**I
AO=AO+D
A1=A1+(I+1)*D
IF(NPRIZ.EQ.1) GO TO 33
DO 3 J=1,8
D1=B(I,J)*ROR**I/TR**(J-1)
A2=A2+(2-J)*D1
3 A3=A3+(J-1)*(2-J)*D1/I
33 CONTINUE
RETURN
END
С Подпрограмма расчета плотности, показателя адиабаты, скорости
С звука
SUBROUTINE TP(ROM)
IMPLICIT REAL*8(A-H,O-Z)
COMMON/BMM/BMM/AI/AO,A1,A2,A3/RM/RM/T/T/TS/RO,PA,W/Z/Z
CALL IDGFU(T,CVOS)
RO=BMM*ROM
R=RM/BMM
A11=1.D0+A1
A21=1.D0+A2
CV=R*(A3+CV0S)
CP=CV+R*A21**2/A11
W=DSQRT(DABS(1.D3*R*T*CP/CV))*DSQRT(DABS(A11))
PA=CP/CV*A11/Z
RETURN
END
С Подпрограмма расчета изохорной теплоемкости в идеально газовом
С состоянии
SUBROUTINE IDGFU(T,CVOS)
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION CPO(8),CVO(8)
COMMON/IDGF/CPC(20,8),TOI(8),MCO(8),MCP(8)/Y/Y(8)/NC/NC
CVOS=0.D0
DO 21 I=1,NC
M=MCP(I)
N=MCO(I)
TAU=T/TOI(I)
S1=0.D0
S2=0.D0
S3=0.D0
S1=CPC(1,I)
IF(M.EQ.0) GO TO 7
DO 9 J=1,M
9 S2=S2+CPC(J+1,I)*TAU**J
7 IF(N.EQ.0)GO TO 11
DO 13 J=1,N
13 S3=S3+CPC(M+J+1,I)/TAU**J
11 CPO(I)=S1+S2+S3
CVO(I)=CPO(I)-l.D0
21 CVOS=CVOS+Y(I)*CVO(I)
RETURN
END
С Подпрограмма расчета вязкости
SUBROUTINE ETAS(ROM)
IMPLICIT REAL*8(A-H,O-Z)
COMMON/ETA/ETA/PARCM/TCM,VCM/BMM/BMM/T/T/PIM/PIM/PCM/PCM
DKSI=TCM**(1D0/6D0)/BMM**.5/PCM**(2D0/3D0)
ROR=VCM*ROM
ТR=Т/TСМ
ETA=78.037D0+3.85612*PIM-29.0053*PIM**2-156.728/TR+145.519/TR**2
*-51.1082/TR**3+6.57895*ROR+(11.7452D0-95.7215*PIM**2/TR)*ROR**2+
*17.1027*ROR**3*PIM+.519623/TR**2*ROR**5
ETA=ETA/DKSI/10.
RETURN
END
BLOCK DATA BDVNIC
IMPLICIT REAL*8(A-H,O-Z)
CHARACTER*26 AR
COMMON/PARCD/VCD(8),TCD(8),PIID(8)/ABIJ/AIJ(10,8),BIJ(10,8)
COMMON/CPCI/CPC1(20,5),CPC2(20,3)/IDGFD/TOID(8),MCOD(8),MCPD(8)
*/AR/AR(25)
DATA TCD/190.67D0,305.57D0,369.96D0,425.4D0,407.96D0,
*125.65D0,304.11D0,373.18D0/
DATA VCD/163.03D0,205.53D0,218.54D0,226.69D0,225.64D0,
*315.36D0,466.74D0,349.37D0/
DATA PIID/0.0006467D0,0.1103D0,0.1764D0,0.2213D0,0.2162D0,
*0.04185D0,0.2203D0,0.042686D0/
DATA AIJ/.6087766D0,-.4596885D0,1.14934D0,-.607501D0,
*-.894094D0,1.144404D0,-.34579D0,-.1235682D0,.1098875D0,
*-.219306D-l,-1.8329l6D0,4.175759D0,-9.404549D0,10.62713D0,
*-3.080591D0,-2.122525D0,1.781466D0,-.4303578D0,-.4963321D-l,
*.347496D-l,1.317145D0,-10.73657D0,23.95808D0,-31.47929D0,
*18.42846D0,-4.092685D0,-.1906595D0,.4015072D0,-.1016264D0,
*-.9129047D-2,-2.837908D0,15.34274D0,-27.71885D0,35.11413D0,
*-23.485D0,7.767802D0,-1.677977D0,.3157961D0,.4008579D-2,0.D0,
*2.606878D0,-11.06722D0,12.79987D0,-12.11554D0,7.580666D0,
*-1.894086D0,4*0.D0,
*-1.15575D0,3.601316D0,-.7326041D0,-1.151685D0,.5403439D0,
*5*0.D0,.9060572D-l,-.5151915D0,.7622076D-l,7*0.D0,
*.4507142D-1,9*0.D0/
DATA BIJ/-.7187864D0,10.67179D0,-25.7687D0,17.13395D0,
*16.17303D0,-24.38953D0,7.156029D0,3.350294D0,-2.806204D0,
*.5728541D0,6.057018D0,-79.47685D0,216.7887D0,-244.732D0,
*78.04753D0,48.70601D0,-41.92715D0,10.00706D0,1.237872D0,
*-.8610273D0,-12.95347D0,220.839D0,-586.4596D0,744.4021D0,
*-447.0704D0,99.6537D0,5.136013D0,-9.5769D0,2.41965D0,
*.2275036D0,15.71955D0,-302.0599D0,684.5968D0,-828.1484D0,
*560.0892D0,-185.9581D0,39.91057D0,-7.567516D0,-.1062596D0,
*0.D0,-13.75957D0,205.541D0,-325.2751D0,284.6518D0,
*-180.8168D0,46.05637D0,4*0.D0,
*6.466081D0,-57.3922D0,36.94793D0,20.77675D0,-12.56783D0,
*5*0.D0,-.9775244D0,2.612338D0,-.4059629D0,7*0.D0,
*-.2298833D0,9*0.D0/
DATA CPCI/1.46696186D+02,-6.56744186D+01,2.02698132D+01,
*-4.20931845D0,6.06743008D-01,-6.l2623969D-02,4.30969226D-03,
*-2.06597572D-04,6.426l581D-06,-1.1680563D-07,9.4095893D-10,
*-2.09233731D+02,2.06925203D+02,-1.35704831D+02,5.64368924D+01,
*-1.34496111D+01,1.39664152D0,3*0.D0,
*6.8120976D+01,-3.0634058D+01,9.5275029D0,-1.6947102D0,
*1.7630585D-01,-9.9545402D-3,2.353643D-4,-8.7407084D+l,
*7.8481374D+I,-4.4865859D+1,1.4654346D+1,-2.0518393D0,8*0.D0,
*-9.209726737D+l,3.070930782D+l,-4.924017995D0,5.045358836D-l,
*-3.140446759D-2,1.076680079D-3,-1.556890669D-5,1.74867128D+2,
*-1.756054503D+2,8.874920732D+l,-1.720610207D+1,9*0.D0,
*-2.096096482D+2,6.877783535D+l,-1.228650555D+l,1.413691547D0,
*-1.002920638D-l,3.985571861D-3,-6.78646087D-5,4.05527285D+2,
*-4.457015773D+2,2.74366735D+2,-8.643867287D+l,1.070428636D+l,
*8*0.D0,
*-3.871419306D+1,4.711104578D+1,-1.758225423D+1,4.183494309D0,
*-5.520042474D-l,3.034658409D-2,2.17160145D+l,-4.4926032D0,
*12*0.D0/
DATA CPC2/0.113129D+2,-0.21596D+1,0.352761D0,-0.321705D-l,
*0.16769D-2,-0.467965D-4,0.542603D-6,-0.174654D+2,0.246205D+2,
*-0.217731D+2,0.116418D+2,-0.342122D+l,0.422296D0,7*0.D0,
*-9.508041394D-1,7.008743711D0.-3.50580167D0, 1.096778D0,
*-2.016835088D-l,1.971024237D-2,-7.860765734D-4,1.087462263D0,
*-7.976765747D-2,-2.837014896D-3,1.479612229D-4,9*0.D0,
*3.91355D0,-6.84851D-2,5.64424D-2,-4.83745D-3,1.71782D-4,
*-2.27537D-6,2*0.D0,1.18658D0,-1.90747D0.8.2852D-1,9*0.D0/
DATA MCOD/6,5,4,5,2,6,4,5/
DATA MCPD/10,6,6,6,5,6,6,5/
DATA TOID/4*100D0,300D0,100D0,300D0,100D0/
DATA AR/' метана (СН4)',' этана (С2Н6)','пропана (С3Н8)',
*' н-бутана (Н-С4Н10)',' и-бутана (и-С4Н10)',' азота (N2)',
*' диоксида углерода (СО2)',' сероводорода (H2S)',
*' ацетилена (С2Н2)',' этилена (С2Н4)',' пропилена (С3Н6)',
*' н-пентана (Н-С5Н12)',' и-пентана (и-С5Н12)',
*' нео-пентана (нео-С5Н12)',' н-гексана (Н-С6Н14)',
*' бензола (С6Н6)',' н-гептана (н-С7Н16)',' толуола (С7Н8)',
*' н-октана (Н-С8Н18)',' н-нонана (н-С9Н20)',
*' н-декана (Н-С10Н22)',' гелия (Не)',' водорода (H2)',
*' моноксида углерода (СО)',' кислорода (O2)'/
END