способ определения показателей физических свойств природного газа и устройство для его автоматического осуществления

Формула изобретения

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

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

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

Описание изобретения к патенту

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

Известен способ определения показателей физических свойств природного газа (плотности при стандартных условиях, динамической вязкости и коэффициента сжимаемости при заданных температуре и давлении, показателя адиабаты, удельной теплоты сгорания газа) по компонентному составу газа (ГОСТ 30319.1-96. Газ природный. Методы расчета свойств природного газа, его компонентов и продуктов переработки). Способ состоит в том, что контролируемый газ подготавливают, определяют компонентный состав, как правило, при помощи хроматографа, и вычисляют перечисленные выше показатели физических свойств газа с помощью вычислительного устройства с программным обеспечением.

Недостатком указанного способа является его сложность, определяемая сложностью определения компонентного состава природного газа, и низкая оперативность.

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

Недостаток известного устройства - его сложность и низкая оперативность получения результатов определения показателей физических свойств природного газа.

Задача, на решение которой направлено предлагаемое изобретение, состоит в том, чтобы создать такое техническое решение, при использовании которого упрощалась бы процедура и повышалась оперативность определения показателей физических свойств природного газа.

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

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

Определение показателей физических свойств природного газа заключается в том, что: контролируемый газ после его подготовки (фильтрации и, при необходимости, осушки) пропускают через последовательно установленные турбулентный и ламинарный дроссели; измеряют абсолютное давление газа перед турбулентным дросселем, в междроссельном участке и после ламинарного дросселя; измеряют температуру газа перед турбулентным дросселем и в междроссельном участке; по измеренным значениям перечисленных параметров вычисляют показатели физических свойств природного газа (плотность газа при стандартных условиях; коэффициент сжимаемости и динамическую вязкость газа при заданных значениях давления и температуры, показатель адиабаты, удельную теплоту сгорания газа), при этом перед турбулентным дросселем поддерживают абсолютное давление, обеспечивающее ламинарный режим течения газа в ламинарном дросселе при максимально-возможном числе Рейнольдса (Re=2300). Плотность природного газа при стандартных условиях вычисляют из формулы (вывод формулы приведен в Приложении 1):

где

_st - плотность газа при стандартных условиях, кг/м ³;

А - постоянный коэффициент, определяемый по экспериментальным данным;

Р₁ - абсолютное давление газа перед турбулентным дросселем, Па;

Р₂ - абсолютное давление газа в междроссельном участке, Па;

Р₃ - абсолютное давление газа в после ламинарного дросселя, Па;

T₁ - температура газа перед турбулентным дросселем, К;

Т₂ - температура газа в междроссельном участке, К;

- динамическая вязкость газа, Па·с;

K_z1 =Z₁/Z_st - коэффициент сжимаемости газа перед турбулентным дросселем;

K_z2=Z₂ /Z_st - коэффициент сжимаемости газа в междроссельном участке;

Z₁, Z₂, Z_st - фактор сжимаемости газа, соответственно, при рабочих и при стандартных условиях.

k - показатель адиабаты.

Значения коэффициентов сжимаемости K_z1 и K_z2 определяют по алгоритму, приведенному в ГОСТ 30319.2-96 (ГОСТ 30319.2-96. Газ природный. Методы расчета физических свойств. Определение коэффициента сжимаемости) в зависимости от плотности газа при стандартных условиях, абсолютного давления и температуры газа, т.е.

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

где - динамическая вязкость газа, мкПа·с;

Х_а - объемная доля азота (N₂);

Х_у - объемная доля диоксида углерода (CO₂).

Значение показателя адиабаты k определяют по формуле, приведенной в ГОСТ 30319.1-96 и имеющей вид

Значения высшей Н_с.в и низшей Н_с.н удельной теплоты сгорания определяют по формулам, приведенным в ГОСТ 30319.1-96 и имеющих вид

Значения Х_а (объемная доля азота (N₂ )) и Х_у (объемная доля диоксида углерода (CO₂ )) определяют по результатам хроматографического анализа.

Значение коэффициента "А" определяют по экспериментальным данным до начала реализации предлагаемого способа (процедура описана в Приложении 1) из формулы

где _st - плотность газа, измеренная одним из известных устройств (например, пикнометром).

Устройство для осуществления предложенного способа приведено на чертеже.

Устройство состоит из: линии контролируемого газа, включающей входной участок 1, междроссельный участок 2 и выходной участок 3, подключенной входом к источнику контролируемого газа (например, к газопроводу, сепаратору и др.), а выходом - к линии 4 утилизации газа; исполнительного механизма (регулирующего клапана) 5, установленного на входном участке 1 линии контролируемого газа; фильтра 6; турбулентного дросселя 7; ламинарного дросселя 8; датчиков абсолютного давления 9, 10, 11, установленных на линии контролируемого газа, соответственно, перед турбулентным дросселем 7, в междроссельном участке, после ламинарного дросселя; датчиков температуры 12, 13, установленных на входном участке перед турбулентным дросселем и в междроссельном участке перед ламинарным дросселем; вычислительного устройства 14, ко входам которого подключены перечисленные датчики, а к выходам - исполнительный механизм и устройство отображения информации (например, жидкокристаллический дисплей) 15. Вычислительное устройство содержит цифровой регулятор давления газа перед турбулентным дросселем 16, программный блок вычисления числа Рейнольдса 17, программный блок (цифровой регулятор) 18, программный блок вычисления _st, , K_z2 19.

Устройство работает следующим образом.

Контролируемый газ из, например, газопровода проходит через исполнительный механизм (регулирующий клапан) 5, фильтр 6, турбулентный 7 и ламинарный 8 дроссели в линию 4 утилизации газа. При помощи регулирующего клапана 5 во входном участке 1 автоматически поддерживается необходимое давление. Фильтр 6 осуществляет очистку газа от механических примесей и жидкости, чем предотвращается засорение дросселей. При прохождении газа через турбулентный дроссель давление Р₂ в междроссельном участке 2 понижается. Величина давления Р₂ зависит от плотности _st газа, а также от давления Р₁ перед турбулентным дросселем, давления Р₃ после ламинарного дросселя, температуры газа T₁, перед турбулентным дросселем, температуры Т₂ газа в междроссельном участке 2, а также от геометрических размеров обоих дросселей и коэффициента расхода турбулентного дросселя. Влияние геометрических размеров дросселей и коэффициента расхода турбулентного дросселя на значение давления P₂ определяется значением коэффициента "А", входящим в вышеприведенное уравнение для определения плотности газа при стандартных условиях. Расчетная формула для определения коэффициента "А" (вывод см. в Приложении 1) имеет вид

где - коэффициент расхода турбулентного дросселя;

F - площадь проходного сечения турбулентного дросселя, м²;

L - длина ламинарного дросселя (капилляра), м;

D _L - диаметр проходного канала ламинарного дросселя, м.

Однако значение коэффициента "А" определяют по экспериментальным данным (как описано в Приложении 1) и тем самым исключают необходимость определения геометрических размеров обоих дросселей и коэффициента расхода турбулентного дросселя для определения значения коэффициента "А" расчетным путем. Практически невозможно получить высокую точность определения коэффициента "А" расчетным путем из-за больших погрешностей определения геометрических размеров обоих дросселей и коэффициента расхода турбулентного дросселя. Поэтому значение коэффициента "А" определяют по экспериментальным данным для каждого устройства. Другими словами - требуется индивидуальная идентификация коэффициента "А" для каждого устройства. Это - неудобно, но другой альтернативы нет.

Давление P ₂ в междроссельном участке 2 является основным информативным параметром, зависящим от плотности газа. Другие параметры (Р ₁, Т₁, Т₂, Р₃) от плотности газа не зависят, но их неучет проводит к недопустимым погрешностям определения плотности газа. Если бы эти параметры (Р₁ , Т₁, Т₂, Р₃) были бы застабилизированы, то плотность газа можно было бы определять только по значению давления P₂. Однако стабилизация параметров с высокой точностью усложнила бы устройство. Проще их измерить и результаты измерения учесть при определении плотности газа. Поэтому в предлагаемом устройстве все параметры измеряются при помощи датчиков 9, 10, 11, 12, 13 и сигналы от этих датчиков подаются на соответствующие входы вычислительного устройства (контроллера) 14. Вычислительное устройство 14 решает две задачи.

Первая задача. Вычислительное устройство 14 по измеренным значениям параметров Р₁ , Р₂, Р₃, Т₁, Т₂ и заранее идентифицированному коэффициенту "А" вычисляет значение плотности _st газа при стандартных условиях и выводит его на устройство отображения информации 15 (например, на жидкокристаллический дисплей). Зная плотность газа, вычислительное устройство вычисляет также другие показатели физических свойств газа. Вычисление плотности и других вышеназванных показателей физических свойств газа осуществляет программный блок 19 вычислительного устройства 14. Вычисления показателей физических свойств газа проводятся из уравнений (1)-(7).

Вторая задача. Вычислительное устройство 14 совместно с исполнительным механизмом (регулирующим клапаном) 5 реализует функцию каскадной автоматической системы регулирования числа Рейнольдса для ламинарного дросселя. Функцию датчика числа Рейнольдса реализует программный блок 17, который вычисляет число Рейнольдса из уравнения (вывод см. в Приложении 1, формула (П.18)).

