ГОСТ 30319-2-96 Газ природный Методы расчета физических свойств Определение коэффициента сжимаемости

ГАЗ ПРИРОДНЫЙ

МЕТОДЫ РАСЧЕТА
ФИЗИЧЕСКИХ СВОЙСТВ

ОПРЕДЕЛЕНИЕ КОЭФФИЦИЕНТА СЖИМАЕМОСТИ

МЕЖГОСУДАРСТВЕННЫЙ СОВЕТ
ПО СТАНДАРТИЗАЦИИ, МЕТРОЛОГИИ И СЕРТИФИКАЦИИ

Минск

Предисловие

1 РАЗРАБОТАН Всероссийским научно-исследовательским центром стандартизации, информации и сертификации сырья, материалов и веществ (ВНИЦ СМВ) Госстандарта России; фирмой «Газприборавтоматика» акционерного общества «Газавтоматика» РАО «Газпром»

ВНЕСЕН Госстандартом Российской Федерации

2 ПРИНЯТ Межгосударственным Советом по стандартизации, метрологии и сертификации (протокол № 9-96 от 12 апреля 1996 г.)

За принятие проголосовали:

3 ПОСТАНОВЛЕНИЕМ Государственного комитета Российской Федерации по стандартизации, метрологии и сертификации от 30 декабря 1996 г. № 723 межгосударственный стандарт ГОСТ 30319.2-96 введен в действие непосредственно в качестве государственного стандарта Российской Федерации с 1 июля 1997 г.

4 ВВЕДЕН ВПЕРВЫЕ

5 ПЕРЕИЗДАНИЕ

СОДЕРЖАНИЕ

Газ природный

МЕТОДЫ РАСЧЕТА ФИЗИЧЕСКИХ СВОЙСТВ

Определение коэффициента сжимаемости

Дата введения 1997-07-01

1 Назначение и область применения

Настоящий стандарт устанавливает четыре метода определения коэффициента сжимаемости природного газа: при неизвестном полном компонентном составе природного газа (два метода) и известном компонентном составе.

Стандарт устанавливает предпочтительные области применения каждого метода по измеряемым параметрам (давление, температура, плотность природного газа при стандартных условиях и компонентный состав природного газа), однако не запрещает использование любого из методов и в других областях.

Допускается применять любые другие методы расчета коэффициента сжимаемости, однако погрешность расчета коэффициента сжимаемости по этим методам не должна превышать погрешностей, приведенных в настоящем стандарте (см. ).

Используемые в настоящем стандарте определения и обозначения приведены в соответствующих разделах ГОСТ 30319.0.

2 Нормативные ссылки

В настоящем стандарте использованы ссылки на следующие стандарты:

ГОСТ 30319.0-96 Газ природный. Методы расчета физических свойств. Общие положения

ГОСТ 30319.1-96 Газ природный. Методы расчета физических свойств. Определение физических свойств природного газа, его компонентов и продуктов его переработки

3 Определение коэффициента сжимаемости

3.1 Общие положения

Коэффициент сжимаемости вычисляют по формуле

К = z/z_c,                                                         (1)

где    z и z_c — фактор сжимаемости соответственно при рабочих и стандартных условиях.

Рабочие условия характеризуются такими давлениями и температурами, которые определяются измерениями в процессе добычи, переработки и транспортирования природного газа. Давление p_c и температура T_c при стандартных условиях приведены в ГОСТ 30319.0.

3.2 Методы расчета коэффициента сжимаемости

3.2.1 Пределы применимости методов расчета и погрешности расчета коэффициента сжимаемости

В приведены общие результаты апробации методов расчета и область их применения. Апробация проведена на обширном массиве высокоточных экспериментальных данных о факторе сжимаемости природного газа [].

Погрешность данных не превышает 0,1 %.

Таблица 1 — Результаты апробации и область применения методов расчета коэффициента сжимаемости природного газа

 

Для расчета коэффициента сжимаемости природного газа при определении его расхода и количества рекомендуется применять:

1) модифицированный метод NX19 мод.- при распределении газа потребителям;

2) модифицированное уравнение состояния (УС) GERG-91 мод. [, ] и УС AGA8-92DC [] — при транспортировании газа по магистральным газопроводам;

3) уравнение состояния ВНИЦСМВ – при добыче и переработке газа». Погрешность расчета коэффициента сжимаемости d приведена в таблице 1 без учета погрешности исходных данных.

Метод NX19 мод. и уравнение состояния GERG-91 мод. могут быть использованы при неизвестном полном компонентном составе природного газа, расчет по этим методам не требует применения ЭВМ.

Расчет по уравнениям состояния AGA8-92DC и ВНИЦ СМВ может быть осуществлен только при наличии ЭВМ и известном полном компонентном составе природного газа, при этом должны быть выдержаны следующие диапазоны концентраций компонентов (в мол. %):

метан                65 — 100              этан                           £ 15

пропан                    £ 3,5             бутаны                      £ 1,5

азот                         £ 15              диоксид углерода    £ 15

сероводород           £ 30              (УС ВНИЦ СМВ) и £ 0,02 (УС AGA8-92DC)

остальные               £ 1

