| |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
ГОСТ 30319.2-96 МЕЖГОСУДАРСТВЕННЫЙ СТАНДАРТ
ГАЗ ПРИРОДНЫЙ МЕТОДЫ РАСЧЕТА ОПРЕДЕЛЕНИЕ КОЭФФИЦИЕНТА СЖИМАЕМОСТИ
МЕЖГОСУДАРСТВЕННЫЙ
СОВЕТ Минск
Предисловие 1 РАЗРАБОТАН Всероссийским научно-исследовательским центром стандартизации, информации и сертификации сырья, материалов и веществ (ВНИЦ СМВ) Госстандарта России; фирмой «Газприборавтоматика» акционерного общества «Газавтоматика» РАО «Газпром» ВНЕСЕН Госстандартом Российской Федерации 2 ПРИНЯТ Межгосударственным Советом по стандартизации, метрологии и сертификации (протокол № 9-96 от 12 апреля 1996 г.) За принятие проголосовали:
3 ПОСТАНОВЛЕНИЕМ Государственного комитета Российской Федерации по стандартизации, метрологии и сертификации от 30 декабря 1996 г. № 723 межгосударственный стандарт ГОСТ 30319.2-96 введен в действие непосредственно в качестве государственного стандарта Российской Федерации с 1 июля 1997 г. 4 ВВЕДЕН ВПЕРВЫЕ 5 ПЕРЕИЗДАНИЕ
СОДЕРЖАНИЕ ГОСТ 30319.2-96 МЕЖГОСУДАРСТВЕННЫЙ СТАНДАРТ Газ природный МЕТОДЫ РАСЧЕТА ФИЗИЧЕСКИХ СВОЙСТВ Определение коэффициента сжимаемости Natural gas. Methods of calculation of
physical properties. Дата введения 1997-07-01 1 Назначение и область примененияНастоящий стандарт устанавливает четыре метода определения коэффициента сжимаемости природного газа: при неизвестном полном компонентном составе природного газа (два метода) и известном компонентном составе. Стандарт устанавливает предпочтительные области применения каждого метода по измеряемым параметрам (давление, температура, плотность природного газа при стандартных условиях и компонентный состав природного газа), однако не запрещает использование любого из методов и в других областях. Допускается применять любые другие методы расчета коэффициента сжимаемости, однако погрешность расчета коэффициента сжимаемости по этим методам не должна превышать погрешностей, приведенных в настоящем стандарте (см. 3.2.1). Используемые в настоящем стандарте определения и обозначения приведены в соответствующих разделах ГОСТ 30319.0. 2 Нормативные ссылкиВ настоящем стандарте использованы ссылки на следующие стандарты: ГОСТ 30319.0-96 Газ природный. Методы расчета физических свойств. Общие положения ГОСТ 30319.1-96 Газ природный. Методы расчета физических свойств. Определение физических свойств природного газа, его компонентов и продуктов его переработки 3 Определение коэффициента сжимаемости3.1 Общие положенияКоэффициент сжимаемости вычисляют по формуле К = z/zc, (1) где z и zc - фактор сжимаемости соответственно при рабочих и стандартных условиях. Рабочие условия характеризуются такими давлениями и температурами, которые определяются измерениями в процессе добычи, переработки и транспортирования природного газа. Давление pc и температура Tc при стандартных условиях приведены в ГОСТ 30319.0. 3.2 Методы расчета коэффициента сжимаемости3.2.1 Пределы применимости методов расчета и область их применения и погрешности расчета коэффициента сжимаемостиВ таблице 1 приведены общие результаты апробации методов расчета. Апробация проведена на обширном массиве высокоточных экспериментальных данных о факторе сжимаемости природного газа [1-12]. Погрешность данных не превышает 0,1 %. Таблица 1 - Результаты апробации и область применения методов расчета коэффициента сжимаемости природного газа
Примечания: 1 При использовании методов расчета NХ19 мод. и УС GERG-91 мод. высшую удельную теплоту сгорания (Нс.в) вычисляют по формуле (52) ГОСТ 30319.1. 2 При использовании методов расчета УС AGA8-92DС и УС ВНИЦ СМВ плотность газа при стандартных условиях (rс) вычисляют по формуле (16) ГОСТ 30319.1, а высшую удельную теплоту сгорания (Нс.в) - по 7.2 ГОСТ 30319.1 (допускается вычислять Нс.в по формуле (52) ГОСТ 30319.1). Новая редакция, (Изм. № 1). Для расчета коэффициента сжимаемости природного газа при определении его расхода и количества рекомендуется применять: 1) модифицированный метод NХ19 мод. - при распределении газа потребителям; 2) модифицированное уравнение состояния (УС) GERG-91 мод. [13, 14] и УСАGА8-92DС [15] - при транспортировании газа по магистральным газопроводам; 3) уравнение состояния ВНИЦСМВ - при добыче и переработке газа. Метод 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 приняты следующие обозначения: 1) dсист - систематическое отклонение от экспериментальных данных
2) diмакс - максимальное отклонение в i-й точке экспериментальных данных
где Красч и Кэксп - соответственно расчетный и экспериментальный коэффициенты сжимаемости; 3) d - погрешность расчета коэффициента сжимаемости по ИСО 5168 [16]
где dст - стандартное отклонение, которое вычисляется из выражения
dэксп - погрешность экспериментальных данных (0,1 %). Погрешность расчета коэффициента сжимаемости d приведена в таблице 1 без учета погрешности исходных данных. Измененная редакция, Изм. № 1. 3.2.2 Модифицированный метод NX19 мод.В соответствии с требованиями стандарта Германии [17] расчет фактора сжимаемости по модифицированному методу NX19 мод. основан на использовании уравнения следующего вида где
Корректирующий множитель F в зависимости от интервалов параметров ра и DТа вычисляют по формулам: при 0 £ ра £ 2 и 0 £ DТа £ 0,3
при 0 £ ра< 1,3 и -0,25 £ DТа < 0
при 1,3 £ pа < 2 и -0,21 £ DTa < 0
где DTa = Ta - 1,09. Параметры рa и Тa определяются по следующим соотношениям:
где рпк и Тпк - псевдокритические значения давления и температуры, определяемые по формулам (48) и (49) ГОСТ 30319.1, а именно: В формулах (17), (18) вместо молярных долей диоксида углерода и азота допускается применять их объемные доли (ry и ra). Коэффициент сжимаемости природного газа вычисляют по формуле (1), при этом фактор сжимаемости при рабочих условиях рассчитывают по формулам (6) - (18) настоящего стандарта, а фактор сжимаемости при стандартных условиях - по формуле (24) ГОСТ 30319.1 Измененная редакция, Изм. № 1. 3.2.3 Модифицированное уравнение состояния GERG-91 мод.Европейская группа газовых исследований на базе экспериментальных данных, собранных в [12], и уравнения состояния вириального типа [18], разработала и опубликовала в [13, 14] УС где Вm и Сm - коэффициенты УС; rм - молярная плотность, кмоль/м3. Коэффициенты уравнения состояния определяют из следующих выражений:
где хэ - молярная доля эквивалентного углеводорода
В формулах (23), (27) Н рассчитывают по выражению
где Мэ - молярная масса эквивалентного углеводорода, значение которой определяется из выражения В выражении (35) молярную долю эквивалентного углеводорода (xэ) рассчитывают с использованием формулы (22), а фактор сжимаемости при стандартных условиях (zс) рассчитывают по формуле (24) ГОСТ 30319.1, а именно После определения коэффициентов уравнения состояния (19) Вm и Сm рассчитывают фактор сжимаемости при заданных давлении (р, МПа) и температуре (Т, К) по формуле где
С0 = b2Cm, (42) Коэффициент сжимаемости природного газа рассчитывают по формуле (1), а именно
Фактор сжимаемости при стандартных условиях zс рассчитывают по формуле (36). Измененная редакция, Изм. № 1. 3.2.4 Уравнение состояния AGA8-92DCВ проекте стандарта ISO/TC 193 SC1 № 62 [15] Американской Газовой Ассоциацией для расчета фактора сжимаемости предложено использовать уравнение состояния где В и Сn* - коэффициенты УС; rм - молярная плотность, кмоль/м3. Константы {bn, cn, kn} УС (45) приведены в таблице A.1. Если состав газа задан в объемных долях, то молярные доли рассчитываются по формуле (12) ГОСТ 30319.1. Приведенную плотность определяют по формуле
Параметр Кт вычисляют по формуле (53). Коэффициенты УС рассчитывают из следующих соотношений: где N - количество компонентов в природном газе. Константы {an, un, gn, qn, fn} и характерные параметры компонентов {Еi, Кi, Gi, Qi, Fi} в формулах (47), (48) приведены соответственно в таблицах А.1 и А.2. Бинарные параметры {Eij, Gij} и параметры {U, G, Km, Q, F} рассчитывают с использованием следующих уравнений:
(i ¹ j)
(i ¹ j)
где {Eij*, Gij*, Uij*, Kij*} - параметры бинарного взаимодействия, которые даны в таблице А.3. Параметры бинарного взаимодействия, которые не приведены в этой таблице, а также при i = j, равны единице. Для расчета фактора сжимаемости по уравнению состояния (45) необходимо определить плотность rм при заданных давлении (р, МПа) и температуре (Т, К). Плотность rм из УС (45) определяют по методу Ньютона в следующем итерационном процессе: 1) начальную плотность определяют по формуле
где приведенное давление вычисляют из выражения
2) плотность на k-м итерационном шаге определяют из выражений
где z(k-1) - рассчитывают из УС (45) при плотности на итерационном шаге (k-1), т.е. при rм(k-1), а безразмерный комплекс А1 определяют из выражения
при этом rп = Кт3rм(k-1); 4) критерий завершения итерационного процесса если критерий (61) не выполняется, то необходимо продолжить итерационный процесс, начиная с пункта 2) алгоритма. После определения фактора сжимаемости при рабочих и стандартных условиях по формуле (1) рассчитывают коэффициент сжимаемости. Измененная редакция, Изм. № 1. 3.2.5 Уравнение состояния ВНИЦ СМВВо Всероссийском научно-исследовательском центре стандартизации, информации и сертификации сырья, материалов и веществ (ВНИЦ СМВ) для расчета фактора сжимаемости природного газа разработано уравнение состояния где ckl - коэффициенты УС; rп = rм/rпк - приведенная плотность; Тп = Т/Тпк - приведенная температура; rм - молярная плотность, кмоль/м3; rпк и Тпк - псевдокритические параметры природного газа. Коэффициенты УС определяют по формуле
где {akl, bkl} - обобщенные коэффициенты УС, которые приведены в таблице Б.1. Псевдокритические параметры природного газа и его фактор Питцера вычисляют по формулам: - псевдокритическую плотность где ( - псевдокритическую температуру где
( - фактор Питцера где В соотношениях (64) - (70) N - число основных компонентов природного газа (метана, этана, пропана, н-бутана, и-бутана, азота, диоксида углерода, сероводорода). Критические параметры компонентов {rкi, rкj, Ткj, Ткj}, их молярная масса {Мi, Мj,} и факторы Питцера {Wi, Wj} приведены в таблице Б.2, а параметры бинарного взаимодействия {xij, lij} - в таблицах Б.3 и Б.4. Если заданный компонентный состав природного газа включает, кроме основных, другие компоненты (но не более 1 % в сумме), то молярные доли этих компонентов прибавляют к соответствующим долям основных компонентов следующим образом: - ацетилен и этилен к этану; - пропилен к пропану; - углеводороды от н-пентана и выше к н-бутану; - прочие компоненты к азоту. Если состав газа задан в объемных долях, то молярные доли рассчитывают по формуле (12) ГОСТ 30319.1. Для расчета фактора сжимаемости по уравнению состояния (62) необходимо определить плотность rм при заданных давлении (р, МПа) и температуре (Т, К). Плотность rм из УС (62) определяют по методу Ньютона в следующем итерационном процессе: 1) начальную плотность определяют по формуле
где приведенное давление вычисляют из выражений
а псевдокритические плотность (rпк), температуру (Тпк) и фактор Питцера (W) рассчитывают по формулам (64), (66) и (69); 2) плотность на k-м итерационном шаге определяется из выражений
где z(k-1) рассчитывают из УС (62) при плотности на итерационном шаге (k-1), т.е. при rм(k-1), a безразмерный комплекс A1 определяют из выражения
4) критерий завершения итерационного процесса. если критерий (81) не выполняется, то необходимо продолжить итерационный процесс, начиная с пункта 2) алгоритма. После определения фактора сжимаемости при рабочих и стандартных условиях по формуле (1) рассчитывают коэффициент сжимаемости. Измененная редакция, Изм. № 1. 4 Влияние погрешности исходных данных на погрешность расчета коэффициента сжимаемостиПри измерении расхода и количества природного газа, транспортируемого в газопроводах, давление (р), температуру (T), плотность при стандартных условиях (rc) и состав (хi) измеряют с определенной погрешностью. Перечисленные параметры являются исходными данными для расчета коэффициента сжимаемости. В соответствии с рекомендациями ИСО 5168 [16] погрешность расчета коэффициента сжимаемости, которая появляется в связи с погрешностью измерения исходных данных, определяют по формуле где dи.д - погрешность расчета коэффициента сжимаемости, связанная с погрешностью измерения исходных данных; dqk - погрешность измерения параметра исходных данных;
qk - условное обозначение k-го параметра исходных данных (р. Т, rc, хi,); `qk - среднее значение k-го параметра в определенный промежуток времени (сутки, месяц, год и т.д.); qkмакс и qkмин - максимальное и минимальное значения k-го параметра в определенный промежуток времени; Nq - количество параметров исходных данных. При
вычислении частных производных по формуле (83) коэффициенты сжимаемости Кqk+ и Кqk-
рассчитывают при средних параметрах Коэффициент сжимаемости `К (среднее значение) рассчитывают по выбранному рекомендуемому методу расчета при средних параметрах qk. Для методов: 1) NX 19 мод. и УС GERG-91 мод. - Nq = 5 и параметрами исходных данных являются давление, температура, плотность при стандартных условиях, молярные доли азота и диоксида углерода; 2) УС AGA8-92DC и УС ВНИЦ СМВ - Nq = 2 + N (N - количество компонентов) и параметрами исходных данных являются давление, температура и молярные доли компонентов природного газа, причем для УС ВНИЦ СМВ учитываются молярные доли только основных компонентов газа. Общую погрешность расчета коэффициента сжимаемости определяют по формуле
где d - погрешность расчета коэффициента сжимаемости, которая для каждого метода приведена в 3.2.1. Для методов NX19 мод. и УС GERG-91 мод. допускается рассчитывать погрешность dи.д по формуле где dТ, dp, drc, dxa и dxy - погрешности измеряемых параметров, соответственно, температуры, давления, плотности природного газа при стандартных условиях, содержания азота и диоксида углерода в нем. Коэффициенты КT, Кр, Кrc, Кxa и Кxу в зависимости от метода, используемого для расчета коэффициента сжимаемости K, определяются по следующим выражениям (см. формулы (34) - (38) или (39) - (43) ГОСТ 30319.1): - при расчете К по методу NX19 мод.
- при расчете К по методу GERG-91
Измененная редакция, Изм. № 1. 5 Программная и техническая реализация расчета коэффициента сжимаемостиРасчет коэффициента сжимаемости природного газа по указанным в стандарте методам реализован на ПЭВМ, совместимых с IBM PC/AT/XT, на языке программирования ФОРТРАН-77. Листинг программы приведен в приложении В. В приложениях Г и Д приведены примеры расчета соответственно коэффициента сжимаемости и погрешности вычисления коэффициента сжимаемости, которая вызвана погрешностью определения исходных данных. ПРИЛОЖЕНИЕ АТаблицы констант и параметров уравнения состояния AGA8-92DCТаблица А.1 - Константы уравнения состояния AGA8-92DC
Таблица А.2 - Характерные параметры компонентов
Таблица А.3 - Параметры бинарного взаимодействия
ПРИЛОЖЕНИЕ БТаблицы коэффициентов и параметров уравнения состояния ВНИЦ СМВТаблица Б.1 - Обобщенные коэффициенты уравнения состояния ВНИЦ СМВ
Таблица Б.2 - Физические свойства компонентов природного газа, используемые в уравнении состояния ВНИЦ СМВ
Таблица Б.3 - Параметры бинарного взаимодействия xij
Таблица Б.4 - Параметры бинарного взаимодействия lij
ПРИЛОЖЕНИЕ В(рекомендуемое)Листинг программы расчета коэффициента сжимаемости природного газа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 кг/м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
Коэффициент сжимаемости (среднее значение) - 0,9520 Погрешность расчета: по формуле (82) - 0,09 %; по формуле (86) - 0,07 %. Д.2 Уравнение состояния GERG-91
Коэффициент сжимаемости (среднее значение) - 0,9521 Погрешность расчета: по формуле (82) - 0,09 %; по формуле (86) - 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 [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
Ключевые слова: природный газ, методы расчета коэффициента сжимаемости, давление, температура, плотность при стандартных условиях, компонентный состав, молярные и объемные доли, коэффициент сжимаемости, фактор сжимаемости, плотность, погрешность, уравнение состояния, итерационный процесс, листинг программы
| |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
© 2013 Ёшкин Кот :-) |