УДК 519.6
Н. Н. БЕЛЯЕВ, Л. Ф. ДОЛИНА, Е. Ю. ГУНЬКО (ДИИТ)
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ПРОЦЕССА ЗАТЕКАНИЯ ТОКСИЧНОГО ГАЗА В ПОМЕЩЕНИЕ ЧЕРЕЗ СИСТЕМУ ВЕНТИЛЯЦИИ ПРИ ЧРЕЗВЫЧАЙНОЙ СИТУАЦИИ НА ТЕРРИТОРИИ ХИМИЧЕСКИ ОПАСНОГО ОБЪЕКТА
Розроблено тривимiрну чисельну модель розрахунку процесу поширення токсичного газу на промпло-щадщ та запканш токсичного газу у примщення при аваршному викидi на пiдприeмствi. Модель базусться на чисельному iнтегруваннi рiвняння конвективно-дифузшного переносу домiшки та моделi течи нестисло! рвдини. Наводяться результати обчислювального експерименту.
Ключовi слова: тривимiрну чисельна модель, токсичний газ, аварiйний викид, нестисла рiдина
Разработана трехмерная численная модель расчета процесса распространения токсичного газа на промп-лощадке и затекания токсичного газа в помещение при аварийном выбросе на предприятии. Модель основывается на численном интегрировании уравнений конвективно-диффузионного переноса примеси и невязкой несжимаемой жидкости. Приводятся результаты вычислительного эксперимента.
Ключевые слова: трехмерная численная модель, токсичный газ, аварийный выброс, несжимаемая жидкость
The 3D-numerical model to simulate the toxic gas dispersion on industrial sites and inflow of toxic gas into the industrial room after accident ejections 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: 3D-numerical model, toxic gas, accident ejections, potential flow
Введение
При возможных авариях на территории химически опасных промышленных объектов крайне важно адекватно оценить масштаб угрозы с целью прогноза риска поражения людей на объекте и разработки эффективных мер защиты. В рамках этой комплексной задачи следует выделить задачу о затекании загрязненного атмосферного воздуха (ЗАВ) через систему вентиляции в помещения зданий, которые расположены на промышленной площадке и последующее загрязнение воздушной среды в этих помещениях. Такой прогноз особенно важен при решении задач промышленной безопасности, т.к. внутри помещений может находиться значительное количество людей и надо оценить уровень угрозы их поражения на рабочих местах. То есть, задачу о прогнозе риска поражения людей на промышленном объекте условно можно разбить на две задачи - поражение людей на открытом пространстве и поражение людей внутри помещений. Так, в принципе, и делается в нормативной методике [11], но риск поражения - т. е. величина людских потерь, в данной методике просто «назначается» без опоры на какие - либо расчеты. С точки зрения промышленной безопасности, здесь крайне важно, при разработке ПЛАСА (план ликвидации аварийной ситуации), рассчитать и научно обосновать величину возможных потерь. Эти потери в различных местах предприятия, в
различных помещениях - то же будут различны. Поэтому необходимо адекватно определить динамику загрязнения воздушной среды, как на открытом пространстве, так и в помещениях, при возможном затекании в них токсичных веществ и, на основе этой информации, оценить риск поражения людей, находящихся в помещениях, возможность их безопасной эвакуации.
Традиционно для решения данной задачи используется нуль-мерная модель прогноза уровня загрязнения воздушной среды в помещениях [8,11] вида
VdCn0м = ЬСЛ - ¿С^ (1),
где V - объем помещения;
Спом - концентрация токсичного вещества (ТВ) в выходящем из помещения воздухе;
¿ - воздухообмен;
С - концентрация ТВ во втекающем воздухе.
Необходимо отметить, что концентрация ТВ во втекающем воздухе (т. е. концентрация ТВ на месте расположения воздухозаборника на здании) постоянно изменяется со временем и определяется из решения «внешней» задачи -т.е. задачи о рассеивании токсичного вещества в атмосфере между зданиями. Таким образом, чтобы решать рассматриваемую задачу - ЗАВ, предварительно, необходимо рассчитать процесс рассеивания ТВ на промышленном объекте. С этой точки зрения решение задачи о затекании загрязненного атмосферного воздуха че-
© Беляев Н. Н., Долина Л. Ф., Гунько Е. Ю., 2011
рез систему вентиляции в помещения становится далеко не тривиальной.
Модель (1) широко использовалась в расчетах для прогноза загрязнения воздушной среды в помещениях химических производств при авариях, в метрополитенах и т.д. [8, 11]. Но следует отметить следующие ее недостатки:
- модель определяет концентрацию токсичного вещества не в помещении, а в удаляемом воздухе. Это не дает возможности оценить реальную степень загрязнения воздушной среды как в самом помещении, так в различных местах внутри него, например, там, где есть риск появления вторичной аварии;
- модель не учитывает влияние положения приточных и вытяжных отверстий вентиляции на процесс переноса токсичного вещества внутри помещения (т. е. не учитывает аэродинамику воздушных потоков в помещении);
- модель не учитывает влияние различных препятствий (технологического оборудования внутри помещения, стеллажей и т. д.) на процесс переноса токсичного вещества и формирование застойных, плохо проветриваемых зон, где возможно накопление загрязнителя.
Но самым главным недостатком модели (1) является то, что «прогнозное» значение концентрации токсичного вещества осредняется по объему помещения, а значит исследователь получает заниженное значение концентрации и тем самым имеет недостоверную информацию о масштабе угрозы - что крайне опасно.
Для ликвидации указанных недостатков следует использовать пространственные (3Б) прогнозные модели. В этой связи важным вопросом является построение эффективных (не только с точки зрения затрат времени на расчет, но и простоты реализации) СББ 3Б моделей для расчета процесса загрязнения воздушной среды в помещении при затекании в них токсичного газа. Для практики крайне важно, чтобы такие модели считали быстро (время расчета задачи составляло порядка нескольких минут) и при этом позволяли учесть аэродинамику воздушных потоков в помещении, положение, форму технологического оборудования в помещении.
Аналогично, крайне важно для практики иметь эффективные СББ 3Б модели для прогноза уровня загрязнения атмосферы в условиях промышленной застройки при аварийных выбросах токсичных веществ. Характерной особенностью задач данного класса является перенос примеси при наличии зданий (препятствий). По этой причине применение аналити-
ческих моделей прогноза для задач данного класса - исключено. Использование на практике прогнозных моделей на базе уравнений На-вье-Стокса и различных моделей турбулентности (k — г модель, LES моделирование и т.п.) для решения задач о переносе загрязняющих веществ в условиях промышленной площадки или внутри помещений является проблемным в настоящее время. Это связано с необходимостью использования для данных моделей крайне мелкой сетки [1, 9, 12, 13, 15 - 18] , необходимостью решения не только уравнений гидродинамики и переноса примеси, но и дополнительных уравнений, применяемых для описания турбулентности. Все это предъявляет высокие требования к ресурсам компьютера. Реализации численных моделей, основанных на уравнениях Навье-Стокса приводит к значительным затратам машинного времени, что является неприемлемым при решении задач, связанных с аварийным загрязнением воздушной среды.
Применение метода физического моделирования для решения рассматриваемого класса задач имеет свои особенности:
1) необходимость использования дорогостоящего оборудования (система Laser Doppler, система Ultra Sonic Turbulence Measurement и т. д.);
2) значительное время на постановку эксперимента, его проведение, обработку данных;
3) проведение физического эксперимента по моделированию процесса переноса токсичных веществ (ТВ) в реальных условиях застройки , при реальных масштабах выброса - опасно, и, по сути, практически неосуществимо;
4) при проведении физического эксперимента в аэродинамических трубах возникает проблема неэквивалентности экспериментальных и реальных условий (например, не соблюдение критерия подобия по числу Рейнольдса).
Применение современных коммерческих кодов для прогноза уровня загрязнения атмосферы, например, таких как WRF, CMAQ, MUSKAT, SILAM, LOTUS-EUROS, MACMOD, WRF-MUSKAT [13 - 15] требует значительного количества метеорологической информации, что при рассмотрении задач локального прогноза - порядка нескольких сотен метров получить невозможно. Реализация этих кодов осуществляется на сетках порядка 10^10 км и более. В этой связи, важным вопросом является создание математических моделей локального прогноза уровня загрязнения атмосферы при аварийных выбросах токсичных веществ на промышленных площадках, позволяющих
учесть основные физические факторы: различный тип выброса, наличие зданий и обеспечить эффективное решение задачи о затекании токсичного газа в помещения.
Целью настоящей работы является построение численной модели и кода для решения задачи прогноза загрязнения воздушной среды в помещениях при затекании в них токсичного газа. Такую задачу можно назвать сочетанием двух задач - «внешней» (перенос примеси в условиях застройки) и «внутренней» (перенос примеси в помещении). Достоинством предложенной модели является возможность учета основных физических факторов, влияющих на процесс переноса токсичного газа в условиях застройки и внутри помещений, и при этом -небольшие затраты машинного времени при практической реализации модели. Расчет гидродинамики воздушного потока в трехмерной постановке, как для «внешней», так и для «внутренней» задач осуществляется в рамках модели невязкой несжимаемой жидкости [3, 4, 9].
Математическая модель загрязнения атмосферы
Для моделирования процесса переноса загрязняющего вещества на промышленной площадке будем использовать трехмерное уравнение миграции примеси [2-4]
дС диС СуС + ОС. _
дt дх ду дг
дх ду
-АГ Сс +_
дх I х дх "
дС И у^
ду У ду
-Си (и.- §)+! я,(! нг -
(2)
где С - концентрация ТВ; и, V, ^ - компоненты вектора скорости воздушной среды; wS -скорость оседания примеси; и = (их, Иу, и.) -коэффициент турбулентной диффузии; Q -интенсивность выброса ТВ; 5(т - г ) - дельта-функция Дирака; г — (х1, yi, ) - координаты
источника выброса.
В построенной численной модели метеорологические параметры - профиль скорости ветра на входе в расчетную область, коэффициенты атмосферной диффузии рассчитываются по зависимостям [3, 4]:
и — и
(I ^
У Ъ )
И. — 0,11- ; Цу = 0,2 и (х, у, г); ц = Цу ,
где щ - значение скорости ветра на высоте 2Х = 10 м; п = 0,15.
Постановка краевых условий для уравнения (1) рассмотрена в работах [6, 7].
Математическая модель загрязнения воздушной среды в помещении
Загрязненный атмосферный воздух попадает в помещение через систему вентиляции.
Для расчета трехмерного процесса загрязнения воздушной среды в помещении при втекании в него ТВ, с учетом геометрической формы помещения, наличия в нем технологического оборудования, с учетом положения отверстий приточно-вытяжной вентиляции используется уравнение переноса примеси (2).
На границе втекания (приточное отверстие вентиляции) загрязненного воздуха в помещение задается концентрация ТВ на месте расположения воздухозаборника вентиляции, которая зависит от времени и определяется в процессе решения «внешней» задачи - т.е. задачи об определении нестационарного поля концентрации ТВ в атмосфере при рассеивания его на промплощадке. На твердых границах - пол, стены, потолок, оборудование и т.п. - ставится граничное условие равенства нулю потока примеси. На границе выхода потока из помещения (выходное отверстие вентиляции) ставится «мягкое» граничное условие [7]. В начальный момент времени полагается, что концентрация ТВ С = 0 внутри помещения.
Модель гидродинамики
Расчет нестационарного процесса загрязнения воздушной среды, как на промышленной площадке, так и внутри помещений при затекании в них ТВ основывается на предварительном решении гидродинамической задачи - определении поля скорости воздушного потока при наличии препятствий (здания, оборудование в помещениях и т.д.). Для расчета поля скорости воздушного потока при обтекании зданий на промышленной площадке или при движении потока внутри помещения делается допущение, что движение воздушной среды -потенциальное [3, 9], тогда компоненты скорости воздушной среды определяются соотношениями:
дР дР дР и ——; V — —; ^ — —, дх ду дг
где Р - потенциал скорости.
Уравнение для определения потенциала скорости имеет вид
д2 Р д2Р д2Р
- + -
- + -
dz2
= 0.
(3)
дх2 ду2
Для уравнения (3) ставятся следующие граничные условия:
• на твердых стенках
дР_ дп
= 0,
где n - единичный вектор внешней нормали;
• на входной границе (граница втекания воздушного потока - входное отверстие вентиляции или ветровой поток)
- = V,
dn п
где Vn - известное значение скорости;
• на выходной границе
P = P (x = ranst, y) + const (условия Дирихле).
Метод решения
Численное интегрирование уравнения (2) осуществляется с использованием неявной попеременно-треугольной разностной схемы расщепления [6]. Для численного интегрирования уравнения (3) используется метод А. А. Самарского [10]. Для формирования вида расчетной области в модели (как при решении внешней, так и внутренней задач) используется метод «porosity technique» [13]. Применение данного метода дает возможность быстро формировать вид расчетной области и «выделять» в численной модели особенности задачи - положение оборудования в помещениях, форму зданий, положение отверстий вентиляции и т.д.
Рис. 1. Схема положения места утечки токсичного
вещества на промышленной площадке и места расположения воздухозаборника на торце здания
Практическое применение численной модели
На основе построенной численной модели создан код CMAP (Computer Modelling of Air Pollution). Данный код был применен для решения следующей задачи. На промышленной площадке расположено здание (рис. 1). Возле него происходит аварийный выброс сероводорода.
Параметры задачи: размеры расчетной области; 60^60x60 м; интенсивность выброса -400 г/с. В расчетах принималось: и1 = 5 м/с -значение скорости ветра на высоте Z1 = 10 м; n = 0,15, ws = 0. При исследовании, ставится задача определить концентрацию токсичного газа в производственном помещении. Воздух в это помещение поступает из воздухозаборника, расположенного в торце здания (рис. 1, место поступления загрязненного атмосферного воздуха в помещение - т.е. положение воздухозаборника условно показано стрелкой). Воздухообмен в помещении (объемом 223 м3) составляет 1 м3/с.
Результаты решения «внешней» задачи показаны на рис. 2 - 4, где представлена зона загрязнения атмосферы на промышленной площадке для различных моментов времени после аварии. Из данных рисунков отчетливо видно формирование обширной зоны загрязнения, охватывающей здание со всех сторон, т.е., в зону поражения выходят все выходы из здания.
Рис. 2. Зона загрязнения атмосферы, t = 5 с (сечение y =27,5 м)
Рис. 3. Зона загрязнения для момента времени / = 8 с (сечение у = 27,5 м)
Рис. 4. Зона загрязнения атмосферы (вид сверху) для момента времени t = 45 с (сечение г = 7,5 м)
Ниже представлены результаты расчета загрязнения воздушной среды в помещении при затекании в него загрязненного атмосферного воздуха и выполненные на базе 3Б-модели переноса примеси (2). Схема производственного помещения показана на рис. 5. В помещении находится технологическое оборудование, которое в численной модели представлено в виде двух параллелепипедов. Поступление воздуха в помещении осуществляется через отверстие в стене, а выходное отверстие вентиляции находится на противоположной стене.
/1 /
1 ¡У / / / /
А > /
х
Рис. 5. Схема производственного помещения, куда втекает загрязненный воздух через воздухозаборник системы вентиляции
На последующих рис. 6 - 9 представлена динамика загрязнения воздушной среды в помещении для различных моментов времени (еще раз отметим, что моменту времени t = 0 соответствует начало выброса токсичного газа на промышленной площадке). Из данных рисунков видно, как происходит дифракция фронта загрязненного воздуха над технологическим оборудованием в помещении. Вблизи входного отверстия формируется зона загрязнения с большим градиентом концентрации.
При проведении расчетов по распространению загрязненного воздуха в помещении принималось: коэффициент турбулентной диффу-
Рис. 6. Зона загрязнения в помещении для момента времени t = 10 с (сечение у =3,6 м)
Рис. 7. Зона загрязнения в помещении для момента времени t = 16 с (сечение у=3,6 м)
1.192Е+01 □
1.100Е»01 г
1.167Е+01 а
1.1Б5Е+01 1
1.143Е+01 л
1.130Е+01 а
1.118Е+01 1
1.105Е+01 в 1.930Е+00
Ш] ж
и
Рис. 8. Зона загрязнения в помещении, t = 26 с (сечение у = 3,6 м)
1.191Е+В1 □ 1.179Е+01 г 1.167Е+В1 4 1.154Е+В1 1 1.142Е+01 л 1.130Е+В1 а 1.117Е+В1 I ).105Е+В1 е 1.925Е+ВВ 1.В02Е+00 1.В7ЯЕ+0В у 1.555Е+0В 1.432Е+В0 ).308Е+В0 1.1Б5Е+ВВ 1.617Е-В1
шш
зии Их — И у —И- — 0,2 м
Рис. 9. Зона загрязнения в помещении для момента времени t = 32 с (сечение у =3,6 м)
Хорошо заметно (рис. 8, 9) как вытягивается фронт зоны загрязнения в направлении выходного отверстия вентиляции. Отчетливо видно,
как формируются подзоны загрязнения в плохо проветриваемых областях: перед оборудованием, за оборудованием и между оборудованием.
На базе 3Б-численной модели имеется возможность прогноза величины концентрации в любой интересующей точке в помещении. Например, концентрация сероводорода перед первым оборудованием (точка А, рис. 5) изменяется с течение времени следующим образом: г = 20 с С = 0,092 г/м3 г = 24 с С = 0,148 г/м3 г = 32 с С = 0,334 г/м3
Принимая во внимание, что ПДК для сероводорода составляет 10 мг/м3, то как видно из представленных данных концентрация загрязнителя в помещении очень быстро превысит это пороговое значение. Таким образом, за такой промежуток времени после аварии эвакуация людей из помещения может быть не осуществлена, и при отсутствии индивидуальных средств защиты следует ожидать токсичное поражение людей внутри помещения.
Отметим, что для решения рассмотренной сопряженной задачи прогноза динамики загрязнения воздушной среды как на промышленной площадке, так и внутри помещения, и при использовании разработанной 3Б численной модели потребовалось около 1 мин. времени работы ПК.
Выводы
В работе рассмотрено численное решение сопряженной задачи - перенос токсичного газа в атмосфере при аварии на химически опасном объекте и процесс загрязнения воздуха в помещении при затекании в него токсичного газа через систему вентиляции. Вычислительный эксперимент проведен на базе построенной численной модели. Проведенный вычислительный эксперимент показал эффективность построенной численной модели - оперативное получение прогнозных данных. При этом в модели осуществляется учет основных физических факторов, влияющих на процесс переноса загрязнителя. Построенная численная модель позволяет на новом уровне осуществлять оценку риска поражения людей на химически опасных объектах при авариях и терактах, что является важным определении уровня безопасности промышленного объекта. Модель может использоваться при разработке ПЛАСА для получения научно-обоснованной информации относительно масштаба поражения.
Дальнейшее совершенствование данного направления необходимо вести по созданию
численной модели для расчета рассеивания
бактериологических загрязнителей на промышленных объектах.
БИБЛИОГРАФИЧЕСКИЙ СПИСОК
1. Андерсон, Д. Вычислительная гидромеханика и теплообмен [Текст] / Д. Андерсон, Дж. Р. Тан-нехилл Плетчер. - М.: Мир, 1990 - 725 с.
2. Беляев, Н. Н. Расчет распространения загрязняющих веществ в условиях застройки [Текст] / Н. Н. Беляев, Е. Ю. Гунько // Вкник Дшпро-петр. нац. ун-ту залiзн. трансп. iм. акад. В. Ла-заряна. - 2007. - Вип. 18. - Д.: Вид-во ДНУЗТ, 2007. - С. 21-24.
3. Берлянд, М. Е. Прогноз и регулирование загрязнения атмосферы [Текст] / М. Е. Берлянд. - Л.: Гидрометеоиздат, 1985. - 273 с.
4. Бесчастнов, М. В. Предупреждение аварий в химических производствах [Текст] / М. В. Бес-частнов, В. М. Соколов. - М.: Химия, 1979. -392 с.
5. Бесчастнов, М. В. Аварии в химических производствах и меры их предупреждения [Текст] / М. В. Бесчастнов, В. М. Соколов, Н. И. Кац. -М.: Химия, 1976. - 368 с.
6. Бруяцкий, Е. В. Теория атмосферной диффузии радиоактивных выбросов [Текст] / Е. В. Бруяц-кий. - К : Ин-т гидромеханики НАН Украины, 2000. - 443 с.
7. Гунько, Е. Ю. Оценка риска токсичного поражения людей при аварийном выбросе химически опасного вещества [Текст] / Е. Ю. Гунько // Вюник Дншропетр. нац.. ун-ту залiзн. трансп. iм. акад. В. Лазаряна. - 2008. - Вип. 20. - Д.: Вид-во ДНУЗТ, 2008. - С. 87-90.
8. Численное моделирование распространения загрязнения в окружающей среде [Текст] / М. З. Згуровский [и др.]. - К.: Наук. думка, 1997. - 368 с.
9. Марчук, Г. И. Математическое моделирование в проблеме окружающей среды [Текст] / Г. И. Марчук. - М.: Наука, 1982. - 320 с.
10. Меньшиков, В. В. Опасные химические объекты и техногенный риск [Текст] / В. В. Меньшиков, А. А. Швыряев. - М.: Изд-во МГУ, 2003. -245 с.
11. Методика прогнозування наслщшв виливу (ви-киду) небезпечних хiмiчних речовин при аварь ях на промислових об'ектах 1 транспорта [Текст]. - К., 2001. - 33 с.
12. Логачев, И. Н. Аэродинамические основы аспирации [Текст] / И. Н. Логачев, К. И. Логачев. -СПб.: Химиздат, 2005. - 659 с.
13. Самарский, А. А. Теория разностных схем [Текст] / А. А. Самарский. - 2-е изд., испр. - М.: Наука, 1983. - 616 с.
14. Эльтерман, В. М. Вентиляция химических производств [Текст] / В. М. Эльтерман. - 3-е изд., перераб. - М.: Химия, 1980. - 288 с.
15. 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.
16. Development of building resolving atmospheric CFD code taking into account atmospheric radiation in complex geometric [Text] / Y. Qu [et al.] // Conf. Abstracts of 31st NATO / SPS Int'l Technical Meeting on Air Pollution Modelling and it's Application (27 Sept. - 01 Oct., Torino, Italy, 2010). - № P1.5.
17. 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 Modelling and it's Application (27 Sept. - 01 Oct., Torino, Italy, 2010). - № 2.17.
18. The use of MM5-CMAQ air pollution modeling system for real time and forecasting air quality impact of industrial emissions / R. San Jose [et al.] //
Advances in Air Pollution Modeling for Environment Security. NATO Science Series, - Springer, 2004. - Vol. 54. - P. 326-335.
19. Murakami, S. Overview of turbulence models applied in CWE - 1997 [Text] / S. Murakami // J. of Wind Engineering and Industrial Aerodynamics. -74-76 (1998). - Р. 1-24.
20. Murakami, S. Past, present and future of CWE: The view from 1999 [Text] / S. Murakami, A. Mochida // Wind Engineering into the 21st Century. - Rotterdam, Balkena, 1999. - P. 91-97.
21. Murakami, S. Comparison of "k-e" Model, ASM and LES with wind tunnel test for flow field around cubic model [Text] / S. Murakami, A. Mochida, H. Yoshihiko // 8th Int'l. Conf. on Wind Engineering (Western Ontario, July 8-11, 1991). -№ 12.
Поступила в редколлегию 11.05.2011.
Принята к печати 13.05.2011.