В области давлений (12 — 30) МПа и температур (260 — 340) К для расчета коэффициента сжимаемости допускается применять уравнения состояния GERG-91 мод. и AGA8-92DC. Погрешность расчета коэффициента сжимаемости природного газа в указанной области давлений и температур составляет: для уравнения GERG-91 мод. — 3,0 % [], для уравнения AGA8-92DC — 0,5 % [].

Выбор конкретного метода расчета коэффициента сжимаемости допускается определять в контракте между потребителем природного газа и его поставщиком с учетом требований настоящего стандарта.

В приняты следующие обозначения:

1) d_сист — систематическое отклонение от экспериментальных данных

,                                                          (2)

2) d_i^макс — максимальное отклонение в i-й точке экспериментальных данных

,                                           (3)

где    К_расч и К_эксп — соответственно расчетный и экспериментальный коэффициенты сжимаемости;

3) d — погрешность расчета коэффициента сжимаемости по ИСО 5168 []

,                                             (4)

где    d_ст — стандартное отклонение, которое вычисляется из выражения

,                                           (5)

d_эксп — погрешность экспериментальных данных (0,1 %).

(Измененная редакция, Изм. № 1).

3.2.2 Модифицированный метод NX19 мод.

В соответствии с требованиями стандарта Германии [] расчет фактора сжимаемости по модифицированному методу NX19 мод. основан на использовании уравнения следующего вида

,                                                 (6)

где ,                                          (7)

,                                         (8)

,                                                       (9)

,                           (10)

,                            (11)

Корректирующий множитель F в зависимости от интервалов параметров р_а и DТ_а вычисляют по формулам:

при 0 £ р_а £ 2 и 0 £ DТ_а £ 0,3

,                (12)

при 0 £ р_а< 1,3 и -0,25 £ DТ_а < 0

,                              (13)

при 1,3 £ р_а< 2 и -0,21 £ DТ_а < 0

,       (14)

где    DT_a = T_a — 1,09.

Параметры р_a и Т_a определяются по следующим соотношениям:

,                                                  (15)

,                                                    (16)

где    р_пк и Т_пк — псевдокритические значения давления и температуры, определяемые по формулам (48) и (49) ГОСТ 30319.1, а именно:

,                               (17)

.                                  (18)

В ), ) вместо молярных долей диоксида углерода и азота допускается применять их объемные доли (r_y и r_a).

Коэффициент сжимаемости природного газа вычисляют по ), при этом фактор сжимаемости при рабочих условиях рассчитывают по ) – ) настоящего стандарта, а фактор сжимаемости при стандартных условиях – по формуле (24) ГОСТ 30319.1

(Измененная редакция, Изм. № 1).

3.2.3 Модифицированное уравнение состояния GERG-91 мод.

Европейская группа газовых исследований на базе экспериментальных данных, собранных в [], и уравнения состояния вириального типа [], разработала и опубликовала в [, ] УС

,                                                    (19)

где    В_m и С_m — коэффициенты УС;

r_м — молярная плотность, кмоль/м³.

Коэффициенты уравнения состояния определяют из следующих выражений:

,              (20)

,    (21)

где х_э — молярная доля эквивалентного углеводорода

х_э = 1 — х_а — х_у,                                                           (22)

,                                                                                                                                                (23)

,                             (24)

,                       (25)

,                              (26)

, (27)

,                          (28)

,                          (29)

,                     (30)

,                    (31)

,                                          (32)

.                                              (33)

В ), Н рассчитывают по выражению

,                                               (34)

где    М_э — молярная масса эквивалентного углеводорода, значение которой определяется из выражения

,                         (35)

В ) молярную долю эквивалентного углеводорода (x_э) рассчитывают с использованием ), а фактор сжимаемости при стандартных условиях (z_с) рассчитывают по формуле (24) ГОСТ 30319.1, а именно

,                          (36)

После определения коэффициентов уравнения состояния В_m и С_m рассчитывают фактор сжимаемости при заданных давлении (р, МПа) и температуре (Т, К) по формуле

,                                             (37)

где

,                                              (38)

,                                                    (39)

,                                                           (40)

,                                                             (41)

С₀ = b²C_m,                                                             (42)

,                                                    (43)

Коэффициент сжимаемости природного газа рассчитывают по ), а именно

,                                                            (44)

Фактор сжимаемости при стандартных условиях z_с рассчитывают по ).

(Измененная редакция, Изм. № 1).

3.2.4 Уравнение состояния AGA8-92DC

В проекте стандарта ISO/TC 193 SC1 № 62 [] Американской Газовой Ассоциацией для расчета фактора сжимаемости предложено использовать уравнение состояния

,                       (45)

где    В и С_n^* — коэффициенты УС;

r_м — молярная плотность, кмоль/м³.

Константы {b_n, c_n, k_n} УС приведены в .

Если состав газа задан в объемных долях, то молярные доли рассчитываются по формуле (12) ГОСТ 30319.1.

Приведенную плотность определяют по формуле

,                                                            (46)

Параметр К_т вычисляют по ).

Коэффициенты УС рассчитывают из следующих соотношений:

,                      (47)

,                      (48)

