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

БЕЗОПАСНОСТЬ В ЧРЕЗВЫЧАЙНЫХ СИТУАЦИЯХ

УДК 614

А.В. Рыбаков, К.Ю. Замана

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

В работе описывается модель разлета осколков поврежденного оборудования при взрыве на пожаровзрывоопасном объекте. Модель представляет собой систему дифференциальных уравнений. Приведены расчеты дальности, высоты и времени полета осколков в конкретных сценариях взрыва, изложены результаты исследования на корректность применения разработанной модели. Алгоритм расчета основан на применении метода Рунге - Кутта четвертого порядка.

Ключевые слова: чрезвычайная ситуация техногенного характера, взрыв, поврежденное оборудование, разлет осколков.

A. Rybakov, K. Zamana

ASSESSMENT MODEL ZONE OF POSSIBLE DAMAGE FROM THE SEPARATION OF FRAGMENTS OF EQUIPMENT DURING MANMADE EMERGENCIES ON FIRE EXPLOSION

OBJECTS

The paper describes a model of separation offragments of damaged equipment in an explosion at fire explosion facility. The model is a system of differential equations. The calculations of distance, altitude and time offlight of the fragments in specific scenarios explosion presents the results of studies on the validity of the application of the developed model. The algorithm of calculation is based on the method of Runge-Kutta fourth order.

Keywords: man-made emergencies, explosion, property damage, scattering debris.

Чрезвычайные ситуации (ЧС) техногенного характера на пожаровзрывоопасных объектах (ПВОО) связаны с возникновением пожаров (диффузионное горение) или взрывов (детонация или дефлаграционное горение). К основным поражающим факторам взрыва относится воздушная ударная волна (избыточное давление в ее фронте) и осколочное поле (количество осколков, их кинетическая энергия и радиус разлета), которое в свою очередь образует возможную зону поражения людей. Как показывают произошедшие ЧС [1], поражающие действия от разлета осколков приводят к большому количеству пострадавшего населения. Так, например, в августе 2015 года в Тяньцзине (север Китая) на севере страны прогремели два мощных взрыва на складе химических веществ, в результате которых погибли, по меньшей мере, 50 человек и пострадали свыше 700 человек. Кроме этого, более 3,5 тысяч человек были вынуждены покинуть свои дома и ночевать под открытым небом во временном укрытии, так как в жилых домах, находящихся наиболее близко к взрыву, были выбиты почти все окна, внутри домов мебель разбросана, словно подхваченная вихрем, сорванные с петель двери между комнат лежали под углом [1].

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

Целью статьи является создание модели разлета осколков и обоснование ее применения.

1. Описание модели.

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

Будем считать, что кинетическая энергия, приобретенная оболочкой разрушенной трубы или оборудования, составляет ^ часть от потенциальной энергии сжатого газа и определяется соотношением

тоб^02 ... т

—б— = пМгАг, (!)

где тоб = 2пЯ0роб8Ь - масса оболочки трубы, кг (здесь мы предполагаем, что толщина стенки трубы пренебрежимо мала по сравнению с радиусом трубы);

Мг = р^пЯ^Ь - масса газа, кг;

рг, роб - плотность газа в трубе и плотность материала оболочки трубы соответственно;

I - длина разрушенного участка;

В = - параметр, характеризующий соотношение масс оболочки и сжатого газа;

™об

Аг - удельная потенциальная энергия (работа расширения) газа для адиабатического процесса;

Уо - начальные скорости фрагментов оболочки трубы, приобретаемые в момент разрыва трубы.

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

Р I (Ро\~

лг=ш-т)\1-(т) I' (2)

где Р - давление сжатого газа в трубе; Р0 - атмосферное давление;

к - показатель адиабаты. Для двухатомных газов (кислород, водород, азот) к = 1,4. От общей энергии взрыва на разлет фрагментов трубы расходуется значительная часть энергии взрыва, в расчетах следует принимать значение ц = 0,5 [4]. Таким образом, зная давление газа в трубе и его

плотность, можно найти Аг, а затем, пользуясь формулой (1), вычислить начальную скорость осколков:

Vo = 4mAY =

РовЧк-Щ VP

'г

м

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

^оск

движутся под действием силы тяжести и силы сопротивления воздуха. Сила сопротивления вычисляется по стандартной формуле:

1 2'

4пр = -~SCpCxpBvV, (4)

где 5ср - площадь миделя осколка;

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

Характерные значения площади осколка по трем ортогональным сечениям

I I 51 =-5, Б2 = пИ5, Б3 =-пИ, (5)

поск поск

где И = 2(Я0 + 5) « 2Я0 - диаметр трубы.

С учетом произвольного вращения фрагмента в воздухе в качестве характерной площади миделя фрагмента можно принять следующее значение:

5 + Б2 + Б3 пЮ

5ср =--т-3 « о-. (6)

ср 3 3п

оск

Здесь мы пренебрегли площадями 51 и 52 в силу малости 5.

Таким образом, уравнения движения осколка в проекциях на вертикальное и горизонтальное направления имеют вид:

dy

dv 1

SepCxPb^v2 + w2v - mog;

(7)

mo^= -2SepCxPu^v2 + w2v - mog;

dx ~dt=W;

dw 1 /—=-г

mo^ü= ~^SepCxPB\V + WW,

v

где y, v - вертикальные значения расстояния и скорости; x, w - горизонтальные значения расстояния и скорости; m0 - масса осколка. Начальные условия: Хо=Уо = 0,

v0 = V0 sin a, w0 = V0 cos a,

где a - угол разлета осколка.

^ср

Пусть т0 = роб5ср—; тогда уравнения движения для скорости преобразуются к виду

поб

к2 = [[£п + 2'Уп + 2к1

(ё.у СхрВПоск / ? ' 7 — = ---^у2 + ж2у - д;

& 2Ьроб ^ СхРвпоск -7

— =--—-^у2 +

М 2ьроб

Найти аналитическое решение системы (7 - 8) дифференциальных уравнений затруднительно, поэтому целесообразно использовать численные методы. В данной модели предлагается использовать метод Рунге - Кутта четвертого порядка. Вычисления следующего шага производятся по следующим итерационным формулам:

У(к) = У0

к =/&п'Уп)

н н 2'Уп + 2

Н Н \ (9)

к3 = /[1п + 2'Уп + 2к2) к4 = КЬп + Н'Уп + Нкз)

н

Уп+1 =Уп+-^(к-1 + 2к2 + 2кз + к4) здесь 0 - начальное время, н - шаг итерирования, кроме того

(Л ^ Ч (-С^У2 + ж2У- д\ (уД _ СхРвпоск пт

У = (у]'/а'У) = 1 1'У0 = 1о)С=(10)

\ -С^у2 + ж2ж ) \Ж0/

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

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

Ключевым моментом в данной модели является использование явного метода Рунге - Кутта 4-ого порядка для численного решения системы дифференциальных уравнений движения осколков, аналитическое решение которой не представляется возможным. Как известно, явные численные методы не являются абсолютно устойчивыми [5], а на соответствующие числа Куранта накладываются ограничения (следовательно, накладывается ограничение и на выбор шага итерирования).

2. Исследование метода Рунге - Кутта на устойчивость.

Исследуем на устойчивость метод Рунге-Кутта четвертой степени для решения системы дифференциальных уравнений, полученной в процессе решения задачи о разлете осколков. Исследование будем проводить методом экспоненциальной подгонки [8]. Суть метода состоит в

получении по методу Рунге - Кутта явной итерационной формулы вида уп+1 = д(а)уп для решения модельного дифференциального уравнения у' = Лу и решении неравенства |д(а)| < 1, которое и определяет область устойчивости. Здесь Л - собственное число матрицы Якоби исходного отображения, задающего систему дифференциальных уравнений; а = ЛН- число Куранта.

Вычисления удобнее проводить с помощью матричных обозначений и таблицы Бутчера [5]. Последняя для метода Рунге - Кутта четвертого порядка имеет вид

А В

Б

0 0 0 0 0

1/2 1/2 0 0 0

1/2 0 1/2 0 0

1 0 0 1 0

1/6 1/3 1/3 1/6

Здесь введены обозначения А, В, И для соответствующих столбца, матрицы и строки таблицы Бутчера.

Для модельного уравнения у' = Лу отображение имеет вид /(Ь, у) = Лу, поэтому

к^ = / (п + ЩН уп + к^ЬцЩ ) = Л (уп + Н^Ь^к-) Л = 1 ...4 (11)

Если ввести столбцы к = (к1 к2 к3 к4)т,е = (111 1)т, то 4 равенства (11) можно записать в виде одного матричного:

к = Лупе + аВк, (12)

или

(I - аВ)к = Лупе. (13)

В нашем случае матрица

1 0 0 0\

С-аВ) = 1-Г -1/2 1 0} (14)

0 0 - а 1 имеет обратную, равную

1 0 0 0\

а/2 1 0 0

(1-аВ)-1 = \а2/4 а/2 1 о). (15)

\а3/4 а2/2 а 1) Значит, существует единственный к, определяемый из формулы (13):

к = Луп(1-аВ)-1е. (16)

Подставляя его в формулу для вычисления следующего шага

4

уп+1 = уп + Н^^^к^ = уп + Нйк, (17)

1=1

получаем

