Математическое моделирование производственных процессов АГНКС Текст научной статьи по специальности «Механика и машиностроение»



10 лет

журналу

Математическое моделирование производственных процессов АГНКС

1А.А. Евстифеев, начальник лаборатории ООО «Газпром ВНИИГАЗ»,

доцент Национального исследовательского ядерного университета «МИФИ», к.т.н.

43

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

__Ключевые слова:

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

дним из применяемых при проектировании методов выбора оптимальной конструкции АГНКС является принцип объемного газового баланса, заключающийся в формировании единого непрерывного газового канала с постоянной или переменной производительностью.

В ряде работ и патентов, опубликованных научными коллективами проектных и научно-исследовательских организаций [1-28], предлагается применять концепцию заправки передвижных автомобильных газовых заправщиков (ПАГЗ) и крупных потребителей на АГНКС по остаточному принципу. Эта концепция сформулирована в работах [15-18]. Суть ее заключается в применении на АГНКС новых методов и средств организации процесса выдачи сжатого газа на заправочные колонки. Наиболее подробно такой подход описан в работе [16]. Исходя из необходимости соблюдения уравнения баланса природного газа в развитие этого подхода были предложены новая схема коммутации газораздаточных линий на АГНКС, алгоритмы управления компрессорным оборудованием станции и технические решения по организации динамического секционного аккумулятора газа с регулировкой по объему.

В процессе принятия решения по наиболее эффективной конструкции АГНКС необходимо выполнить следующие задачи:

• сформировать эквивалентную обобщенную структурную схему производственного объекта с выделением ключевых элементов и линий коммутации;

• построить интегральную математическую модель объекта производства компримированного природного газа (КПГ), используя формализованные

\m\ww\Y

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

• провести моделирование режимов работы для предлагаемых конструктивных решений станции;

• выполнить процедуру принятия решения, используя принципы эффективности.

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

В общем виде схема большинства действующих АГНКС представлена на рисунке.

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

Математическая модель рабочего процесса линии нагнетания с использованием многоступенчатых поршневых компрессоров

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

Основы системы дифференциальных уравнений, характеризующих изменение параметров газа в отдельных полостях ступеней, описаны и проанализированы в работе Б.С. Фотина [1]. В данной работе при рассмотрении процессов, связанных с линией нагнетания, из-за изменения параметров в межступенчатых коммуникациях, соединяющих следующие друг за другом ступени, и сложности процессов, происходящих в межступенчатых коммуникациях, будут приняты следующие допущения:

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

• движение газа в трубопроводах и аппаратах коммуникаций одномерное и дозвуковое;

• коэффициенты трения газа о стенки и местных сопротивлений, определенные в стационарном потоке, справедливы и для нестационарного течения;

Обобщенная модель классической АГНКС с общей линией нагнетания и общей линией заправки

• коэффициенты теплоотдачи, определяемые из критериальных уравнений, полученных для стационарного потока, справедливы и для нестационарного течения;

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

компрессоре приведены ниже.

1. Уравнение неразрывности потока сжимаемой жидкости:

Эр Эс Эр с ЭА Л

—+ р — + с — + р--= 0,

Эт Эх Эх А Эх

(1)

где р - плотность газа; т - время; с - скорость; х - координата вдоль оси трубопровода; А - площадь поперечного сечения трубопровода.

2. Уравнение изменения количества движения с учетом гидравлического трения и местных сопротивлений:

Эс Эс 1 Эр _

— + с— +--- = -Ф,

Эт Эх р Эх

(2)

( ° Л

где Ф= А, + £- (с | с 112D); X - коэффициент гидравлического трения;

V А*)

£ - коэффициент местного сопротивления, Д - гидравлический диаметр.

Потери давления на трение Ар = Хр1с2 I 2D и местное сопротивление Др5 = £р/с21 2.

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

3. Уравнение энергии:

Аэ = £ DАx.

_ du d

д + с Ф =-+ р-

d т d т

т

V р J

(3)

Здесь и - удельная внутренняя энергия газа; сФ - подвод тепла трения к газу; д - удельный тепловой поток, который определяется по формуле

= аПие^х (Гст - Т) = аПпер (Тст - Т) рAdx рА