где    N — количество компонентов в природном газе.

Константы {a_n, u_n, g_n, q_n, f_n} и характерные параметры компонентов {Е_i, К_i, G_i, Q_i, F_i} в ), ) приведены соответственно в и .

Бинарные параметры {E_ij, G_ij} и параметры {U, G, K_m, Q, F} рассчитывают с использованием следующих уравнений:

,                                                  (49)

(i ¹ j)

,                                               (50)

(i ¹ j)

 ,                             (51)

,                                   (52)

,                             (53)

,                                                         (54)

,                                                          (55)

где    {E_ij^*, G_ij^*, U_ij^*, K_ij^*} — параметры бинарного взаимодействия, которые даны в .

Параметры бинарного взаимодействия, которые не приведены в этой таблице, а также при i=j, равны единице.

Для расчета фактора сжимаемости по уравнению состояния необходимо определить плотность r_м при заданных давлении (р, МПа) и температуре (Т, К).

Плотность r_м из УС определяют по методу Ньютона в следующем итерационном процессе:

1) начальную плотность определяют по формуле

,                                          (56)

где приведенное давление вычисляют из выражения

,                                                                 (57)

2) плотность на k-м итерационном шаге определяют из выражений

.                                 (58)

,                                                      (59)

где z^(k-1) — рассчитывают из УС ) при плотности на итерационном шаге (k-1), т.е. при r_м^(k-1), а безразмерный комплекс А₁ определяют из выражения

,(60)

при этом r_п = К_т³r_м^(k-1);

4) критерий завершения итерационного процесса

,                                                    (61)

если критерий ) не выполняется, то необходимо продолжить итерационный процесс, начиная с пункта 2) алгоритма.

После определения фактора сжимаемости при рабочих и стандартных условиях по ) рассчитывают коэффициент сжимаемости.

(Измененная редакция, Изм. № 1).

3.2.5 Уравнение состояния ВНИЦ СМВ

Во Всероссийском научно-исследовательском центре стандартизации, информации и сертификации сырья, материалов и веществ (ВНИЦ СМВ) для расчета фактора сжимаемости природного газа разработано уравнение состояния

,                                                  (62)

где c_kl — коэффициенты УС;

r_п = r_м/r_пк — приведенная плотность;

Т_п = Т/Т_пк — приведенная температура;

r_м — молярная плотность, кмоль/м³;

r_пк и Т_пк — псевдокритические параметры природного газа.

Коэффициенты УС определяют по формуле

,                                                   (63)

где {a_kl, b_kl} — обобщенные коэффициенты УС, которые приведены в .

Псевдокритические параметры природного газа и его фактор Питцера вычисляют по формулам:

— псевдокритическую плотность

,                                           (64)

где ,                     (65)

(; )                                                 

— псевдокритическую температуру

,                                               (66)

где ,                                     (67)

;                                      (68)

(; )                                             

— фактор Питцера

,                                     (69)

где ,              (70)

В N — число основных компонентов природного газа (метана, этана, пропана, н-бутана, и-бутана, азота, диоксида углерода, сероводорода).

Критические параметры компонентов {r_кi, r_кj, Т_кj, Т_кj}, их молярная масса {М_i, М_j,} и факторы Питцера {W_i, W_j} приведены в , а параметры бинарного взаимодействия {x_ij, l_ij} — в и .

Если заданный компонентный состав природного газа включает, кроме основных, другие компоненты (но не более 1 % в сумме), то молярные доли этих компонентов прибавляют к соответствующим долям основных компонентов следующим образом:

— ацетилен и этилен к этану;

— пропилен к пропану;

— углеводороды от н-пентана и выше к н-бутану;

— прочие компоненты к азоту.

Если состав газа задан в объемных долях, то молярные доли рассчитывают по формуле (12) ГОСТ 30319.1.

Для расчета фактора сжимаемости по уравнению состояния ) необходимо определить плотность r_м при заданных давлении (р, МПа) и температуре (Т, К).

Плотность r_м из УС ) определяют по методу Ньютона в следующем итерационном процессе:

1) начальную плотность определяют по формуле

,                                          (75)

где приведенное давление вычисляют из выражений

,                                    (76)

,                                                             (77)

а псевдокритические плотность (r_пк), температуру (Т_пк) и фактор Питцера (W) рассчитывают по ), ) и );

2) плотность на k-м итерационном шаге определяется из выражений

,                                  (78)

,                                                     (79)

где    z^(k-1) рассчитывают из УС ) при плотности на итерационном шаге (k-1), т.е. при r_м^(k-1), a безразмерный комплекс A₁ определяют из выражения

,                                              (80)

4) критерий завершения итерационного процесса.

,                                                        (81)

если критерий ) не выполняется, то необходимо продолжить итерационный процесс, начиная с пункта 2) алгоритма.

После определения фактора сжимаемости при рабочих и стандартных условиях по ) рассчитывают коэффициент сжимаемости.

(Измененная редакция, Изм. № 1).

4 Влияние погрешности исходных данных на погрешность расчета коэффициента сжимаемости

