"ГОСТ 30319.2-96. Межгосударственный стандарт. Газ природный. Методы расчета физических свойств. Определение коэффициента сжимаемости" (введен в действие Постановлением Госстандарта России от 30.12.1996 N 723) (ред. от 10.03.2004)

Главная // Актуальные документы // ГОСТ (Государственный стандарт)

СПРАВКА

Источник публикации

В данном виде документ опубликован не был.

Первоначальный текст документа опубликован в издании

М.: ИПК Издательство стандартов, 2002.

Информацию о публикации документов, создающих данную редакцию, см. в справке к этим документам.

Примечание к документу

Документ утратил силу с 1 января 2017 года в связи с изданием Приказа Росстандарта от 10.11.2015 N 1744-ст. Взамен введен в действие ГОСТ 30319.2-2015.

С 1 июля 2003 года до вступления в силу технических регламентов акты федеральных органов исполнительной власти в сфере технического регулирования носят рекомендательный характер и подлежат обязательному исполнению только в части, соответствующей целям, указанным в пункте 1 статьи 46 Федерального закона от 27.12.2002 N 184-ФЗ.

Изменение N 1 введено в действие Постановлением Госстандарта России от 10.03.2004 N 167-ст с 01.06.2004.

Название документа

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

(введен в действие Постановлением Госстандарта России от 30.12.1996 N 723)

(ред. от 10.03.2004)

"ГОСТ 30319.2-96. Межгосударственный стандарт. Газ природный. Методы расчета физических свойств. Определение коэффициента сжимаемости"
(введен в действие Постановлением Госстандарта России от 30.12.1996 N 723)
(ред. от 10.03.2004)

Содержание

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

Предисловие

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Приложение Е. Библиография

Введен в действие

Постановлением

Государственного комитета

Российской Федерации

по стандартизации,

метрологии и сертификации

от 30 декабря 1996 г. N 723

МЕЖГОСУДАРСТВЕННЫЙ СТАНДАРТ

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

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

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

Natural gas. Methods of calculation of physical properties.

Definition of compressibility coefficient

ГОСТ 30319.2-96

Список изменяющих документов

(в ред. Изменения N 1, введенного в действие Постановлением

Госстандарта России от 10.03.2004 N 167-ст)

Группа Б19

МКС 75.060

ОКСТУ 0203

Дата введения

1 июля 1997 года

Предисловие

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

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

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

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

Наименование государства	Наименование национального органа по стандартизации
Азербайджанская Республика	Азгосстандарт
Республика Армения	Армгосстандарт
Республика Беларусь	Госстандарт Беларуси
Республика Грузия	Грузстандарт
Республика Казахстан	Госстандарт Республики Казахстан
Киргизская Республика	Киргизстандарт
Республика Молдова	Молдовастандарт
Российская Федерация	Госстандарт России
Республика Таджикистан	Таджикский государственный центр по стандартизации, метрологии и сертификации
Туркменистан	Главгосинспекция Туркменистана
Украина	Госстандарт Украины

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

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

5. ПЕРЕИЗДАНИЕ. Январь 2002 г.

1. НАЗНАЧЕНИЕ И ОБЛАСТЬ ПРИМЕНЕНИЯ

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

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

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

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

2. НОРМАТИВНЫЕ ССЫЛКИ

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

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

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

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

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

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

K = z/z_c, (1)

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

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

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

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

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

(в ред. Изменения N 1, введенного в действие Постановлением Госстандарта России от 10.03.2004 N 167-ст)

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

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

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

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

3) уравнение состояния ВНИИЦСМВ - при добыче и переработке газа.

(в ред. Изменения N 1, введенного в действие Постановлением Госстандарта России от 10.03.2004 N 167-ст)

Таблица 1

Результаты апробации и область применения методов расчета

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

(таблица 1 в ред. Изменения N 1, введенного в действие

Постановлением Госстандарта России от 10.03.2004 N 167-ст)

Метод расчета	Область применения и погрешность метода расчета				Отклонения от экспериментальных данных
Метод расчета	Область применения	, кг/м³	p, МПа	Погрешность , %	, %	, %
NX19 мод.	32 <= H_с.в, МДж/м³ <= 40	< 0,70	< 3	0,12	-0,02	+0,07	-0,09
			3 - 7	0,18	-0,01	+0,37	-0,10
			> 7	0,41	0,17	+0,59	-0,08
	0,66 <= , кг/м³ <= 1,05 0 <= x_а, мол.% <= 15 0 <= x_y, мол.% <= 15 250 <= T, К <= 340 0,1 <= p, МПа <= 12,0	0,70 - 0,75	< 3	0,13	0,01	+0,14	-0,13
			3 - 7	0,29	0,12	+0,46	-0,15
			> 7	0,42	0,27	+0,66	-0,12
		> 0,75	< 3	0,20	0,05	+0,41	-0,13
			3 - 7	0,57	0,24	+1,06	-0,25
			> 7	1,09	0,34	+1,65	-0,40
		0,74 - 1,00 (смеси с H₂S)	0,1 - 11	0,15	-0,02	+0,09	-0,10
