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

УДК 519.6

Н. Н. БЕЛЯЕВ, Е. Ю. ГУНЬКО, Т. П. РЕШЕТНЯК (ДИИТ)

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

Розроблено тривимiрну чисельну модель розрахунку процесу поширення токсичного газу на проммай-данчику та затiканнi токсичного газу у примщення при аваршному викидi на пiдприeмствi. Модель базуегь-ся на чисельному iнтегруваннi рiвняння конвективно-дифузiйного переносу домiшки та моделi течи нестис-ло! рiдини. Наводяться результати обчислювального експерименту.

Ключовi слова: аваршний викид токсичного газу на проммайданчику та його затiкання у примщення, рь вняння конвективно-дифузiйного переносу домшки, модель течи нестисло! рвдини, тривимiрна чисельна модель

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

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

The 3D numerical model to simulate the toxic gas dispersion at an industrial site and inflow of toxic gas into the industrial room after an accident ejection was developed. The model is based on the K-gradient transport model and equation of potential flow. The results of numerical experiment are presented.

Keywords: accident ejection of toxic gas at industrial site and its inflow into industrial room, K-gradient transport model, equation of potential flow, 3D numerical model

Введение

В настоящее время большую актуальность приобретает вопрос оценки риска поражения людей на территории химически опасных промышленных объектов в случае аварий, терактов [1, 2, 5, 8]. Особенно важно адекватно оценить масштаб угрозы вблизи места аварии - на промплощадках (локальный прогноз), т.к. промышленный объект - это, как правило, место концентрации людей. Распространение опасных веществ в атмосфере на промплощадках имеет свои особенности - это перенос загрязнителя в условиях застройки. Наличие зданий способно существенно изменить картину переноса. Поэтому при решении задачи прогноза риска поражения людей на промышленном объекте крайне важно учесть влияние зданий на процесс рассеивания опасных веществ в атмосфере. Постановка физических экспериментов для решения таких задач имеет свои особенности:

1. необходимость использования дорогостоящего оборудования (система Laser Doppler, система Ultra Sonic Turbulence Measurement и т.д);

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

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

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

Создание прогнозных моделей на базе уравнений Навье-Стокса и различных моделей турбулентности (k-г модель, LES модель и т.п.) для решения задач о переносе загрязняющих веществ в условиях промплощадки является проблемным в настоящее время. Это связано, в первую очередь, с необходимостью использования мелкой сетки [12, 13], больших ресурсов компьютера и приводит к значительным затратам машинного времени. Например, при практической реализации коммерческого кода

© Беляев Н. Н., Гунько Е. Ю., Решетняк Т. П., 2011

SERENE [13], реализующего 3D RANS модель, совместно с « к-г » моделью, потребовалось 3,5 суток на получение прогнозного результата по расчету гидродинамики обтекания нескольких контейнеров, поставленных в ряд (рис. 1) и задачи теплопереноса. Отметим, что расчет выполнялся во Франции с использованием 9 процессоров, и при этом длина расчетной области была менее 50 м.

Безусловно, что такие большие затраты компьютерного времени на получение прогнозного результата не позволяют еще считать моделирование на базе уравнений Навье-Стокса каждодневным инструментом решения задач рассматриваемого класса - прогноза аварийного загрязнения атмосферы на промпло-щадках при разработке ПЛАСа (план ликвидации аварийной ситуации). Реализация других, современных коммерческих кодов для прогноза уровня загрязнения атмосферы, например таких - WRF, CMAQ, MUSKAT, SILAM, LOTUSEUROS, MACMOD, WRF-MUSKAT [14, 15] -требует значительного количества метеорологической информации, что при рассмотрении задач локального прогноза - порядка нескольких сотен метров получить практически крайне трудно. Реализация этих кодов осуществляется на сетках порядка 10 х 10 км и более. В этой связи, актуальным является создание математических моделей локального прогноза уровня загрязнения атмосферы при аварийных выбросах токсичных веществ на промышленных площадках, позволяющих учесть основные физические факторы: различный тип выброса, наличие зданий и т.п.

->

Щ- т

ais ."■.'J.

. v

Рис. 1. Схема расположений контейнеров в работе [13]

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