При измерении расхода и количества природного газа, транспортируемого в газопроводах, давление (р), температуру (T), плотность при стандартных условиях (r_c) и состав (х_i) измеряют с определенной погрешностью. Перечисленные параметры являются исходными данными для расчета коэффициента сжимаемости.

В соответствии с рекомендациями ИСО 5168 [] погрешность расчета коэффициента сжимаемости, которая появляется в связи с погрешностью измерения исходных данных, определяют по формуле

,                                (82)

где    d_ид — погрешность расчета коэффициента сжимаемости, связанная с погрешностью измерения исходных данных;

d_qk — погрешность измерения параметра исходных данных;

,                                              (83)

,                                                  (84)

В :

q_k — условное обозначение k-го параметра исходных данных (р. Т, r_c, х_i,);

`q_k — среднее значение k-го параметра в определенный промежуток времени (сутки, месяц, год и т.д.);

q_k^макс и q_k^мин — максимальное и минимальное значения k-го параметра в определенный промежуток времени;

N_q — количество параметров исходных данных.

При вычислении частных производных по коэффициенты сжимаемости К_qk⁺ и К_qk^— – рассчитывают при средних параметрах и параметрах и соответственно. Рекомендуется выбирать

Коэффициент сжимаемости `К (среднее значение) рассчитывают по выбранному рекомендуемому методу расчета при средних параметрах q_k.

Для методов:

1) NX 19 мод. и УС GERG-91 мод. — N_q = 5 и параметрами исходных данных являются давление, температура, плотность при стандартных условиях, молярные доли азота и диоксида углерода;

2) УС AGA8-92DC и УС ВНИЦ СМВ — N_q = 2 + N (N — количество компонентов) и параметрами исходных данных являются давление, температура и молярные доли компонентов природного газа, причем для УС ВНИЦ СМВ учитываются молярные доли только основных компонентов газа.

Общую погрешность расчета коэффициента сжимаемости определяют по формуле

,                                                        (85)

где d — погрешность расчета коэффициента сжимаемости, которая для каждого метода приведена в .

Для методов NX19 мод. и УС GERG-91 мод. допускается рассчитывать погрешность d_и.д по формуле

, (86)

где    d_Т, d_p, d_rc, d_xa и d_xy — погрешности измеряемых параметров, соответственно, температуры, давления, плотности природного газа при стандартных условиях, содержания азота и диоксида углерода в нем.

Коэффициенты К_T, К_р, К_rc, К_xa и К_xу в зависимости от метода, используемого для расчета коэффициента сжимаемости K, определяются по следующим выражениям (см. формулы (34) — (38) или (39) — (43) ГОСТ 30319.1):

— при расчете К по методу NX19 мод.

,                                                 (87)

,                                                    (88)

,                                                       (89)

,                                                      (90)

,                                                       (91)

— при расчете К по методу GERG-91

,                                                  (92)

,                                                   (93)

,                                                              (94)

,                                                      (95)

.                                                      (96)

(Измененная редакция, Изм. № 1).

5 Программная и техническая реализация расчета коэффициента сжимаемости

Расчет коэффициента сжимаемости природного газа по указанным в стандарте методам реализован на ПЭВМ, совместимых с IBM PC/AT/XT, на языке программирования ФОРТРАН-77. Листинг программы приведен в .

В и приведены примеры расчета соответственно коэффициента сжимаемости и погрешности вычисления коэффициента сжимаемости, которая вызвана погрешностью определения исходных данных.

ПРИЛОЖЕНИЕ А

(обязательное)

Таблицы констант и параметров уравнения состояния AGA8-92DC

Таблица А.1 — Константы уравнения состояния AGA8-92DC

Таблица А.2 — Характерные параметры компонентов

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

ПРИЛОЖЕНИЕ Б

(обязательное)

Таблицы коэффициентов и параметров уравнения состояния ВНИЦ СМВ

Таблица Б.1 — Обобщенные коэффициенты уравнения состояния ВНИЦ СМВ

Таблица Б.2 — Физические свойства компонентов природного газа, используемые в уравнении состояния ВНИЦ СМВ

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

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

ПРИЛОЖЕНИЕ В

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

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

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)*

*EIJM(I,I)"UN(N)*KI(I)*KI(I)*KI(I)

DO 11 I=1,NC-1

DO 11 J=I+1,NC

11     В1(N)=B1(N)+2.*Y(I)*Y(J)(GIJМ(I,J)+1D0-GN(N))**GN(N)*