При этом значения плотности _st, динамической вязкости ₂ и коэффициента сжимаемости K_z2 газа подаются в программный блок 17 из программного блока 19, а значения Р ₂, Р₃, T₂ - вводятся от датчиков соответственно 10, 11 и 13. Значения D_L и L вводят вручную. Вычисленное значение числа Рейнольдса подается на вход цифрового регулятора (программного блока) 18 числа Рейнольдса. На задающий вход этого регулятора вводится (пользователем) заданное значение числа Рейнольдса (Re_z=2300). Цифровой регулятор 18 числа Рейнольдса совместно с датчиком 17 числа Рейнольдса образуют внешний контур каскадной автоматической системы регулирования числа Рейнольдса. Регулирующим воздействием этого контура регулирования является заданное значение давления газа P_1z перед турбулентным дросселем 7. При реализации ПИД-закона регулирования регулирующее воздействие вычисляется по рекуррентному выражению (см. Приложение 1)

где Re(i)=Re(i)-Re _z(i) - ошибка регулирования числа Рейнольдса на i-м шаге;

Re(i) - текущее значение числа Рейнольдса на i-м шаге;

Re_z(i) - заданное значение числа Рейнольдса на i-м шаге;

P_1z(i) - задание регулятору давления на i-м шаге (регулирующее воздействие регулятора внешнего контура регулирования);

q₀, q₁, q₂ - параметры настройки цифрового регулятора.

При малых тактах квантования параметры настройки q₀, q₁ , q₂ можно вычислять, используя параметры настройки К, T₁ и T_D аналогового ПИД-регулятора, по формулам:

q₀=K*(1+T₀/(2*T ₁)+T_D/T₀);

q₁=-K*(1+2*T _D/T₀-T₀/(2*T₁));

q₂=K*T_D/T₀;

где К - коэффициент передачи аналогового ПИД-регулятора:

T ₁ - постоянная интегрирования аналогового ПИД-регулятора;

Т_D - постоянная дифференцирования аналогового ПИД-регулятора;

Т₀ - такт квантования, с.

Из выражения (11) следует, что для вычисления заданного значения давления на i-м шаге P_1z(i) в памяти вычислительного устройства 14 (конкретнее, в программном блоке 18) необходимо хранить значение заданного значения давления P_1z(i-1) на предыдущем шаге и значения ошибки регулирования на i-1 и i-2 шагах.

Заданное значение P_1z(i) на i-м шаге подается на задающий вход цифрового регулятора (программного блока) 16 давления Р₁ газа. На другой вход этого регулятора подается текущее значение давления P₁ от датчика 9 давления. Цифровой регулятор (программный блок) 16 совместно с датчиком давления 10 и исполнительным механизмом (регулирующим клапаном) 5 образуют внутренний контур каскадной автоматической системы регулирования числа Рейнольдса. Цифровой регулятор 16 вычисляет разность между текущим и заданным значениями давления Р₁, т.е. ошибку регулирования, и в зависимости от ее величины и знака воздействует на исполнительный механизм 5 до тех пор, пока давление P₁ не станет равным заданному.

Наличие в каскадной системе автоматического регулирования внутреннего контура регулирования давления P₁ дает возможность заранее погасить внешние возмущения, например, возникшие при изменении давления газа в газопроводе, и тем самым улучшить качество процесса регулирования числа Рейнольдса.

Использование предлагаемого технического решения дает возможность упростить процедуру и повысить оперативность (выполнять в режиме реального времени) определения показателей физических свойств природного газа.

Приложение 1.

ОБОСНОВАНИЕ СПОСОБА ОПРЕДЕЛЕНИЯ ПОКАЗАТЕЛЕЙ ФИЗИЧЕСКИХ СВОЙСТВ ПРИРОДНОГО ГАЗА ПО ПАРАМЕТРАМ ИСТЕЧЕНИЯ ЧЕРЕЗ ТУРБУЛЕНТНЫЙ И ЛАМИНАРНЫЙ ДРОССЕЛИ

1. Связь между параметрами истечения и плотностью газа

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

По характеру течения газа в каналах дроссели делят на ламинарные и турбулентные [1]. Турбулентные дроссели характеризуются малыми отношениями длины канала дросселя к его диаметру и турбулентным режимом течения газа. Обычно в турбулентных дросселях отношение длины канала к его диаметру не превышает 10. Так как канал имеет малую длину, а скорость течения велика, то протекающий по дросселю газ не успевает обменяться теплом со стенками канала, и термодинамический процесс в дросселях такого типа считают адиабатическим [1].

Течение газа в турбулентных дросселях может происходить как с дозвуковыми (докритический режим), так и со звуковыми скоростями (надкритический режим). Режим истечения через турбулентный дроссель определяется величиной отношения давлений Р₁ и P₂ до и после него. Отношение давлений, при котором происходит переход от дозвуковой скорости к звуковой, называют критическим и обозначают (Р ₂/Р₁)_кр. Известно [1], что

где k - показатель адиабаты.

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

Рассмотрим режим докритического течения газа через турбулентный дроссель.

Массовый расход G_тд газа (в кг/с) через турбулентный дроссель рассчитывают по формуле [1]

где

- коэффициент расхода;

F - площадь отверстия дросселя, м²;

P₁ - абсолютное давление газа перед дросселем, Па;

Р₂ - абсолютное давление газа после дросселя, Па;

R - газовая постоянная, Дж/(кг К);

T₁ - температура газа перед дросселем, К;

K_z1=Z₁/Z_st - коэффициент сжимаемости газа;

Z₁, Z_st - фактор сжимаемости газа соответственно при рабочих и при стандартных условиях;

k - показатель адиабаты.

