CC BY
45
Т. Р. Тюменев, С. И. Поникаров, В. С. Гасилов
ПРОГНОЗИРОВАНИЕ РАСПРОСТРАНЕНИЯ ОБЛАКОВ ЛЕГКИХ И НЕЙТРАЛЬНЫХ ОХВ В УСЛОВИЯХ УСТОЙЧИВОЙ АТМОСФЕРЫ ПРИ ПОМОЩИ ЧИСЛЕННОГО МОДЕЛИРОВАНИЯ
Исследуется процессы распространения облаков легких и нейтральных газов при помощи численного моделирования турбулентного движения идеального несжимаемого газа.
Прогнозирование масштабов заражения при авариях на предприятиях с оборотом опасных химических веществ (ОХВ) на сегодняшний день, как и прежде, немаловажная задача. Неспособность существующих методик в отдельных случаях качественно описать процесс распространения облака газа с одной стороны и наличие достаточно мощных вычислительных ресурсов (компьютеров и программ) с другой, приводят нас к выводу об актуальности численного решения задачи.
Определения. Легкий газ - газ, отношение плотности которого к плотности воздуха не превышает 0,8.
Нейтральный газ - газ, отношение плотности которого к плотности воздуха лежит в пределах 0,8.. .1,2.
Описание расчетной модели. Выбрана стандартная к—є модель - математическая модель турбулентного несжимаемого движения газа в поле ветра, состоящая из уравнения неразрывности, осредненных по Рейнольдсу уравнений движения, уравнения энергии для температуры, двух уравнений переноса турбулентных характеристик к, є и соотношения, связывающего турбулентную вязкость с турбулентными характеристиками.
Данная модель имеет вид [1]:
V V = 0; (1)
йУ
" У(йт-УУт) = 0;
йї
р- — = р- д-УР + м■V2V + д
йї
дх,
Мг-
ди, ди,
—L +—1-дх, дх,
Р' й = йР + ІЯ +Л'У Т + -УТ) + Ф + р-± 1-Г-йт-УУт ]
т=1
йк д 1 М + Мт дк + Мт ■ диі ди, ~] +
йї дх, |_Р Рг, дх, 2-Р 1дх; дхі )
йє д 1 ц + мт дє С2 + — МТ є д V ди, '1 + —
йї дх, Р РГє дх, ] 2 ■ к ■ р кдх1 дхі )
Сз-є2
к
(2)
(3)
(4)
(5)
(6)
2
Мт =
С-p- к
= -М - 2. 1 = CP-Mr .
-----------. AT =-----------
S
PrT
Ф =
2
du dui
—L + —L
ydxj
dx,
(7)
где p - плотность;
P - давление;
T - температура;
V- скорость с компонентами U\, и2, U3;
Dm - коэффициент бинарной диффузии в газовой фазе;
Ym - массовая концентрация m-ой компоненты;
V - оператор Гамильтона;
D( ') = д^"-) + V ■ V(...) - субстанциональная производная от скалярной функции;
g - ускорение свободного падения; t - время; h - энтальпия;
X - коэффициент теплопроводности;
Хт - коэффициент турбулентной теплопроводности;
Л - коэффициент динамической вязкости;
Лт - коэффициент турбулентной вязкости;
dQ/dt - заданная объемная мощность источников тепла;
Ср - теплоемкость при постоянном давлении;
Cv - теплоемкость при постоянном объеме;
Xi, x2, x3 - координаты радиус-вектора точки;
Sj - символ Кронекера; к - кинетическая энергия турбулентности;
е - скорость диссипации кинетической энергии турбулентности;
Cj, Сл - заданные константы;
РГт - турбулентное число Прандтля;
РГк и Рге - турбулентные числа Прандтля для к и е соответственно.
Для учета влияния сил плавучести на турбулентность в к- и е-уравнения включаются члены, описывающие силы плавучести [2].
Данная система уравнений решается численно с применением пакета FLUENT. В этом пакете для получения дискретных аналогов решаемых уравнений используется метод контрольного объема. Основная идея метода контрольного объема легко понятна и поддается прямой физической интерпретации. Расчетную область разбивают на некоторое число непересекающихся контрольных объемов таким образом, что каждая узловая точка содержится в одном контрольном объеме. Дифференциальное уравнение интегрируют по каждому контрольному объему. Для вычисления интегралов используют кусочные профили, которые описывают изменение переменной (компонента скорости, давление, температура, турбулентные характеристики) между узловыми точками. В результате находят дискретный аналог дифференциального уравнения, в который входят значения переменной в нескольких узловых точках.
Полученный подобным образом дискретный аналог выражает закон сохранения переменной для конечного контрольного объема точно так же, как дифференциальное уравнение выражает закон сохранения для бесконечно малого контрольного объема. Одним из важных свойств метода контрольного объема является то, что в нем заложено точное инте-
гральное сохранение таких величин, как масса, количество движения и энергия на любой группе контрольных объемов и, следовательно, на всей расчетной области. Это свойство проявляется при любом числе узловых точек, а не только в предельном случае очень большого их числа. Таким образом, даже решение на грубой сетке удовлетворяет точным интегральным балансам.
Проиллюстрируем вышесказанное на примере решения уравнения переноса для переменной ф:
рУУф = йу V, (8)
где V - вектор скорости; йф - коэффициент диффузии для величины ф; р - плотность.
Воспользовавшись формулой Остроградского-Гаусса, проинтегрируем уравнение (1) по контрольному объему:
$ p(pVdS = $ йрУрСЗ , (9)
где З - площадь поверхности, ограничивающей контрольный объем. Полученное уравнение отражает выполнение закона сохранения применительно к контрольному объему (ячейке) вычислительной области.
С учетом дискретизации можем записать:
N
2ЧрА = Ё й (ур)А, (10)
/
где N - число граней, замыкающих ячейку; З/ - площадь грани; (V ф)п - градиент представленный в виде конечно-разностного отношения.
Для нахождения конвективных слагаемых использовалась схема против потока первого порядка точности, суть которой в том, что значение ф на грани контрольного объема равно значению в соседней узловой точке с подветренной стороны грани.
Диффузионные члены аппроксимировались центральными разностями второго порядка точности.
В данной работе для решения уравнений сохранения используется численный алгоритм, который называется подходом с сегрегированием: одно уравнение решается отдельно от других, причем по очереди для каждой переменной.
В данном подходе решения каждое дискретизированное уравнение линеаризуется неявно относительно переменной этого уравнения, то есть для данной переменной ее значение в каждой ячейке вычислено, используя соотношение, которое включает и известные, и неизвестные величины соседних ячеек. Таким образом, каждая неизвестная величина появляется в более чем одном уравнении в системе и эти уравнения должны быть решены одновременно, чтобы получить значения этих величин.
Линейная форма дискретизированного уравнения (3) имеет вид:
аР = Ё 3пьРпЬ+Ь , (11)
пЬ
где ф и фпь - значения переменных в центре рассматриваемой и соседних с ней ячейках соответственно; Эр и апь - линеаризующие коэффициенты для ф и фпь соответственно.
Подобные уравнения могут быть записаны для каждой ячейки в сетке. В результате получаем систему линейных алгебраических уравнений со сложной матрицей коэффициентов, причем для каждой ячейки - одно уравнение. Полученная система решается с помощью алгебраического многосеточного метода в соединении с итерационным методом Гаусса-Зейделя.
Для определения давления из уравнения неразрывности применяется алгоритм SIMPLE (Semi-Implicit Method for Pressure-Linked Equations). Для решения уравнения неразрывности следует интерполировать значения скорости от центра ячейки к грани. Линейная интерполяция в данном случае приводит к нефизическому результату; вместо этого используется взвешенное усреднение, основанное на использовании коэффициентов уравнения импульса. Отсюда, поток через грань равен:
J, = J, + d, p - pc,), (12)
где Рсо и рс1 - значения давления в соседних с гранью ячейках; Jf - значение потока, определяемого скоростью;
. р-S2
df =-=-L, (13)
aP
где ap - усредненный коэффициент.
При решении уравнения сохранения импульса с предполагаемым полем давления поток через грань контрольного объема равен:
j; = j; + d, (p; - p;,). (14)
Это значение потока не удовлетворяет уравнению неразрывности, поэтому необходимо ввести такую поправку Jf, чтобы скорректированное значение потока
Jf = J; + J'f (15)
удовлетворяло уравнению непрерывности.
Алгоритм SIMPLE постулирует, что
J',= df (P'cО - P'c ) . (16)
где p' - поправка давления.
Подставляя уравнение коррекции потоков в дискретизированное уравнение неразрывности, получаем уравнение для нахождения поправки давления p':
ap = Y> а„Ль + b, (17)
nb v 7
N
где b = ^ J; .
f
В результате решения определим давление в ячейке и исправленное значение потока через грань:
p = p" + ap - p', (18)
j, = Jf + df (pcо - pc,), (19)
где aP - параметр релаксации.
Заменяя предыдущие промежуточные значения давления и компонент скорости необходимо повторить приведенную процедуру, пока решение не сойдется.
Для дискретизации нестационарных членов применялась полностью неявная схема, которая, безусловно, устойчива относительно шага времени. Принципом полностью неявной схемы является то, что в пределах всего шага по времени переменная принимается равной новому значению.
В данной работе рассматривается поток с высоким числом Рейнольдса, поэтому вязкий подслой не моделируется. Вместо этого используется набор полуэмпирических формул и функций, которые связывают соответствующие величины на поверхности и в
смежных с ней ячейках. К этим функциям относится “закон стенки” и формулы для турбулентных характеристик вблизи поверхности.
При построении сетки необходимо, чтобы первый от стенки расчетный узел располагался в логарифмической области течения. Граничные условия для всех описываемых уравнениями переноса неизвестных определяются в этом узле на основе функций стенки.
“Закон стенки” для поля средней скорости имеет вид:
где х - постоянная Кармана (х = 0,41);
Е - эмпирическая константа (Е = 9,81);
ир - средняя скорость потока в точке Р;
кр - турбулентная кинетическая энергия в точке Р;
Ур - расстояние от точки Р до поверхности;
тт - касательное напряжение на стенке;
Л - динамическая вязкость жидкости.
Генерация турбулентной кинетической энергии Ок и её скорость диссипации в смежных с поверхностью ячейках вычислены на основе локальной гипотезы равновесия. Согласно этому предположению Ок и е принимаются равными по всему контрольному
Для моделирования шероховатости применяется предположение, что распределение средней скорости (изображенной в полулогарифмическом масштабе), имеет тот же самый наклон (1/к), но различную точку пересечения, то есть, включается в «закон стенки» дополнительная константа. Таким образом, модифицированный «закон стенки» будет выглядеть следующим образом:
Функция шероховатости для полностью шероховатого режима имеет следующий
вид:
(20)
объему.
ди
(21)
(22)
(23)
лв = 1Ц1+Ок, ■ к;),
К; Р- К • и * к С й
где К3 =------------; Кэ - высота шероховатости; Скэ - константа, величина которой за-
М
висит от вида шероховатости (ОКэ = 0,5).
Принимались следующие граничные условия. На твердой поверхности (стенке) накладывались условия прилипания и нулевого градиента концентрации. Температура стенки задавалась на 3 К ниже температуры потока воздуха - моделирование инверсии. Шероховатость моделировалась заданием значений высоты (0,01 м) и константы шероховатости. На входе области задавалась скорость ветрового потока (1 м/с) и его направление, масштаб и интенсивность турбулентности (0,2 м и 0,1 %). Последние значения уточнялись в процессе вычисления. На верхнюю границу накладывалось условие симметрии (градиент всех переменных через границу равен нулю).
Использовалась структурированная неравномерная сетка. Для сходимости численного решения и получения адекватных результатов, но в то же время разумного использования ресурсов при построении сетки производилось сгущение узлов вблизи твердой поверхности (стенки), источника и в зонах излома рельефа, т.е. в местах с потенциальным значительным изменением концентрации. По мере удаления от этих зон, в местах с незначительными концентрациями и градиентом концентраций ячейки сетки увеличиваются в размерах и могут иметь искаженную (отличную от прямоугольника) форму, т.к. это мало влияет на количественные и еще в меньшей степени на качественные результаты.
Начальный шаг по времени Д£ выбирался из расчета, чтобы он не превышал ¡/V, где I -наименьший линейный размер ячейки в направлении ветра; V - скорость ветра. Если сходимость в процессе расчета достигалась за меньшее, чем 10 число итераций, шаг по времени увеличивался.
Результаты расчетов
Расхождение между распределениями объемных концентраций полученных экспериментальным путем [3] и рассчитанных с использованием аналогичной модели при изо-термии составляет не более 20 % [4]. Данный факт подтверждает, что численный метод достаточно точно отображает реальность.
В работе исследовалось влияние угла наклона подстилающей поверхности и атмосферной устойчивости на процесс распространения облаков легких и нейтральных газов. А именно метана, аммиака, этилена и этана. Рассматривались углы 0, 15 и 30°; движение облака «в горку» и «с горки»; состояние атмосферы - инверсия. Время наблюдения - 1800 с (30 мин). Результаты показали интересные качественные закономерности.
Легкие газы. Аммиак и метан стелились по земле. При движении облака аммиака в горизонте и в горку 15 и 30° оно стелилось по земле. При движении облака аммиака с горки 15 и 30° оно оторвалось от склона и продолжало распространяться параллельно горизонту!
Распространение облака метана проходило не так устойчиво, как облако аммиака. На начальном этапе при движении в горизонте оно стелилось по земле, но вскоре оторвалось от нее и стало двигаться вверх под углом к горизонту ~ 40°. При движении в горку, при угле 15° оно оторвалось земли, не дойдя до склона, но двигалось в целом параллельно ему. При угле 30° отрыва от земли не произошло. При движении с горки об-
лако метана отрывалось от земли в месте перелома рельефа и двигалось вверх под углом к горизонту 40-45°.
Нейтральные газы. Облака этана (рис. 1) при движении в горку 15 и 30° и этилена при 15° (в меньшей степени) движутся против ветра в горизонтальной плоскости, и лишь на достаточном удалении от источника начинают двигаться по ветру в сторону склона. Это явление возникает оттого, что в сложившихся условиях (удар ветрового потока о склон, инверсия, свойства газа) возникает обратный (рециркуляционный) поток воздуха в «тонком» слое над поверхностью земли. Поэтому это явление проявляется у более тяжелых газов - этана и этилена. Этим же явлением объясним отрыв от земли легкого метана при движении в горку 15°.
При распространении облака этилена с горки 15° оно стекало вниз по склону, верхняя граница облака оставалась горизонтальной. При распространении этилена с горки 30° облако разделилось - основная масса продолжила движение параллельно горизонту, часть стекала по склону.
Вывод
Ввиду того, что существующие методики не могут дать полученные нами результаты, численное моделирование представляется нам перспективным инструментом прогнозирования масштабов поражения при авариях.
Рис. 1 - Распространение этана «в горку» при инверсии и угле подстилающей поверхности 15 (а) и 30° (б)
Литература
1. Селезнев В.Е., Клишин Г.С., Алешин В.В. Математический анализ газовой опасности при выбросах природного газа // Инженерная экология. 2000. № 5. С. 29-36.
2. FLUENT 6.0 User’s Guide. - Fluent Inc., 2001.
3. Colenbrander G.V., Puttock J.S. Atmospheric dispersion of heavy gases and small particles // IUTAM symposium Netherlands. 1983. Р. 277-295.
4. Галеев А.Д., Гасилов В.С., Поникаров С.И., Чепегин И.В. Численное моделирование процесса рассеивания «тяжелых» примесей в условиях сложного рельефа местности // Безопасность жизнедеятельности. 2004. № 6.
5. ОНД-86. Методика расчета концентраций в атмосферном воздухе вредных веществ, содержащихся в выбросах предприятий.
6. Методика оценки последствий химических аварий (методика ТОКСИ-2). М.: НТЦ «Промышленная безопасность», 1999.
© Т. Р. Тюменев - асп. каф. машин и аппаратов химических производств КГТУ; С. И. Поникаров -д-р техн. наук, проф., зав. каф. машин и аппаратов химических производств КГТУ; В. С. Гасилов -кан. тех. наук, доц. той же кафедры.