*(QI(I)*QI(J)+1D0-QN(N))**QN(N)*((FI(I)*FI(J))**.5+

1D0-FN(N))**FN(N)*EIJM(I,J)**UN(N)*(KI(I)*KI(J)**1.5

7       CONTINUE

DO 13 N=8,53

13     C1(N)=AN(N)*(GM+1D0-GN(N))**GN(N)*(QM**2+1D0-QN(N))**

*QN(N)*(FM+1D0-FN(N))**FN(N)*UM**UN(N)

RETURN

END

SUBROUTINE PHASE1

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

COMMON/Z/Z/RM/RM/T/T/P/P/AI1/AO,A1/AN/AN(53)

*/COEF1/B1(13),C1(53)/COEF2/B,C(53)/UN/UN(53)

CALL PCONT1(P,T)

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

B=0D0

DO 1 N=1,13

1       B=B+AN(N)/T**UN(N)*B1(N)

DO 3 N=8,53

3       C(N)=C1(N)/T**UN(N)

PR=P/5.

RO=9D3*P/(RM*T*(1.1*PR+0.7D0))

CALL FUN1(RO)

Z=1D0+AO

134   RETURN

END

С      Подпрограмма, реализующая итерационный процесс определения

С      плотности из уравнения состояния (метод Ньютона)

SUBROUTINE FUN1(X)

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

COMMON/P/P/RM/RM/T/T/AI1/AO,A1

ITER=1

1       CONTINUE

CALL COMPL1(X)

Z=1.D0+AO

FX= 1 .D6*(P-(1.D-3*RM*T*Z*X))

F= 1 .D3*RM*(1.D0+A1)

DR=FX/F

X=X+DR

IF(ITER.GT.10) GO TO 4

ITER=ITER+1

IF(DABS(DR/X).GT.1.D-6) GO TO 1

4       CALL COMPL1(X)

RETURN

END

SUBROUTINE PCONT1(P,T)

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

COMMON/Z/Z

Z=-1D0

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

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

RETURN

END

SUBROUTINE COMPL1(RO)

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

REAL*8 KM

INTEGER BN,CN

COMMON/KM/KM/COEF2/B,C(53)/BCKN/BN(53),CN(53),KN(53)/AI1/AO,A1

ROR=KM*RO

S1=0D0

S2=0D0

S3=0D0

DO 1 N=8,53

EXP=DEXP(-CN(N)*ROR**KN(N))

IF(N.LE.13) S1=S1+C(N)

S2=S2+C(N)*(BN(N)-CN(N)*KN(N)*ROR**KN(N))*ROR**BN(N)*EXP 1

S3=S3+C(N)*(-CN(N)*KN(N)**2*KM*ROR**(KN(N)-1)*ROR**BN(N)*

*EXP+(BN(N)-CN(N)*KN(N)*ROR**KN(N))*BN(N)*KM*ROR**(BN(N)-1)*

*EXP-(BN(N)-CN(N)*KN(N)*ROR**KN(N))*ROR**BN(N)*EXP*CN(N)*KN(N)*

*KM*ROR**(KN(N)-1))AO1=B*RO-ROR*S1

AO=AO1+S2

A1=AO+AO1+RO*S3

RETURN

END

BLOCK DATA DCAGA8

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

REAL*8 KD

INTEGER BN,CN,GN,QN,FN

COMMON/EFD/ED(19),KD(19),GD(19),QD(19),FD(19)

*/BCKN/BN(53),CN(53),KN(53)/UN/UN(53)

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

DATA ED/151.3183D0,244.1667D0,298.1183D0,337.6389D0,324.0689D0,

*99.73778D0,241.9606D0,296.355D0,370.6823D0,365.5999D0,

*402.8429D0,427.5391D0,450.6472D0,472.1194D0:488.7633D0,

*2.610111D0,26.95794D0,105.5348D0,122.7667D0/

DATA KD/.4619255D0,.5279209D0,.583749D0,.6341423D0,.6406937D0,

*.4479153D0,.4557489D0,.4618263D0,.6798307D0,.6738577D0,

*.7139987D0,.7503628D0,.7851933D0,.8157596D0,.8389542D0,

*.3589888D0,.3514916D0,.4533894D0,.4186954D0/

DATA GD/0D0,.0793D0,.141239D0,.281835D0,.256692D0,.027815D0,

*.189065D0,.0885D0,.366911D0,.332267D0,.432254D0,.512507D0,

*.576242D0,.648601D0,.716574D0,0D0,.034369D0,.038953D0,.021D0/

DATA QD/6*0D0,.69D0,12*0D0/,FD/16*0D0,1D0,2*0D0/

DATA AN/.1538326D0,1.341953D0,-2.998583D0,-.04831228D0,

*.3757965D0,-1.589575D0,-.05358847D0,2.29129D-9,1576724D0,

*-.4363864D0,-.04408159D0,-.003433888D0,.03205905D0,.02487355D0, —

*.07332279D0,-.001600573D0,.6424706D0,-.4162601D0,-.06689957D0,

*.2791795D0,-.6966051D0,-.002860589D0,-.008098836D0,3.150547D0,

*.007224479D0,-.7057529D0,.5349792D0,-.07931491D0-1.418465D0,

*-5.99905D-17,.1058402D0,.03431729D0,-.007022847D0,.02495587D0,

*.04296818D0,.7465453D0,-.2919613D0,7.294616D0,-9.936757D0,

*-.005399808D0,-.2432567D0,.04987016D0,.003733797D0,1.874951D0,

*.002168144D0,-.6587164D0,.000205518D0,.009776195D0,-.02048708D0,

*.01557322D0,.006862415D0,-.001226752D0,.002850906D0/

DATA BN/13*1,9*2,10*3,7*4,5*5,2*6,2*7,3*8,2*9/

DATA CN/7*0,6*1,2*0,7*1,0,9*1,2*0,5*1,0,4*1,0,1,0,6*1/

DATA KN/7*0,3,3*2,2*4,2*0,3*2,4*4,0,2*1,2*2,2*3,3*4,2*0,3*2,

*2*4,0,2*2,2*4,0,2,0,2,1,4*2/

DATA UN/0D0,.5D0,1D0,3.5D0,-.5D0,4.5D0,.5D0,-6D0,2D0,3D0,2*2D0,

*11D0,-.5D0,.5D0,0D0,4D0,6D0,21D0,23D0,22D0,-1D0,-.5D0,7D0,-1D0,

*6D0,4D0,1D0,9D0,-13D0,21D0,8D0,-.5D0,0D0,2D0,7D0,9D0,22D0,23D0,

*1D0,9D0,3D0,8D0,23D0,1.5D0,5D0,-.5D0,4D0,7D0,3D0,0D0,1D0,0D0/

DATA GN/4*0,2*1,13*0,1,3*0,1,2*0,3*1,16*0,1,2*0,1,0,1,2*0/

DATA QN/6*0,1,3*0,1,9*0,1,0,1,8*0,1,4*0,1,4*0,1,0,1,2*0,1,5*0,1/

DATA FN/7*0,1,13*0,1,2*0,1,4*0,1,23*0/

END

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

С      *                                                                                                                  *

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

С      * газа по уравнению состояния ВНИЦ СМВ.                                        *

С      *                                                                                                                  *

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

SUBROUTINE VNIC(ICALC)

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

REAL*8 LIJ(8,8)

DIMENSION VC(8),TC(8),PII(8),DIJ(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/Z/Z

RM=8.31451DO

IF(ICALC.NE.1) GO TO 1

CALL COMPON

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

CALL DDIJ(DIJ,LIJ)

DO 75 I=1,NC

TC(I)=TCD(NI(I))

VC(I)=BM(I)/VCD(NI(I))

PII(I)=PIID(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,PII,PIM)

DO 27 I=1,10

DO 27 J=l,8

27     B(I,J)=AIJ(I,J)+BIJ(I,J)*PIM

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/NPR/NPR

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/

DATAR0I/0.6682D0,1.2601D0,1.8641D0,2.4956D0,2.488D0,

*1.1649D0,1.8393D0,1.4311D0/

DO 100 I=1,25

100   YI(I)=YC(I)

IF(NPR.EQ.1) GO TO 333

BMM=0D0

DO 3333 I=1,25

3333 ВММ=ВММ+YI(I)*ВМI(I)

333   YS=0D0

DO 55 I=9,25

YS=YS+YI(I)

55     CONTINUE

YS1=0D0

DO 67 I=12,21

67     YS1=YS1+YI(I)

YS2=0D0

DO 69 1=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(NPR.EQ.1.AND.YI(5).LT.0.01D0.AND.YS3.LT.0.03D0) YI(4)=YS3

IF(NPR.EQ.1.AND.YI(5).LT.0.01D0.AND.YS3.LT.0.03D0) YI(5)=0D0

IF(NPR.EQ.0.AND.Y1(5).LT.0.01D0.AND.YS3.LE.0.03D0) YI(4)=YS3

IF(NPR.EQ.0.AND.YI(5).LT.0.01D0.AND.YS3.LE.0.03D0)YI(5)=0D0

YI(6)=YI(6)+YS2

IF(NPR.EQ.0) 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 1=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(NC)=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(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.YI(8).GT.0.3D0)Z=0D0

RETURN

END

SUBROUTINE DDIJ(DIJ,LIJ)

IMPLICIT REAL-8(A-H,O-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(1,2)=0.036D0

DIJ(1,3)=0.076D0

DIJ(1,4)=0.121D0

DIJ(1,5)=0.129D0

DIJ(1,6)=0.06D0

DIJ(1,7)=0.074D0

DIJ(2,6)=0.106D0

DIJ(2,7)=0.093D0

DIJ(6,7)=0.022D0

DIJ(1,8)=0.089D0

DIJ(2,8)=0.079D0

DU(6,8)=0.211D0

DIJ(7,8)=0.089D0

LIJ(1,2)=-0.074D0

LIJ(1,3)=-0.146D0

LIJ(1,4)=-0.258D0

LIJ(1,5)=-0.222D0

LIJ(1,6)=-0.023D0

LIJ(1,7)=-0.086D0

LIJ(6,7)=-0.064D0

LIJ(7,8)=-0.062D0

RETURN

END

SUBROUTINE PARMIX(DIJ,LIJ,TC,VC,PII,PIM)

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

DO 1 I=1,NC

1       V13(I)=VC(I)**(1.DO/3.DO)

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)=(1.DO-LIJ(I,J))*((V13(I)+V13(J))/2.)**3

PIIJ(I,J)=(VC(I)*PII(I)+VC(J)*PII(J))/(VC(I)+VC(J))

TCU(I,J)=(1.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

TCM=(TCM/VCM)**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/AI/AO,A1

IF(T.LT.250D0.OR.T.GT.340D0.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=1.D0+AO

134   RETURN

END

С      Подпрограмма, реализующая итерационный процесс определения

С      плотности из уравнения состояния (метод Ньютона)

SUBROUTINE FUN(X)

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

COMMON/P/P/RM/RM/T/T/AI/AO,A1

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.1.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

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

IF(NPRIZ.NE.0) GO TO 7

TR=T/TCM

DO 1 I=1,10

BK(I)=0

DO 1 J=1,8

1       BK(I)=BK(I)+B(I,J)/TR**(J-1)

7       ROR=RO*VCM

AO=0.D0

A1=0.D0

DO 33 I=1,10

D=BK(I)*ROR**I

AO=AO+D

33     A1=A1+(I+1)*D

RETURN

END

BLOCK DATA BDVNIC

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

COMMON/PARCD/VCD(8),TCD(8),PIID(8)/ABIJ/AIJ(10,8),BIJ(10,8)

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-1,-1.832916D0,4.175759D0,-9.404549D0,10.62713D0,

*-3.080591D0,-2.122525D0,1.781466D0,-.4303578D0,-.4963321D-1,

*.347496D-1,1.317145D0,-10.73657D0,23.95808 D0,-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-1,-.5151915D0,.7622076D-1,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/

END

ПРИЛОЖЕНИЕ Г

(обязательное)

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

Г.1 Модифицированный метод NX19

Плотность при 0,101325 МПа и 293,15 К: 0,6799 кг/м³

Содержание:

азота ………………………………………………………………………………………….. 0,8858 мол. %

диоксида углерода ……………………………………………………………………… 0,0668 мол. %

Давление ……………………………………………………………………………………. 2,001 МПа

Температура ………………………………………………………………………………. 270,00 К

Коэффициент сжимаемости ……………………………………………………….. 0,9520

Давление ……………………………………………………………………………………. 2,494 МПа

Температура ………………………………………………………………………………. 280,00 К

Коэффициент сжимаемости ……………………………………………………….. 0,9473

Давление ……………………………………………………………………………………. 0,900 МПа

Температура ………………………………………………………………………………. 290,00 К

Коэффициент сжимаемости ……………………………………………………….. 0,9844

Г.2 Уравнение состояния GERG-91

Плотность при 0,101325 МПа и 293,15 К: 0,6799 кг/м³

Содержание:

азота ………………………………………………………………………………………….. 0,8858 мол. %

диоксида углерода ……………………………………………………………………… 0,0668 мол. %

Давление ……………………………………………………………………………………. 2,001 МПа

Температура ………………………………………………………………………………. 270,00 К

Коэффициент сжимаемости ……………………………………………………….. 0,9521

Давление ……………………………………………………………………………………. 3,997 МПа

Температура ………………………………………………………………………………. 290,00 К

Коэффициент сжимаемости ……………………………………………………….. 0,9262

Давление ……………………………………………………………………………………. 7,503 МПа

Температура ………………………………………………………………………………. 330,00 К

Коэффициент сжимаемости ……………………………………………………….. 0,9244

Г.3 Уравнение состояния AGA8-92DC

Состав природного газа в молярных процентах:

метан …………………………………………………………………………………………. 98,2722

этан …………………………………………………………………………………………… 0,5159

пропан ……………………………………………………………………………………….. 0,1607

н-бутан ………………………………………………………………………………………. 0,0592

азот ……………………………………………………………………………………………. 0,8858

диоксид углерода ……………………………………………………………………….. 0,0668

н-пентан …………………………………………………………………………………….. 0,0157

н-гексан …………………………………………………………………………………….. 0,0055

н-гептан …………………………………………………………………………………….. 0,0016

н-октан ………………………………………………………………………………………. 0,0009

гелий …………………………………………………………………………………………. 0,0157

Плотность при 0,101325 МПа и 293,15 К: 0,6799 кг/м³

Давление ……………………………………………………………………………………. 2,001 МПа

Температура ………………………………………………………………………………. 270,00 К

Коэффициент сжимаемости ……………………………………………………….. 0,9520

Давление ……………………………………………………………………………………. 3,997 МПа

Температура ………………………………………………………………………………. 290,00 К

Коэффициент сжимаемости ……………………………………………………….. 0,9262

Давление ……………………………………………………………………………………. 7,503 МПа

Температура ………………………………………………………………………………. 330,00 К

Коэффициент сжимаемости ……………………………………………………….. 0,9246

Г.4 Уравнение состояния ВНИЦ СМВ

Состав природного газа в молярных процентах:

метан …………………………………………………………………………………………. 89,2700

этан …………………………………………………………………………………………… 2,2600

пропан ……………………………………………………………………………………….. 1,0600

и-бутан ………………………………………………………………………………………. 0,0100

азот ……………………………………………………………………………………………. 0,0400

диоксид углерода ……………………………………………………………………….. 4,3000

сероводород ……………………………………………………………………………….. 3,0500

пропилен …………………………………………………………………………………… 0,0100

Плотность при 0,101325 МПа и 293,15 К: 0,7675 кг/м3

Давление ……………………………………………………………………………………. 1,081 МПа

Температура ………………………………………………………………………………. 323,15 К

Коэффициент сжимаемости ……………………………………………………….. 0,9853

Давление ……………………………………………………………………………………. 4,869 МПа

Температура ………………………………………………………………………………. 323,15 К

Коэффициент сжимаемости ……………………………………………………….. 0,9302

Давление ……………………………………………………………………………………. 9,950 МПа

Температура ………………………………………………………………………………. 323,15 К

Коэффициент сжимаемости ……………………………………………………….. 0,8709

ПРИЛОЖЕНИЕ Д

(обязательное)

Влияние погрешности исходных данных на погрешность расчета коэффициента сжимаемости природного газа (примеры расчета)

Д.1 Модифицированный метод NX19

Коэффициент сжимаемости (среднее значение) — 0,9520

Погрешность расчета: по ) — 0,09 %; по ) — 0,07 %.

Д.2 Уравнение состояния GERG-91

Коэффициент сжимаемости (среднее значение) — 0,9521

Погрешность расчета: по ) — 0,09 %; по ) — 0,09 %.

