Моделирование загрязнения атмосферы при испарении жидкости из грунта Текст научной статьи по специальности «Математика»

УДК 519.6

Е. Ю. ГУНЬКО (ДИИТ)

МОДЕЛИРОВАНИЕ ЗАГРЯЗНЕНИЯ АТМОСФЕРЫ ПРИ ИСПАРЕНИИ ЖИДКОСТИ ИЗ ГРУНТА

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

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

The 3D model to simulate the pollutant transfer near building after accidents was developed. The model is based on the transport gradient model and the model of invisid imcompressible fluid. The results of numerical experiments are presented.

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

Q = (5,83 + 4,1Г)Рну[м ,

где V - скорость ветра; Рн - давление насыщенных паров; М - молекулярная масса вещества.

Важно отметить, что данная зависимость была получена для случая испарения жидкости со свободной поверхности, но ввиду отсутствия модели, позволяющей учесть процесс испарения из грунта, эта же формула применяется на практике при оценке количества испарившегося вещества из грунта. Целью данной работы является разработка математической модели для моделирования процесса испарения жидкости из грунта с учетом тех основных физических факторов, которые влияют на данный процесс. Более конкретно - в работе рассматривается решение сопряженной задачи - «испарение из грунта + перенос паров загрязнителя в атмосфере», с учетом деформации зоны загрязнения в атмосфере, встречающей на своем пути препятствия - здания (рис. 1). С этой целью, к математической модели испарения жидкости из грунта добавляется математическая модель переноса паров загрязнителя, выделившихся из грунта, в атмосфере, с учетом зданий на пути миграции паров. Анализ научных публикаций, связанных в той или иной мере с рассматриваемой задачей, показал отсутствие научных разработок, посвященных решению данной задачи в той сложной постановке, которая рассматривается в данной работе, т. е. в формате «испарение из грунта + перенос паров загрязнителя в атмосфере с учетом влияния на этот перенос зданий».

Математическая модель испарения жидкости из грунта. Будем считать, что в результате аварийного разлива в грунте образовалась зона глубиной к, пропитанная жидкостью.

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

дС

д 2С

п"дГ=

(1)

где ц - коэффициент диффузии паров в подзоне 1; С - концентрация паров в этой подзоне; z0 - координата; п - объемная пористость грунта; t - время; пе - поверхностная пористость.

Ось г0 направлена вертикально вниз, величине = 0 соответствует положение поверхности земли, а величине г0 = к - глубина загрязнения грунта разлитой жидкостью. Граничные условия для данного уравнения имеют вид: на нижней границе (граница раздела подзоны 1 и подзоны 2) ставится условие С = Сн, где Сн - концентрация насыщенных паров разлитой жидкости при заданной температуре Т0. Величина Сн рассчитывается по зависимости

[5]:

16МР

Сн =-

Т

где М - молекулярный вес вещества; Рн - давление насыщенных паров; Т - температура.

На границе г0 = 0 (граница раздела «грунт-

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

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

Для экспресс-прогноза можно воспользоваться моделью Поясова Н. П. [2]:

Ц = Г|(Х-Х 0 V 0,

где п - коэффициент извилистости; X - общая пористость; X 0 - объем пор, не участвующих в процессе диффузии; ц 0 - коэффициент диффузии паров загрязнителя в атмосфере.

С другой стороны, при экспресс-прогнозе можно считать, что коэффициент ц примерно в 3.. .5 раз меньше коэффициента ц 0 [2].

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

2 = ж—,

* дг

где 5" - площадь поверхности зоны разлива; Л -промежуток времени; г - вертикальная ось в системе декартовых координат х, у, г, которая применяется для моделирования процесса переноса паров в атмосфере после их «выхода» из грунта.

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

Для численного интегрирования уравнения (1) будем использовать неявную схему условной аппроксимации [5]. В этой схеме на каждом шаге расщепления расчет неизвестной величины С осуществляется по явной схеме бегущего счета.

Математическая модель миграции паров загрязнителя в атмосфере. Для расчета процесса загрязнения атмосферы при миграции в ней паров жидкости, выделившихся вследствие испарения жидкости из грунта, и в условиях застройки будем использовать трехмерное уравнение миграции примеси (модель градиентного типа) [2-4]:

дС диС дуС д( - ws )С =

дt дх

ду

дС

= -(Ц

дх ^ х дх

ц

дг

дС

ду [ ду

1 ^ z

C

X Q (t )d(r - r),

(2)

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

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

u =-

dP_

дх

dP

dP

dy ' dz

д 2P д 2P д 2P

■ + -

dx2 dy2

dz2

= 0

(3)

Постановка краевых условий для уравнения (2) рассмотрена в работе [3]. Для уравнения (3) ставятся следующие граничные условия: • на твердых стенках:

дР=0,

дп

где п - единичный вектор внешней нормали; на входной границе (границы втекания воздушного потока в помещение):

dP_ дп

= Vn

u = u

(z у (z \m

, Hz = kl

(zl

( zi У

где и1 - скорость ветра на высоте ; к = 0,2 ; п = 0,16 ; т « 1.

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

Для численного интегрирования уравнения (3) используется идея установления решения по времени, т. е. интегрируется уравнение вида:

дР д 2Р д 2Р д 2Р • + ——+ -

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

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

дт дх2

dy2

&2

где Vn - известное значение скорости; • на выходной границе:

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

Процессами теплообмена будем пренебрегать.

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

где т - фиктивное время.

Численное интегрирование данного уравнения проводится с использованием неявной разностной схемы условной аппроксимации. В разработанном пакете программ осуществляется также численное интегрирование уравнения (3) с помощью метода Либмана и метода условной аппроксимации. Расчет поля потенциала проводится одновременно с использованием указанным методов с целью контроля результатов расчета.

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

Рассматривается следующая постановка задачи. Размеры расчетной области 60^60^30 м; скорость ветра 5,5 м/с (на высоте 10 м). Профиль ветра - неравномерный; вертикальный коэффициент диффузии изменяется с высотой. Площадь зоны разлива составляет 50 м2; глубина смачивания грунта нитробензолом - 0,35 м. Температура воздуха 200С; принимаем, что такая же температура у разлитого нитробензола. На площадке расположено 2 здания: одно здание прямоугольной формы, а другое здание Г-образной формы (рис. 1.), высотой 7,5 м. Объемная пористость грунта - 0,2; коэффициент

диффузии паров принимается равным = 0,005 м2/с.

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

г = 50 с....................................................О = 8,7 г

г = 70 с...................................................О = 12,1 г

г = 100 с................................................О = 17,4 г

г = 120 с................................................О = 20,8 г

г = 4960 с............................................О = 863,2 г

г = 29760 с......................................О = 5182,3 г

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

г

ТБ = | Сйг,

0

где С - концентрация токсичного газа в точке расположения рецептора (человека) на маршруте эвакуации; г - время.

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

Перейдем к рассмотрению результатов проведенного вычислительного эксперимента. Рассмотрим сначала вопрос о токсичном поражении людей при их эвакуации. Будем считать, что персонал передвигается по маршруту эвакуации со скоростью 2 м/с.

Динамика изменения величины токсодозы для рассматриваемого сценария выглядит следующим образом:

г = 50 с.............................ТБ = 0,00097 мг-мин/л

г = 70 с...............................ТБ = 0,0012 мг-мин/л

г = 100 с.............................ТБ = 0,0015 мг-мин/л

г = 120 с...............................ТБ = 0,017 мг-мин/л

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

На последующих рисунках представлена динамика загрязнения воздушной среды на промплощадке после аварии. На рис. 3-5 отчетливо видно как шлейф, формирующийся над зоной разлива, огибает промышленные здания. Из этих рисунков также видно, что при этом произойдет попадание загрязняющего вещества через окна и вентиляционные отверстия в помещения. На рис. 6 показано как загрязняющее вещество «обходит» здания сбоку и полностью «накрывает» маршрут эвакуации, проходящий возле зданий (зона А, см. также рис. 1). На рис. 7 показана зона загрязнения, формирующаяся на промплощадке при разливе на грунт, но уже в другом месте - между зданиями. Из этого рисунка хорошо видно, что формирующийся шлейф также накроет часть маршрута эвакуации.

Рис. 1. Схема размещения объектов (план): 1 - зона разлива (на грунт); 2 - здание; 3 - маршрут эвакуации людей

20=0

20=Ь

Рис. 2. Схема зоны загрязнения грунта

2

Рис. 3. Зона загрязнения для момента времени t = 10 с (сечение у = 60 м)

Рис. 4. Зона загрязнения для момента времени t = 20 с (сечение у = 60 м)

Рис. 5. Зона загрязнения для момента времени t = 40 с (сечение у = 60 м)

Рис. 6. Зона загрязнения для момента времени t = 60 с (сечение г = 10 м)

Отметим, для решения каждой задачи требуется около 10 с машинного времени ПК.

Выводы. В работе разработана трехмерная численная модель, на основе которой построен

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

1,142Е+02 о 1.133Е+В2 г 1.123Е+02 а 1.114Е+В2 1 1.105Е+В2 п 1.960Е+01 а 1.868Е+01 ± 1.777Е+01 е I 68БЕ+01 1.594Е+01 1.503Е+01 у 1.411Е+01 1,320Е+В1 1.229Е+01 1.137Е+01 1.457Е+00

0.125Е*01

—>

соопНпа'Ье

Рис. 7. Зона загрязнения для момента времени t = 180 с (сечение г = 10 м)

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

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

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

2. Голубев В. С., Гарибянц А. А. Гетерогенные процессы геохимической миграции. - М.: Недра, 1984 - 262 с.

3. Численное моделирование распространения загрязнения в окружающей среде / М. З. Згуров-ский, В. В. Скопецкий, В. К. Хрущ, Н. Н. Беляев. - К.: Наук. Думка, 1997 - 368 с.

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

5. Мацак В. Г., Хоцянов Л. К. Гигиеническое значение скорости испарения и давления пара токсических веществ, применяемых в производстве. - М.: Медгиз, 1959 - 231 с.

6. Самарский А. А. Теория разностных схем. -М.: Наука, 1983 - 616 с.

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