где а - коэффициент теплоотдачи; Ппер - периметр поперечного сечения трубы; Гст - температуры стенки (определяется по эмпирическим формулам: для трубопроводов всасывания Гст = 1,1 Гн; для трубопроводов нагнетания Тст = 0,9 Гн; по длине холодильника температура изменяется по линейному закону). Значение а определяется из критериальных уравнений типа

10 лет

журналу

45

№д = 0 ,1 3 Яе 33 Ог° '1 (при ламинарном течении), N11 = 0,0 18 Яе 0 '8 (при турбулентном течении), где Яе - определяется по мгновенным скоростям; Яе, Ог - критерии Нуссельта, Рейнольдса, Грасгофа.

Числа Яе и величины а в нестационарном потоке отличаются от квазистационарного из-за пульсации давлений и скоростей газа. Однако для практических расчетов

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

4. Уравнение состояния:

р = рЯТ . (4)

Система уравнений (1)...(4) включает в себя нелинейные дифференциальные уравнения с частными производными. Для ее решения использовался конечно-разностный метод Куранта, Изаксона и Риса, или, как его еще называют, метод «пространственных разностей вперед и назад». Разности во времени берутся вперед, а разности по пространству для каждого уравнения выбираются вперед и назад согласно наклону соответствующей характеристики. Этот метод требует приведения уравнений к характеристическому виду. Система уравнений в характеристической форме запишется в следующем виде:

дх

для — = с + а

дг

др др

— + (с + а)— + ра дт дх

да

для — = с - а

дг

др др -=Г + (с - а)^--ра дт дх

дх

для — = с

дг

дс дс

— + (с + а)— дт дх

дс дс

— + (с + а)— дт дх

= (д + сФ)р(к -1) - раФ - ра:

= (д + сФ)р(к -1) + раФ - ра:

сёЛ ; Лек'

сёЛ ; Лех

^ + с Зр-а2 Г ^ + с ^ дт дх ^ дт дх

= (д + сФ)р(к -1);

где р = р / ЯТ.

Если система дифференциальных уравнений в частных производных принадлежит к гиперболическому типу и записана в характеристическом виде

с + и^х) = в,

то она аппроксимируется разностной схемой.

В конечных разностях уравнение колебаний в коммуникациях примет вид

Рк,1+1 рк ,1 , / , ч рк ,1 рк-1,1 , „