Д.3 Уравнение состояния AGA8-92DC

Коэффициент сжимаемости (среднее значение) — 0,9520

Погрешность расчета — 0,08 %

Д.4 Уравнение состояния ВНИЦ СМВ

Коэффициент сжимаемости (среднее значение) — 0,9853

Погрешность расчета — 0,03 %

ПРИЛОЖЕНИЕ Е

(справочное)

Библиография

[1] Сычев В.В. и др. Термодинамические свойства метана. — М., Изд-во стандартов, 1979, 348 с

[2] Kleinrahm R., Duschek W., Wagner W. Measurement and correlation of the (pressure, density, temperature) relation of methane in the temperature range from 273.15 К to 323.15 К at pressures up to 8 MPa. — J. Chem. Thermodynamics, 1988, v.20, p.621-631

[3] Robinson R.L., Jacoby R.H. Better compressibility factors. — Hydrocarbon Processing, 1965,v.44,No.4,p.141-145

[4] Achtermann H.-J., Klobasa F.,Rogener H. Realgasfaktoren von Erdgasen. Teil I: Bestimmung von Realgasfaktoren aus Brechungsindex-Messungen. — Brennstoff-Warme-Kraft, 1982, Bd.34, No.5, s.266-271

[5] Achtermann H.-J., Klobasa F.,Rogener H. Realgasfaktoren von Erdgasen. Teil II: Bestimmung von Realgasfaktoren mit eener Burnett-Apparatur. — Brennstoff-Warme-Kraft, 1982, Bd.34, No.6, s.311-314