УС GERG-91 мод.	20 <= H_с.в, МДж/м³ <= 48	< 0,70	< 3	0,11	0,01	+0,13	-0,04
			3 - 7	0,15	0,02	+0,51	-0,06
			> 7	0,20	0,03	+0,63	-0,06
	0,66 <= , кг/м³ <= 1,05 0 <= x_a, мол.% <= 15 0 <= x_y, мол.% <= 15 250 <= T, К <= 340 0,1 <= p, МПа <= 12,0	0,70 - 0,75	< 3	0,12	-0,01	+0,08	-0,17
			3 - 7	0,15	-0,02	+0,11	-0,43
			> 7	0,19	0,02	+0,16	-0,34
		> 0,75	< 3	0,13	0,01	+0,26	-0,12
			3 - 7	0,15	-0,01	+0,15	-0,30
			> 7	0,19	0,01	+0,65	-0,31
		0,74 - 1,00 (смеси с H₂S)	0,1 - 11	2,10	-0,66	+0,06	-3,10
УС AGA8-92DC	20 <= H_с.в, МДж/м³ <= 48	< 0,70	< 3	0,10	-0,01	+0,03	-0,06
			3 - 7	0,11	-0,01	+0,15	-0,06
			> 7	0,12	0,02	+0,19	-0,04
	0,66 <= , кг/м³ <= 1,05 0 <= x_a, мол.% <= 15 0 <= x_y, мол.% <= 15 250 <= T, К <= 340 0,1 <= p, МПа <= 12,0	0,70 - 0,75	< 3	0,12	-0,01	+0,08	-0,18
			3 - 7	0,15	-0,03	+0,11	-0,43
			> 7	0,19	0,01	+0,16	-0,37
		> 0,75	< 3	0,12	0,01	+0,25	-0,11
			3 - 7	0,15	-0,02	+0,24	-0,24
			> 7	0,17	0,01	+0,31	-0,17
		0,74 - 1,00 (смеси с H₂S)	0,1 - 11	1,30	-0,38	+0,06	-1,88
УС ВНИЦ СМВ	20 <= H_с.в, МДж/м³ <= 48	< 0,70	< 3	0,11	-0,04	+0,01	-0,10
			3 - 7	0,12	-0,04	+0,05	-0,11
			> 7	0,12	-0,01	+0,06	-0,14
	0,66 <= , кг/м³ <= 1,05 0 <= x_a, мол.% <= 15 0 <= x_y, мол.% <= 15 250 <= T, К <= 340 0,1 <= p, МПа <= 12,0	0,70 - 0,75	< 3	0,12	-0,03	+0,08	-0,17
			3 - 7	0,15	-0,02	+0,11	-0,33
			> 7	0,18	0,02	+0,13	-0,27
		> 0,75	< 3	0,13	-0,01	+0,25	-0,11
			3 - 7	0,15	-0,01	+0,18	-0,25
			> 7	0,24	-0,01	+0,28	-0,33
		0,74 - 1,00 (смеси с H₂S)	0,1 - 11	0,36	0,10	+0,54	-0,24