Эту задачу условно можно считать сочетанием двух задач: «внешней» задачи - т.е. расчета рассеивания ТВ в атмосфере с учетом застройки и «внутренней» задачи - расчет динамики изменения концентрации загрязнения в помещениях. Решение «внешней » задачи позволяет определить значение концентрации ТВ в любом месте, в том числе там, где располагаются воздухозаборники, обеспечивающие подачу воздуха в помещения. Решение «внутренней» задачи позволяет выявить риск токсичного поражения людей в помещениях и возможность возникновения вторичной аварии. Целью настоящей работы является построение численной модели для решений перечисленных «внешней» и «внутренней» задач. Достоинством предложенной модели является возможность учета основных физических факторов, влияющих на процесс переноса токсичного газа в условиях застройки и внутри помещений, и при этом - небольшие затраты машинного времени при практической реализации модели. Расчет гидродинамики воздушного потока в трехмерной постановке осуществляется в рамках модели невязкой несжимаемой жидкости.

Математическая модель

загрязнения атмосферы

Для моделирования процесса переноса загрязняющего вещества на промплощадке будем использовать трехмерное уравнение миграции примеси [2 - 4]

дС диС ЗуС д(w - ws)С — +-+-+ ^-=

дt

дх

= 1^ х

dz

Í dC ^ » y —

ду

С )+_£

дх ) ду У у ду у

§)+£атг-г), №

где С - концентрация загрязняющего вещества; и, V, w - компоненты вектора скорости воздушной среды; wS - скорость оседания примеси; ц = (ц х, ц у, ц г) - коэффициент турбулентной диффузии; Q - интенсивность выброса токсичного вещества; 5( г - г,) - дельта-функция Дирака; г, =( х^, у{, 21) - координаты источника выброса.

Метеорологические параметры - профиль скорости ветра, коэффициенты атмосферной

диффузии в построенной численной модели рассчитываются по зависимостям [3 ,4]

u = щ

f z Vn

v Z1 у

ц z = 0,11z;

Цу = 0,2 и(х, у, г); ц = Цу ,

где щ - значение скорости ветра на высоте = = 10 м; п = 0,15.

Постановка краевых условий для уравнения (1) рассмотрена в работах [6, 7].

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

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

Для экспресс-прогноза используется нульмерная модель [8, 11] вида

VdC = LCdt - LC

(2)

где V - объем помещения; Спом - концентрация ТВ в выходящем из помещения воздухе; Ь - воздухообмен; С - концентрация ТВ во втекающем воздухе.

Необходимио отметить, что концентрация ТВ во втекающем воздухе (т.е. концентрация ТВ на месте расположения воздухозаборника на здании) постоянно изменяется со временем и определяется из решения «внешней» задачи -т.е. задачи о рассеивании токсичного вещества в атмосфере между зданиями.

Для решения данного уравнения используется метод Эйлера [10].

Относительно модели (2) необходимо отметить следующее:

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

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

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

Для расчета 3Б процесса загрязнения воздушной среды в помещении с учетом его геометрической формы, наличия в нем оборудования, с учетом положения отверстий приточно-вытяжной вентиляции используется вторая модель - уравнение переноса примеси (1).

На границе втекания (приточное отверстие вентиляции) загрязненного воздуха в помещение задается концентрация ТВ. На твердых границах - пол, стены, потолок - ставится граничное условие равенства нулю потока примеси, на границе выхода потока из помещения (выходное отверстие вентиляции) ставится «мягкое» граничное условие [9]. В начальный момент времени полагается, что концентрация ТВ С = 0 внутри помещения.

Модель гидродинамики

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

дР дР дР и =—; V = —; ^ = —, дх ду дг

где Р - потенциал скорости.

Уравнение для определения потенциала скорости имеет вид

д P д P д P п —Т + —Т +—Т = 0.

дх2 ду2 dz2

(3)

Для уравнения (3) ставятся следующие граничные условия:

• на твердых стенках

^ = 0, дп

где п - единичный вектор внешней нормали;

• на входной границе (граница втекания воздушного потока)

- = vn,

дп

где Vn - известное значение скорости;

• на выходной границе

P = P (х = ranst, у) + const (условия Дирихле).

Метод решения

Численное интегрирование уравнения (1) осуществляется с использованием неявной попеременно-треугольной разностной схемы расщепления [6]. Для численного интегрирования уравнения (3) используется метод Либмана [9, 10]. Для формирования вида расчетной области в модели используется метод «porosity technique» [12].

Практическое применение численной модели

Разработанная численная модель и созданный на ее основе код были применены для решения следующей задачи. На промплощадке расположены два здания, со смещением относительно друг друга (рис. 2). Первое здание имеет П-образную форму, и возле него происходит аварийный выброс сероводорода. Второе здание имеет Г-образную форму.

Рис. 2. Схема размещение зданий на промышленной площадке

Параметры задачи: размеры расчетной области: 60 х 60 х 60 м; интенсивность выброса -400 г/с. В расчетах принималось: щ = 3 м/с -значение скорости ветра на высоте = 10 м; п = 0,15, ws = 0. При исследовании ставится задача определить концентрацию токсичного газа в производственном помещении, расположенном в торце второго здания (рис. 2, место поступления загрязненного атмосферного воздуха в помещение - т. е. положение воздухозаборника, условно показано стрелкой); координаты воздухозаборника: х = 42,5 м; у = 17,5 м; г = 7,5 м). Воздухообмен в помещении (объемом 223 м3) составляет 0,37 м3/с.