[6] Eubank Ph.T., Scheloske J., Hall K.R., Holste J.C. Densities and mixture virial coefficients for wet natural gas mixtures. — Journal of Chemical and Engineering Data, 1987, v.32, No.2, p.230-233

[7] Jaeschke М., Julicher H.P. Realgasfaktoren von Erdgasen. Bestimmung von Realgasfaktoren nach der Expansionsmethode. — Brennstoff-Warme-Kraft, 1984, Bd.36, No.11, s.445-451

[8] Jaeschke М. Realgasverhalten Einheitliche Berechnungsmoglichkeiten von Erdgas L und H. — Gas und Wasserfach. Gas/Erdgas, 1988, v.129, No.l, s.30-37

[9] Blanke W., Weiss R. pvT-Eigenschaften und Adsorptions- verhalten von Erdgas bei Temperaturen zwischen 260 К und 330 К mit Drucken bis 3 MPa. — Erdol-Erdgas-Kohle, 1988, Bd.104, H.10, s.412-417

[10] Samirendra N.B. et al Compressibility Isotherms of Simulated Natural Gases. — J. Chem. Eng. Data, 1990, v.35, No.l, p.35-38

[11] Fitzgerald M.P., Sutton C.M. Measurements of Kapuni and Maui natural gas compressibility factors and comparison with calculated values. — New Zealand Journal of Technology, 1987, v.3, No.4, p.215-218