уп+1 = уп(1 + аО(1-аВ)-1е). (18)

Откуда

234 а2 а3 а4

ц(а) = 1 + аИ(1 - аВ) е = 1 + а + — + — + —. (19)

2 6 24

Точки области устойчивости являются решениями неравенства

234 а2 а3 а4

1 + а + — + — + — 2 6 24

< 1. (20)

Теперь на время оставим это неравенство и вернемся к исходной задаче. Метод Рунге - Кутта применялся для численного решения следующей системы дифференциальных уравнений:

dy dt

dv --

— = —Cуv2 + w2v — д dt

dx dt

= ж

dw

—— = —Су V2 + { dt у

Вводя столбец у = (у V х ж)т, систему (21) можно записать в виде

V

у = !(у) =

—С^2 + w2v —

д

ж

Матрица Якоби отображения / имеет вид

—С У V2 +

(22)

0 0

](П =

— С

1

2v2 + ж2 V V2 + ж2

— С

vw

— С

vw

ж2

— С

ж2

1

V2 + 2ж2

V V2 + ж2)

(23)

Собственные числа этой матрицы находятся из равенства определителя матрицы ((/) — Л1) нулю. Раскладывая этот определитель по 1-му и 3-му столбцу, приходим к уравнению 2v2 + ж2 vw

Л2

— С

— Л — С

ж2

ж2

— С

vw

V2 + 2ж2 —С —Л

ж2

ж2

= Л2(л2 + 3СЛуй2+^2 + 2С2&2 + ж2)) = 0 (24)

Решая его, получаем собственные числа

Л12 = 0, Л3 = —С^й^Г^2, Л4 = —2С^Ут+^2. (25)

Заметим попутно, что VV2 + ж2 есть абсолютная величина скорости осколка. Обозначим ее V, тогда собственные числа (25) запишутся в виде

Л121 = 0, Л3 = , Л4 = —2^ (26)

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

а + — + — + —. Вторая ее производная

(27)

(28)

а2 1 , 1 ч"(а) = 1 + а + - = -(1 + а)2+->0

положительна, поэтому первая производная

а2 а3

д'(а) = 1 + а + — + —

1

строго возрастает. При этом д (0) = 1 > 0,ц (—2) = — - < 0, поэтому на интервале (—2; 0) первая

производная имеет нуль. Значение ао соответствует точке д'(а) = 0. В силу своего строгого возрастания функция ц'(а^) отрицательна при а < ао и положительна при а > ао. Значит, сама функция д(а) убывает при а < а о и возрастает при а > ао, имея при этом глобальный минимум в точке а о. Замечая, что

а

д(а) = Ч'(а)+—,

(29)

имеем

ао4

Ч(ао) =а24>°,

(30)

0

0

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

неравенстве |ц( а)| < 1 можно снять и рассматривать неравенство ц(а) < 1. Кроме того, поскольку

2

—2 < а0 < 0, то ц(а0) < - < 1, значит уравнение ц(а) = 1 имеет два решения (одно из них получается при пересечении графика функции ц(а) с прямой Ц = 1 слева от точки а0, а второе -справа), а тогда решением неравенства ц(а) < 1 является отрезок с концами в этих решениях. Правый конец этого отрезка, очевидно, есть ноль ( ц(0) = 1), а левый является единственным корнем кубического уравнения

2 3

а а2 а3

1+- + — + — = 0. (31)

2 6 24

Отметим, что уравнение (31), как известно, разрешимо в радикалах, а также допускает решение по стандартным численным методам (метод бисекции, метод касательных и т.д.) Мы не будем здесь останавливаться на решении этого уравнения, а сразу приведем приближенное значение его корня

а = —К^ —2,7853. (32)

Итак, искомым решением неравенства ц(а) < 1 является отрезок [—К;0]. Теперь осталось лишь проверить, что числа Куранта, соответствующие найденным ранее собственным числам, попадают в этот отрезок. Все эти собственные числа неположительны, а наибольшим по модулю является Л4, поэтому достаточно рассмотреть только его. Условие попадания числа Куранта а4 = Л4к в отрезок [—К; 0] приводит к неравенству

2СУк < К, (33)

из которого находим соответствующее ограничение на выбор шага итерирования к:

К

к<--(34)

2СУ

Согласно физическому смыслу задачи, скорость осколка в течение полета не может превзойти его начальной скорости по абсолютному значению в силу наличия силы сопротивления, поэтому неравенство (34) можно усилить, если вместо V подставить начальную скорость Уо- Тогда окончательно условие устойчивости имеет вид

К

h <

2CVn

(35)

3. Пример расчета.

Сценарий С(ГНВ): «Гильотинный разрыв в галерее нагнетателей надземного нагнетательного газопровода ближайшего к главному щиту управления газоперекачивающего агрегата, с образованием воздушной ударной волны (ВУВ) и пожаром».