Примечания: 1 При использовании методов расчета NX19 мод. и УС GERG-91 мод. высшую удельную теплоту сгорания (H_с.в) вычисляют по формуле (52) ГОСТ 30319.1. 2 При использовании методов расчета УС AGA8-92DC и УС ВНИЦ СМВ плотность газа при стандартных условиях вычисляют по формуле (16) ГОСТ 30319.1, а высшую удельную теплоту сгорания (H_с.в) - по 7.2 ГОСТ 30319.1 (допускается вычислять по (H_с.в) формуле (52) ГОСТ 30319.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% [14], для уравнения AGA8-92DC - 0,5% [15].

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

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

- систематическое отклонение от экспериментальных данных

; (2)

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

, (3)

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

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

, (4)

где

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

, (5)

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

Погрешность расчета коэффициента сжимаемости

приведена в таблице 1 без учета погрешности исходных данных.

(абзац введен Изменением N 1, введенным в действие Постановлением Госстандарта России от 10.03.2004 N 167-ст)

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

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

, (6)

где

, (7)

, (8)

, (9)

, (10)

. (11)

Корректирующий множитель F в зависимости от интервалов параметров p_а и

вычисляют по формулам:

при 0 <= p_а <= 2 и 0 <=

<= 0,3

, (12)

при 0 <= p_а < 1,3 и -0,25 <=

< 0

, (13)

при 1,3 <= p_а < 2 и -0,21 <=

< 0

(в ред. Изменения N 1, введенного в действие Постановлением

Госстандарта России от 10.03.2004 N 167-ст)

, (14)

где

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

p_а = 0,6714(p/p_пк) + 0,0147, (15)

T_а = 0,71892(T/T_пк) + 0,0007, (16)

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

, (17)

. (18)

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

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

(в ред. Изменения N 1, введенного в действие Постановлением Госстандарта России от 10.03.2004 N 167-ст)

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

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

, (19)

где B_m и C_m - коэффициенты УС;

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

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

, (20)

, (21)

где x_э - молярная доля эквивалентного углеводорода

x_э = 1 - x_а - x_у, (22)

B₁ = -0,425468 + 2,865·10^-3T - 4,62073·10^-6T² +

+ (8,77118·10^-4 - 5,56281·10^-6T + 8,81514·10^-9T²)H +

+ (-8,24747·10^-7 + 4,31436·10^-9T - 6,08319·10^-12T²)x H², (23)

(в ред. Изменения N 1, введенного в действие Постановлением

Госстандарта России от 10.03.2004 N 167-ст)

B₂ = -0,1446 + 7,4091·10^-4T - 9,1195·10^-7T², (24)

B₂₃ = -0,339693 + 1,61176·10^-3T - 2,04429·10^-6T², (25)

B₃ = -0,86834 + 4,0376·10^-3T - 5,1657·10^-6T², (26)

C₁ = -0,302488 + 1,95861·10^-3T - 3,16302·10^-6T² +

+ (6,46422·10^-4 - 4,22876·10^-6T + 6,88157·10^-9T²)H +

+ (-3,32805·10^-7 + 2,2316·10^-9T - 3,67713·10^-12T²) x H², (27)

C₂ = 7,8498·10^-3 - 3,9895·10^-5T + 6,1187·10^-8T², (28)

C₃ = 2,0513·10^-3 + 3,4888·10^-5T - 8,3703·10^-8T², (29)

C₂₂₃ = 5,52066·10^-3 - 1,68609·10^-5T + 1,57169·10^-8T², (30)

C₂₃₃ = 3,58783·10^-3 + 8,06674·10^-6T - 3,25789·10^-8T², (31)

B^* = 0,72 + 1,875·10^-5(320 - T)², (32)

C^* = 0,92 + 0,0013(T - 270). (33)

В формулах (23), (27) H рассчитывают по выражению

H = 128,64 + 47,479M_э, (34)

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

. (35)

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

. (36)

После определения коэффициентов уравнения состояния (19) B_m и C_m рассчитывают фактор сжимаемости при заданных давлении (p, МПа) и температуре (T, К) по формуле

, (37)

где

, (38)

A₀ = 1 + 1,5(B₀ + C₀), (39)

A₁ = 1 + B₀, (40)

B₀ = bB_m, (41)

C₀ = b²C_m, (42)

b = 10³p/(2,7715T). (43)

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

. (44)

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

(в ред. Изменения N 1, введенного в действие Постановлением Госстандарта России от 10.03.2004 N 167-ст)

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

В проекте стандарта ИСО/ТС 193 SC1 N 62 [15] Американской Газовой Ассоциацией для расчета фактора сжимаемости предложено использовать уравнение состояния

, (45)

где B и

- коэффициенты УС;

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

Константы

УС (45) приведены в таблице А.1.

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

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

. (46)

Параметр K_m вычисляют по формуле (53).

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

, (47)

, (48)

(в ред. Изменения N 1, введенного в действие Постановлением

Госстандарта России от 10.03.2004 N 167-ст)

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

Константы

и характерные параметры компонентов

в формулах (47), (48) приведены соответственно в таблицах А.1 и А.2.

Бинарные параметры

и параметры

рассчитывают с использованием следующих уравнений:

, (49)

, (50)

, (51)

, (52)

, (53)

, (54)

, (55)

(в ред. Изменения N 1, введенного в действие Постановлением

Госстандарта России от 10.03.2004 N 167-ст)

где

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

(в ред. Изменения N 1, введенного в действие Постановлением Госстандарта России от 10.03.2004 N 167-ст)

Для расчета фактора сжимаемости по уравнению состояния (45) необходимо определить плотность

при заданных давлении (p, МПа) и температуре (T, К).

Плотность

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

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

, (56)

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

p_п = p/5; (57)

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

, (58)

, (59)

где z⁽^k^-1) рассчитывают из УС (45) при плотности на итерационном шаге (k - 1), т.е. при

, а безразмерный комплекс A₁ определяют из выражения

, (60)

(в ред. Изменения N 1, введенного в действие Постановлением

Госстандарта России от 10.03.2004 N 167-ст)

при этом

;

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

, (61)

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

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

(в ред. Изменения N 1, введенного в действие Постановлением Госстандарта России от 10.03.2004 N 167-ст)

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

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

, (62)

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

- приведенная плотность;

T_п = T/T_пк - приведенная температура;

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

и T_пк - псевдокритические параметры природного газа.

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

, (63)

где

- обобщенные коэффициенты УС, которые приведены в таблице Б.1.

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

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

, (64)

где

; (65)

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

, (66)

где

, (67)

; (68)

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

, (69)

где

. (70)

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

Критические параметры компонентов

, их молярная масса

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

приведены в таблице Б.2, а параметры бинарного взаимодействия

- в таблицах Б.3 и Б.4.

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

(в ред. Изменения N 1, введенного в действие Постановлением Госстандарта России от 10.03.2004 N 167-ст)

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

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

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

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

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

(в ред. Изменения N 1, введенного в действие Постановлением Госстандарта России от 10.03.2004 N 167-ст)

Формулы (71) - (74) исключены с 1 июня 2004 года. - Изменение N 1, введенное в действие Постановлением Госстандарта России от 10.03.2004 N 167-ст).

Для расчета фактора сжимаемости по уравнению состояния (62) необходимо определить плотность

при заданных давлении (p, МПа) и температуре (T, К).

Плотность

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

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

, (75)

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

, (76)

p_п = p/p_пк, (77)

а псевдокритические плотность

, температуру (T_пк) и фактор Питцера

рассчитывают по формулам (64), (66) и (69);

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

, (78)

, (79)

где

рассчитывают из УС (62) при плотности на итерационном шаге (k - 1), т.е. при

, а безразмерный комплекс A₁ определяют из выражения