Результаты решения задачи показаны на рис. 3 - 5, где представлена зона загрязнения атмосферы на промплощадке для момента времени ^ = 150 с после аварии. Из данных рисунков отчетливо видно формирование зон загрязнения в «застойных» областях. Видно, что в различных сечениях на промплощадке форма зоны загрязнения - различна. Так, в сечении у = = 37,5 м зона загрязнения между зданиями значительно больше, чем в сечении у = 27,5 м.

Рис. 3. Зона загрязнения атмосферы, / = 150 с (вид сверху, сечение г = 7,5 м)

Рис. 4. Зона загрязнения атмосферы, / = 150 с (вид сбоку, сечение у = 27,5 м)

Рис. 5. Зона загрязнения атмосферы, / = 150 с (вид сбоку, сечение у = 37,5 м)

Динамика изменения уровня загрязнения воздушной среды в помещении при затекании в него загрязненного атмосферного воздуха после аварии (расчет на базе нуль-мерной модели) выглядит следующим образом:

t = 38 с С = 4,87 мг/м3

t = 48 с С -- = 8,05 мг/м3

t = 56 с С = 10,68 мг/м"

t = 60 с С = 12,01 мг/м

t = 65 с С = 13,33 мг/м

Отметим, что моменту времени ^ = 0 соответствует начало выброса ТВ на промплощад-ке. Принимая во внимание, что ПДК для сероводорода составляет 10 мг/м3, то примерно через 1 мин после аварии концентрация загрязнителя в помещении превысит это пороговое

значение. Таким образом, данный промежуток времени определяет время эвакуации людей из помещения.

Ниже представлены результаты расчета загрязнения воздушной среды в помещении при затекании в него загрязненного атмосферного воздуха и выполненные на базе второй модели - 3Б модели переноса примеси (1). Схема производственного помещения показана на рис. 6. В помещении находится технологическое оборудование в виде параллелепипеда. Поступление воздуха в помещении осуществляется снизу, а выходное отверстие вентиляции находится на середине противоположной стены.

Рис. 6. Схема производственного помещения

Рис. 7. Зона загрязнения в помещении, t = 5 с (вид сбоку, сечение у = 3,63 м)

). 176Е+00 о 164Е+00 г ). 153Е+00 й ). 142Е+00 1 ). 130Е+00 л !.119Е*00 1 ). 108Е+00 \ ¡.963Е-01 е ¡.В50Е-01 ¡.737Е-01 ).623Е-01 у ¡.510Е-01 ¡.397Е-01 ¡.283Е-01 ).170Е-01 ¡.567Е-02

Рис. 9. Зона загрязнения в помещении, t = 38 с (вид сбоку, сечение у = 3,63 м)

При проведении расчетов по распространению загрязненного воздуха в помещении принималось: коэффициент турбулентной диффузии цх = цу = цг = 0,2 м2/с . Скорость втекания

загрязненного воздуха в помещении - 1,5 м/с; площадь входного отверстия системы вентиляции - 0,66 м2.

На последующих рисунках предоставлена динамика загрязнения воздушной среды в помещении для различных моментов времени (еще раз отметим, что моменту времени t = 0 соответствует начало выброса токсичного газа на промплощадке). Из данных рисунков видно, как происходит дифракция фронта загрязненного воздуха над технологическим оборудованием в помещении. Вблизи входного отверстия формируется зона загрязнения с большим градиентом концентрации. Отчетливо видно движение загрязненного воздуха к потолку помещения.

На базе 3Б численной модели имеется возможность прогноза величины концентрации в любой интересующей точке в помещении. Например, концентрация сероводорода над оборудованием (точка А, рис. 7) изменяется с течением времени следующим образом:

t = 28 с С = 0,012 г/м3

t = 34 с t = 45с t = 54 с

С = 0,019 г/м3 С = 0,032 г/м3 С = 0,044 г/м3

Рис. 8. Зона загрязнения в помещении, t = 26 с (вид сбоку, сечение у = 3,63 м)

Обратим внимание на то, что первая модель, нуль-мерная, дает заниженное значение концентрации загрязнителя в помещении, чем 3Б численная модель. Это объясняется тем, что в нуль-мерной модели осуществляется осреднение прогнозного значения концентрации по объему помещения. Совершенно очевидно, что