Зависимости для расчета факторов сжимаемости Z=Z( _st, P_st, T_st) и коэффициента адиабаты k=k( _st, P_st, T_st), где _st(кг/м³) - плотность газа при стандартных условиях, приведены соответственно в ГОСТ 30319.2-96 [2] и ГОСТ 30319.1-96 [3]. Как видим, коэффициент сжимаемости и показатель адиабаты зависят от плотности газа при стандартных условиях (Р _st=101325 Па, T_st=293,15 К, ГОСТ 30319.0-96.

В формулу (П.2) входит газовая постоянная, значение которой вычисляют по известному выражению

где R_u=8314,31 Дж/(кг К) - универсальная газовая постоянная;

М - относительная молекулярная масса газа, кг/кмоль.

При стандартных условиях плотность газа (кг/м ³) связана с относительной молекулярной массой газа известным соотношением

где V_st=24,0546 м³/кмоль - удельный мольный объем газа при стандартных условиях.

Из формул (П.4) и (П.5) получаем

где C_r=R_u/V_st=345,64 Дж/(м³·K).

Подставив в формулу (П.2) выражение (П.6), получим уравнение массового расхода газа через турбулентный дроссель

где K_z1=Kz( _st,P₁,T₁) и k=k( _st,P₁,T₁) - функции от плотности газа при стандартных условиях, рабочих давлений и температуры газа.

Перейдем к рассмотрению течения газа через ламинарный дроссель.

Массовый расход G_лд газа (в кг/с) через ламинарный дроссель(капилляр) определяют по формуле Пуазейля [1, стр.37]

где =3,14;

D_L - диаметр проходного сечения капилляра, м;

Р₂ - абсолютное давление газа перед ламинарным дросселем, Па;

Р₃ - абсолютное давление газа после ламинарного дросселя, Па;

- динамическая вязкость, Па·с;

L - длина ламинарного дросселя, м;

R - газовая постоянная, Дж/(кг·град);

Т₂ - температура газа, К;

К_z2 - коэффициент сжимаемости газа при P₂ и Т₂ .

Подставим выражение (П.6) в (П.8). Получим

или

Мы принимаем, что давление газа после турбулентного дросселя равно давлению перед ламинарным дросселем. Другими словами: мы принимаем, что потери давления в междроссельном участке линии контролируемого газа отсутствуют(или ими можно пренебречь в силу их исчезающе малой величины). В установившемся режиме массовый расход газа через турбулентный дроссель равен массовому расходу газа через ламинарный дроссель. Приравнивая левые части уравнений (П.7) и (П.10), получаем

где

Значения коэффициентов сжимаемости K_z1 и К _z2 определяют по алгоритму, приведенному в ГОСТ 30319.2-96 [2], в зависимости от плотности газа при стандартных условиях, абсолютного давления и температуры газа, т.е.

Значение показателя адиабаты k определяют по формуле, приведенной в ГОСТ 30319.1-96 [1] и имеющей вид

Динамическая вязкость природного газа и других газовых смесей, как показано в [3] и [5], однозначно зависит от плотности газа при стандартных условиях. В ГОСТ 30319.1-96 [3] приведена зависимость динамической вязкости природного газа от его плотности при стандартных условиях. При принятых выше условных обозначениях она имеет вид

где - динамическая вязкость газа, мкПа·с;

Х_а - объемная доля азота (N₂);

Х_y - объемная доля диоксида углерода (СО₂).

Для перевода динамической вязкости из мкПа·с, вычисленной по уравнению (П.17), в Па·с, необходимо ее значение разделить на 10 ⁶.

Согласно ГОСТ 30319.1-96 [3] формула (П.17) применима при давлениях газа до 0.5 МПа и в диапазоне температур 240-360 К. Погрешность определения вязкости газа в этом диапазоне не превышает для метана 1.0%. Природный газ, как правило, содержит в основном метан (95-100%). Поэтому ГОСТ 30319.1-96 рекомендует формулу (П.17) для определения динамической вязкости природного газа.

Уравнение (П1.11) устанавливает однозначную связь между плотностью _st газа при стандартных условиях и параметрами течения газа P₁, T₁, P₂, T₂ , P₃. В явном виде зависимость _st от указанных параметров выразить невозможно из-за сложности зависимостей коэффициентов сжимаемости, показателя адиабаты и динамической вязкости от плотности газа _st. Однако численным методом решить уравнение (П.11) относительно _st можно. Для этого необходимо знать значения Р ₁, Т₁, Р₂, Т₂, Р₃ и А.

Значения параметров Р₁, Т₁ , Р₂, Т₂ и Р₃ могут быть измерены инструментальными средствами, а значение коэффициента «А» (будем называть этот коэффициент постоянной устройства) может и должно быть определено по экспериментальным данным. Для этого проводят эксперимент, состоящий в том, что через последовательно соединенные турбулентный и ламинарный дроссели пропускают газ, измеряют параметры P₁, T₁, P₂, T₂, P ₃, определяют плотность _st газа при стандартных условиях (пикнометром или по компонентному составу газа) и из уравнения (П.11) определяют значение «А» постоянной устройства (например, численным методом половинного деления, на ПЭВМ). Эксперимент проводят несколько раз и принимают среднее значение коэффициента А. Вычисление постоянной устройства по формуле (П1.12) приводит к недопустимой погрешности из-за большой погрешности определения , F, D_L, L.

Из приведенного выше анализа вытекает способ определения плотности газа при стандартных условиях. Он состоит в следующем. Контролируемый газ пропускают через турбулентный и ламинарный дроссели, установленные последовательно на линии контролируемого газа; измеряют абсолютное давление Р₁ и температуру T₁ газа перед турбулентным дросселем, абсолютное давление Р₂ и температуру Т₂ газа в междроссельной линии, абсолютное давление Р₃ газа после ламинарного дросселя и по измеренным параметрам вычисляют плотность _st газа при стандартных условиях из уравнения (П.11), при этом используют значение коэффициента «А», полученное заранее.

2. Математическое моделирование процесса истечения газа через турбулентный и ламинарный дроссели

Основным информативным параметром, зависящим от плотности газа, является давление Р ₂ в междроссельном участке линии контролируемого газа. Значения таких параметров, как давление Р₁ перед турбулентным дросселем и давление Р₃ после ламинарного дросселя, предопределяются выбором пользователя и не являются информативными в смысле их зависимости от плотности газа. Это - возмущающие параметры. Их значения должны учитываться в процессе определения плотности газа. То же касается температур T₁ и Т ₂ газа. Их значения предопределяются выбором пользователя.

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

Для реализации этой цели была разработана математическая модель процесса течения газа через два последовательно установленные турбулентный и ламинарные дроссели и программа ROTULA 1.BAS, реализующая расчет зависимости давления газа в междроссельной камере от плотности газа. (Текст программы прилагается, Приложение 2 ).

Программа написана на алгоритмическом языке TURBOBASIC. Программа состоит из главной программы и одной подпрограммы. Главная программа реализует алгоритм интерактивного ввода исходных данных и вывода результатов расчета по одному из трех, задаваемых пользователем, направлений: на монитор, на печать, в файл. Главная программа осуществляет перевод параметров в нужные размерности, вычисляет (в цикле) значения плотности газа, значения давления в междроссельном участке линии контролируемого газа, значения массового расхода газа и критерия Рейнольдса для ламинарного дросселя. Подпрограмма вычисляет значения коэффициента сжимаемости газа (по алгоритму, приведенному в ГОСТ 30319.1-96 [2], показателя адиабаты и динамической вязкости газа, а также удельную теплоту сгорания газа.

В качестве исходных данных используются значения D_t, D_L, L, T ₁, T₂, P₁, P₃, а также начальное значение плотности газа _n и приращение плотности газа d. Различные значения формируются в цикле с заданным приращением плотности газа. Исходные данные вводят в ПЭВМ в интерактивном режиме (по запросу ПЭВМ). В результате работы программы для каждого значения плотности газа при стандартных условиях вычисляются значения давления Р ₂, массового расхода газа через турбулентный и ламинарный дроссели, динамической вязкости, коэффициентов сжимаемости газа, показателя адиабаты и числа Рейнольдса для ламинарного дросселя, а также коэффициента f=×F, что в формуле (П.7).

Число (критерий) Рейнольдса для ламинарного дросселя определяется по формуле

Эта формула получена следующим образом.

Известно [6], что

где v - скорость газа, м/с;

d_L - внутренний диаметр капилляра, м;

- плотность газа при рабочих условиях, кг/м³;

₂ - динамическая вязкость газа, Па·с.

Скорость газа в ламинарном дросселе

где F=d _L ²/4;

Q - объемный расход газа через ламинарный дроссель, м³/с;

G_лд - массовый расход газа через ламинарный дроссель, кг/с.

Плотность газа при рабочих условиях определяют по формуле

Подставив в (П.20) уравнения (П.8) и (П.21), получим

Подставив выражения (П.22) и (П.21) в уравнение (П.19), после несложных преобразований получим формулу (П.18) для расчета критерия Рейнольдса для ламинарного дросселя. При моделировании принимали: Избыточное давление перед турбулентным дросселем Р ₁=10000 мм в.ст. (Абсолютное давление вычисляли в программе). Абсолютное давление после ламинарного дросселя Р ₃=100000 Па (барометрическое давление).

Диаметр турбулентного D_t=0.2 мм, диаметр ламинарного дросселя D_L =0.3 мм.

Длина ламинарного дросселя L=100 мм.

Начальная плотность газа при стандартных условиях _n=0.45 кг/м³.

Приращение плотности газа d=0.05 кг/м³.

Количество шагов по плотности газа N=12.

Варьировали значение Р₁. Выводили на печать: для каждого значения плотности газа - давление Р₂ в междроссельной линии, массовый расход газа через турбулентный и ламинарный дроссели, число Рейнольдса для ламинарного дросселя.

Результаты моделирования показали, что давление P ₂ в междроссельном участке существенно зависит от плотности газа. При этом, чем больше давление Р₁ перед турбулентным дросселем, тем больше чувствительность. На фиг.2 приведены зависимости приращения избыточного давления dР₂ от плотности газа при различных давлениях перед турбулентным дросселем. Под приращением давления dР₂ понимали разность между давлением P ₂ при =0.5 кг/м³ и давлением P₂ при =1.05 кг/м³. Из фиг.2 следует, что зависимость приращения давления Р₂ от плотности газа в общем случае нелинейная. При изменении плотности газа от 0,55 кг/м³ до 1,05 кг/м³, т.е. на d=0.55 кг/м³ приращение давления dP₂ зависит от давления Р₁ газа перед турбулентным дросселем. Значения dP₂ при d=0.55 кг/м³ и различных P₁ приведены в таблице 1. Там же приведены значения чувствительности, под которой понимается отношение Kr=dP₂/d.

Таблица 1 Зависимость чувствительности Kr от давления P1 перед турбулентным дросселем
P1, мм вод.ст.	4000	6000	8000	10000
dP 2, мм вод.ст.	767.3	1144.3	1493.6	1811.3
Kr, мм вод.ст./кг/м 3	1395.1	2080.5	2715.6	3293.3

Из таблицы 1 следует, что для повышения чувствительности, а следовательно, и точности определения плотности газа давление перед турбулентным дросселем следует повышать. Однако на повышение давления Р₁ накладываются ограничения, связанные с необходимостью обеспечить ламинарный режим движения газа через ламинарный дроссель. На фиг.3, построенном на основании результатов моделирования, приведены зависимости числа Рейнольдса от плотности газа при различных давлениях P₁. Как видно, ламинарный режим (Re<=2300) движения газа в ламинарном дросселе при P ₁=10000 мм вод.ст. может быть обеспечен только при плотности газа меньше 0,62 кг/м³. При давлении Р ₁=4000 мм вод.ст. ламинарный режим движения газа может быть обеспечен во всем диапазоне изменения плотности газа (от 0,5 до 1,1 кг/м³). Однако при Р₁=4000 мм вод.ст. имеем самую низкую чувствительность. Отсюда следует важный вывод: с целью достижения максимальной чувствительности контроля плотности газа давление перед турбулентным дросселем необходимо поддерживать таким, чтобы число Рейнольдса для ламинарного дросселя было равно Re=2300.

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

Первый вариант: реализуют одноконтурную АСР с виртуальным датчиком числа Рейнольдса, ПИД-регулятором и исполнительным механизмом на линии контролируемого газа. Виртуальный датчик числа Рейнольдса представляет собой вычислительное устройство (контроллер), которое в реальном масштабе времени вычисляет число Рейнольдса по измеренным параметрам контролируемого газа. Вычисления осуществляют по уравнению (П.18), при этом при помощи соответствующих датчиков измеряют Р₂, Р ₃ и Т₂, a _st определяют по предлагаемому способу. Выходной сигнал виртуального датчика числа Рейнольдса вводят, например, в ПИД-регулятор как регулируемую величину. На задающий вход регулятора вводят заданное значение числа Рейнольдса (Re=2300). Рассогласование этих сигналов преобразуется по принятому закону (например, по ПИД-закону) в регулирующее воздействие на исполнительный механизм. В результате автоматически поддерживается заданное число Рейнольдса (Re=2300). Понятно, что вычислительные операции, связанные с реализацией виртуального датчика числа Рейнольдса и, например, ПИД-закона регулирования, могут выполняться одним вычислительным устройством (контроллером).

Второй вариант: реализуют каскадную АСР, имеющую внутренний и внешний контуры регулирования. Внутренний контур - обычная одноконтурная АСР давления Р₁, содержащая датчик давления газа (Р₁), автоматический регулятор (например, ПИД-регулятор) и упомянутый выше исполнительный механизм на линии контролируемого газа. Задание этой АСР (заданное значение давления P₁) вводится внешним контуром, содержащим виртуальный датчик числа Рейнольдса и, например, ПИД-регулятор числа Рейнольдса. Заданное значение числа Рейнольдса вводит пользователь. При отклонении текущего значения числа Рейнольдса от заданного значения внешний контур формирует выходной сигнал, который подается в качестве заданного значения на задающий вход регулятора давления. Последний поддерживает такое давление Р₁, при котором текущее значение числа Рейнольдса равно заданному (Re=2300). Второй вариант АСР обеспечивает лучшее качество регулирования, т.к. возникающие отклонения давления Р₁ (от внешних возмущений) заранее устраняются внутренним контуром регулирования.

Алгоритм работы цифрового ПИД-регулятора [7, стр.82] числа Рейнольдса имеет вид

где Re(i)=Re(i)-Rez(i) - ошибка регулирования числа Рейнольдса на i-м шаге;

Re(i) - текущее значение числа Рейнольдса на i-м шаге;

Rez(i) - заданное значение числа Рейнольдса на i-м шаге;

P ₁z(i) - задание регулятору давления на i-м шаге (регулирующее воздействие регулятора внешнего контура регулирования).

При малых тактах квантования параметры настройки q₀ , q₁, q₂ можно вычислять, используя параметры настройки К, T₁ и Т_D аналогового ПИД-регулятора [7], по формулам

q₀=K*(1+T₀/(2*T ₁)+T_D/T₀);

q₁=-K*(1+2*T _D/T₀-T₀/(2*T₁));

q₂=K*T_D/T₀;

где К - коэффициент передачи аналогового ПИД-регулятора:

T ₁ - постоянная интегрирования аналогового ПИД-регулятора;

Т_D - постоянная дифференцирования аналогового ПИД-регулятора;

Т₀ - такт квантования, с.

Аналогично, алгоритм работы цифрового ПИД-регулятора давления имеет вид

где U(i) - регулирующее воздействие регулятора давления на i-м шаге;

P ₁(i)=P₁(i)-P₁z(i) - ошибка регулирования давления Р₁ на i-м шаге;

P₁(i) - текущее значение давления Р₁ на i-м шаге;

P ₁z(i) - заданное значение давления Р₁ на i-м шаге.

3. Идентификация коэффициента "А"

Идентификацию коэффициента "А" проводят по экспериментальным данным.

Технология эксперимента описана выше. Из уравнения (П.11) следует, что

Значения коэффициентов сжимаемости K_z1 и K _z2 определяют по алгоритму, приведенному в ГОСТ 30319.2-96 [2], в зависимости от плотности газа при стандартных условиях, абсолютного давления и температуры газа.

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

Динамическую вязкость природного газа определяют по формуле, приведенной в ГОСТ 30319.1-96 [1], в зависимости от плотности газа при стандартных условиях и температуры газа.

Для вычисления коэффициента "А" на ПЭВМ была разработана программа ROTULA2.BAS. Программа написана на алгоритмическом языке TURBOBASIC. Текст программы прилагается, Приложение 2. Программа состоит из основной программы и подпрограммы. Основная программа осуществляет: ввод исходных данных в интерактивном режиме, выбор направления вывода результатов расчета (на дисплей, на принтер или в файл), распечатку исходных данных, вычисление абсолютного давления по избыточному давлению, вычисление коэффициента "А", вывод значения коэффициента "А" по одному из заданных пользователем направлений. Подпрограмма осуществляет вычисление: коэффициентов сжимаемости, показателя адиабаты и динамической вязкости газа.

Работа программы была проверена на конкретных данных. Рассмотрены два примера. В первом примере в качестве исходных данных принимали: Р₁=4000 мм Н₂ О, Р₃=100000 Па, T₁=300 К, Т₂ =300 К, =0.7, D_t=0.20 мм, D_L=0.3 мм, L=100 мм, X _a=0, X_y=0. Во втором: Р₁=4000 мм Н₂О, Р₃=100000 Па, Т₁=300 К, Т₂=300 К, =0.7, D_t=0.20 мм, D_L=0.3 мм, L=100 мм, Х _а=0, Х_y=0.

В первом примере плотность газа при стандартных условиях составляла =0.7 кг/м³, а избыточное давление в междроссельном участке - Р₂=2505.9 мм вод.ст., во втором примере: =0.8 кг/м³, а Р₂=3136.1 мм вод.ст. В первом примере получено значение коэффициента А=0.5827×10⁹ , во втором - А=0.5816×10⁹, т.е. практически одинаковые значения. Вычисление коэффициента А по формуле (П.12), которая была использована при моделировании процесса истечения газа, дало А=0.5817×10⁹. Таким образом, можно утверждать, что программа идентификации коэффициента "А" дает хорошие результаты.

4. Расчет плотности газа при стандартных условиях по параметрам истечения

Вычисление плотности газа при стандартных условиях, как указывалось выше, осуществляют из формулы (П.11)

где

При этом коэффициент А определяют по экспериментальным данным, как указывалось в предыдущем разделе.

Расчет _st из формулы (П.11) вручную практически невозможен из-за громоздкости вычисления K_z1, К_z2, k, и P_ot. Поэтому была разработана программа ROTULA3.BAS. Программа разработана на алгоритмическом языке TURBOBASIC. Текст программы прилагается, Приложение 2. Программа состоит из главной программы и одной подпрограммы. Главная программа осуществляет: ввод и распечатку исходных данных в интерактивном режиме, выбор направления вывода значений исходных данных и результатов расчета (на дисплей, принтер или в файл), пересчет вводимых значений избыточного давления в абсолютные давления, расчет и вывод на печать плотности газа при стандартных условиях, числа Рейнольдса и массовый расход газа через ламинарный дроссель. Подпрограмма осуществляет: расчет коэффициентов сжимаемости газа, вязкости газа в ламинарном дросселе, показателя адиабаты, высшей и низшей теплотворной способности газа. Все эти величины выводятся на печать. Плотность газа при стандартных условиях вычисляют методом последовательных приближений (методом половинного деления). Алгоритм вычисления плотности газа имеет вид

Ron=0.2

Rok=1.8

3 Ro=0.5*(Ron+Rok)

Call Zet(P1m, T1, Ro, Xa, Xy, Kz1, Mu1, Kad1, Hh1, Hl1)

Call Zet(P2m, T2, Ro, Xa, Xy, Kz2, Mu2, Kad2, Hh2, Hl2)

'Определение расчетного коэффициента Ar'

Pot=(P2a/P1a)^(2/Kad1)-(P2a/P1a)^((Kad1+1)/Kad1)

B1=P1a*(Mu2/1000000)*Kz2*T2/((P2a-P3)*(P2a+P3))

B2=SQR(Kad1*Pot/(Kz1*(Kad1-1)*T1))

B3=B1*B2

Ar=SQR(Ro)/B3

EA=Ar-A

If Abs(EA)<0.001*A then goto 1

End If

If EA<0 then

Ron=Ro goto 3 else Rok=Ro goto 3

End If

1 Дальнейший расчет и распечатка результатов расчета.

Алгоритм включает:

1. Расчет среднего (половинного) значения плотности газа по принятым начальному Ron и конечному Rok значениям плотности. (Принимали: на первом, шаге Ron=0.2, Rok=1.8, что перекрывает возможный диапазон изменения плотности природного газа).

2. Расчет Kz1, Kz2, Mu2, k по подпрограмме.

3. Расчет коэффициента Ar (Расчетное значение коэффициента "А") по формуле (П.28)

4. Определение невязки расчетного значения коэффициента Ar с фактическим А значением по формуле ЕА=Ar-А.

5. Проверку условия: Если абсолютное значение невязки ЕА меньше 0.001*А, то вычисленное в п.1 значение плотности газа принимается за решение задачи. В противном случае проверяется знак невязки. Если невязка отрицательная, то принимают Ron=Ro, иначе - Rok=Ro, и возращаются к п.1.

6. Процесс расчета продолжают до тех пор, пока невязка не станет меньше 0.001 А. Точность расчета плотности газа Ro можно сколь угодно повысить, приняв соответствующую величину допустимой невязки.

Программа ROTULA3.BAS была апробирована на двух конкретных примерах. В первом примере в качестве исходных принимали для =0.7 кг/м³ Р₁=4000 мм вод.ст., Р₂ =2505.9 мм вод.ст., T₁=Т₂=300 К, D _L=0.3 мм, L=100 мм и барометрическое давление Pb=10125 Па. Коэффициент "А" принимали по результатам идентификации (см. предыдущий раздел): А=0.5817х10⁹. В результате расчета получили плотность газа =0.6984 кг/м³, против 0.7 кг/м³, как должно быть. Такая сходимость результата говорит о приличной работоспособности программы. Кроме плотности газа были вычислены: динамическая вязкость =0.1111×10 ² мкПа·с газа в ламинарном дросселе, высшая Hh=38,55 МДж/м³ и низшая Hl=34,77 МДж/м³ теплотворные способности газа, а также показатель адиабаты Kad=1.2937 и коэффициенты сжимаемости газа K_z1=1.0521, К_z2=1.0524.

Во втором примере в качестве исходных принимали данные: Р₁=6000 мм Н₂О, Р2=3136 мм Н₂ О, Р₃=100000 Па, Т₁=Т₂=300 К, А=0.5817×10⁹, D_L=0.30 мм, L=100 мм для =0.8 кг/м³. В результате решения задачи получили =0.8 кг/м³, т.е. абсолютно точное значение. Таким образом, разработанная программа ROTULA3.BAS работоспособна и может быть использована для расчета плотности газа при стандартных условиях по параметрам истечения газа через последовательно установленные турбулентный и ламинарный дроссели.

5. Влияние содержания азота и диоксида углерода на точность определения плотности природного газа

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

В качестве примера рассматривали случай истечения газа через два дросселя при исходных данных: Р₁=6000 мм вод.ст., Р₃=100000 Па, Т ₁=Т₂=300 К, А=0.5817×10⁹, D _L=0.3 мм, L=100 мм, Pb=101325 Па, Х_a=Х_y =0%. В этом случае давление газа в междроссельном участке Р₂=3136.1 мм вод.ст и получили плотность газа _st=0,8 кг/м³. При тех же исходных данных плотность газа определяли при различных значениях Х_a и Х_y. Результаты решения задачи приведены в таблице

Ха, %	0	2	4	6	8
Х y, %	0	2	4	6	8
st, кг/м3	0.8000	0.8219	0.8438	0.8687	0.8922

Как видно, неучет содержания азота и диоксида углерода при определении плотности газа может привести к существенной погрешности. Эта погрешность проявится в том случае, когда идентификацию коэффициента "А" проводили на природном газе, не содержавшем азота и диоксида углерода, а определяли плотность газа, содержащего эти компоненты. Отсюда следуют правила: 1. Идентификацию коэффициента "А" следует проводить на том газе, плотность которого будет определяться. 2. При определении плотности газа значения Х_а и Х_y следует вводить в вычислительное устройство по результатам, например, хроматографического анализа.

Источники информации

1. Дмитриев В.Н., Градецкий В.Г. Основы пневмоавтоматики. М.: Машиностроение, 1973, 360 с.

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

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

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

5. Голубев И.Ф., Гнездилов Н.Е. Вязкость газовых смесей. М.: Изд-во "Стандартов". - 1971. 327 с.

6. Кремлевский П.П. Расходомеры и счетчики количества. Изд.3-е, перераб. и доп. Л.: Машиностроение (Ленинградское отделение), 1975. 776 с. с ил.

7. Изерман Р. Цифровые системы управления. Пер.с англ. - М.: Мир, 1984. - 541 с., ил.

Приложение 2.

ПРОГРАММЫ:

1. ROTULA1.BAS: "РАСЧЕТ ДАВЛЕНИЯ ГАЗА В МЕЖДРОССЕЛЬНОЙ КАМЕРЕ В ЗАВИСИМОСТИ ОТ ПЛОТНОСТИ ПРИРОДНОГО ГАЗА. ДРОССЕЛИ: ТУРБУЛЕНТНЫЙ-ЛАМИНАРНЫЙ"

2. ROTULA2.BAS: "ИДЕНТИФИКАЦИЯ КОЭФФИЦИЕНТА "А" ПО ПАРАМЕТРАМ ИСТЕЧЕНИЯ ГАЗА ЧЕРЕЗ ТУРБУЛЕНТНЫЙ И ЛАМИНАРНЫЙ ДРОССЕЛИ"

3. ROTULA3.BAS: "РАСЧЕТ ПЛОТНОСТИ ГАЗА ПО ПАРАМЕТРАМ ИСТЕЧЕНИЯ ЧЕРЕЗ ТУРБУЛЕНТНЫЙ И ЛАМИНАРНЫЙ ДРОССЕЛИ"

'ПРОГРАММА ROTULA1.BAS'

'РАСЧЕТ ДАВЛЕНИЯ ГАЗА В МЕЖДРОССЕЛЬНОЙ КАМЕРЕ'

'В ЗАВИСИМОСТИ ОТ ПЛОТНОСТИ ПРИРОДНОГО ГАЗА'

'ДРОССЕЛИ: ТУРБУЛЕНТНЫЙ - ЛАМИНАРНЫЙ'

DEFSNG P, T, E, V, R, D, G, Z, K, A, B, C, M, X

DEFINT i, j, n

77 cls

print" Введите избыточное давление"

print" перед турбулентным дросселем (мм вод.ст.)"

print" P1=";

input P1

print" Введите абсолютное давление после лам.др. (Па)"

print" P3=";

input Р3

print" Введите температуру газа перед турб.дрос. (К)"

print" T1=";

input T1

print" Введите температуру газа перед ламин.дрос. (К)"

print" T2=";

input T2

print" Введите диаметр турбулентного дросселя (мм)"

print" Dt=";

input Dt

print" Введите коэффициент расхода турбулентного дросселя (мм)"

print" E=";

input E

print" Введите диаметр ламинарного дросселя(мм)"

print" DL=";

input DL

print" Введите длину ламинарного дросселя (мм)"

print" L=";

input L

print" Введите начальную плотность газа"

print" при стандартных условиях(кг/м3)"

print" Ron=";

input Ron

print" Введите шаг изменения плотности газа"

print" DRo=";

input DRo

print" Введите количество шагов изменения плотности газа"

print" N=";

input N

print" Введите мольную концентрацию CO2 в газе(%)"

print" Ху=";

input Xy

print" Введите мольную концентрацию азота в газе(%)"

print" Xa=";

input Xa

'ВЫБОР НАПРАВЛЕНИЯ ВЫВОДА РЕЗУЛЬТАТОВ'

print "Задайте направление вывода результатов расчета"

print "1 - дисплей"

print "2 - печать"

print "3 - файл REZTULA.TXT"

print "NAP=";

input NAP

print "Если будете корректировать исходные данные,"

print "задайте kor=1, иначе - kor=0"

print "kor⁼";

input kor

if kor=0 goto 78

goto 77

78 if NAP = 1 then viw$ = "CON"

if NAP = 2 then viw$ = "PRN"

if NAP = 3 then viw$ = "REZTULA.TXT"

open viw$ for output as #2

print #2,

print #2," ПРОГРАММА ROTULA 1.BAS"

print #2,

print #2," РАСЧЕТ ЗАВИСИМОСТИ ДАВЛЕНИЯ ГАЗА В"

print #2," МЕЖДРОССЕЛЬНОЙ КАМЕРЕ ОТ ЕГО ПЛОТНОСТИ"

print #2," (ДРОССЕЛИ: ТУРБУЛЕНТНЫЙ - ЛАМИНАРНЫЙ)"

print #2,

print #2," ИСХОДНЫЕ ДАННЫЕ"

print #2,

a1$="P1=##### мм Н2O P3=###### Па T1=###.# К Т2=###.# E==#.##"

print #2, using a1$; P1, P3, T1, T2, E

a2$="Dt=#.## мм DL=#.## мм L=###.# мм Xa=#.#### мол.%"

print #2, using a2$; Dt, D1, L, Xa

a3$="Ron=#.#### кг/м3 DRo=#.#### кг/м3 N=## Xy=#.#### мол.%"

print #2, using a3$; Ron, DRo, N, Xy

print #2,

print #2," РЕЗУЛЬТАТЫ РАСЧЕТА"

print #2,

print #2,

"Расчет коэффициента истечения газа через дроссель'

Ef=E*3.14*(Dt*0.001)^2/4

Р1=Р1*9.80665 'Перевод мм Н2O в Па'

Р1а=Р1+101325

DL=0.001*DL 'Перевод мм в м'

L=0.001*L 'Перевод мм в м'

'Печать результатов расчета'

a4$="Et=+.####^^^^

print #2, using a4$; Ef

print #2,

print #2," i Ro	P2	Gt	GL	Re"
print #2," - кг/м3 мм Н2O		г/с	г/с	-"

print #2,

P1m=P1a/(10^6) 'Перевод Па в МПа'

i=0

17 i=i+1

Ro=Ron+i*DRo

Р2n=Р1a

P2k=P3

1 P2=0.5*(P2n+P2k)

P2m=P2/(10^6) 'Перевод Па в МПа'

Call Zet(P1m, Т1, Ro, Xa, Xy, Kz1, Mu1, Kad1, Hh1, Hl1)

Call Zet(P2m, T2, Ro, Xa, Xy, Kz2, Mu2, Kad2, Hh2, Hl2)

'Определение критерия Рейнольдса(ламинарный дроссель)'

Re1=DL^Ro*(P2^2-P3^2)

Re2=22121*Kz2*T2*(Mu2/1000000^2*L

Re=Re1/Re2

'Определение расхода газа через турбулентный дроссель'

Pot=(P2/P1a)^(2/Kad1)-(P2/P1a)^((Kad1+1)/Kad1)

В1=SQR(2*Kad1*Ro*Pot/(345.64*Kz1*T1*(Kad1-1)))

Gt=Ef*Р1а*В1

'Определение расхода газа через ламинарный дроссель'

В5=3.14*(DL^4)*Ro*(P2^2-P3^2)

B6=88483*Mu2*L*Kz2*T2/1000000

GL=B5/B6

EpsG=Gt-GL

If Abs(EpsG)<0.00000001 then goto 5

else

goto 6

End if

6 If EpsG<0 then

P2n=P2 goto 1 else P2'k=P2 goto 1

End if

5 Gt=(10^3)Gt 'Перевод кг/с в г/с'

GL=(10^3)*GL

P21=(P2-101325)/9.80665 'Перевод Па в мм вод.ст.'

'Печать результатов расчета'

а5$="## #.### #####.# +.####^^^^ +.####^^^^ #####"

print #2, using a5$; i, Ro, P21, Gt, GL, Re

If i<N goto 17

12 print #2,

print #2, "Расчет закончен"

End

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

'по модифицированному уравнению состояния GERG-91 мод,'

ГОСТ 30319.2-96, а также динамической вязкости,'

'показателя адиабаты и теплотворной способности газа'

'по ГОСТ 30319.1-96'

Sub Zet(P, T, Ro, Xa, Xy, K, Mu, Kad, Hh, H1)

"Размерности: Р-МПа, Т-град.К, Ro-кг/м3, Ха, Ху-мол.%'

'Расчет фактора сжимаемости при стандартных условиях'

Zc=1-(0.0741*Ro-0.006*Xa-0.0575*Xy)

'Расчет молярной доли эквивалентного углеводорода'

Хе=1-Ха-Ху

'Расчет молярной массы эквивалентного углеводорода'

Me=(24.05525*Zc*Ro-28.0135*Xa-44.01*Xy)/Xe

'Расчет показателя Н'

Н=128.64+47.479*Ме

'Расчет коэффициентов уравнения состояния'

В11=(8.77118/10^-5.56281*T/10^6+8.8151*T^2/10^9)*Н В12=(-

8.24747/10^7+4.31436*Т/10^9-6.08319*T^2/10^12)*H^2 В1=-0.425468+2.865 *Т/10^3-4.62073*T^2/10^6+B11+B12

B2=-0.1446+7.4091*Т/10^4-9.1195*T^2/10^7 В23=-0.339693+1.61176*T/10^3-2.04429*T^2/10^6 В3=-0.86834+4.0376*Т/10^3-5.1657*T^2/10^6

С11=(6.46422/10^4-4.22876*T/10^6+6.88157*T^2/10^9)*Н С12=(-

3.32805/10^7+2.2316*T/10^9-3.67713*T^2/10^12)*H^2 С1=-0.302488+1.95861*T/10^3-3.16302*T^2/10^6+C11+С12

C2=7.8498/10^3-3.9895*T/10^5+6.1187*T^2/10^8 C3=2.0513/10^3+3.4888*T/10^5-8.3703*T^2/10^8-C223=5.52066/10^3-1.68609*T/10^5+1.57169*T^2/10^8

C233=3.58783/10^3+8.06674*T/10^6-3.25798*T^2/10^8

Bz=0.72+1.875*(320-T)^2/10^5

Cz=0.92+0.0013*(T-320)

Bm1=Xe^2*B1+Xe*Xa*Bz*(B1+B2)-1.73*Xe*Xy*(B1*B3)^0.5

Bm2=Xa^2*B2+2*Xa*Xy*B23+Xy^2*B3

Bm=Bm1+Bm2

Cm1=Xe^3*C1+3*Xe^2*Xa*Cz*(C1^2*C2)^(1/3)

Cm2=2.76*Xe^2*Xy*(C1^2*C3)^(1/3)+3*Xe*Xa^2*Cz*(C1*C2^2)^(1/3)

Cm3=6.6*Xe*Xa*Xy*(C1*C2*C3)^(1/3)+2.76*Xe*Xy^2*(C1*C3^)^(1/3)

Cm4=Xa^3*C2+3*Xa^2*Xy*C223+3*Xa*Xy^2*C233+Xy^3*C3

Cm=Cm1+Cm2+Cm3+Cm4

'Расчет коэффициентов уравнения для вычисления Z'

В=1000*Р/(2.7715*Т)

C0=B^2*Cm

В0=В*Bm

А1=1+В0

А0=1+1.5*(В0+С0)

A2=(A0-(A0^2-A1^3)^0.5)^1/3)

'Расчет фактора сжимаемости'

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

'Расчет коэффициента сжимаемости'

K=Z/Zc

'Расчет динамической вязкости'

'Расчет псевдокритического давления и температуры'

Ppk=2.9585*(1.608-0.05994*Ro+Xy-0.392*Xa)Tpk=88.25*(0.9915+1.759*Ro-Xy-1.681*Ха)

'Расчет приведенного давления и температуры'

Ppr=P/Ppk

Tpr=T/Tpk

'Расчет динамической вязкости при давлении до 0.5 МПа'

Mu1=T^0.5+1.37-9.09*Ro^0.125

Mu2=Ro^0.5+2.08-1.5*(Ха+Ху)

Mu=3.24*Mu1/Mu2 'Размерность - мкПа.с'

'Расчет динамической вязкости при давлении 0.5 МПа и больше'

If P>0.5 then

Cmu=1+Ppr^2/30*(Tpr-1) goto 9

else

goto 2

End if

9 Mu=Cmu*Mu 'Размерность - мкПа.с'

'Расчет удельной объемной теплоты сгорания' '(теплотворной способности) природного газа'

'Размерность - МДж/куб.м'

2 Hh=92.819*(0.51447*Ro+0.05603-0.65689*Xa-Xy) 'высшая'

H1=85.453*(0.5219*Ro+0.04242-0.65197*Xa-Xy)'низшая'

'Расчет показателя адиабаты при Т=240-360 К и Р до 10 МПа'

Kad1=1.556*(1+0.074*Ха)-3.9*Т*(1-0.68*Ха)/10^4-0.208*Ro

Kad2=(p/T)^1.43*(384*(1-Xa)*(p/T)^0.8+26.4*Xa)

Kad=Kad1+Kad2

End sub

'ПРОГРАММА ROTULA2.BAS'

'ИДЕНТИФИКАЦИЯ КОЭФФИЦИЕНТА "А" ПО ПАРАМЕТРАМ ИСТЕЧЕНИЯ'

'ГАЗА ЧЕРЕЗ ТУРБУЛЕНТНЫЙ И ЛАМИНАРНЫЙ ДРОССЕЛИ'

DEFSNG P, T, E, V, R, D, G, Z, K, A, B, C, M, X

DEFINT i, j, n

77 cls

print" Введите избыточное давление"

print" перед турбулентным дросселем (мм вод.ст.)"

print" P1=":

input P1

print" Введите избыточное давление"

print" между дросселями (мм вод.ст.)"

print" P2=";

input P2

print" Введите абсолютное давление после лам.дрос.(Па)"

print" P3=":

input Р3

print" Введите барометрическое давление(Па)"

print" Pb=";

input Pb

print" Введите температуру газа перед турб.дрос.(К)"

print" T1=";

input T1

print" Введите температуру газа перед лам.дрос.(К)"

print" Т2=";

input T2

print" Введите плотность газа при ст.условиях Ro (кг/м3)"

print" Ro=";

input Ro

print" Введите мольную концентрацию CO2 в газе(%)"

print" Ху=";

input Xy

print" Введите мольную концентрацию азота в газе(%)"

print" Ха=";

input Xa

'ВЫБОР НАПРАВЛЕНИЯ ВЫВОДА РЕЗУЛЬТАТОВ'

print "Задайте направление вывода результатов расчета"

print "1 - дисплей"

print "2 - печать"

print "3 - файл REZTULA2.TXT"

print "NAP=";

input NAP

print "Если будете корректировать исходные данные,"

print "задайте kor=1, иначе - kor=0"

print "kor=";

input kor

if kor = 0 goto 78

goto 77

78 if NAP = 1 then viw$ = "CON"

if NAP = 2 then viw$ = "PRN"

if NAP - 3 then viw$ = "REZTULA2.TXT"

open viw$ for output as #2

print #2,

print #2," ПРОГРАММА ROTULA2.BAS"

print #2,

print #2," ИДЕНТИФИКАЦИЯ КОЭФФИЦИЕНТА А ПО ПАРАМЕТРАМ"

print #2," ИСТЕЧЕНИЯ ГАЗА ЧЕРЕЗ ТУРБУЛЕНТНЫЙ И ЛАМИНАРНЫЙ ДРОССЕЛИ"

print #2,

print #2," ИСХОДНЫЕ ДАННЫЕ"

print #2,

a1$="P1=#### мм H2O P2=#### мм H2O P3=###### Па T1=###.# K"

print #2, using a1$; P1, P2, P3, T1

a2$="T2=###.# Xa=#.#### мол.% Ху=#.#### мол.%"

print #2, using a2$; T2, Xa, Xy

a3$="Ro=#.#### Pb=######.# Па"

print #2, using a3$; Ro, Pb

print #2,

print #2, " РЕЗУЛЬТАТЫ РАСЧЕТА"

print #2,

P2=P2*9.80665 'Перевод мм Н2O в Па'

Р2а=Р2+Pb 'Вычисл.абсолютного давления'

Р1=Р1*9.80665 'Перевод мм Н2O в Па'

P1a=P1+Pb 'Вычисл.абсолютного давления'

P1m=P1a/(10^6) 'Перевод Па в МПа'

P2m=P2a/(10^6)

Call Zet(P1m, Т1, Ro, Xa, Xy, Kz1, Mu1, Kad1, Hh1, Hl1)

Call Zet(P2m, T2, Ro, Xa, Xy, Kz2, Mu2, Kad2, Hh2, Hl2)

'Определение коэффициента А'

Pot=(P2a/P1a)^(2/Kad1)-(P2a/P1a)^((Kad1+1)/Kad1)

В1=P1a*(Mu2/1000000)*Kz2*T2/((P2a-P3)*(P2a+P3))

B2=SQR(Kad1*Pot/(Kz1*(Kad1-1)*T1))

B3=B1*B2

A=SQR(Ro)/B3

a5$=" Коэффициент А=+.####^^^^"

print #2, using a5$; A

print #2,

print #2, "Расчет закончен"

End

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

'по модифицированному уравнению состояния GERG-91 мод,'

'ГОСТ 30319.2-96, а также динамической вязкости,'

'показателя адиабаты и теплотворной способности газа'

'по ГОСТ 30319.1-96'

Sub Zet(P, T, Ro, Xa, Xy, K, Mu, Kad, Hh, Hl)

'Размерности: Р-МПа, Т-град.К, Ro-кг/м3, Ха, Ху-мол.%'

'Расчет фактора сжимаемости при стандартных условиях'

Zc=1-(0.0741 *Ro-0.006*Xa-0.0575*Xy)

'Расчет молярной доли эквивалентного углеводорода'

Хе=1-Ха-Ху

'Расчет молярной массы эквивалентного углеводорода'

Me=(24.05525*Zc*Ro-28.0135*Xa-44.01*Xy)/Xe

'Расчет показателя Н'

Н=128.64+47.479*Ме

'Расчет коэффициентов уравнения состояния'

В11=(8.77118/10^4-5.56281*T/10^6+8.8151*T^2/10^9)*H

В12=(-8.24747/10^7+4.31436*T/10^9-6.08319*T^2/10^12)*H^2

В1=-0.425468+2.865 *Т/10^3-4.62073*T^2/10^6+B11+В12

В2=-0.1446+7.4091*Т/10^4-9.1195*T^2/10^7

В23=-0.339693+1.61176*Т/10^3-2.04429*T^2/10^6

В3=-0.86834+4.0376*Т/10^3-5.1657*T^2/10^6

С11=(6.46422/10^4-4.22876*T/10^6+6.88157*T^2/10^9)*H

C12=(-3.32805/10^7+2.2316*T/10^9-3.67713*T^2/10^12)*H^2

С1=-0.302488+1.95861*T/10^3-3.16302*T^2/10^6+C11+C12

С2=7.8498/10^3-3.9895*T/10^5+6.1187*T^2/10^8

C3=2.0513/10^3+3.4888*T/10^5-8.3703*T^2/10^8

C223=5.52066/10^3-1.68609*T/10^5+1.57169*T^2/10^8

C233=3.58783/10^3+8.06674*T/10^6-3.25798*T^2/10^8

Bz=0.72+1.875*(320-T)^2/10^5

Cz=0.92+0.0013*(T-320)

Bm1=Xe^2*B1+Xe*Xa*Bz*(B1+B2)-1.73*Xe*Xy*(B1*B3)^0.5

Bm2=Xa^2*B2+2*Xa*Xy*B23+Xy^2*B3

Bm=Bm1+Bm2

Cm1=Xe^3*C1+3*Xe^2*Xa*Cz*(C1^2*C2)^(1/3)

Cm2=2.76*Xe^2*Xy*(C1^2*C3)^(1/3)+3*Xe*Xa^2*Cz*(C1*C2^2)^(1/3)

Cm3=6.6*Xe*Xa*Xy*(C1*C2*C3)^(1/3)+2.76*Xe*Xy^2*(C1*С3^2)^(1/3)

Cm4=Xa^3*C2+3*Xa^2*Xy*C223+3*Xa*Xy^2*C233+Xy^3*C3

Cm=Cm1+Cm2+Cm3+Cm4

'Расчет коэффициентов уравнения для вычисления Z'

В=1000*Р/(2.7715*Т)

C0=B^2*Cm

В0=В*Bm

А1=1+В0

А0=1+1.5*(В0+С0)

A2=(A0-(A0^2-A1^3)^0.5)^(1/3)

"Расчет фактора сжимаемости'

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

'Расчет коэффициента сжимаемости'

K=Z/Zc

'Расчет динамической вязкости'

'Расчет псевдокритического давления и температуры'

Ppk=2.9585*(1.608-0.05994*Ro+Xy-0.392*Xa)

Tpk=88.25*(0.9915+1.759*Ro-Xy-1.681*Xa)

'Расчет приведенного давления и температуры'

Ppr=P/Ppk

Tpr=T/Tpk

'Расчет динамической вязкости при давлении до 0.5 МПа'

Mu1=T^0.5+1.37-9.09*Ro^0.125

Mu2=Ro^0.5-2.08-1.5*(Ха+Ху)

Mu=3.24*Mu1/Mu2 'Размерность - мкПа.с'

'Расчет динамической вязкости при давлении 0.5 МПа и больше'

If P>0.5 then

Cmu=1+Ppr^2/30*(Tpr-1) goto 9

else

goto 2

End if

9 Mu=Cmu*Mu 'Размерность - мкПа.с'

'Расчет удельной объемной теплоты сгорания'

'(теплотворной способности) природного газа'

'Размерность - МДж/куб.м'

2 Hh=92.819*(0.51447*Ro+0.05603-0.65689*Xa-Xy) 'высшая'

Hl=85.453*(0.5219*Ro+0.04242-0.65197*Xa-Xy) 'низшая'

'Расчет показателя адиабаты при Т=240-360 К и Р до 10 МПа'

Kad1=1.556*(1+0.074*Xa)-3.9*T*(1-0.68*Xa)/10^4-0.208*Ro

Kad2=(p/T)^1.43*(384*(1-Xa)*(p/T)^0.8+26.4*Xa)

Kad=Kad1+Kad2

End sub

'ПРОГРАММА ROTULA3.BAS '

'РАСЧЕТ ПЛОТНОСТИ ГАЗА ПО ПАРАМЕТРАМ ИСТЕЧЕНИЯ'

'ЧЕРЕЗ ТУРБУЛЕНТНЫЙ И ЛАМИНАРНЫЙ ДРОССЕЛИ'

DEFSNG P, T, E, V, R, D, G, Z, K, A, B, C, M, X

DEFINT i, j, n

77 cls

print" Введите избыточное давление"

print" перед турбулентным дросселем(мм вод.ст.)"

print" P1=";

input P1

print" Введите избыточное давление"

print" между дросселями (мм вод.ст.)"

print" P2=";

input P2

print" Введите абсолютное давление после лам.дрос.(Па)"

print" P3=";

input Р3

print" Введите барометрическое давление(Па)"

print" Pb=";

input Pb

print" Введите температуру газа перед турб.дрос.(К)"

print" T1=":

input T1

print" Введите температуру газа перед лам.дрос.(К)"

print" T2=":

input T2

print" Введите коэффициент А"

print" A=";

input A

print" Введите мольную концентрацию CO2 в газе(%)"

print" Ху=":

input Xy

print" Введите мольную концентрацию азота в газе(%)"

print" Xa=":

input Xa

print" Введите диаметр ламинарного дросселя(мм)"

print" DL=":

input DL

print" Введите длину ламинарного дросселя(мм)"

print" L=";

input L

'ВЫБОР НАПРАВЛЕНИЯ ВЫВОДА РЕЗУЛЬТАТОВ'

print "Задайте направление вывода результатов расчета"

print "1 - дисплей"

print "2 - печать"

print "3 - файл REZTULA3.TXT"

print "NAP=";

input NAP

print "Если будете корректировать исходные данные,"

print "задайте kor=1, иначе - kor=0"

print "kor=":

input kor

if kor=0 goto 78

goto 77

78 if NAP = 1 then viw$ = "CON"

if NAP = 2 then viw$ = "PRN"

if NAP = 3 then viw$ = "REZTULA3.TXT"

open viw$ for output as#2

print #2,

print #2," ПРОГРАММА ROTULA3.BAS"

print #2,

print #2," РАСЧЕТ ПЛОТНОСТИ ГАЗА ПО ПАРАМЕТРАМ ИСТЕЧЕНИЯ "

print #2," ЧЕРЕЗ ТУРБУЛЕНТНЫЙ И ЛАМИНАРНЫЙ ДРОССЕЛИ"

print #2,

print #2," ИСХОДНЫЕ ДАННЫЕ"

print #2,

a1$="P1=#### мм H2О P2=#### мм H2O P3=###### Па Т1=###.# К"

print #2, using a1$; P1, P2, P3, T1

a2$="T2=###.# A=+.####^^^^ Xa=#.#### мол.% Ху=#.#### мол.%"

print #2, using a2$; T2, A, Xa, Xy

a3$="DL=#.### mm L=###.## мм Pb=######.# Па"

print #2, using a3$; DL, L, Pb

print #2,

print #2," РЕЗУЛЬТАТЫ РАСЧЕТА"

print #2,

P2=P2*9.80665 'Перевод мм Н2O в Па'

Р2а=Р2+Pb 'Вычисл.абсолютного давления'

Р1=Р1*9.80665 'Перевод мм Н2O в Па'

Р1а=Р1+Pb 'Вычисл.абсолютного давления'

Р1m=P1a/(10^6) 'Перевод Па в МПа'

P2m=P2a/(10^6)

Ron=0.2

Rok=1.8

3 Ro=0.5*(Ron+Rok)

Call Zet(P1m, T1, Ro, Xa, Xy, Kz1, Mul, Kad1, Hh1, Hl1)

Call Zet(P2m, T2, Ro, Xa, Xy, Kz2, Mu2, Kad2, Hh2, Hl2)

'Определение расчетного коэффициента Ar'

Pot=(P2a/P1a)^(2/Kad1)-(P2a/P1a)^((Kad1+1)/Kad1)

В1=Р1a*(Mu2/1000000)*Kz2*T2/((P2a-P3)*(P2a+P3))

B2=SQR(Kad1*Pot/(Kz1*(Kad1-1)*T1))

B3=B1*B2

Ar=SQR(Ro)/B3

EA=Ar-A

If Abs(EA)<0.001*A then goto 1

End If

If EA<0 then

Ron=Ro goto 3 else Rok=Ro goto 3

End If

'Определение критерия Рейнольдса(ламинарный дроссель)'

1 DL=0.001*DL

L=0.001*L

Re1=DL^3*(P2a^2-P3^2)*Ro

Re2=22121 *Kz2*T2*(Mu2/1000000)^2*L

Re=Re1/Re2

'Определение расхода газа через ламинарный дроссель'

B5=3.14*(DL^4)*Ro*(P2a^2-P3^2)

B6=88483.84*Mu2*L*Kz2*T2/1000000

GL=B5/B6

GL=1000*GL 'Перевод кг/с в г/с'

а5$=" Плотность газа Ro =#.#### кг/м3 Re=##### GL=+.####^^^^ г/с" print #2, using a5$; Ro, Re, GL

а8$=" Вязкость газа Mu=+.####^^^^ мкПа·с"

print #2, using a8$; Mu2

а7$=" Теплотворная способность газа Hh=###.## H1=###.## МДж/м3" print #2, using a7$; Hh1, Hl2

a6$="Kz1=#.#### Kz2=#.#### Kad=#.####"

print #2, using a6$; Kz1, Kz2, Kad1

print #2,

print #2, "Расчет закончен"

End

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

'по модифицированному уравнению состояния GERG-91 мод,'

'ГОСТ 30319.2-96, а также динамической вязкости,'

'показателя адиабаты и теплотворной способности газа'

'по ГОСТ 30319.1-96'

Sub Zet(P, T, Ro, Xa, Xy, K, Mu, Kad, Hh, Hl)

'Размерности: Р-МПа, Т-град.К, Ro-кг/м3, Ха, Ху-мол.%'

'Расчет фактора сжимаемости при стандартных условиях'

Zc=1-(0.0741*Ro-0.006*Xa-0.0575*Xy)

'Расчет молярной доли эквивалентного углеводорода'

Хе=1-Ха-Ху

'Расчет молярной массы эквивалентного углеводорода'

Me=(24.05525*Zc*Ro-28.0135*Ха-44.01*Ху)/Хе

'Расчет покачателя Н'

Н=128.64+47.479*Ме

'Расчет коэффициентов уравнения состояния'

В11=(8.771 18/10^4-5.56281*Т/10^6+8.8151*-Т^2/10^9)*Н В12=(-

8.24747/10^+4.31436*Т/10^9-6.08319*T^2/10^12)*H^2 В1=-0.425468+2.865*Т/10^3-

4.62073*T^2/10^6+B11+B12

В2=-0.1446+7.4091*T/10^4-9.1195*T^2/10^7 В23=-0.339693+1.61176*Т/10^3-

2.04429*T^2/10^6 В3=-0.86834+4.0376^T/10^-5.1657*T^2/10^6

C11=(6.46422/10^4-4.22876*T/10^6+6.88157*T^2/10^9)*H C12=(-

3.32805/10^7+2.2316*T/10^9-3.67713*7^2/10^12)*H^2 С1=-0.302488+1.95861*T/10^3-

3.16302*T^2/10^6+C11+C12

С2=7.8498/10^3-3.9895*T/10^5+6.1187*T^2/10^8 С3=2.0513/10^+3.4888*T/10^5-

8.3703*T^2/10^8 C223=5.52066/10^3-1.68609*T/10^5+1.57169*T^2/10^8

C233=3.58783/10^3+8.06674*T/10^6-3.25798*T^2/10^8

Bz=0.72+1.875*(320-T)^2/10^5

Cz=0.92+0.0013*(T-320)

Bm1=Xe^2*B1+Xe*Xa*Bz*(B1+B2)-1.73*Xe*Xy*(B1*B3)^0.5

Bm2=Xa^2*B2+2*Xa*Xy*B23+Xy^2*B3

Bm=Bm1+Bm2

Cm1=Xe^3*C1+3*Xe^2*Xa*Cz*(C1^2*C2)^(1/3)

Cm2=2.76*Xe^2*Xy*(C1^2*C3)^(1/3)+3*Xe*Xa^2*Cz*(C1*C2^2)^(1/3)

Cm3=6.6*Xe*Xa*Xy*(C1*C2*C3)^(1/3)+2.76*Xe*Xy^2*(C1*C3^2)^(1/3)

Cm4=Xa^3*C2+3*Xa^2*Xy*C223+3*Xa*Xy^2*C233+Xy^3*C3

Cm=Cm1+Cm2+Cm3+Cm4

'Расчет коэффициентов уравнения для вычисления Z'

В=1000*Р/(2.7715*Т)

C0=B^2*Cm

B0=B*Bm

А1=1+В0

А0=1+1.5*(В0+С0)

A2=(A0-(A0^2-A1^3)^0.5)^1/3)

"Расчет фактора сжимаемости'

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

'Расчет коэффициента сжимаемости'

K=Z/Zc

'Расчет динамической вязкости'

'Расчет псевдокритическою давления и температуры'

Ppk=2.9585*(1.608-0.05994*Ro+Xy-0.392*Xa)Tpk=88.25*(0.9915+1.759*Ro-Xy-1.681*Ха)

'Расчет приведенного давления и температуры'

Ppr=P/Ppk

Tpr=T/Tpk

'Расчет динамической вязкости при давлении до 0.5 МПа'

Mu1=T^0.5+1.37-9.09*Ro^0.125

Mu2=Ro^0.5+2.08-1.5*(Ха+Ху)

Mu=3.24*Mu1/Mu2 'Размерность - мкПа.с'

'Расчет динамической вязкости при давлении 0.5 МПа и больше'

If P>0.5 then

Cmu=1+Ppr^2/30*(Tpr-1) goto 9

else

goto 2

End if

9 Mu=Cmu*Mu 'Размерность - мкПа.с'

'Расчет удельной объемной теплоты сгорания'

'(теплотворной способности) природного газа'

'Размерность - МДж/куб.м'

2 Hh=92.819*(0.51447*Ro+0.05603-0.65689*Xa-Xy) 'высшая'

H1=85.453*(0.5219*Ro+0.04242-0.65197*Xa-Xy)'низшая'

'Расчет показателя адиабаты при Т=240-360 К и Р до 10 МПа'

Kad1=1.556*(1+0.074*Ха)-3.9*Т*(1-0.68*Ха)/10^4-0.208*Ro

Kad2=(p/T)^1.43*(384*(1-Xa)*(p/T)^0.8+26.4*Xa)

Kad^Kad1+Kad2

End sub