, (80)

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

, (81)

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

(в ред. Изменения N 1, введенного в действие Постановлением Госстандарта России от 10.03.2004 N 167-ст)

4. ВЛИЯНИЕ ПОГРЕШНОСТИ ИСХОДНЫХ ДАННЫХ НА ПОГРЕШНОСТЬ РАСЧЕТА КОЭФФИЦИЕНТА СЖИМАЕМОСТИ

При измерении расхода и количества природного газа, транспортируемого в газопроводах, давление (p), температуру (T), плотность при стандартных условиях

и состав

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

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

, (82)

(в ред. Изменения N 1, введенного в действие Постановлением

Госстандарта России от 10.03.2004 N 167-ст)

где

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

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

; (83)

(в ред. Изменения N 1, введенного в действие Постановлением

Госстандарта России от 10.03.2004 N 167-ст)

. (84)

В формулах (82) - (84):

q_k - условное обозначение k-го параметра исходных данных

;

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

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

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

При вычислении частных производных по формуле (83) коэффициенты сжимаемости K_qk₊ и K_qk_- рассчитывают при средних параметрах

и параметрах

соответственно. Рекомендуется выбирать

(в ред. Изменения N 1, введенного в действие Постановлением Госстандарта России от 10.03.2004 N 167-ст)

Коэффициент сжимаемости

(среднее значение) рассчитывают по выбранному рекомендуемому методу расчета при средних параметрах

Для методов:

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

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

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

, (85)

где

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

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

по формуле

, (86)

где

- погрешности измеряемых параметров, соответственно, температуры, давления, плотности природного газа при стандартных условиях, содержания азота и диоксида углерода в нем.

Коэффициенты K_T, K_p,

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

- при расчете K по методу NX19 мод.

, (87)

, (88)

, (89)

, (90)

; (91)

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

, (92)

, (93)

, (94)

, (95)

. (96)

5. ПРОГРАММНАЯ И ТЕХНИЧЕСКАЯ РЕАЛИЗАЦИЯ РАСЧЕТА КОЭФФИЦИЕНТА СЖИМАЕМОСТИ

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

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

ПРИЛОЖЕНИЕ А

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

ТАБЛИЦЫ КОНСТАНТ И ПАРАМЕТРОВ УРАВНЕНИЯ СОСТОЯНИЯ

AGA8-92DC

Таблица А.1

Константы уравнения состояния AGA8-92DC

n	a_n	b_n	c_n	k_n	u_n	g_n	q_n	f_n
1	0,153832600	1	0	0	0,0	0	0	0
2	1,341953000	1	0	0	0,5	0	0	0
3	-2,998583000	1	0	0	1,0	0	0	0
4	-0,048312280	1	0	0	3,5	0	0	0
5	0,375796500	1	0	0	-0,5	1	0	0
6	-1,589575000	1	0	0	4,5	1	0	0
7	-0,053588470	1	0	0	0,5	0	1	0
8	2,29129E-9	1	1	3	-6,0	0	0	1
9	0,157672400	1	1	2	2,0	0	0	0
10	-0,436386400	1	1	2	3,0	0	0	0
11	-0,044081590	1	1	2	2,0	0	1	0
12	-0,003433888	1	1	4	2,0	0	0	0
13	0,032059050	1	1	4	11,0	0	0	0
14	0,024873550	2	0	0	-0,5	0	0	0
15	0,073322790	2	0	0	0,5	0	0	0
16	-0,001600573	2	1	2	0,0	0	0	0
17	0,642470600	2	1	2	4,0	0	0	0
18	-0,416260100	2	1	2	6,0	0	0	0
19	-0,066899570	2	1	4	21,0	0	0	0
20	0,279179500	2	1	4	23,0	1	0	0
21	-0,696605100	2	1	4	22,0	0	1	0
22	-0,002860589	2	1	4	-1,0	0	0	1
23	-0,008098836	3	0	0	-0,5	0	1	0
24	3,150547000	3	1	1	7,0	1	0	0
25	0,007224479	3	1	1	-1,0	0	0	1
26	-0,705752900	3	1	2	6,0	0	0	0
27	0,534979200	3	1	2	4,0	1	0	0
28	-0,079314910	3	1	3	1,0	1	0	0
29	-1,418465000	3	1	3	9,0	1	0	0
30	-5,99905Е-17	3	1	4	-13,0	0	0	1
31	0,105840200	3	1	4	21,0	0	0	0
32	0,034317290	3	1	4	8,0	0	1	0
33	-0,007022847	4	0	0	-0,5	0	0	0
34	0,024955870	4	0	0	0,0	0	0	0
35	0,042968180	4	1	2	2,0	0	0	0
36	0,746545300	4	1	2	7,0	0	0	0
37	-0,291961300	4	1	2	9,0	0	1	0
38	7,294616000	4	1	4	22,0	0	0	0
39	-9,936757000	4	1	4	23,0	0	0	0
40	-0,005399808	5	0	0	1,0	0	0	0
41	-0,243256700	5	1	2	9,0	0	0	0
42	0,049870160	5	1	2	3,0	0	1	0
43	0,003733797	5	1	4	8,0	0	0	0
44	1,874951000	5	1	4	23,0	0	1	0
45	0,002168144	6	0	0	1,5	0	0	0
46	-0,658716400	6	1	2	5,0	1	0	0
47	0,000205518	7	0	0	-0,5	0	1	0
48	0,009776195	7	1	2	4,0	0	0	0
49	-0,020487080	8	1	1	7,0	1	0	0
50	0,015573220	8	1	2	3,0	0	0	0
51	0,006862415	8	1	2	0,0	1	0	0
52	-0,001226752	9	1	2	1,0	0	0	0
53	0,002850906	9	1	2	0,0	0	1	0