такая «прогнозная» информация может существенно дезориентировать проектировщика.

Очевидно, что применение второй модели дает существенное преимущество по сравнению с нуль-мерной, т.к. имеется возможность выполнить 3D прогноз развития зоны загрязнения в помещении за 10...20 секунд с учетом его формы, наличия оборудования, места втекания воздуха и т.д.

В заключение отметим, что для решения сопряженной задачи прогноза динамики загрязнения воздушной среды как на промплощадке, так и помещении, и при использовании 3D численной модели потребовалось около 1 мин. времени работы ПК.

Выводы

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

БИБЛИОГРАФИЧЕСКИЙ СПИСОК

1. Аварии и катастрофы. Предупреждение и ликвидация последствий [Текст] : учеб. пособие в 6-ти кн. / под ред. В. А. Котляревского и А. В. Забегаева. - М.: Изд-во АСВ, 2001 - 200 с.

2. Беляев, Н. Н. Расчет распространения загрязняющих веществ в условиях застройки [Текст] / Н. Н. Беляев, Е. Ю. Гунько // Вюник Дшпро-петр. нац. ун-ту залiзн. трансп. iм. акад. В. Ла-заряна. - 2007. - Вип. 18. - Д.: Вид-во ДНУЗТ, 2007. - С. 21-24.

3. Берлянд, М. Е. Прогноз и регулирование загрязнения атмосферы [Текст] / М. Е. Берлянд. - Л.: Гидрометеоиздат, 1985. - 273 с.

4. Бруяцкий, Е. В. Теория атмосферной диффузии радиоактивных выбросов [Текст] / Е. В. Бруяц-

кий. - К.: Ин-т гидромеханики НАН Украины, 2000. - 443 с.

5. Гунько, Е. Ю. Оценка риска токсичного поражения людей при аварийном выбросе химически опасного вещества [Текст] / Е. Ю. Гунько // Вюник Дншропетр. нац. ун-ту залiзн. трансп. iм. акад. В. Лазаряна. - 2008. - Вип. 20. - Д.: Вид-во ДНУЗТ, 2008. - С. 87-90.

6. Численное моделирование распространения загрязнения в окружающей среде [Текст] / М. З. Згуровский [и др.]. - К.: Наук. думка, 1997. - 368 с.

7. Марчук, Г. И. Математическое моделирование в проблеме окружающей среды [Текст] / Г. И. Марчук. - М.: Наука, 1982. - 320 с.

8. Меньшиков, В. В. Опасные химические объекты и техногенный риск [Текст] / В. В. Меньшиков, А. А. Швыряев. - М.: Изд-во МГУ, 2003. -245 с.

9. Роуч, П. Вычислительная гидродинамика [Текст] / П. Роуч. - М.: Мир, 1980. - 616 с.

10. Самарский, А. А. Теория разностных схем [Текст] / А. А. Самарский. - 2-е изд., испр. - М.: Наука, 1983. - 616 с.

11. Эльтерман, В. М. Вентиляция химических производств [Текст] / В. М. Эльтерман. - 3-е изд., перераб. - М.: Химия, 1980. - 288 с.

12. Mahfound, K. Computation of Wind Flow and Air Pollution for Regions Having a Complex Topography [Text] // Kadja Mahfound, Anagnostopoulos [et al.] / Proc. of 3rd European & African Conf. on Wind Engineering (Eindhoven Univ. of Technology, Netherlands, July 2 - 6, 2001). - P. 355-358.

13. Development of building resolving atmospheric CFD code taking into account atmospheric radiation in complex geometric [Text] / Qu Y [et al.] // Conf. Abstracts of 31st NATO / SPS Int'l Technical Meeting on Air Pollution Modelling and it's Application (27 Sept. - 01 Oct. 2010, Torino, Italy). -№ P1.5.

14. Rakitin, A. Modeling the impact of urban emissions in Russia on air quality in northern Europe [Text] / A. Rakitin, M. Makarova, D. Ionov // Conf. Abstracts of 31st NATO / SPS Int'l Technical Meeting on Air Pollution Modeling and it's Application (27 Sept. - 01 Oct. 2010, Torino, Italy). - № 2.17.

15. The use of MM5 - CMAQ air pollution modeling system for real time and forecasting air quality impact of industrial emissions [Text] / R. San Jose [et al.] // Advances in Air Pollution Modeling for Environment Security. - NATO Science Series, Springer. 2004. - Vol. 54. - P. 326-335.

Поступила в редколлегию 04.11.2010.

Принята к печати 17.11.2010.