CC BY
58
УДК 519.6:504.3.054
П. Б. МАШИХИНА (ДИИТ)
МОДЕЛИРОВАНИЕ ЗАГРЯЗНЕНИЯ АТМОСФЕРЫ С УЧЕТОМ РЕЛЬЕФА МЕСТНОСТИ
На 6a3i тривимiрноï чисельно1 моделi виконано розрахунок процесу поширення токсично1 речовини в aтмосферi з урахуванням рельефу. Модель базуеться на чисельному iнтегрувaннi рiвняння конвективно-дифузiйного переносу домiшки та моделi потенцiйного руху. Наводяться результати обчислювального екс-перименту.
На базе трехмерной численной модели выполнен расчет распространения токсичного газа в атмосфере с учетом рельефа. Модель основывается на численном интегрировании уравнения конвективно-диффузионного переноса примеси и модели потенциального течения. Приводятся результаты вычислительного эксперимента.
The 3D numerical model was used to simulate the toxic gas dispersion over complex terrain. The model is based on the ^-gradient transport model and the model of potential flow. The results of numerical experiment are presented.
Введение
В настоящее время актуальным является вопрос разработки эффективных методов прогноза уровня загрязнения атмосферы при авариях, связанных с эмиссией токсичных веществ [1]. Нормативная методика, используемая для прогноза последствий при авариях с опасными веществами, дает возможность только оценить размеры (площадь) зоны загрязнения и не позволяет рассчитать концентрацию токсичного газа на различном расстоянии от места аварии [5]. Повышение требований к качеству прогнозной информации, с целью получения научно-обоснованных расчетных данных, ставит задачу создания моделей, учитывающих такой важный фактор как рельеф местности, поскольку он оказывает самое существенное влияние на процесс рассеивания токсичных веществ в атмосфере. В этом случае рассматриваемые задачи переноса примеси в атмосфере характеризуются такими особенностями:
• сложность реализации в модели граничных условий на поверхности рельефа;
• необходимость учета деформации поля скорости при обтекании рельефа;
• изменение с высотой скорости ветра, величины коэффициента турбулентной диффузии;
• необходимость быстрой перенастройки модели на решение другой задачи (иная конфигурация рельефа местности, формы облака и т.д.).
Именно поэтому применение аналитических моделей для решения задач переноса токсич-
ных газов в атмосфере с учетом рельефа исключено. Расчет рассеивания загрязнителя в атмосфере с учетом реального рельефа местности (немодельные задачи) может быть выполнен только с помощью численных моделей (СББ модели) [8 - 11]. При этом для задач рассматриваемого класса немаловажным является вопрос получения прогнозной информации за минимальное расчетное время (порядка нескольких минут). Целью настоящей работы является применение разработанной автором численной модели и созданного на её основе кода для решения задачи прогноза загрязнения атмосферы при миграции облака токсичного газа над поверхностью со сложным рельефом. Модель позволяет рассчитать процесс переноса загрязнителя в течение нескольких секунд на компьютерах средней мощности.
Математическая модель
Для моделирования процесса переноса токсичного газа в атмосфере будем использовать трехмерное уравнение переноса примеси [2, 3]
дС диС дуС д(^ - wS)С — +-+-+ —-—— =
дt дх ду дг
д , дС. д , дС. д , дС. = — (ц* —)+—(ц у —)+—(ц2 —) +
дх дх ду ду дг дг
+
I й ( t Щ r - r ).
(1)
где С - концентрация примеси, попадающего в атмосферу при аварии; и, V, ^ - компоненты вектора скорости воздушной среды; - скорость оседания примеси (в данной работе при-
© Машихина П. Б., 2010
нята равной нулю); ц = (ц , цу , ц) - коэффициент турбулентной диффузии; Q - интенсивность выброса токсичного вещества; 5( г - гг) -
дельта-функция Дирака; гг= (хг, у , гг) - координаты источника выброса.
Если рассматривается задача о миграции облака, то в модели задается форма облака и концентрация загрязнителя в нем.
В разработанной математической модели используется неравномерный профиль скорости ветра и вертикального коэффициента диффузии:
ду+С ]+1,,]-1,к
д у Ау
^ } +1,ксг,}+1,к ~хг]ксг]к
д у А у
д м>+с ^ },к+1сг}к - ™1]ксг, },к-1
д А.г
д м>~с , ,к+1сг,},к+1 - ™1]ксг,},к _
д г
А г
=гс
= гсп
п+1
и = и
( г >
V г1 У
; vz = 0,11*,
где и1 - скорость ветра на высоте г1 = 10 м; п = 0,15.
Постановка краевых условий для данного уравнения рассмотрена в работах [4, 2].
Метод решения
Выполним следующую аппроксимацию производных. Заменим производную по времени разделённой разностью «назад»:
дс ск - сп
ы
"У'к ук
А!
Конвективные производные представим в виде:
Компоненты вектора скорости будут рассчитываться, при решении гидродинамической задачи, на гранях разностной ячейки (контрольного объема). Вторые производные аппроксимируем следующими выражениями:
д( дС ^ - Сп1з, к - с"+к
д х
д х
V х
V /
С"+1 — с"+1
к ^1-1,}, к
А х2
Ах2
= Ы~сп+1 + М+хсп+1;
д ( дсл
V
ду V у ду у
V у
V у
с "+1 — с "+1
и,к и-1,к = М—Усп+1 + М+П,сп+1'
с п+1 - с"+1
с1,}+1,к с1}к
Ау2
А у2
уу
уу
диС ди+С ди С
дх =--1 дх ь-: дх
дуС ду+С ду~С
ду ду ду
дм>С С дм>~ С -
дг дг дг и + |и| и - |и| + У + Ы
где и =■
2
V - IV
; и =■
V = ■
; =
2
м +М
V = -
2
м -Ы
ди+С и1+\,},к С"}к - ичк С"+
д х
,}, к г }к ¡}к г-1,}, к
А х
= Г с
п+1 .
ди С иг+1,},к Сг"-l\},к иг}к С}
д х
А х
= Г С
п+1
д д г
г
( дс >
V
сп+1 сп+1
',},к ~ г},к-1 = М-гсп+1 + М+С+1
А2 2 22 22
V,
сп+1 - сп+1
}, к+1 г }к
А?
В данных выражениях Гх, Ьх, Гу, Ьу, Ь+, Ьу ,
М+,М-х и т.д. - обозначения разностных операторов. С учетом этих обозначений разностный аналог трехмерного уравнения переноса примеси будет иметь вид:
2 2 2 Для аппроксимации конвективных производных используем выражения:
Сг}к - Сг}к
А(
ГСп+1 + Г Сп+1 + Г, С
п+1
+ ГС п+1 + Г С п+1 + ГС п+1 + о Сп+ =
у г г г} к
= (М+ С п+1 + М-С п+1 + М+ Сп+1 +
\ хх хх уу
+ М- С п+1 + М + С п+1 + М- с п+1).
уу
Расщепим решение данного разностного уравнения при интегрировании на временном интервале с1г так:
, 1
- на первом шаге к = — :
С п+к _С п 1
11 _ 11 + - (( Ск +1+ Ск +1+ Ск) + Аг 2V у г '
+ | Ск]к = - (М+ Ск + М_ Сп + МУу Ск + + МУу Сп + М+ С + М_ Сп ); (3)
7 1 1
- на втором шаге к = п + ^ с = п + ~ :
Ск1к _СГ/к + 1 (( Ск + Ц Ск + Ц Ск)+ ° Ск, = Аг 2\ х у г ' 4 11
=- (МХхСк + М+Сс + М_у Ск + + МУу Сс + М_ Ск + М+ Сс);
д2Р д2 Р д2 Р
(4)
/ 3 1
- на третьем шаге к = п + 4; с = п + 2 используется формула (4);
7 1 3
- на четвертом шаге к = п +1; с = п + 4 используется формула (3).
На пятом шаге расщепления учитывается влияние источников на изменение концентрации примеси, и расчетная зависимость на этом шаге имеет вид:
5 п+1 5 п
С1,1, к_ С1,1, к =у 4 (г"+1/2) 5
А г А х Ау А г 1
В дискретном виде дельта-функция Дирака «размазывается» по объему разностной ячейки с учётом сбережения суммарного количества выбрасываемого загрязнения. Функции тождественно равняются нулю, кроме ячеек, где расположен 1-й источник загрязнения. На каждом шаге расщепления расчет неизвестной концентрации осуществляется по явной формуле бегущего счета.
Для расчета поля скорости воздушного потока при обтекании сложного рельефа используется модель потенциального течения. Применение данной модели позволяет, в отличие от гидродинамических моделей, учитывающих вязкость, рассчитать поле скорости воздушной среды в течение нескольких секунд. Уравнение для определения потенциала скорости имеет вид [3]:
дх2 ду2 dz2
Компоненты вектора скорости рассчитываются по зависимостям
дР дР u = —; v = — и т.д.
д х ду
Постановка граничных условий для данного уравнения рассмотрена в работе [3]. Для численного интегрирования этого уравнения используется идея установления решения по времени, Численное интегрирование проводится с помощью неявной схемы условной аппроксимации [7].
В построенной численной модели используется метод маркирования, позволяющий выделить разностные ячейки, которые относятся к области течения и к твердым границам (рельеф). Данный метод называют за рубежом «porosity technique» [8]. Применение этого метода дает возможность построить эффективный алгоритм расчета, при котором можно задавать любую форму рельефа местности.
Тестовые расчеты
Для тестирования блока модели, осуществляющего расчет поля потенциала скорости, было найдено численное решение уравнения теплопроводности вида:
дТ_
Г д 2т
■ = a
д2T д T
2
дх2 ду2 дz
которое
области 0 < х <<, 0 < у <<х>, 0 < г << при граничных условиях Т = 0 на границах х = 0, у = 0, г = 0 и начальном условии Т = Т0 имеет аналитическое решение вида [6]:
T = To erf
erf
2*J~at) v 2*J~at
у
erf
2\Jat ,
Для численного решения задачи необходимо ограничить область расчета, поэтому при численном интегрировании этого уравнения с помощью неявной схемы условной аппроксимации использовались такие данные: длина области в направлении оси х - 9,6 м; в направлении оси у - 8 м; в направлении оси г - 8 м; а = 1, координаты расчетной точки: х = 1 м; у = 1 м; г =1 м. Т0 = 100 °С - начальная температура.
Результаты расчета температуры в заданной точке и полученные на основе численного мо-
делирования и с помощью аналитическом модели представлены в табл. 1.
Таблица 1 Динамика изменения температуры
г, с Аналитическое решение Численное решение, сетка 24*21*21 узел Численное решение, сетка 48*42*42 узла
9,76 0,573 0,593 0,587
12,76 0,386 0,405 0,398
17,30 0,240 0,254 0,248
зоны: г = 1.6 м (точка х =297 м, у = 146 м). Принимая во внимание, что для данного токсичного вещества ПДК составляет 20 мг/м3, то очевидно, что в случае аварии, для рассматриваемого сценария, последствия будут крайне тяжелыми, со смертельным исходом. Максимальная концентрация токсичного газа будет значительно превосходить ПДК. Наиболее интенсивное загрязнение селитебной зоны будет примерно через 50... 60 с после выброса. Совершенно очевидно, за это время люди не успеют уйти из зоны поражения.
В данной таблице представлены результаты, полученные на двух сетках: 24*21*21 узел и 48x42x42 узла. Представленные данные позволяют отметить, что разработанная численная модель дает удовлетворительное согласование данных моделирования с данными, полученными аналитически.
Практическая реализация модели
Рассматривается сценарий аварии возле вокзала ст. Днепропетровск-Южный (рис. 1), в результате которой в атмосферу произошел выброс токсичного газа - аммиака. Место аварии показано на рис. 1 «кружком». Вокзал расположен на склоне, внизу склона - жилая зона. Считаем, что на месте аварии образовалось облако сложной формы (рис. 2), которое начинает двигаться, под действием ветра, вниз по склону, в направлении жилмассива Победа-1. Средняя высота склона составляет порядка 10 м. Скорость ветра составляет 5,3 м/с (на высоте 10 м).
Численное моделирование выполнено при следующих параметрах: размеры расчетной области 420*320*40 м, коэффициент диффузии в направлении осей х, у - 0,5 м2/с. В облаке находится около 2 т токсичного вещества.
На рис. 2 - 4 показана зона загрязнения атмосферы для различных моментов времени. Как видно из данных рисунков, облако, образовавшееся на месте аварии, быстро размывается, теряя первоначальную форму, и вытягивается вдоль склона.
Для практики, с целью прогноза масштаба опасности аварии, представляет интерес оценка концентрации токсичного газа возле селитебной зоны.
В табл. 2 представлены прогнозные значения концентрации токсичного газа внизу склона, перед застройкой, на уровне дыхательной
Рис. 1. Место аварии на ст. Днепропетровск-Южный
Рис. 2. Форма зоны загрязнения атмосферы, г = 0,0001 с (сечение у = 146 м)
Рис. 3. Зона загрязнения атмосферы г = 3 с (сечение у = 146 м)
В табл. 3 представлены данные о динамике изменения концентрации загрязнителя в расчетной точке для варианта, когда масса токсичного газа в облаке составляла порядка 670 кг. Учитывая большие объемы перевозок опасных
грузов железной дорогой, можно считать такой выброс «малым».
0.167Е+01
0.175Е+02 coordinate х 0.403Е+03
Рис. 4. Зона загрязнения атмосферы t =7 с (сечение y =146 м)
Таблица 2
Динамика изменения концентрации токсичного газа вблизи селитебной зоны (мощность выброса 2 т)
t, c 24 50 60
С, г/м3 4,32 13,41 11,82
Таблица 3
Динамика изменения концентрации токсичного газа вблизи селитебной зоны (мощность выброса 670 кг)
t, c 24 50 60
С, г/м3 1,38 4,30 3,80
Как видно из табл. 3, даже при такой мощности эмиссии концентрация токсичного газа в селитебной зоне будет превышать величину ПДК. Следовательно, жилая застройка у склона, прилегающего к железной дороге, находится в опасной зоне.
Отметим, что для решения данной задачи потребовалось 10 с.
Выводы
В работе на основе разработанной численной модели и созданного на ее основе кода выполнен прогноз уровня загрязнения атмосферы при аварии на железнодорожной станции. Расчеты показали, что в случае аварийной ситуации произойдет крайне опасное загрязнение воздушной среды в жилой зоне, прилегающей к месту аварии. Дальнейшее совершенствование рассмотренной в работе модели необходимо
проводить в направлении ее адаптации к моделированию процесса нейтрализации токсичного
газа.
БИБЛИОГРАФИЧЕСКИЙ СПИСОК
1. Аварии и катастрофы. Предупреждение и ликвидация последствий [Текст] : учеб. пособие в 5-ти кн. / под ред. В. А. Котляревского и А. В. Забегаева. - М.: Изд-во АСВ, 2001. -200 с.
2. Численное моделирование распространения загрязнения в окружающей среде [Текст] / М. З. Згуровский и др. - К.: Наук. думка, 1997. - 368 с.
3. Лойцянский, Л. Г. Механика жидкости и газа [Текст] / Л. Г. Лойцянский. - М.: Наука, 1978. -735 с.
4. Марчук, Г. И. Математическое моделирование в проблеме окружающей среды [Текст] / Г. И. Марчук. - М.: Наука, 1982. - 320 с.
5. Методика прогнозування насладив впливу (ви-киду) небезпечних х1м1чних речовин при аварь ях на промислових об'ектах i транспорт! [Текст]. - К., 2001. - 33 с.
6. Справочник по точным решениям уравнений тепло- и массопереноса [Текст] / А. Д. Полянин и др. - М.: Факториал, 1998. - 368 с.
7. Самарский, А. А. Теория разностных схем [Текст] / А. А. Самарский. - М.: Наука, 1983. -616 с.
8. Castino, F. Parameterization of Convective and Stable Interval Boundary Layers into Mass Consistent Models [Text] / F. Castino, M. Tombrou // 2 EACWE (Genova, Italy, 1997). - P. 309-314.
9. Computation of Wind Flow and Air Pollution for Regions Having a Complex Topography [Text] / Mahfound Kadja, Anagnostopoulos et al. // Proc. of 3rd European & African Conf. on Wind Engineering (Eindhoven Univ. of Technology, Netherlands, July 2-6, 2001). - P. 355-358.
10. Maurizi, A. Numerical Simulation of the Atmospheric Flow in a Mountainous Region of North Portugal [Text] / A. Maurizi, J. M. L. M. Palma, F. A. Castro // Proc. of 2nd EACWE (Genova, Italy, 1997). - Vol. 1. - P. 293-300.
11. Montavon, C. Validation of Non-hydrostatic Numerical Model to Simulate Stratified Wind Fields over Complex Topography [Text] / C. Montavon // 2 EACWE (Genova, Italy, 1997). - P. 357-362.
Поступила в редколлегию 14.01.2010.
Принята к печати 20.01.2010.