Таблица А.2

Характерные параметры компонентов

Компонент	Молярная масса	Характерные параметры
Компонент	Молярная масса	E, К	K, м³/кмоль	G	Q	F
Метан	16,0430	151,3183	0,4619255	0,0	0,0	0,0
Этан	30,0700	244,1667	0,5279209	0,079300	0,0	0,0
Пропан	44,0970	298,1183	0,5837490	0,141239	0,0	0,0
н-Бутан	58,1230	337,6389	0,6341423	0,281835	0,0	0,0
и-Бутан	58,1230	324,0689	0,6406937	0,256692	0,0	0,0
Азот	28,0135	99,73778	0,4479153	0,027815	0,0	0,0
Диоксид углерода	44,0100	241,9606	0,4557489	0,189065	0,69	0,0
Сероводород	34,0820	296,3550	0,4618263	0,088500	0,0	0,0
н-Пентан	72,1500	370,6823	0,6798307	0,366911	0,0	0,0
и-Пентан	72,1500	365,5999	0,6738577	0,332267	0,0	0,0
н-Гексан	86,1770	402,8429	0,7139987	0,432254	0,0	0,0
н-Гептан	100,2040	427,5391	0,7503628	0,512507	0,0	0,0
н-Октан	114,2310	450,6472	0,7851933	0,576242	0,0	0,0
Гелий	4,0026	2,610111	0,3589888	0,0	0,0	0,0
Моноксид углерода	28,0100	105,5348	0,4533894	0,038953	0,0	0,0
Кислород	31,9988	122,7667	0,4186954	0,021000	0,0	0,0
Аргон	39,9480	119,6299	0,4216551	0,0	0,0	0,0
Вода	18,0153	514,0156	0,3825868	0,332500	0,0	0,0

Таблица А.3

Параметры бинарного взаимодействия

Компоненты		Параметры бинарного взаимодействия
i	j	E_ij^*	U_ij	K_ij	G_ij^*
Метан	Азот	0,971640	0,886106	1,003630
	Диоксид углерода	0,960644	0,963827	0,995933	0,807653
	Пропан	0,996050	1,023960
	Моноксид углерода	0,990126
	и-Бутан	1,019530
	н-Бутан	0,995474	1,021280
	и-Пентан	1,002350
	н-Пентан	1,003050
	н-Гексан	1,012930
	н-Гептан	0,999758
	н-Октан	0,988563
Азот	Диоксид углевода	1,022740	0,835058	0,982361	0,982746
	Этан	0,970120	0,816431	1,007960
	Пропан	0,945939	0,915502
	Моноксид углерода	1,005710
	и-Бутан	0,946914
	н-Бутан	0,973384	0,993556
	и-Пентан	0,959340
	н-Пентан	0,945520
	н-Гексан	0,937880
	н-Гептан	0,935977
	н-Октан	0,933269
Диоксид углерода	Этан	0,925053	0,969870	1,008510	0,370296
	Пропан	0,960237
	Моноксид углерода	1,500000	0,900000
	и-Бутан	0,906849
	н-Бутан	0,897362
	и-Пентан	0,726255
	н-Пентан	0,859764
	н-Гексан	0,766923
	н-Гептан	0,782718
	н-Октан	0,805823
Этан	Пропан	1,035020	1,080500	1,000460
	и-Бутан		1,250000
	н-Бутан	1,013060	1,250000
	и-Пентан		1,250000
	н-Пентан	1,005320	1,250000
Пропан	н-Бутан	1,004900

ПРИЛОЖЕНИЕ Б

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

ТАБЛИЦЫ КОЭФФИЦИЕНТОВ И ПАРАМЕТРОВ

УРАВНЕНИЯ СОСТОЯНИЯ ВНИЦ СМВ

Таблица Б.1

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

