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

УДК 502.55:519.711.3

МОДЕЛИРОВАНИЕ РАСПРОСТРАНЕНИЯ ВЫБРОСОВ ОПАСНЫХ ВЕЩЕСТВ С ОБЛАКАМИ ГОРЯЧЕГО ГАЗА В УСЛОВИЯХ ПРОМЫШЛЕННОЙ И ГОРОДСКОЙ ЗАСТРОЙКИ

А.В. Аксаков

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

Модель подъема облака нагретого газа

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

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

Расчет плотности температуры и скоростей потоков газа проводился численным интегрированием по времени полных уравнений Навье-Стокса [1]

дЛ х '

+ V (руу ) + Ур = р§ + Уа + Р -р/ [п х у ]

+ У(рЕ + р )у = р§у + У(упУ/ + ау ) + Ру

где Е = / + (1 + рр/р)(уу/2 - §г) - удельная полная энергия системы;

J = р![(у - 1)р] = суТ - внутренняя энергия системы;

Т - температура; р - давление;

V - скорость среды; р - плотность среды;

р - плотность примеси;

он - тензор вязких напряжений;

П - вязкость воздуха;

у = ср /су « 1,4 - показатель адиабаты;

« 719,89 - удельная изохорическая теплоемкость воздуха; ср « 1007,85 - удельная изобарическая теплоемкость воздуха; § - ускорение свободного падения;

Р - объемные силы, действующие со стороны частиц примеси; П - нормаль к поверхности земли; / - параметр Кориолиса.

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

Для корректного описания турбулентного рассеяния энергии и импульса в уравнениях газовой динамики коэффициент молекулярной вязкости п заменялся на полуэмпирический коэффициент турбулентной вязкости пэфф , рассчитанный из полуэмпирической К£-модели турбулентности [1]. Сущность К£-модели состоит в добавлении к исходной системе соответствующих уравнений для кинетической энергии турбулентности К и скорости ее диссипации £. Запишем эти соотношения в следующем виде:

д(рК )

Ы д(р8)

Ы

+ У(рКГ ) = V + У(р8Г ) = V

(

п

эфф

V ак

пэфф

Л

VK + В-р8

)

+ -|(С1 В -с2р8) К

Ту = Пэфф

дх, ,

3 (рк + ), В

^дvг ^

vдX )

пэфф = Сцр '

К2

где В - член, характеризующий производство турбулентной энергии;

т.. - тензор напряжений.

При численном решении эмпирические константы К£-модели имели следующие значения [1]: сц_ 0,09 ; с1 _ 1,45; с2 = 1,9 ; аК _ 1; ае _ 1,3 .

На нижней границе расчетной области скорость газа принималась постоянной и равной скорости воздуха на поверхности земли. Также принимались нулевыми потоки тепла, массы, кинетической энергии турбулентности и скорость диссипации турбулентной энергии. На остальной части границы расчетной области выбирались условия протекания с экспериментальным распределением атмосферных параметров по высоте [2]. Распределение кинетической энергии турбулентности и скорости диссипации турбулентной энергии по высоте находилось из условия

баланса между процессами генерации и диссипации турбулентной энергии в рамках К£-модели

турбулентности.

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

В общем случае перенос примеси описывается уравнением турбулентной диффузии [4]:

дС д(у,-С)_^ ^ - дС

дг г_! дх. .._! дх.

(

VК дх, , V 1J

(3)

где V. - компоненты скорости среды;

К.. - коэффициент турбулентной диффузии.

У

Для компоненты вектора смещения отдельной лагранжевой частицы это уравнение переписывается в виде [4-6]

Ах. _ |у (X, г) + £ Кт^} ^+ ^ *"12К1 (х, г)<* , (4)

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

Аг - шаг по времени;

1 - случайная нормально распределенная величина с единичной дисперсией и нулевым

математическим ожиданием.

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

К (^ 0 = ((фф (x, t)/р(x, 0) = ^ (x, О5*, (5)

где Vэфф (х, /) - кинематическая турбулентная вязкость.

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

дvэфф . . 0/2

ах = Vxdtа/+(2vэффёг) $х;

дх дv

ау = Vy а/ ^+(2v эффа^ ) $ у;

(6)

дV э

ск = ( - ^ )+^ ^+(2v эффа/) $ ,

где - скорость гравитационного оседания примеси.

Модель атмосферного переноса примеси

В общем случае перенос примеси в атмосфере описывается уравнением (3). Так же как и в случае полного газодинамического расчета, данное уравнение решалось методом случайных блужданий (4)-(6). Скорости течения газа и коэффициенты турбулентной диффузии будем получать, восстанавливая профили атмосферных параметров по результатам наземных наблюдений [2]. Особенностью данного типа течений является зависимость атмосферных параметров только от вертикальной координаты. Следовательно, мы можем переписать трехмерное стохастическое дифференциальное уравнение Ито (4) в следующем виде [6]:

дш / / ч\ дш / / ч\ дш „ / / ч\ д2ш / / \\ д2ш

аГ ^ ( »&-( %+к" (»ай+К"( %

дК„

^ = ( -^)+—^а/+(к» $г д2

ш(х,у,/о) = 5(х-хо)5(у-Уо), Нш ш(х,у,/) = 0

(х2+у2 )<»

(7)

где ш(х, у, /) - плотность вероятности перехода отдельной лагранжевой частицы из точки

у о, / о) в точку (х, у, /). Уравнение для плотности вероятности перехода допускает полуаналитическое решение:

ах ^) = х0 + | ^ ( ^)) аУ (t) = У0 +1 ^у ( ^))ds

t0 t0

о x (t) = 1 Кх (z (s )) а у ^ ) = } ^ (z (s ))

ш

(х у, t ) =

4л>/аД0а"(7)

exp

(х ^)- ax ^))2 (у (t)- ay ())2