[12] Jaeschke М., Humphreys A.E. The GERG Databank of High Accuracy Compressibility Factor Measurements. GERG TM4 1990. — GERG Technical Monograph, 1990, 477 p

[13] Jaeschke М., Humphreys A.E. Standard GERG Virial Equation for Field Use. Simplification of the Input Data Requirements for the GERG Virial Equation — an Alternative Means of Compressibility Factor Calculation for Natural Gases and Similar Mixtures. GERG TM5 1991. — GERG Technical Monograph, 1991, 173 p

[14] ISO/TC 193 SC1 № 63. Natural gas — calculation of compression factor. Part 3: Calculation using measured physical properties

[15] ISO/TC 193 SC1 № 62. Natural gas — calculation of compression factor. Part 2: Calculation using a molar composition analysis

[16] ISO 5168:1978 International Standard. Measurement of fluid flow — Estimation of uncertainty of a flow-rate measurement

[17] VDI/VDE 2040, part 2, 1987. Calculation principles for measurement of fluid flow using orifice plates, nozzles and venturi tubes. Equations and formulas

[18] Jaeschke М. et al. High Accuracy Compressibility Factor Calculation for Natural Gases and Similar Mixtures by Use of a Truncated Virial Equation. GERG TM2 1988. — GERG Technical Monograph, 1988, 163 p