■ ■ -+(ск ^+ак ^—^+рк а,,

Ат

А х

ск,¡+1 ск,1 + ( + ) ск,1 ск-1,1 "■Ч6 к л^ик ,1!'

Ат

СиЛ, . - Л,

Ах

= [д+ск ,фк,, ] рк,, (к -1) - рк а ,фк - - рк— ~—^;

Л, . А х

к ,1

рк ,1+1 рк ,1 , / ч рк+1,1 рк ,1 -+ (ск 1 - ак 1)-:--рк 1ак;

Ат

А х

Ат

- + (ск 1 - ак 1)-

Ах

=(к -1) р, , [д+с, , Ф к, ] - рк а 1 ф к , - рк а. -г Лк-1А Лк' ■ (5)

Лк; А х

Для уравнения (5) будет два аналога, так как коэффициент Н может менять свой знак в зависимости от величины ск1 :

1=1

ск, +1 - ск,

^, > о

при с-

Рк ,1+1 — Рк ,1

+ С,

Рк,1 Рк 1,1

--/

Ат ' А х = (Як 1 + Ск,ф к,) Рк (к-1);

рк ,1+1 рк Ат

+С

рк ,1+1 рк Ах

(6)

10 лет

журналу

47

с,„< 0

при Ск

рк ,1 +1 — рк ,1 Ат

+С

рк +1,1 рк ,1 2

Ах

-а,,

рк ,1+1 Рк,1 . рк +1,1 рк,1

, , ■ + С, , ,

Ат

Ах

= (Як,1 + Ск,1ф к,1) Рк ,1(к -1).

(7)

Разностные уравнения относительно рь+1, скм, рк1+1 записаны в неявной форме. Для удобства дальнейшего расчета решим их относительно искомых величин:

рк ,1 +1 = рк ,1

Ск ,1 + ак ,1

[ркрк-1, + Рк, гак,1 (Ск,1 - Ск-1,1)]--

Аt Ск1 - а

1 к, 1

А х

х

х[к+и - рк,, -Рк,1ак, 1 (ск+1, 1 - Ск;+ (Як,, + Ск, 1фк,1) Рк, 1(к - 1)А -

Ах

2 Ск ,1 ААк ,1 Л, Рк ,ак,,~т—г^ Аt;

" , Ак1 А х

(8)

Ск ,1+1 = Ск ,1

Ск ,1 + ак ,1

х

рк +1,1 рк ,1 Рк, ,ак,1

рк ,1 рк -1,1 V Рк, А ,1

Л

+ Ск ,1 Ск-1,1

Аt + Ск ,1 - ак ,1 х

--1 х

Ах 2

+ Ск +1,1 Ск ,1

Аt л А ^-Ф к ,1 Аt.

Ах

(9)

Математическая модель движения газа в линии высокого давления

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

( рабс рабс К )

I

■ = 1,4 -1 0-

'к, 1922^ ^4 ( е У

-1+-1 дг,

^ е.

где рабс н - абсолютное давление в начале газопровода, МПа; рабс к - абсолютное давление в конце газопровода, МПа; I - длина рассчитываемого участка газопровода, м.

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

к

2

2

расчетной длины газопровода на 5.. .10 %. Расчетная длина наружных надземных и внутренних газопроводов определяется по формуле

1 = 'д + X &

где 1д - действительная длина газопровода, м; ££ - сумма коэффициентов местных сопротивлений участка газопровода длиной /д; 1Э - условная эквивалентная длина (м) прямолинейного участка газопровода, потери давления на котором равны потерям давления в местном сопротивлении со значениями коэффициента £=1.

Эквивалентная длина газопровода в зависимости от режима движения газа в нем определяется по формулам, приведенным в таблице.

Режим движения газа Эквивалентная длина, м

Ламинарный 1э = 5,5 -10-6 Q V

Критический , /1794 • d4v '' = i Q

Турбулентный I = d 3 / ч0.25 11 Г k + 1922vd ^ I d QH J

В работе [3] приведены таблицы с коэффициентами местных сопротивлений при турбулентном движении газа и номограммы получения коэффициентов для различных диаметров и давлений газа в газопроводе. В проектах АГНКС в линии подачи газа используются трубы с внутренним диаметром 40 и 14 мм.

Математическая модель заправки транспортных средств КПГ

1. Уравнение работы сжатия при заправке

Минимальная работа сжатия подаваемого на заправочное устройство природного газа в ходе заполнения баллона транспортного средства (ТС) при наличии блока охлаждения КПГ, обеспечивающего поддержание постоянной температуры (Тг = const), и постоянном геометрическом объеме заполняемого баллона ТС (^бал= const), может быть вычислена по формуле

k-1

k ' ' 4

Lmin =-

к-1

f Pp V к

-1

ВС у

где Ь^ - минимальная возможная работа сжатия, необходимая для заправки транспортного средства; к - показатель адиабаты сжимаемого газа; Я - универсальная газовая постоянная; Твс - температура на входе компрессорной установки; г - коэффициент сжимаемости природного газа при температуре Тг и давлениирр;рвс ирр - давления всасывания компрессора и рабочее газобаллонной установки.

z

При использовании аккумуляторного блока АГНКС работа, затрачиваемая на заполнение аккумуляторной емкости с давлением рак, выражается формулой

—

к RT„r

к -1

—

ВС у

Поскольку, как показано в работе [1], использование общего аккумулятора приводит к непроизводительному расходованию 50 % затрачиваемой энергии, используется двух- или трехступенчатая заправка транспортных средств. При этом давление в аккумуляторных емкостях поддерживается на меньшем уровне в сравнении с рабочим давлением транспортного средства. На практике широкое распространение получила многоступенчатая заправка, предполагающая наличие нескольких аккумуляторных блоков с разным объемом и рабочим давлением. Заправка сначала производится из аккумулятора большей емкости с меньшим давлением, а затем из аккумулятора меньшей емкости и большим давлением.

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

к

к -1

(

-1

ВС у

где Емс - минимальная возможная работа сжатия при многоступенчатой заправке, необходимая для заправки ТС; п - число ступеней заправки; р1 - давление в аккумуляторной емкости 2-й ступени заправки.

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

2. Система уравнений газовой динамики при заправке транспортного средства

Несмотря на внешнюю простоту процесса заправки ТС он характеризуется резко выраженной нестабильностью, значительными перепадами давления и температуры потока газа, высокими скоростями, достигающими скорости звука. Таким образом, процесс заправки относится к области газовой динамики, в том числе газовой динамики высоких скоростей. В процессе заправки транспортного средства можно рассматривать выпуск газа из сосуда большого объема с постоянным давлением р0 через сужающийся трубопровод в баллоны ТС с давлением ра. В данном случае при выполнении условия УШ>>У1С скорость течения газа определяется по формуле Сен-Венана - Венцеля:

V —

2К \к-1

gRTо

—

о V

где К - показатель адиабаты; g - ускорение свободного падения; Т0 - температура газа в сосуде.

Для перепада давлений в 10 МПа при температурах от -50 до 50 °С и режиме выпуска сжатого газа из сосуда через трубопровод в случае адиабатического течения

10 лет

журналу

49

к -1

к

г

к -1

к

г

1—1

к-1

к

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

Из раздела физики по газовой динамике известно, что при критической скорости удельный массовый расход газа через трубопровод максимален и определяется выражением V = pv.

Скорость звука для природного газа характеризуется сжимаемостью и определяется выражением

c = \kgRT -

г - pT

- дг -

др

_2к_ k + 1

gRT-

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

Сущность процесса заправки на начальном этапе соответствует адиабатическому течению с критической скоростью, что описывается уравнениями состояния для критического давления, температуры и плотности:

кр ад

2

к +1

к к-1

т =-

крад к + 1

т

-1 п

к+1

1

к -1

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

С

к

__р

С„

Для природного газа в интервале давлений 0,1...25 МПа и температур от -50 до 60 °С к изменяется в интервале от 1,3 до 3,0.

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

о

4р = ™п{Е<6,; 4СОтс; 4СОзк+зУ}: 1=1

где О - число баллонов, установленных на ТС; 4сб, - диаметр проходного сечения баллонного вентиля; 4СОтс - минимальный диаметр проходного сечения линии высокого давления в газовой части топливной системы газобаллонного ТС; 4Созк+зУ - минимальный диаметр проходного сечения арматуры заправочной колонки и заправочного устройства.

2

г

0

2

Р

Р

0

Правила построения интегральной математической модели объекта производства КПГ

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

• соответствие модели схеме коммутации элементов производства сжатого газа и анализ наличия разветвлений в обвязке газопровода высокого давления;

• соответствие производительности и вместимости моделируемых вариантов основного и вспомогательного оборудования существующим образцам, произво-

димым промышленностью.

Литература

1. Фотин Б.С., Пирумов И.Б. Поршневые компрессоры. - М.: Машиностроение, 1987.

2. Евстифеев А.А., Люгай С.В. Основы логико-вероятностного анализа безопасности транспортных средств на газовом топливе. - М.: ВНИИГАЗ, 2017.

3. Евстифеев А.А., Ермолаев А.Е. Влияние холостых пробегов газовых городских автобусов на показатели производственно-хозяйственной деятельности // Транспорт на альтернативном топливе. - 2016. - № 4 (52). - С. 23-30.

4. Люгай С.В., Балашов М.Л., Евстифеев А.А. Оценка времени ожидания заправки транспортного средства на АГНКС // Транспорт на альтернативном топливе. - 2016. - № 6 (54). - С. 50-54.

5. Евстифеев А.А., Заева М.А., Сергеев М.С. Метод обеспечения работоустой-чивости системы управления питанием газового транспортного средства // Транспорт на альтернативном топливе. - 2016. - № 3 (51). - С. 51-60.

6. Евстифеев А.А., Заева М.А., Хетагуров Я.А. Методы и средства оптимизации размещения объектов производственно-сбытовой инфраструктуры / В книге: Научная сессия НИЯУ МИФИ-2015. Аннотации докладов: в 3-х томах. Отв. ред. О.Н. Голотюк. - 2015. - С. 74.

7. Попов М.А., Егорова А.Н., Евстифеев А.А. Моделирование и оптимизация мест размещения объектов газовой заправки транспорта / В книге: Научная сессия НИЯУ МИФИ-2015. Аннотации докладов: в 3-х томах. Отв. ред. О.Н. Голотюк. -2015. - С. 97a.

8. Евстифеев А.А., Дрыгина Ю.Н., Ермолаев А.Е. Моделирование и оптимизация процесса развития производственно-сбытовой сети газовых заправочных станций // Газовая промышленность. - 2015. - № S3 (728). - С. 30-33.

9. Евстифеев А.А. Математическая модель определения численности и производительности заправочных колонок на АГНКС // Газовая промышленность. -2015. - № 8 (726). - С. 95-97.

10. Евстифеев А.А. Размещение объектов инфраструктуры газовой заправки транспорта // Транспорт на альтернативном топливе. - 2015. - № 6 (48). - С. 26-39.

11. Евстифеев А.А. Анализ эффективности производственного процесса на АГНКС // Транспорт на альтернативном топливе. - 2015. - № 5 (47). - С. 27-33.

12. Евстифеев А.А. Метод формирования адекватной стоимости газового моторного топлива // Транспорт на альтернативном топливе. - 2015. - № 2 (44). -С. 41-46.

13. Евстифеев А.А. Современные подходы к выбору производительности основного технологического оборудования объектов заправки природным газом //

10 лет

журналу

51

Транспорт на альтернативном топливе. - 2015. - № 4 (46). - С. 48-54.

14. Evstifeev A., Zaeva M., Krasnikova S., Shuvalov V. Multi-criteria equipment control in complex engineering systems // Asian Journal of Applied Sciences. - 2015. -Т. 8. - № 1. - P. 86-91.

15. Евстифеев А.А. Структурный синтез и алгоритмы решения для математической модели системы газовой заправки транспорта и газоснабжения автономных потребителей // Вести газовой науки. - 2015. - № 1 (21). - С. 79-85.

16. Евстифеев А.А. Обеспечение муниципальных автобусных парков мегаполиса газовым моторным топливом // Газовая промышленность. - 2014. - № 2 (702).

- С. 86-89.

17. Евстифеев А.А. Математическая модель процесса заправки транспортных средств КПГ на АГНКС // Транспорт на альтернативном топливе. - 2014. - № 1 (37).

- С. 24-31.

18. Евстифеев А.А. Многокритериальное управление технологическим оборудованием сложных технических систем с использованием методов получения нечетких выводов // Транспорт на альтернативном топливе. - 2014. - № 5 (41). -С. 44-48.

19. Евстифеев А.А. Методология рационального построения и непрерывного совершенствования региональной сети АГНКС // Транспорт на альтернативном топливе. - 2014. - № 3 (39). - С. 53-60.

20. Евстифеев А.А. Математическая модель анализа потребности в КПГ и СПГ на вновь газифицируемых территориях // Газовая промышленность. - 2013. -№ 1 (685). - С. 87-88.

21. Люгай С.В., Евстифеев А.А., Тимофеев В.В., Балашов М.Л., Дрыгина Ю.Н. Сравнение экономических показателей при использовании жидкого моторного и газомоторного топлив // Транспорт на альтернативном топливе. - 2013. -№ 5 (35). - С. 14-19.

22. Евстифеев А.А. Модель прогнозирования потребления газового моторного топлива в населенном пункте // Транспорт на альтернативном топливе. - 2013. -№ 3 (33). - С. 43-47.

23. Евстифеев А.А., Балашов М.Л. Методика определения границы экономической эффективности перехода на природный газ в качестве моторного топлива // Транспорт на альтернативном топливе. - 2013. - № 2 (32). - С. 4-5.

24. Евстифеев А.А. Расчет надежности системы поставок газового моторного топлива потребителям // Транспорт на альтернативном топливе. - 2013. - № 4 (34).

- С. 61-65.

25. Евстифеев А.А., Заева М.А., Хетагуров Я.А. Применение математического моделирования при испытаниях и отработке сложных технических систем // Вестник Национального исследовательского ядерного университета МИФИ. -2013. - Т. 2. - № 1. - С. 115.

26. Евстифеев А.А., Заева М.А. Автоматизированная система единого государственного реестра газобаллонного оборудования // В книге: Научная сессия НИЯУ МИФИ-2012. Аннотации докладов: в 3-х томах. - 2012. - С. 285.

27. Дедков В.К., Евстифеев А.А. Методика оценивания надежности технической системы по результатам испытаний // Вопросы теории безопасности и устойчивости систем. - 2010. - № 12. - С. 215-221.

28. Евстифеев A.A., Северцев Н.А. Модели минимизации направленного ущерба транспортной системы при отсутствии информации // Вопросы теории безопасности и устойчивости систем. - 2009. - № 11. - С. 137-145.