Идентификация фактической производительности центробежных нагнетателей газоперекачивающих агрегатов в реальном времени Текст научной статьи по специальности «Математика»

УДК 681.327.12.001.362; 519.7; 519.81

ИДЕНТИФИКАЦИЯ ФАКТИЧЕСКОЙ ПРОИЗВОДИТЕЛЬНОСТИ ЦЕНТРОБЕЖНЫХ НАГНЕТАТЕЛЕЙ ГАЗОПЕРЕКАЧИВАЮЩИХ АГРЕГАТОВ В РЕАЛЬНОМ ВРЕМЕНИ

АДАМЕНКОВ.А., АДАМЕНКО А.В., ТЕВЯШЕВА О.А.

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

1. Введение

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

2. Математическая модель ЦБН

Согласно [1], основными термодинамическими характеристиками ЦБН являются: зависимости степе-

ни сжатия є = —-, политропического коэффициента

Рн

полезного действия ц и относительной приведенной

( N 1

внутренней мощности

н Упр

ной производительности Qnp

от приведенной объем-(м3/мин) и приведен-

ного относительного числа оборотов

( n ^

n0

, где Рн ,

пр

Рк — давления газа на входе и на выходе ЦБН (кгс/см2 ); no — номинальное число оборотов на ЦБН (об/мин); n— обороты привода на ЦБН

(об/мин). Степень сжатия є = 1

( n '

Q

( n ^

пр

V no7

при

n0

v - ~'пР7 можно представить в виде многочлена

пр

второй степени [2]:

(Q пр д)_ ao + aiQ пр + a2Q

пр

2

(1)

Функции ц^пр) и I — I (qпр) хорошо приближа-

пр

ются многочленами третьей степени [2]:

p(Qпр) = do + diQпр + d2Qпр2 + d3Qпр3 , (2)

^пр) _ с0 + c1Qпр + c2Qпр + c3Q

пр . (3)

пр

Формулы приведения имеют вид [2]:

( n ^

V n0 7

пр

n

n0

Z пр R пр Тпр ' Z н R н Тн

' А1 J Д0

Д н J пр I n Jr н ,

Q пр = ^ у 0

ZH R н Тн

Рн•104

q • 106 1440 ,

(4)

кгс • м

где Zпр , Rпр (-к), Тпр (К)— приведенные зна-

кг • к

чения коэффициента сжимаемости, газовой посто-

кгс • м

янной и температуры газа; Zн , Rн (--—), Тн

кг • K

(К)— коэффициент сжимаемости газа, газовая постоянная и температура на входе ЦБН; у н — удельный вес газа перед ЦБН (кгс/м3 ); N— внутренняя мощность, потребляемая ЦБН (кВт); q—коммерческий расход газа на ГПА (млн. м3/сут); у0 — удельный

вес газа в нормальных условиях (кгс/м3 ).

В выражениях (1)—(3) a0 , a1, а2, С0 , q , С2 , С3, d0,

d1, d2 , d3 — коэффициенты аппроксимации соответствующих функций.

Таким образом, основными термодинамическими характеристиками ЦБН являются приведенные характеристики (1)—(3).

Формула для пересчета є

Q

пр

n0

пр 7

при измене-

нии числа оборотов на ЦБН и физических параметров газа имеет вид [2]:

є =

1 +

( n ^

2 / m-1

V n0 J

пр

-1

m

m-1

(5)

где m— показатель политропы.

Важной характеристикой течения газа в нагнетателе является также температура газа. Абсолютная температура газа на выходе нагнетателя Тк (К) определяется через температуру на его входе Тн (К) [2]:

n

m

0

44

РИ, 1999, № 4

m -1

Tk = ТнЄ• (6)

При идентификации производительности ЦБН ГПА необходимо учитывать технологические ограничения, которые задаются в виде следующих неравенств:

Qnp.min — Qпр — Qпр.тах , (7)

n ■ < n < n P. < pmax T < Tmax (8)

где Qnp.min, Qnp.max — минимально и максимально допустимые значения приведенной объемной производительности (м3/мин); nmin, nmax — минимальная и максимальная частота вращения вала нагнетателя

(об/мин); pmax—максимальное давление нагнетателя, определяемое прочностью труб (кгс/см2 ); Tmax — ограничение сверху на температуру газа на выходе ЦБН, зависящее от свойств изоляционного покрытия (К).

Для постановки и решения задачи идентификации производительности ЦБН ГПА в качестве основных уравнений и неравенств модели ЦБН будем использовать выражения (5)—(8). В целях расширения области изменения переменных модели ЦБН, обеспечения монотонности зависимости по каждой из этих переменных, повышения устойчивости модели и обеспечения условий ее разрешимости эти уравнения и неравенства необходимо модифицировать таким образом, чтобы в области изменения реальных значений переменных модели исходные и модифицированные уравнения и неравенства модели совпадали. Модифицированную систему уравнений и неравенств (5)—(8) модели ЦБН представим в виде:

MГПА - { є -

1 +

( n ^

V no7

2 / m-1 Л

р m _ і є0 1

пр

m

m-1

, (9)

f D Л

m-1

Tk = Th

p

Vін у

Qnp.min — Qnp — Qnp.max ,

(10)

(11)

m

k

k _ k0 . 1 | dCp k0 -1

k -1 k0 -1 I R k0 ,

k — показатель адиабаты;

1

ZCp • (1 + X p), (13)

k0

k0 -1

= 5.15+

(5.65 + 0.017 • tcp) A 1.987

k0 — показатель “изоэнтропы” газа в идеальном состоянии; Zcp = —(Z h + Zk) — средний приведенный коэффициент сжимаемости газа; Zh = z(ph ,Th)—

коэффициент сжимаемости газа на входе ЦБН; Zk = z(Pk ,Tk) — коэффициент сжимаемости газа на

выходе ЦБН; А = ^н — относительная плотность ^ ’ 1.206

газа по воздуху; рн — плотность сухого газа в нормальном состоянии (кг / м3);

Qnp = — У 0

n Zh • Rн • InTmax(Tj q „

n0„ н н Tmin' h' q m2

1440

102

P

н

Tmin, Tmax — минимально и максимально допустимые значения температуры газа (К);

z(p,t)=1 -

(InP'"ax(p) - ^1 ^

' pmin 'I 102

0.345 Д 0.446

10

Л

+ 0.015

: 1.3 - 0.0144 InTmax(T- 283.2

' ' Tmin 4 7 "

— коэффициент сжимаемости газа; Pmin , Pmax — минимально и максимально допустимые значения

давления газа (кгс/см2 );

x, xmin < x < xmax,

In xmax xmin

W=\

xmin, x — xmin, xmax, x — xmax,

— функция проецирования точки на область;

nmin * n < nmax , Pk * Pmx , Tk < Tmax } , (12)

p, 2

где Є = —k ; Є0 = a0 + a1Qnp + a2Qnp ;

m =

k p

k (p -1) +1 ;

2 3

p = d0 + d1Q np + d2Q np + d3Q np ;

2

V n0 7

np

ZnpR npTnp ( n У

■^пр^пр Апр

Z н • r н ■ -nTmaxfr.)

v n0 7

если Zh InTm-fe) > ZminTm

2

V n0 у

, в ост. случаях.

— приведенные обороты ЦБН;

2

н

РИ, 1999, № 4

45

Zmin — минимально допустимое значение коэффициента сжимаемости;

X =

(

1.23 + 0.12 P,

пр.ср

T

V пРсР

2

Л

0.061

P

пр. ср

TZ

* пр.ср cp

Измеренные значения давлений, температур, числа оборотов представим в виде:

Рн = Рнист +^Рн , Pk = Ркист +^Pk , Тн = Тнист +^Тн , Тк = Ткист +^Тк ,

— коэффициент изобарического сжатия газа;

n = nист n , (14)

dCp

"Г"

Рпр.ср ' (246 + 0.12 Рпр.ср)

Т

3

пр.ср

— средний приведенный коэффициент теплоемкости газа;

где Рнист , Р^, Тист, T^, nист — истинные значения давлений, температур и числа оборотов, а |рн ,

^Рк , ^Тн , ^Тк , |п — ошибки измерений соответствующих переменных.

Рпр .ср - 2 (Рпр (Рн ) + Рпр (Рк )) — среднее приведенное давление;

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

Тпр .ср - 2 (Тпр (Тн ) + Тпр (Тк ))

1пр Цк J

- средняя приведенная температура;

^Рн ~N(0,оРн), ^Рк~Ц0,а2Рк ),

^Тн ~N(0,4н), ^Tk~N(0,4к ), |n~N(0,стП)

Рпр (Р)

Р +1,033 Рпк

— приведенное давление;

Тпр (Т)

Т

—— — приведенная температура;

Тпк

Рпк = 30.168 X

х (0.05993 (26.831 -рн) + NCO2 -0.392 NN2 )

Принято полагать, что ошибки измерений давлений, температур, числа оборотов статистически независимы друг от друга [3]. В этом случае совместная функция плотности распределения ошибок будет равна произведению функций плотности распределения ошибок переменных модели. Функция максимального правдоподобия с учетом статистических свойств ошибок переменных модели ЦБН и с учетом выражения (14) примет вид:

— псевдокритическое давление; NCO2 , NN2 — молярные концентрации углекислого газа и азота в транспортируемом газе в долях единицы;

Тпк = 88.25 х

х(1.7591(0.56364 + рн)-NCO2 -1.681NN2)

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

tcp = |(Тн + Тк) - 273.15

— среднее значение температуры.

3. Метод идентификации фактической производительности ЦБН ГПА

Задача идентификации производительности ЦБН ГПА сводится к решению системы уравнений и неравенств модели ЦБН (9)—(12) относительно неизвестной переменной q. Все остальные переменные модели ЦБН считаются известными. Система (9)— (12) будет переопределенной. Поэтому ее решение можно найти только в статистическом смысле.

Пусть имеются наборы измеренных значений давлений, температур и числа оборотов на ЦБН: Рн , Рк ,

Тн , Тк , П .

6

f = (2*)-2 -ст^н 'стТ1н 'стТк •стй1

Рн Рк Тн Тк

(

х exp

V

(Рн - Рн У

2 ст2

>\ ( • exp

У

V

1

н

(Рк - Рк)2

2 а2

Рк

1

X

х

/

(

х exp

V

1

(Тн - Тн У

2 ст2

>\ ( • exp

У

V

н

(Тк - Тк)2

2 а2

Тк

1

X

х

/

f

х exp

V

1 (п - n> С

2 стП )

(15)

Прологарифмируем выражение (15) и получим логарифмическую функцию максимального правдоподобия:

lnF = -31п(2л)-lnстрн -lnCTPk -lnсттн -1 (Рн - Рн У _ 1 (Рк - РкГ

- ln СТТк - ln Стn —

2 СТР Рн

2 »Рк

ЙТн - Т, 1 _ I (Тк - Тк? _ 1 (п - n)2

2 стП .

2 стТ

Тн

2 стТк

46

РИ, 1999, № 4

Функции F и lnF достигают экстремума при одних и тех же значениях аргументов (следствие монотонно возрастающего характера функции lnF), что позволяет формулировать исходную задачу как задачу максимизации функции lnF или как задачу минимизации функции (-lnF). В этом случае задача идентификации фактической производительности ЦБН ГПА принимает вид:

(рн ~ Рн f , (Pk ~ PkF , (тн - TH У , (тк - Tkf

<j±

а:

Pk

а:

Тн

а:

Тк

н

(и - и)2

Рн,Рк,Тн,Tk,n,q eQ

->min

(16)

где область ограничений q описывается уравнениями и неравенствами модели ЦБН (9)—(12).

Задача (16) является задачей условной минимизации. Переход от задачи условной минимизации к задаче безусловной минимизации можно осуществить, если

выразить переменную Pk из уравнения (9), а переменную Тк из уравнения (10) и подставить полученные аналитические выражения в (16). Но найти

аналитические выражения для Pk и Tk в явном виде нельзя. Однако, если пренебречь зависимостями

показателя адиабаты k от переменных Pk , Tk в выражениях (9)—(10), то тогда можно получить аналитические выражения для Pk из уравнения (9)

и для Tk из уравнения (10). С учетом сказанного, предлагается специальный алгоритм решения задачи условной минимизации вида (16).

Задаются начальные приближения по переменным Pн , Тн, n, равные их измеренным значениям (P<i0) = Pн , Ti0) = Тн , n(0) = n ), задается начальное приближение по переменной q (q(0)) и вычисляется значение для показателя адиабаты k (k(0)) по формуле (13) с учетом значений P^), Ti0), n(0), Pk , Tk . i-я итерация алгоритма решения задачи условной

минимизации состоит в следующем (i = 1, КИ , где КИ— количество итераций):

1. Находятся аналитические выражения для Pk из

уравнения (9) и для Tk из уравнения (10) в предположении, что показатель адиабаты k является постоянной величиной (значение для k получается на предыдущей итерации: k = kO-i)).

решении задачи безусловной минимизации не учитываются, их необходимо проверить после ее решения). Переменными полученной задачи безусловной минимизации являются: Pн , Тн, n, q.

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

полученное на предыдущей итерации:

P(i-1)

-Тн

T(1_1)

тн

q(M)

T(1) т н ,

, n(l 1). Пусть этим решением является: P^), q(1), n(l).

4. По аналитическим выражениям, полученным в п.1, и с учетом значений P^), Тн^, q(l), n(l) вычисляются значения для Pk , Tk : P^(l), т(^ .

5. По найденным значениям

P(1)

Pн

T(1)

тн

q(1), n(l),

Pkl), Tkl) вычисляется значение для показателя адиабаты k по формуле (13): ((і) .

6. Проверяется критерий выхода из условной минимизации. Если он не выполняется— осуществляется переход к п.1. Если он выполняется, то итерационный процесс завершен.

Таким образом, решение задачи условной минимизации вида (16) сводится к решению последовательности задач безусловной минимизации. При этом пренебрегаются зависимости показателя адиабаты k от переменных P( , T( в выражениях (9), (10) на каждой итерации, т.е. показатель адиабаты k считается равным константе на каждой итерации (значение его на каждой итерации берется из предыдущей итерации). Такой подход является вполне обоснованным, поскольку показатель адиабаты k слабо

зависит от Pk , Tk .

После того, как получено решение задачи (16), необходимо провести его анализ. Обозначим решение задачи (16) P*, P(*, Т* , т(, n*, q* .

Если выполняются условия

P

- Pн

<5

Тн) max ,

Pk* - Pk <

5

(Pk)

max ,

It*-T I<8(Tli)

Ан Ан — °max ,

T* - тк <smTx),

n - n <5

(n)

(17)

2. Найденные аналитические выражения для P( и Т( подставляются в функцию цели (16); таким образом, ограничения на равенства (9), (10) исключаются и, соответственно, исключаются переменные P( , Т( . В результате получается задача безусловной минимизации (ограничения на неравенства (11), (12) при

где sm>alx), 8mPakx), Smsix , , 8max — максимальные

погрешности измерений начального и конечного давлений, начальной и конечной температуры, числа оборотов соответственно, то считается, что ошибок в измерениях нет и что модель ЦБН ГПА является адекватной. Если же хотя бы одно из условий не

РИ, 1999, № 4

47

выполняется, то тогда возможны два варианта: либо были ошибки в процессе передачи измерений параметров на вход предлагаемого алгоритма или другие ошибки, которые привели к неадекватности входных данных (тогда нужно проверить правильность ввода исходных данных и, в случае ошибки, попытаться заново решить задачу); либо модель ЦБН ГПА не является адекватной, т.е. присутствуют ошибки модели и, возможно, некоторые параметры модели нуждаются в корректировке (но это уже отдельная задача).

4. Пример решения задачи идентификации фактической производительности ЦБН ГПА

|Тн* -Тн| = 0.706 < S^aX , |Tk* - Tk| = 0.444 < S^aX ,

|n* - n| = 342 -10-7 <5®) .

Условия (17) выполняются. Следовательно, ошибок в процессе передачи измерений параметров на вход предлагаемого алгоритма нет и модель ЦБН ГПА является адекватной.

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

5. Заключение

Приведем результаты решения задачи идентификации фактической производительности ЦБН ГПА. При решении этой задачи были использованы следующие значения переменных модели: Qmin = 150

м3/мин , Qmax = 450 м3/мин , Zmin = 10_1° ,Rпр = 50 кгс • м кгс • м

■, R н = 49

Тпр = 288 К,

Z пр = 0.91;

кг • К ’ н кг • К

У0 =0.70511 кгс/м3 , Tгр = 278 К, рн =0.7236 кг /м3 ,

Nco2 = 0.003 , Nn2 =0.044, a0 = 1.21226 ,

a1 = 0.00084532, a2 = -2.589934-Ш”6 , d0 = 0.4546676,

ф = 0.002372 , d2 =-7.69645-10-7, d3 = -8.18505-10“9 . Измеренные значения давлений, температур, числа

оборотов: рн = 46 кгс/см2 , Pj = 53.1752 кгс/см2 , Tн = 288 K, Tk = 298.051 K, її = 4320 . Дисперсии

давлений, температур, числа оборотов: стрн = 0.09 ,

CTpk = 0.1681, стТн = 0.14 , ст^ = °.°9 , CTn = 857.10-7 . Максимальные погрешности измерений давлений,

Предложен эффективный метод идентификации фактической производительности центробежных нагнетателей газоперекачивающих агрегатов. Модифицированная математическая модель ЦБН с расширенной областью допустимых решений обеспечивает непрерывность и монотонность зависимостей по каждой из переменных модели. Результаты решения задачи идентификации фактической производительности ЦБН ГПА подтверждают возможность использования предлагаемого метода в системах технической диагностики реального времени.

Литература: 1. Альбом характеристик центробежных нагнетателей природного газа. М., 1985. 86 с. 2. Евдокимов А.Г., Тевяшев А.Д. Оперативное управление потокораспределением в инженерных сетях. Харьков: Изд-во при Харьк. ун-те издательского объединения “Вища школа”, 1980. 144с. 3. Тевяшев А., Козыренко С., Адаменко А. Статистически устойчивая идентификация состояния модели стационарного режима транспорта газа в магистральном газопроводе// Транспортування, контроль якості та облік енергоносіїв. Львів: “Львівська політехніка”. 1998. С. 48—57.

Поступила в редколлегию 20.12.99 Рецензент: д-р техн. наук Евдокимов А.Г.

температур, числа оборотов: S^aX = 0.6 кгс/см2 ,

smax=0.82 кгс/см2, sgax*=0.75 к sgax)=0.6 к

5ma)x = 7 • Алгоритм решения задачи условной минимизации вида (16), приведенный в п.2, был программно реализован. Результаты решения задач, возникающих на каждой итерации условной минимизации, приведены в таблице.

Как видно из таблицы, итерационный процесс быстро сходится и результатом решения задачи условной минимизации вида (16) будет:

Адаменко Вера Анатольевна, канд. техн. наук, ассистент кафедры прикладной математики ХТУРЭ. Научные интересы: системный анализ, математическое моделирование, компьютерная графика, методы оптимизации, численные методы, математическое программирование. Увлечения и хобби: спортивный бридж, настольный теннис, компьютерные игры, плавание. Адрес: Украина, 61166, Харьков, пр. Ленина,14, тел. (0572) 40-94-36.

Адаменко Андрей Викторович, научный сотрудник кафедры прикладной математики ХТУРЭ. Научные интересы: системный анализ, математическое моделирование, компьютерная графика, методы оптимизации, численные методы, математическое программирование. Увлечения и хобби: спортивный бридж, настольный теннис, компьютерные игры, шахматы. Адрес: Украина, 61166, Харьков, пр. Ленина, 14, тел. (0572) 40-

Итерационный процесс решения задачи условной минимизации

№ итерации, i Pj(i) , кгс/см2 Т® , К q(i) , млн. м3/сут n(i) P® , кгс/см2 Т® , К

1 45,5811 288,7 20,6817 4320 53,8403 297,611

2 45,5775 288,706 20,6565 4320 53,8458 297,607

3 45,5776 288,706 20,6572 4320 53,8456 297,608

4 45,5776 288,706 20,6572 4320 53,8456 297,608

Рн = 45.5776 кгс/см2 , pJ = 53.8456 кг^см2, Ті* = 288.706 К, Tk* = 297.608 К, n* = 4320 , q* = 20.6572 млн. м3/сут .

Проверяются условия (17):

ІРн* -Рн| = 0.422 < 8rpaX), |Pk* -Pkl = 0.671 <smax ,

94-36.

Тевяшева Ольга Андреевна, студентка 5 курса ХГПУ. Научные интересы: системный анализ и теория оптимальных решений. Увлечения и хобби: горные лыжи, туризм, музыка. Адрес: Украина, 61176, Харьков, ул. Велозаводская, 38, кв.38, тел. (0572) 11-26-73.

48

РИ, 1999, № 4