k	1	a_kl	b_kl
1	0	6,087766·10^-1	-7,187864·10^-1
2	0	-4,596885·10^-1	1,067179·10¹
3	0	1,149340·10⁰	-2,576870·10¹
4	0	-6,075010·10^-1	1,713395·10¹
5	0	-8,940940·10^-1	1,617303·10¹
6	0	1,144404·10⁰	-2,438953·10¹
7	0	-3,457900·10^-1	7,156029·10⁰
8	0	-1,235682·10^-1	3,350294·10⁰
9	0	1,098875·10^-1	-2,806204·10⁰
10	0	-2,193060·10^-2	5,728541·10^-1
1	1	-1,832916·10⁰	6,057018·10⁰
2	1	4,175759·10⁰	-7,947685·10¹
3	1	-9,404549·10⁰	2,167887·10²
4	1	1,062713·10¹	-2,447320·10²
5	1	-3,080591·10⁰	7,804753·10¹
6	1	-2,122525·10⁰	4,870601·10¹
7	1	1,781466·10⁰	-4,192715·10¹
8	1	-4,303578·10^-1	1,000706·10¹
9	1	-4,963321·10^-2	1,237872·10⁰
10	1	3,474960·10^-2	-8,610273·10^-1
1	2	1,317145·10⁰	-1,295347·10¹
2	2	-1,073657·10¹	2,208390·10²
3	2	2,395808·10¹	-5,864596·10²
4	2	-3,147929·10¹	7,444021·10²
5	2	1,842846·10¹	-4,470704·10²
6	2	-4,092685·10⁰	9,965370·10¹
7	2	-1,906595·10^-1	5,136013·10⁰
8	2	4,015072·10^-1	-9,576900·10⁰
9	2	-1,016264·10^-1	2,419650·10⁰
10	2	-9,129047·10^-3	2,275036·10^-1
1	3	-2,837908·10⁰	1,571955·10¹
2	3	1,534274·10¹	-3,020599·10²
3	3	-2,771885·10¹	6,845968·10²
4	3	3,511413·10¹	-8,281484·10²
5	3	-2,348500·10¹	5,600892·10²
6	3	7,767802·10⁰	-1,859581·10²
7	3	-1,677977·10⁰	3,991057·10¹
8	3	3,157961·10^-1	-7,567516·10⁰
9	3	4,008579·10^-3	-1,062596·10^-1
1	4	2,606878·10⁰	-1,375957·10¹
2	4	-1,106722·10¹	2,055410·10²
3	4	1,279987·10¹	-3,252751·10²
4	4	-1,211554·10¹	2,846518·10²
5	4	7,580666·10⁰	-1,808168·10²
6	4	-1,894086·10⁰	4,605637·10¹
1	5	-1,155750·10⁰	6,466081·10⁰
2	5	3,601316·10⁰	-5,739220·10¹
3	5	-7,326041·10^-1	3,694793·10¹
4	5	-1,151685·10⁰	2,077675·10¹
5	5	5,403439·10^-1	-1,256783·10¹
1	6	9,060572·10^-2	-9,775244·10^-1
2	6	-5,151915·10^-1	2,612338·10⁰
3	6	7,622076·10^-2	-4,059629·10^-1
1	7	4,507142·10^-2	-2,298833·10^-1

Таблица Б.2

Физические свойства компонентов природного газа,

используемые в уравнении состояния ВНИЦ СМВ

Компоненты	Химическая формула	Молярная масса M_i	Критические параметры				, кг/м³	Фактор Питцера
Компоненты	Химическая формула	Молярная масса M_i	p_ki, МПа	, кг/м³	T_k_i, К	z_k_i	, кг/м³	Фактор Питцера
Метан	CH₄	16,043	4,5988	163,03	190,67	0,2862	0,6682	0,0006467
Этан	C₂H₆	30,070	4,88	205,53	305,57	0,2822	1,2601	0,1103
Пропан	C₃H₈	44,097	4,25	218,54	369,96	0,2787	1,8641	0,1764
н-Бутан	н-C₄H₁₀	58,123	3,784	226,69	425,40	0,2761	2,4956	0,2213
и-Бутан	и-C₄H₁₀	58,123	3,648	225,64	407,96	0,2769	2,488	0,2162
Азот	N₂	28,0135	3,390	315,36	125,65	0,2850	1,16490	0,04185
Диоксид углерода	CO₂	44,010	7,386	466,74	304,11	0,2744	1,8393	0,2203
Сероводород	H₂S	34,082	8,940	349,37	373,18	0,2810	1,4311	0,042686
Примечания: 1 Плотность (), температура (T_k_i) в критической точке и фактор Питцера () отличаются от литературных данных и применимы только для уравнения состояния ВНИЦ СМВ. 2 - плотность i-го компонента при стандартных условиях

Таблица Б.3

Параметры бинарного взаимодействия

j	i
j	CH₄	C2H₆	C3H₈	н-C₄H₁₀	и-C₄H₁₀	N₂	CO₂	H₂S
CH₄	0,0	0,036	0,076	0,121	0,129	0,060	0,074	0,089
C₂H₆	-	0,0	0,0	0,0	0,0	0,106	0,093	0,079
C₃H₈	-	-	0,0	0,0	0,0	0,0	0,0	0,0
н-C₄H₁₀	-	-	-	0,0	0,0	0,0	0,0	0,0
и-C₄H₁₀	-	-	-	-	0,0	0,0	0,0	0,0
N₂	-	-	-	-	-	0,0	0,022	0,211
CO₂	-	-	-	-	-	-	0,0	0,089
H₂S	-	-	-	-	-	-	-	0,0

Таблица Б.4

Параметры бинарного взаимодействия

j	i
j	CH₄	C₂H₆	C₃H₈	н-C₄H₁₀	и-C₄H₁₀	N₂	CO₂	H₂S
CH₄	0,0	-0,074	-0,146	-0,258	-0,222	-0,023	-0,086	0,0
C₂H₆	-	0,0	0,0	0,0	0,0	0,0	0,0	0,0
C₃H₈	-	-	0,0	0,0	0,0	0,0	0,0	0,0
н-C₄H₁₀	-	-	-	0,0	0,0	0,0	0,0	0,0
и-C₄H₁₀	-	-	-	-	0,0	0,0	0,0	0,0
N₂	-	-	-	-	-	0,0	-0,064	0,0
CO₂	-	-	-	-	-	-	0,0	-0,062
H₂S	-	-	-	-	-	-	-	0,0