Пусть радиус трубы, её длина и толщина её оболочки равны R0 = 510 мм, L = 10 м, S = 21,5 мм соответственно; плотность материала трубы роg = 7800 кг/м3. Давление газа внутри трубы примем равным Р = 7,5 МПа, его показатель адиабаты к = 1,4. При нормальных условиях атмосферное давление воздуха Р0 = 101 кПа, его плотность рв = 1,225 кг/м3. Часть переданной кинетической энергии ц = 0,5, коэффициент сопротивления осколка в воздухе сх = 0,2. Шаг по времени выберем h = 0,01.

По формуле (3) находим значение начальной скорости осколков: Vq = 142 м/с. Затем проводим расчеты по методу Рунге - Кутта (9-10). Шаг итерирования принимаем равным 0,01. Результаты расчетов дальности полета, высоты подъема и времени полета осколков в зависимости от их числа и угла разлета приведены в таблице 1.

Таблица 1

Результаты расчетов дальности, высоты и времени полета осколков в метрах по сценарию С(ГНВ)

^■оск 10 ^■оск 100 ^■оск 1000

L, м H, м t, с L, м H, м t, с L, м H, м t, с

а = 30° 1749 255 14,4 1498 234 13,8 691 144 10,7

а = 45° 2007 507 20,4 1658 456 19,3 688 255 14,3

а = 60° 1737 760 24,9 1416 673 23,4 563 358 17

Посмотрим теперь, какой максимальный шаг интегрирования позволяет использовать условие устойчивости (35) в указанном сценарии. Поскольку коэффициент С увеличивается с увеличением числа осколков (а значит, ограничение на шаг становится строже), то будем проводить расчеты для поск = 1000. В результате С = 1,571 • 10-3, У0 = 142, а значит ктах = 6,243. Предлагаемый в решении задачи шаг к = 0,01 с большим запасом удовлетворяет условиям устойчивости для рассматриваемого сценария.

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

Время, с

Рис.1. Пространственное распределение параметров разлета осколков

В результате (см. рис. 1, «ось дальность полета») получаем, что персонал или население, находящиеся в радиусе 2 км, окажутся в зоне ударных воздействий от разлета осколков. Это влечет необходимость разработки инженерно-технических мероприятий, которые позволили бы снизить возможные последствия поражения населения.

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

Рис. 2 (А, Б) 3Б-модель оценки параметров разлета осколков

Рис. 2 (В, Г). 3Б-модель оценки параметров разлета осколков (продолжение)

4. Выводы.

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

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

Реализация разработанной модели в виде программно-вычислительного комплекса позволяет в оперативном режиме выполнять оценочные и прогнозные расчеты по анализу возможных техногенных факторов опасности и их последствий, что, в свою очередь, способствует своевременной разработке актуальных регламентов по оценке безопасности объектов и проведению инженерно-технических мероприятий гражданской обороны [9], разработку декларации промышленной безопасности.

Литература

1. http://www.bbc.com/russian/international/2015/08/150813_china_blasts_death_toll_rises (дата обращения: 21.10.15).

2. Бесчастнов М.В. Промышленные взрывы. Оценка и предупреждение. - М.: Химия, 1991. 432 с.

3. Зиновьев В.Е., Плотникова Е.С., Гончар А.Э. Оценка последствий аварийных взрывов топливно-воздушных смесей в подземных емкостях нефти, нефтепродуктов // Наука и технологии трубопроводного транспорта нефти и нефтепродуктов. - М.: ООО «ТрансПресс», 2014, № 4(16), С. 74 - 83.

4. Методические указания по проведению анализа риска для опасных производственных объектов газотранспортных предприятий ОАО «Газпром». - М.: ОАО «Газпром», 2009.

5. Калиткин Н. Н., Корякин П. В. Численные методы в двух книгах. Книга 2: Методы математической физики.- М.: Издательский центр «Академия», 2013.

6. Седов Л. И. Механика сплошной среды. Том 2. - М.: Наука, 1970.

7. Глаголев К. В., Морозов А. Н. Физическая термодинамика. - М.: МГТУ, 2004.

8. Холодов А.С., Лобанов А.И., Евдокимов А.В. Разностные схемы для решения жестких обыкновенных дифференциальных уравнений в пространстве неопределенных коэффициентов. - М.: МФТИ, 1985.

9. Шеломенцев, В.Н. Правовые аспекты освоения разработки и защиты Арктики // Научный вестник Московского государственного технического университета гражданской авиации. - 2015. - №216. - С. 113-117.

Рецензент: доктор технических наук, профессор Воскобоев В.Ф.