CC BY
29
УДК 669:658.011.56
Ю.В. ШАРИКОВ, д-р техн. наук, профессор, yvshar@mail. ru И.Н. БЕЛОГЛАЗОВ, д-р техн. наук, профессор
Yu.V. SHARIKOV, Dr. in eng. sc., professor, yvshar@mail.ru I.N. BELOGLAZOV Dr. in eng. sc., professor
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ АВАРИЙНЫХ СИТУАЦИЙ, ВОЗНИКАЮЩИХ ПРИ ХРАНЕНИИ И ОБРАБОТКЕ ПОТЕНЦИАЛЬНО ОПАСНЫХ ХИМИЧЕСКИХ ВЕЩЕСТВ
Рассмотрены возможные причины разгерметизации и разрушения хранилищ и аппаратов для обработки и хранения потенциально опасных материалов. Отмечено, что тепловой взрыв является одной из основных причин разрушения хранилищ и реакторов для обработки таких продуктов. Приведены основные уравнения для описания таких процессов для различной формы аппаратов и различных условий теплообмена на их внешней поверхности; предложена и опробована численная процедура для решения уравнений на основе программного комплекса ^егтЕх.
Изучены также процессы движения и эволюции облаков опасных веществ, образовавшихся при внезапной разгерметизации хранилищ и реакторов. Приведены основные уравнения для математического описания рассмотренных процессов. Предложена и опробована процедура численного моделирования процесса рассеяния веществ из облака на основе программного пакета ReactOp.
Ключевые слова: математическое моделирование, тепловой взрыв, аварийные ситуации, эволюция облаков опасных веществ.
MATHEMATICAL SIMULATION
OF EMERGENCY SITUATIONS ARISING AT STORAGE AND TREATMENT OF POTENTIAL HAZARD CHEMICAL SUBSTANCES
Possible reasons arising destroying of storage and apparatuses for potential hazard materials have been analyzed. Thermal explosion has been indicated as one of most impotent reasons of destroying chemical reactors and storages. Basic equations for describing of this phenomena have been evaluated for different type geometry and conditions of heat exchange on outside surface. Numerical procedures have been proposed and tested for solving and analyzing of evaluated equations with using special program package ThermEx.
Processes of evaluation arising clouds of hazard substances have been considered. Basic mathematical equations have been proposed for analyze of considered processes. The procedure for solving of considered problems has been proposed and tested on base of special program package ReactOp.
Key words: mathematical simulation, thermal explosion, emergency situation, evolution of hazard substances clouds.
Как показал анализ главных опасностей, основными причинами возникновения и развития аварийных ситуаций [3] является нарушение герметичности хранилищ потенциально
опасных химических веществ и аппаратуры для их получения и переработки. Рассмотрим основные, наиболее распространенные причины разгерметизации:
1. Нарушение герметичности вследствие разрушения материала аппарата или хранилища под действием либо коррозии со стороны химического продукта, либо с наружной стороны, под действием атмосферных осадков с примесью агрессивных химических продуктов.
2. Разрыв прокладок, сальниковых уплотнений и других устройств герметизации аппаратуры.
3. Разрушение аппаратуры или трубопроводов вследствие непредусмотренных механических воздействий.
4. Нарушение герметичности аппаратуры вследствие неожиданного роста температуры и давления, вызванного внезапным отклонением хода процесса от регламентных норм.
Такую ситуацию в зарубежной литературе называют «run-away». Явления, возникающие при несанкционированном и неконтролируемом отклонении хода процесса от регламента, чаще всего приводят к так называемому «тепловому взрыву». Тепловой взрыв явился первопричиной очень многих крупных аварий в промышленности.
Математическое моделирование процессов возникновения и развития теплового взрыва. Тепловой взрыв является одной из основных причин возникновения крупных выбросов пожаровзрывоопасных и токсических веществ, приводящих к техногенным катастрофам. Например, причиной известной аварии в Бхопале (Индия), приведшей к выбросу больших количеств изоционата и отравлению огромного числа жителей, был тепловой взрыв в реакторе синтеза полиуретанов. В 20-е гг. ХХ в. в Германии произошел взрыв почти 200 т аммиачной селитры, сопровождавшийся большими разрушениями и человеческими жертвами. Причиной его также был тепловой взрыв. Тепловой взрыв может возникнуть при хранении больших масс химических продуктов, склонных к разложению с выделением тепла. Если реакция разложения является экзотермической, то выделяющееся при реакции тепло может частично отводиться в окружающую среду, а частично расходоваться на повышение температуры разлагающегося продукта. Соотношение этих долей зависит от теплоемкости и теплопро-
водности разлагающегося продукта, и от скорости теплоотвода в окружающую среду
[4].
В каждой точке хранилища или аппарата синтеза скорости тепловыделения и теплоотвода в единице объема соответственно
^gen " . dq^
—— = Ь Wi AHi
dt i=i
dt
= dT_
AF spec j '
dn
где wi и AHi - скорость и тепловой эффект i-й реакции соответственно; - удельная
поверхность единицы объема; X - коэффициент теплопроводности реагирующего
объема;
. dT
dn
градиент температур по нор-
мали к рассматриваемой поверхности.
В качестве хранилищ и емкостей для хранения и транспортировки химических продуктов используют сосуды различной геометрической формы, в последнее время все чаще сферической. Сосуды такой формы наиболее экономичны, так как сфера имеет минимальную наружную поверхность при равном объеме по сравнению с сосудами другой формы, и на изготовление емкостей тратится минимальное количество металла. Однако сосуды цилиндрической формы остаются наиболее распространенными для транспортных емкостей и хранилищ.
В связи с этим рассмотрим вывод уравнений теплового баланса для сосудов сферической и цилиндрической форм при протекании в них процессов терморазложения химических продуктов.
Составив уравнение материального и теплового баланса для сосуда сферической формы, получим дифференциальные уравнения, описывающие процесс тепловыделения при разложении химических продуктов, сопровождающийся тепловыделением:
дТ
X (d2T 2 dT Л
+—г 1 + 2
AH,
дг c pi dr r dr I j=ic p
vr v 7 V
W j
(1)
где AHj и Wj - тепловой эффект и скорость ]-й реакции в единице объема соответственно; су - теплоемкость единицы объема реакционного пространства; р - плотность реакционного пространства.
Это дифференциальное уравнение в частных производных 2-го порядка. Для него нужно задать начальные и граничные условия. Начальные условия задаются в виде заданного профиля поля температур в начальный момент времени
Т (0, г) = Т (г).
Граничные условия задаются следующим образом. Первое граничное условие задается из условий изотропности свойств и симметричного распределения поля температур относительно центра сферы:
dT (т, r) dr
= 0.
(2)
Второе граничное условие задается либо в виде заданной температуры на внешней границе сферы (граничное условие 1 -го рода)
T (т, R) = T
given '
(3)
либо в виде заданного теплового потока на границе (граничное условие 2-го рода)
dT (т, r) dr
qsurf ;
(4)
либо в виде равенства тепловых потоков к поверхности сферы изнутри и с поверхности сферы в окружающую среду (это граничные условия 3-го рода):
dT (т, r)
dr
= -a[T(0,R) - Tsur],
(5)
где qsuгf - заданный тепловой поток на границе с окружающей средой; Т8иг - заданная температура окружающей среды.
Выведя аналогичным образом уравнения материального баланса, получим уравнение для изменения концентрации реагента внутри сферы:
^ = ц
дт '
Гд2 ^
2 д cA
дг r dr
-Z w
j=i
', j ■
где Di - коэффициент диффузии /-го компонента; а - концентрация /-го компонента; ^ у - скорость изменения концентрации /-го компонента в у-й реакции.
Поскольку стенки емкости непроницаемы для потока вещества, граничные условия для этого уравнения задаются в виде
dc¡ (х,0) dc¡ (х, Я)
dr
dr
= 0.
а начальные условия в виде заданного поля концентраций по радиусу сферы
Ci(r,т) = c. . (т).
1 х 7 ' г given v 7
Тогда математическая модель для описания возникновения и развития теплового взрыва в хранилищах сферической формы примет вид
( д2 T 2 dT Л
dT
дт дт
= a
кд r2
+ — r дr
, AHjWi
+ Z——
j=1 Cv P
= D
fd2c1 кд r2
+
2 де^ r дr
+ Z
w
j=1
ci (0, r) = c(r); T(0, r) = T (r); д (т,0) д (т, R)
dr
дr
= 0,
где а - коэффициент температуропроводности.
Для температуры граничные условия задаются в соответствии с уравнениями (2)-(5) в зависимости от условий теплообмена на внешней поверхности хранилища или аппарата.
На основании уравнения теплового и материального баланса для сосуда цилиндрической формы, получим уравнения математической модели процесса в цилиндрических координатах в виде
дт (д2 T 1 дт Л
дr
= a
\д r2
+ — r dr
+Z
AH i
nv p
dT
"дТ ar
(d2 T 1 dT Л d2 Tm AH
2 +-Kdr r dr
+ax^T Z
Л2 Z - -Wj,
dx j=1 nvp
dni t=D
(дП+1 dc
dr2 r dr
\
+ Dx
d 2n
AH,
dx
2 Z^— Wi, j .
1 =1 nv P
Для случая реакции в плоской пластине получим уравнение
dT — = a—2 дт dx2
д 2T AH
- + ■
-w„
cv P
Анализ уравнений для сосудов трех геометрических форм позволяет, введя параметр геометрической формы Г, получить одно обобщенное уравнение:
dT
dz
= a
fd2 T
dr
+ ТдТ r dr
Л
AH nv p
-Wr
* = D dz
( d2 c Г dc
\
dr
r dr
± Wr
J
где Г = 0 для плоской пластины; Г = 1 для цилиндра и Г = 2 для сферы.
Для сосуда цилиндрической формы распределение тепла вдоль оси цилиндра можно учесть эквивалентным дополнительным плоским потоком вдоль оси, которое нужно прибавить к двум другим слагаемым, описывающим изменение по радиусу. В итоге для этого случая получим следующую систему уравнений модели:
dT С — = a dz
v
d2 т i dT
dr2 r dr
Л
+ ax
d2T
AH
H--Wr
dx nV p
dÄ = Dr dz r
fd2n 1 dn dr2 r dr
Л
+ Dx
J
'äX2
Wr
Для этой системы записываются граничные условия 1-го или 2-го рода в зависимости
от условий теплообмена на боковой поверхности и стенках цилиндра. Такими уравнениями можно описывать процессы в горизонтальных емкостях и в железнодорожных цистернах. Однако эти уравнения применимы к емкостям, заполненным полностью. Если заполнение не полное, нарушается осевая симметрия в распространении тепла и массы и в уравнения необходимо добавить производные по угловой координате, учитывающей несимметричность в заполнении цистерны или хранилища, и граничные условия по угловой координате в центре и на горизонтальной поверхности уровня в цистерне. Вышеприведенные уравнения модели пригодны для описания процессов в твердых материалах или в высоковязких жидкостях, когда можно пренебречь переносом тепла и массы за счет движения жидкостей. Для решения приведенных уравнений и анализа условий возникновения и развития теплового взрыва, может быть использован программный комплекс ТЬегшБх [1], который позволяет создать модель процесса для выбранной геометрии хранилища, ввести физико-химические и тепло-физические свойства содержимого хранилища, вид кинетического уравнения и кинетические параметры для реакции разложения.
2.0016
Высота, м
Радиус, м
О 90
Л
Радиус, м
Высота, м
Рис. 1. Эволюция температурного поля при исследовании экзотермической реакции в цилиндрическом хранилище при равномерных (а) и неравномерных (б) граничных условиях
2
С использованием численных методов решения дифференциальных уравнений в частных производных можно провести многократное моделирование процесса и определить условия переработки, хранения и транспортировки продуктов, обеспечивающие безопасные условия.
Из рис.1 видно, что при граничных условиях, соответствующих неравномерным условиям теплообмена, температурное поле значительно деформируется и максимально горячая точка сдвигается в пространстве.
Эволюция облака при его движении. После наступления теплового взрыва происходит выброс содержимого аппарата или хранилища в атмосферу, образуется облако паров взрывоопасных или токсичных компонентов, которое движется в атмосфере, рассеиваясь в окружающей среде. При этом важно знать, за сколько времени и на каком расстоянии от источника выброса произойдет рассеяние и снижение концентрации выброшенных веществ в облаке до безопасного уровня. Этот уровень определяется предельно допустимыми концентрациями для паров токсических веществ либо величинами, соответствующими нижнему пределу взрываемости или воспламенения для пожаровзрывоопасных веществ.
Для описания процесса эволюции облака используются уравнения конвективной диффузии с подвижными границами и разрабатываются методы численного моделирования. Для реализации таких методов исследования создан целый ряд программных комплексов, основанных на различных моделях эволюции облаков. Согласно этим моделям процесс рассеяния паровых облаков зависит, в основном, от плавучести данного газа и от состояния атмосферных условий.
Плавучесть обычно выражается через соотношение плотностей газа и воздуха следующим образом:
Ф = (Ра1Г "Р^)/ Р
Бая '
Атмосферные условия характеризуют вертикальным градиентом температур - изменением температуры с высотой. В качестве стандартной величины принято, что при
увеличении высоты на 1 м температура уменьшается на 0,01 °С. Таким образом эта величина, называемая сухоадиабатическим градиентом Г = "0,01 °С, является мерой устойчивости состояния атмосферы. Чем больше вертикальный градиент температур, тем более неустойчиво состояние атмосферы, тем больше внутренняя циркуляция воздуха и интенсивнее рассеяние паровых облаков.
Рассеяние облаков паров легких газов. Для рассеяния облаков легких газов с так называемой «нулевой плавучестью», т.е. не отличающихся по плотности от воздуха, был разработан метод расчета, основанный на представлении о том, что в горизонтальном и вертикальном направлениях распределение концентраций соответствует распределению Гаусса
I (с) =
1
1л/2п
где I(с) - функция распределения концентрации газа в облаке; я - дисперсия концентрации газа относительно средней концентрации.
При этом горизонтальная составляющая этого распределения расположена по направлению ветра. Эволюция вертикальной и горизонтальной составляющей этого распределения определяется погодными условиями и состоянием атмосферы и может быть выражена через изменение параметров распределения концентрации вещества в движущемся облаке, которое зависит от погодных условий и состояния атмосферы:
я = F (ЛТШг, ушШ), (6)
где ЛТа1Г - градиент температуры по высоте в данной точке; - скорость ветра в
данном месте.
Вид функциональной зависимости определяется в конкретных экспериментах по рассеянию паровых облаков. Такой подход был впервые применен к изучению рассеяния облаков с так называемой нулевой плавучестью, которые образуются, например, при выбросе газообразных радиоактивных отходов.
Рассеяние облаков паров тяжелых газов. Облака паров тяжелых газов возни-
с /2я
кают обычно при разливе криогенных жидкостей или сжиженных газов. В процессе рассеяния таких паровых облаков можно выделить две основные стадии [3]:
1) облако резко опускается под действием силы тяжести. Продолжительность этой стадии около 1 мин, а погодные условия несущественно влияют на продолжительность этой стадии;
2) опустившееся облако рассеивается при смешении его с воздухом. Продолжительность этой стадии может изменяться от десятков до сотен минут и ее длительность зависит от погодных условий.
Количественную оценку скорости движения границ облака можно выполнить с помощью следующих формул:
uv = c^gh (pi- Р2)/ Pi
где иу - скорость горизонтального перемещения вертикальной границы, м/с; g - ускорение свободного падения, м/с2; р1 -плотность тяжелой фазы, кг/м3; р2 - плотность легкой фазы, кг/м3; ^ - высота парового облака, м; с - безразмерная константа, из теоретических предпосылок константа с может быть принята близкой к единице (с = 1).
Позднее была предложена формула, учитывающая объем перемешиваемого облака:
u
Xr) = c 1 rV&[(Po -Ра)! Po]Vo/ % ,
где r - расстояние от точки возникновения парового облака, м; V0 - объем парового облака тяжелого газа в момент его образования, м3; pa - плотность воздуха, кг/м3; p0 - плотность газа в начальный
момент времени, кг/м3.
Вертикальное перемешивание в газовом облаке незначительно в тех случаях, когда uv >> 2u/r, где ufr - скорость турбулентного перемешивания. Таким образом, переход от гравитационного опускания к турбулентному перемешиванию появляется только тогда, когда uv « 2u/r и p (r) = pa,
т.е. плотность облака достигает плотности
156_
ISSN 0135-3500. Записки Горного института. Т.187
воздуха. В этих условиях гравитационное опускание облака замедляется и начинается его турбулентное перемешивание с воздухом.
При проведении крупномасштабных экспериментов по изучению движения и рассеяния облака тяжелого газа и проверки результатов расчета по указанным формулам было обнаружено, что через 80 с после начала выброса скорость и^ = 0,25 м/с.
Полученные результаты позволили сделать следующие выводы:
• модели, основанные на распределении Гаусса, не применимы для описания поведения тяжелого газа без соответствующей поправки на растекание плотной фазы;
• сильный ветер и нагрев солнцем приводит к ослаблению эффекта растекания;
• находиться с наветренной стороны небезопасно, так как растекание может происходить против ветра за счет турбулентного перемешивания (в ночное время растекающееся вещество может охватить достаточно большую площадь).
Из многочисленных работ, посвященных разработке методов расчета процессов образования и рассеяния паровых облаков, можно сделать вывод, что модели, применяемые для расчета поведения облаков с нулевой плавучестью, нельзя использовать для расчета процессов рассеяния тяжелых газов.
Модели, применяемые наиболее успешно для описания процессов рассеяния облаков, можно разделить на два класса [3]: модели К-типа, или модели турбулентной вязкости, и модели слоя или ящика.
Наиболее строгими являются модели турбулентной вязкости (К-типа). В этих моделях используются дифференциальные уравнения, основанные на сохранении массы, импульса и энергии. Эти уравнения интегрируются численно при соответствующих граничных и начальных условиях.
Выбросы при наличии стационарных источников характерны для работы предприятий при номинальных значениях их параметров в режиме нормальной эксплуатации.
В этом случае для описания процессов распространения и рассеяния выбросов используются математические модели, осно-
ванные на уравнениях конвективной диффузии с источниками и с подвижными границами. Положение источника либо строго фиксировано, если на заводе предусмотрен сбор отсосов от местных источников выбросов и сбор их для выброса через общую трубу рассеяния, либо этот источник распределен по территории завода.
Модели К-типа, или модели турбулентной вязкости. Уравнения математической модели для процессов движения и рассеяния облаков вредных и опасных выбросов в общем случае имеют вид
р^ + Р
Gt
dwx , Gwx , Gwx wx— + Wy— + wz — GX Gy GZ
\
-Pgx
Gx
^Gwx . „
P~w+P
Wx
Gwy
Gx
■ + Wy
= Pgy
Gwy
' Gy Gp
Gy '
■ + Wz
Gwy ~dz~
+P
(
Gci _ „ X~/2
— - Di V Ci -
Gt
Gwz , 5wz , dwz
Wx— + Wy— + Wz —
GX Gy Gz
f
\
_Pgz
dci dci GCi
Wx— + Wy— + Wz —
Gx Gy Gz
(7)
dp H '
Л
f
GWx . Gw
- + -
y +,Gwz
v Gx Gy Gz у
+ Ti .
(8)
Система уравнений (7) описывает изменение скорости движения облака, а уравнение (8) описывает изменение концентрации компонентов с, в движущемся облаке. Здесь Wy и - компоненты вектора скорости вдоль координатных осей; р - плотность смеси; [, - коэффициент диффузии компонентов смеси; у2 - оператор Лапласа (сумма вторых производных концентраций компонентов по координатным осям); у, -объемный источник выброса компонентов из места аварии.
На подвижных границах облака происходит рассеяние компонентов облака в со-
ответствии с законами массопередачи следующим образом:
- Dt.d IT1 _ Px (Ci - Cai) ; - Dt.d ^ _ Py (Ci - Cai); Gx Gx
Ос'
- [ы-Г- = Р г (С, - Са,)-ОХ
где - коэффициент турбулентной
диффузии на границе облака с внутренней стороны; РХ, Ру е Рг - коэффициенты
массопередачи на наружной стороне облака по направлению соответствующих осей координат.
Интенсивность рассеяния облака определяется наличием и направлением ветра, а также атмосферными условиями, которые характеризуются так называемым сухо-адиабатическим градиентом температуры, а также скоростью и направлением ветра.
Модель ящика. Модель слоя или ящика учитывает массообмен на границах ящика и эквивалентна модели реактора идеального перемешивания (перемешивание внутри облака считается очень быстрым и не рассматривается):
dC _ р Fsp (Ci - Cai); Fsp _
F
int
V cloud
P_ ; 8bil _ f (Vi.d > Vwind),
8,
(9)
b.l
где - внешняя поверхность облака (общий интерфейс); Кс1оиа - общий объем облака; Р - общий коэффициент массопереда-чи из объема облака в окружающую среду;
- коэффициент турбулентной диффузии; 5б./ -толщина диффузионного пограничного слоя; уы - коэффициент турбулентной вязкости; Уиы - скорость ветра.
Интенсивность рассеяния облака определяется величиной коэффициента массопе-редачи на границе облака и зависит от степени турбулентности и скорости движения облака, а также от разности плотностей воздуха и газа.
Модель турбулентной диффузии более точно описывает эволюцию облака с учетом изменения распределения концентраций веществ внутри облака и массопередачи с границ облака в окружающую среду. Однако решение уравнений этой модели сопряжено с большими математическими трудностями и, хотя в настоящее время имеются пакеты специальных программ для решения
таких задач, их применение требует использования мощных ЭВМ и связано с большими затратами машинного времени. Такие задачи также можно решать с помощью вышеупомянутого пакета ^егтех [1], если ввести в программу вместо коэффициента теплопроводности коэффициент эффективного перемешивания, зависящий от скорости конвективного переноса внутри облака.
а
0 100 200 300 400 500 600 700 800 900 1000
Время, мин
б
Время, мин
Рис.2. Результаты моделирования в программном комплексе ReactOp процесса рассеяния облака при использовании модели идеального перемешивания внутри облака в программном комплексе ReactOp: а - изменение во времени средней концентрации вещества в облаке; б - изменение концентрации вещества
в окружающей среде с учетом разбавления
Модель слоя или ящика по существу эквивалентна модели реактора идеального перемешивания с граничными условиями, отражающими интенсивность массообмена на границах облака. Для расчета процесса рассеяния по этой модели необходимо решение систем обыкновенных дифференциальных уравнений, что значительно упрощает задачу. При этом получаемая информация позволяет оценить эволюцию изменения средней концентрации в облаке и определить условия достижения ее безопасного уровня. Задавшись дисперсией распределения 5 можно, в соответствии с формулой (6), рассчитать и распределение концентраций внутри облака при нормальном законе распределения. Для решения уравнений модели ящика (9) можно использовать программный комплекс ЯеасЮр [1], позволяющий легко создать модель процесса рассеяния, используя нестандартную модель реактора идеального перемешивания в соответствии с системой (9) с граничными условиями, отражающими массобмен на внешней поверхности облака (рис.2). Необходимые для решения модели (9) коэффициенты были рассчитаны по известным литературным данным [2].
Выводы
1. Рассмотрены закономерности возникновения и развития теплового взрыва при хранении и обработке веществ, склонных к тепловому взрыву применительно к аппаратам различной геометрической формы.
2. Предложена процедура численного анализа процессов на основе предложенных моделей с использованием программного комплекса ^егтЕх. Приведены результаты численного моделирования.
3. Проведен анализ эволюции движения и рассеяния облаков опасных веществ, образовавшихся при потере герметичности хранилищ и аппаратов.
4. Показана возможность анализа процессов рассеяния облаков опасных продуктов с использованием модели идеального перемешивания и специализированного программного комплекса ЯеасЮр.
ЛИТЕРАТУРА
1. Бенин А.И. Автоматизированная система исследований теплового взрыва химико-технологических процессов / А.И.Бенин, А.И.Коссой, Ю.В.Шариков // Журнал ВХО им.Д.И.Менделеева. 1990. № 4.
2. Левич В.Г. Физико-химическая гидродинамика. М., 1959.
3. Маршалл В. Основные опасности химических производств. М., 1989.
4. Франк-Каменецкий Д.А. Диффузия и теплопередача в химической кинетике. М., 1987.
REFERENCES
1. Benin V.I., Kossoy A.A., Sharikov Y.V. Auttomated system of investigation of thermal explosion of chemical and technological processes // Journal VChO im. D.I.Mendeleeva. 1990. № 4.
2. Levitch V.G. Physic-chemical Hydrodynamics. Moscow, 1959.
3. Marshall V. Basic hazards of chemical productions. Moscow, 1989.
4. Frank-Kamenetskiy D.A. Diffusion and heat transfer in chemical Kinetics. Moscow, 1987.