ПРИЛОЖЕНИЕ В

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

ЛИСТИНГ ПРОГРАММЫ РАСЧЕТА

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

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

C * *

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

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/' метана (CH4)',' этана (C2H6)',' пропана (C3H8)',

*' н-бутана (н-C4H10)',' и-бутана (и-C4H10)',' азота (N2)',

*' диоксида углерода (CO2)',' сероводорода (H2S)',

*' ацетилена (C2H2)',' этилена (C2H4)',' пропилена (C3H6)',

*' н-пентана (н-C5H12)',' и-пентана (и-C5H12)',

*' нео-пентана (нео-C5H12)',' н-гексана (н-C6H14)',

*' бензола (C6H6)',' н-гептана (н-C7H16)',' толуола (C7H8)',

*' н-октана (н-C8H18)',' н-нонана (н-C9H20)',

*' н-декана (н-C10H22)',' гелия (He)',' водорода (H2)',

*' моноксида углерода (CO)',' кислорода (O2)'/

200 WRITE(*,100)

CALL VAR(NVAR)

IF(NVAR.EQ.5) GO TO 134

WRITE(*,100)

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.1) WRITE(*,33)

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

WRITE(*,'(A\)') ' азота (N2)

READ(*,*)YA

YA = YA/100.

WRITE(*,'(A\)') ' диоксида углерода (CO2) '

READ(*,*)YY

YY = YY/100.

ELSE

WRITE(*,35)

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

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

READ(*,*)NPR

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

IF(NPR.EQ.1) 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)

D0 7 J=1,NT

T=TI(J)

IF(NVAR.EQ.1) 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 ИС=0

DO 11 J=1,NT

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

ИС=ИС+1

11 CONTINUE

IF(ИС.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

TI(I)=TI(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,' Расчет коэффициента сжимаемости природного газа'//

*10X,' ----------------------Метод расчета------------------------- '/

*10X,' '/

*10X,' 1. Модифицированный метод NX19 '/

*10X,' '/

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

*10X,' '/

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

*10X,' '/

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

*10X,' '/

*10X,' ------------------------------------------------------------ '/

WRITE(*,5)

5 FORMAT(/,3X,

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

*\)

READ(*,*)NVAR

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(YA,YY,PI,TI,ZP,NP,NT,NVAR,AR)

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

CHARACTER*26 AR(25),FNAME

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

*AT(6)*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*'----------'/,ALIN2/5*'----------'/,LIN3/6*'----------'/,

*LIN4/'----------'/,A/'----------'/

DATA AT/

*' T, К',' T, К',' T, К',' T, К',

*' T, К',' T, К'/

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))',

*'(3X,F9.6,1X,A6,5(3X,F6.4))','(3X,F9.6,1X,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.1) 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,880)METH(NVAR)

88 FORMAT(

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

*18X,A31/)

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

ELSE

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

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

NW=NW+1

I=1

9 J=I+1

13 CONTINUE

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

WRITE(l,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

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

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.1) 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.1) WRITE(1,21)LIN4

21 FORMAT(' p, МПа ',6AА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

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(NL1.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.1) 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.1) 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.1) PAUSE

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

NW=0

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

IF(NL.EQ.1) 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.1) 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(1)

WRITE(*,7)

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

WRITE(*,66)

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

*', 0 - нет, 1 - да '\)

READ(*,*)NBOLB

IF(NBOLB.EQ.1) GO TO 22

RETURN

END

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

C * *

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

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

C * *

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

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.305131 D0)*T1/TA**4

B1=2.*T1/3.-T0**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 **************************************************************

C * *

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

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

C * *

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,T0,R

DATA BMO/,0838137D0/,BMl/-,00851644D0/,WDO/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

X11=X1*X1

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

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

H=47,479*BMY+128,64D0

RETURN

3 T=T1

TC=T1-T0

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,BH)

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),BRHH2(3),BR22(3),BR23(3),BR33(3)

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

COMMON/XBLOK/Xl,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/Xl,X2,X3,X11,X12,X13,X22,X23,X33

T2=T*T

C111=CR111H0(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.0D0.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.*X1*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

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

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

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

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

B1=1D3*P/2.7715/T

B0=B1*Bm

C0=B1**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-А01**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 CR111HO/-.302488D0,.195861D-2,-.316302D-5/,

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

* CR111 H2/-.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/,

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

END 46

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

C * *

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

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.1) 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)=0D0

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 I=1,25

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

*ИС=ИС+1

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

NC=NC+1 NI(NC)=I-ИС

Y(NC)=YI(I)

YSUM=YSUM+Y(NC)

11 CONTINUE

CALL MOLDO1(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)

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.5D-5) Z=0D0

RETURN

END

SUBROUTINE PARIN1

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

REAL*8 K.IJ

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

EIJ(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

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 B1(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 B1(N)=B1(N)+2.*Y(I)*Y(J)*(GIJM(I,J)+1D0-GN(N))**GN(N)*

*(QI(I)*QI(J)+lD0-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)

*/COEFl/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

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

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

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*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 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.5391 D0,450.6472D0,472.1194D0,488.7633D0,

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

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

*.4479153D0,.4557489 D0,.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 * *

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

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

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),BIJ10,8)