4а x (t)

4а у (t)

(8)

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

Вычисление концентраций и токсонагрузок в приземном слое толщиной 1г :

I

г г í

Б (х, у, t )=1 С (х, у, s )(у, 0 (х ) = |1

К (х у, 0

0 х < 0

1 х > 0

(9)

где тг - масса, сосредоточенная в г-й частице;

С(х, у, t) - концентрация в точке (х, у) на поверхности земли в момент времени t;

Б(х, у, t) - токсонагрузка.

Подобное усложнение численного метода (7)-(9) по сравнению с (6) позволяет существенно сократить общее время расчета.

Примеры расчетов

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

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

облака можно оценить как Б0 =

12 у- 1 Ео

,1/3

, где Е0 - масса взрывчатого вещества в тро-

V П У Ра J

тиловом эквиваленте; ра - атмосферное давление на поверхности земли. Начальное отноше-

1

ние плотности вещества в облаке к плотности окружающего воздуха оценивалось как = 0,126 [8].

Известно, что при подрыве емкостей с жидкостью образуется облако, состоящее из пара и мельчайших капель жидкости. Общую массу выброса при взрыве мощностью в 20 кг тротило-вого эквивалента берем порядка 500 кг. Из них 120 кг - пары люизита и хлористого мышьяка, остальная масса облака состоит из капель диаметром 1,0-0,1 мм, равномерно распределенных по размеру. В начальный момент времени примесь равномерно распределена по облаку. Используя полученное облако в качестве начальных данных для численного решения уравнений (1)-(9), были вычислены концентрации токсичных химических веществ (ТХВ) в атмосфере за время порядка трех часов с момента взрыва.

Для примера на рис. 1 и рис. 2 приведены результаты расчета обтекания облаком взрыва препятствий различной высоты в виде элементов промышленной застройки. В первом случае (см. рис. 1.) представлены расчеты для газообразных продуктов взрыва. Во втором случае (см. рис. 2.) - для аэрозоля.

На рис. 3 показана эволюция облака взрыва при обтекании препятствия.

0 100000

0 0 00010 -"-"-"-I-I-I-1-

0 2000 4000 6000 ВООО 10000 12000 14Л00 16000

Х|т1

Рис. 1 - Зависимость трехчасового прогноза токсической нагрузки (Д кг-с/м3) для газообразных продуктов взрыва от высоты препятствия (к, м).

X - расстояние до эпицентра взрыва. Скорость ветра на высоте десяти метров 4 м/с. Четвертый класс устойчивости стратификации атмосферы (Т-ИЭМ) [7]

и =э И =7

, \ чХ. \

0 2000 4000 6000 0000 10000 12000 14000 16000

Х[п.1

Рис. 2 - Зависимость трехчасового прогноза токсической нагрузки (Д кг-с/м3) для аэрозоля от высоты препятствия (к, м). X - расстояние до эпицентра взрыва. Скорость ветра на высоте десяти метров 4 м/с. Четвертый класс устойчивости стратификации атмосферы (Т-ИЭМ) [7]

Рис. 3 - Обтекание неподвижной стенки (й=7м) облаком взрыва. Скорость ветра на высоте десяти метров 4 м/с.

Четвертый класс устойчивости стратификации атмосферы (Т-ИЭМ) [7]

В целом была проведена серия экспериментов по моделированию переноса примеси с облаком взрыва. Общей для этих опытов является тенденция к увеличению приземной концентрации ТХВ с ростом скорости ветра при неустойчивой стратификации атмосферы.

Работа выполнена при финансовой поддержке гранта РФФИ № 01-01-96444, гранта научно-технической программы «Университеты России» № УР.03.01.015 и гранта МНТЦ № 2065.

ЛИТЕРАТУРА

1. Белоцерковский О.М. Численное моделирование в механике сплошных сред. - М.: Наука, 1984.

2. Мурин А.В. Математическое моделирование на параллельных системах последствий химических аварий. - Дис... канд. физ.-мат. наук. - Ижевск, 2002.

3. W.-Mei Jiang, H. Liu, and H.-Nian Liu.The Numerical Simulation on Atmospheric Transport and Dispersion of the Spray Atomized from Flood Discharging by Hydropower Station over Complex Terrain. - Meteorol. Atmos. Phys. 70, 215-226 (1999).

4. D. Brickman, P.C. Smith. Lagrangian Stochastic Modelling in Coastal Oceanography. - Journal of Atmospheric and Oceanic Technology. - Vol. 19. - № 1, 2002. - Pp. 83-99.

5. V.P. Reshetin, T.S. Zenich, R.V. Arutyunyan, V.V. Belikov, V.P. Kiselev, V.N. Semenov,

D.N. Tokarchuk. Electric Power Plant on the enviroment and human healf - Seventh International Conference On Harmonisation within Atmospheric Dispersion Modelling for Regulatory Purposes Organised by JRC-EI May 28-31. - Belgirate, Italy, 2001. - Pp. 163-167.

6. Кузнецов Д.Ф. Численное моделирование стохастических дифференциальных уравнений и стохастических интегралов. - СПб: Наука, 1999.

7. Атмосфера. Справочник. - Л.:Гидрометеоиздат, 1991.

8. Гончаров Е.А. Модель, описывающая динамику подъема облака неядерного взрыва /

E.А. Гончаров, В.Н. Пискунов, А.И. Харченко, Ф.Дж. Мартин, Х.У. Черч // Вопросы атомной науки и техники. Сер. Теоретическая и прикладная физика. - 1995. - Вып. 3/1. - С. 59-68.