*/B/B(10,8)/RM/RM/Y/Y(8)/BM/BM(8)/NI/NI(8)/NC/NC/Z/Z

RM=8.31451D0

IF(ICALC.NE.l) 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=1,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/

DATA ROI/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 BMM=BMM+YI(I)*BMI(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 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(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.YI(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 Gl(I)=YI(I)*R01(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)=1

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(I).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

DIJ(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.D0/3.D0)

D0 3 I=1,NC

VCIJ(I,I)=VC(I)

PIIJ(I,I)=PII(I)

TCIJ(I,I)=TC(I)

D0 3 J=1,NC

IF(I.GE.J) GO TO 3

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

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

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

SUBROUTINE FUN(X)

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

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

ITER=1

1 CONTINUE

NPRIZ=0

IF(ITER.NE.1) 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.2203 D0,0.042686D0/

DATA AIJ/.6087766D0,-.4596885D0,1.14934D0,-.607501D0,

*-.894094D0,1.144404 D0,-.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.95808D0,-31.47929D0,

*18.42846D0,-4.092685D0,-.1906595D0,.4015072D0,-.1016264D0,

*-.9129047D-2,-2.837908D0,15.34274D0,-27.71885D0,35.11413D0,

*-23.485D0,7.767802D0,-l.677977D0,.3157961 D0,.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.136013 D0,-9.5769D0,2.41965 D0,

*.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 кг/м3

Содержание:

азота 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 кг/м3

Содержание:

азота 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 кг /м3

Давление 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

Исходные данные (заданные параметры)	Значения
Исходные данные (заданные параметры)	минимальное	максимальное	погрешности, %
Давление, МПа	1,991	2,011	1,00
Температура, К	269,50	270,50	0,35
Плотность, кг/м³ (0,101325 МПа, 293,15 К)	0,6790	0,6808	0,25
Содержание, мол.%:
азота (N₂)	0,8769	0,8947	2,00
диоксида углерода (CO₂)	0,0661	0,0675	2,00

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

Погрешность расчета: по формуле (82) - 0,09%; по формуле (86) - 0,07%.

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

Исходные данные (заданные параметры)	Значения
Исходные данные (заданные параметры)	минимальное	максимальное	погрешности, %
Давление, МПа	1,991	2,011	1,00
Температура, К	269,50	270,50	0,35
Плотность, кг/м³ (0,101325 МПа, 293,15 К)	0,6790	0,6808	0,25
Содержание, мол.%:
азота (N₂)	0,8769	0,8947	2,00
диоксида углерода (CO₂)	0,0661	0,0675	2,00

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

Погрешность расчета: по формуле (82) - 0,09%; по формуле (86) - 0,09%.

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

Исходные данные (заданные параметры)	Значения
Исходные данные (заданные параметры)	минимальное	максимальное	погрешности, %
Давление, МПа	1,991	2,011	1,00
Температура, К	269,50	270,50	0,35
Содержание, мол.%:
метана (CH₄)	97,2722	99,2722	2,00
этана (C₂H₆)	0,5030	0,5288	5,00
пропана (C₃H₈)	0,1607	0,1607	-
н-бутана (н-C₄H₁₀)	0,0592	0,0592	-
азота (N₂)	0,8769	0,8947	2,00
диоксида углерода (CO₂)	0,0661	0,0675	2,00
н-пентана (н-C₂H₁₂)	0,0157	0,0157	-
н-гексана (н-C₆H₁₄)	0,0055	0,0055	-
н-гептана (н-C₇H₁₆)	0,0016	0,0016	-
н-октана (н-C₈H₁₈)	0,0009	0,0009	-
гелия (He)	0,0157	0,0157	-

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

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

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

Исходные данные (заданные параметры)	Значения
Исходные данные (заданные параметры)	минимальное	максимальное	погрешности, %
Давление, МПа	1,076	1,086	1,00
Температура, К	322,65	323,65	0,31
Содержание, мол.%:
метана (CH₄)	88,3700	90,1700	2,00
этана (C₂H₆)	2,2030	2,3170	5,00
пропана (C₃H₈)	1,0600	1,0600	-
н-бутана (н-C₄H₁₀)	0,0100	0,0100	-
азота (N₂)	0,0396	0,0404	2,00
диоксида углерода (CO₂)	4,2570	4,3430	2,00
сероводорода (H₂S)	3,0500	3,0500	-
пропилена (C₃H₆)	0,0100	0,0100	-

Коэффициент сжимаемости (среднее значение) - 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 K 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. - Hvdrocarbon 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, D. 230 - 233.

[7] Jaeschke M., 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 M. Realgasverhalten Einheitliche Berechnungsmoglichkeiten von Erdgas L und H. - Gas und Wasserfach. Gas/Erdgas, 1988, v. 129, No. 1, s. 30 - 37.

[9] Blanke W., Weiss R. pvT-Eigenschaften und Adsorptions- verhalten von Erdgas bei Temperaturen zwischen 260 K und 330 K 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 M., Humphreys A.E. The GERG Databank of High Accuracy Compressibility Factor Measurements. GERG TM4 1990. - GERG Technical Monograph, 1990, 477 p.

[13] Jaeschke M., 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] ICO/TC 193 SC1 N 63. Natural gas - calculation of compression factor. Part 3: Calculation using measured physical properties.

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

[16] ИСО 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 M. 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.