CC BY
61
УДК 004.9:504:519.6
А.В. ХАЛЧЕНКОВ*, И.В. КОВАЛЕЦ*
ИСПОЛЬЗОВАНИЕ ВЫЧИСЛИТЕЛЬНОЙ ГИДРОДИНАМИЧЕСКОЙ МОДЕЛИ OPENFOAM ДЛЯ ОЦЕНКИ ВЛИЯНИЯ ЗАГРЯЗНЕННЫХ ЗДАНИЙ НА ПРИЛЕГАЮЩИЕ ТЕРРИТОРИИ
Институт проблем математических машин и систем НАН Украины, Киев, Украина
Анотаця. У po6omi представленi результати застосування програмного комплексу OpenFoam для оцтки середньорiчних характеристик забруднення навколо будiвель колишнього уранового ви-робництва. Використано «конструктор моделей» OpenFoam, за допомогою якого iснуючi моделi були адаптоваш для опису вивiтрювання забруднення з виробничих будiвель шляхом введення до-даткових члетв у вихiднi рiвняння. Показано високу ефективтсть використаних чисельних мето-дiв. Представлено результати розрахунюв безрозмiрних концентрацт поблизу будiвель. Ключевые слова: обчислювальна гiдродинамiка, атмосферний перенос, оцтка безпеки.
Аннотация. В работе представлены результаты применения программного комплекса OpenFoam для оценки среднегодовых характеристик загрязнения в окрестности зданий бывшего уранового производства. Использован «конструктор моделей» OpenFoam, с помощью которого существующие модели были адаптированы для описания выветривания загрязнения из производственных зданий путем введения дополнительных членов в исходные уравнения. Показана высокая эффективность используемых численных методов. Представлены результаты расчетов безразмерных концентраций вблизи зданий.
Ключевые слова: вычислительная гидродинамика, атмосферный перенос, оценка безопасности.
Abstract. The paper presents the results of OpenFoam software application to assess the average annual pollution characteristics around the buildings of the former uranium industrial plant. It was used OpenFoam "designer of models " with the help of which the existing models were adapted to describe emission of the pollutants from the industrial buildings by introducing additional terms in the original equations. High numerical efficiency of the model was demonstrated. The calculated non-dimensional concentrations near buildings were presented.
Keywords: computational fluid dynamics, atmospheric dispersion, safety assessment. 1. Введение
Загрязненные здания опасных производств могут оказывать неблагоприятное влияние на окружающую среду. Для оценки воздушного загрязнения в окрестности таких зданий необходимо использовать математические модели. Если для оценки загрязнения на расстояниях несколько сот метров можно использовать упрощенные гауссовы модели, то для оценки загрязнения в непосредственной близости от зданий (расстояния от нескольких метров до нескольких десятков метров) требуется использовать модели вычислительной гидродинамики для детального учета влияния геометрии зданий на радиоактивное загрязнение.
В предыдущей работе [1] для этой цели была разработана нестационарная гидродинамическая модель на основе полных уравнений гидродинамики. Эта модель может быть использована для оценки загрязнения при конкретных метеорологических сценариях. Однако для получения среднегодовых значений загрязнения требуется проводить массивные расчеты. Это, в свою очередь, налагает следующие дополнительные требования к математическим моделям:
1) возможность упрощения уравнений модели;
2) возможность расчетов на неструктурированных сетках, что позволяет сгущать сетку вокруг произвольно расположенных объектов;
© Халченков А.В., Ковалец И.В., 2014
ISSN 1028-9763. Математичш машини i системи, 2014, № 2
3) использование алгоритмов параллельных вычислений.
Всем этим требованиям удовлетворяет пакет программ ОрепБоат [2], свободно распространяемый через Интернет. ОрепБоат представляет собой свободораспространяемый в исходных кодах программный комплекс, реализующий различные гидродинамические модели. Здесь пользователю предоставляется так называемый «конструктор моделей», позволяющий модифицировать модели под особенности решаемой задачи.
Целью данной работы является применение программного комплекса ОрепБоат для оценки среднегодовых характеристик загрязнения в окрестности зданий бывшего уранового производства «Приднепровский химический завод» с использованием перечисленных выше возможностей данной модели.
2. Задача оценки среднегодовых характеристик загрязнения вблизи зданий уранового производства
Приднепровский химический завод (ПХЗ) был одним из крупнейших предприятий по переработке урана в Советском Союзе [3]. Это предприятие работало с 1948 по 1991 гг., а после распада СССР ПХЗ был разделен на несколько компаний, и обработка урана на нем прекратилась. На территории ПХЗ присутствуют загрязненные здания и другие объекты, которые были задействованы в производственном технологическом цикле. При этом промплощадка ПХЗ находится в непосредственной близости от жилой зоны г. Днепродзержинск (около 1 км), а на территории бывшего ПХЗ функционируют действующие предприятия.
Одним из наиболее опасных зданий на территории бывшего ПХЗ является здание № 103 (рис. 1), в непосредственной близости от которого, на расстоянии около 10 м с северной стороны, функционирует действующее предприятие в здании № 102.
Для обоснования мероприятий по демонтажу здания № 103 в рамках «Государственной программы по приведению опасных объектов ПО «ПХЗ» в экологически безопасное состояние и обеспечению защиты населения от вредного воздействия ионизирующего излучения на 2005-2014 гг.» требуется оценить возможное влияние этого здания на облучение людей, работающих в здании № 102. В качестве одного из вероятных сценариев такого влияния был выбран сценарий выветривания радионуклидов через разбитые окна в здании № 103. Эта задача может быть решена средствами математического моделирования с применением вычислительной гидродинамической модели, описанной в следующем разделе.
3. Метод решения
3.1. Уравнения модели
Рассмотрим вкратце основные уравнения и их реализацию в ОрепБоат [2], которые использовались в настоящей работе. Расчет поля скорости вокруг зданий осуществлялся путем решения стационарных уравнений гидродинамики несжимаемой жидкости, осреднен-ных по Рейнольдсу, без учета процессов теплообмена. В работе использовался модуль в1тр1еБоат, в котором решается следующий упрощенный вариант уравнений гидродинамики:
Здесь (1) - уравнение неразрывности, (2) - уравнения сохранения импульса, и -вектор скорости, р - возмущение давления относительно произвольного фонового значения, пт - эффективная кинематическая вязкость, определяемая по формуле
У и = 0, и-Уи = У • (пт Ум )-Ур.
(1) (2)
^ k2 vt = с, —,
m e
где к - кинетическая энергия турбулентности, е - скорость диссипации энергии турбулентности.
Для замыкания уравнений используется стандартная к — е модель турбулентности:
dk ди
и -= t —L
J dx J дх
■e + -
д_ дх
f
v + -
s
Л
dk
k У
дх
u de = с e du с e2 + d
J dxJ el k J dxJ e2 k dxJ
v + -
K
'e У
dk_ дх
(4)
(5)
где т. = (ди. /дх. + ди. /дх.). В работе использовались стандартные значения констант: Ст = 0,09; СА = 1,44; С£2 = 1,92; аъ = 1,0; аа = 1,3.
Система (1)-(5) дополняется следующими граничными условиями на верхней границе: м11 = мг1 (Н), и2| = иг2 (Н), и31 = 0, где мг1 (г),иг2 (г) - компоненты скорости ветра в невозмущенном приземном слое атмосферы, Н - высота вычислительной области. На входных боковых границах ставится условие равенства компонент скорости соответствующим значениям в невозмущенном приземном слое, на выходных границах - условие ди/дп = 0, п - нормаль к границе. На твердых границах задается поток импульса за счет турбулентного трения, соответствующий логарифмическому пограничному слою.
Решение уравнений модели турбулентности реализовано в специальной библиотеке моделей турбулентности; конкретная модель турбулентности задается значением входной переменной ЯЛ8Моёе1 (в нашем случае КЛ8Моёе1=кер811оп).
Расчет переноса загрязнения производился путем решения стационарного уравнения переноса с учетом процессов турбулентного перемешивания:
div(uC)-D((D0 +K )С)-qs = 0
(6)
где qs - объемная плотность источников, С - концентрация, D0 = l, 2 10-5 м2 / с - коэффициент молекулярной диффузии, vc - коэффициент турбулентной диффузии, Sc = vT / vc = 0,9 - число Шмидта.
Для решения уравнения (6) использовался модуль ScalarTransport. Модуль ScalarTransport решает уравнение (6) без учета источника. Поэтому с помощью конструктора OpenFoam был реализован стационарный решатель, учитывающий наличие объемных источников выброса. Для этого в исходном файле программы ScalarTransport, в котором производится решение уравнения (6), было добавлено слагаемое, описывающее объемные источники выброса:
fvm::div(phi, w_c) - fvm::laplacian(ADT*(D_ca+nut/0.9), w_c) - qs.
Приведенный отрывок кода на языке С использует функции модуля (namespace) fvm системы OpenFoam, в котором реализованы функции расчета Лапласиана (laplacian) и дивергенции (div) методом конечных объемов на произвольных (неструктурированных) сетках. Первый и второй члены программно реализовывают соответствующие члены уравнения (6), а последний член описывает объемную мощность источника.
Заметим, что, хотя в исходной постановке задачи источником является граница вычислительной области (окна зданий), преобразование поверхностного источника в объем-
ный осуществляется по следующей процедуре. В ячейки, прилегающие к границе области, от которой действует поверхностный источник, вводится объемный источник по формуле
а = — с|, (7)
ч, ¿У г « > V )
где интегрирование производится по поверхности ячейки, дУ - объем ячейки, ] - поток вещества через поверхность ячейки (без учета адвективного и диффузионного потоков).
Вычислительная сетка построена таким образом, что поверхности граничных ячеек совпадают с поверхностью стен зданий. Кроме того, вблизи зданий вычислительная сетка равномерна и состоит из прямоугольных ячеек размеров Нх, Ну, Н2 (см. ниже). Таким образом, формула (7) приводит к соотношению: д, = ] / Нх.
Поток ] через окна здания №103 вычисляется так же, как и в работе [1]:
] = ( щп ((Д • й )2+(Д • й )2 ")+Дуй) с.
Здесь С, - концентрация радионуклидов в воздухе внутри здания № 103. Векторы
коэффициентов эффективности вентиляции, с учетом того, что нормаль к стене здания №103, которая является источником, параллельна оси у (рис. 1), определены следующим образом: Д = ( 0.3,0,0); Ду = ( 0,0.8,0); Д =( 0,0,0.3).
3.2. Вычислительная сетка
Преимуществом программного комплекса ОрепБоаш является возможность решения уравнений гидродинамики и переноса вещества на произвольных сетках, включая неструктурированные сетки. Это позволяет сгущать сетку в областях, где необходимо получить более детальные данные о распределении расчетных параметров. В исследуемой задаче основной интерес представляет распределение концентраций вблизи зданий №№ 102, 103, 104. Поэтому в области, ограничивающей территорию этих зданий (называемая ниже областью застройки) по горизонтали (рис. 1) и по вертикали (нижний слой 28 м), была выбрана равномерная прямоугольная сетка.
По горизонтали за пределами области застройки ячейки имели трапецевидную форму и равномерно увеличивались в обоих направлениях при удалении от области застройки (рис. 1). Выше области застройки в вертикальном сечении сетка оставалась прямоугольной, но вертикальные размеры ячеек увеличивались с коэффициентом Н2. / Н2._ = 1,15.
Построение сетки выше области застройки осуществлялось путем экструдиро-
вания горизонтальной сетки с помощью утилиты ех^цёеМезЬ программного комплекса
Рис. 1. Вычислительная сетка в окрестности зданий №№ 102, 103, 104 ПХЗ
ОрепБоаш. Размеры всей вычислительной области в направлениях х, у, г были равны 200х300х300 м.
3.3. Метод расчета среднегодовых концентраций
Для получения среднегодовых характеристик загрязнения в общем случае требуется проводить непрерывный расчет атмосферного переноса при заданных характеристиках источника за период в несколько лет (например [4]). Но при использовании вычислительных гидродинамических моделей такой подход требует непомерно большого расчетного времени. Поэтому на практике для получения среднегодовых характеристик загрязнения проводят расчеты при различных характеристиках скорости и направления ветра, категории устойчивости и результаты усредняют, используя информацию относительно характеристик повторяемости тех или иных метеорологических условий.
В настоящей работе расчеты проводились без учета влияния стратификации на характеристики турбулентности. Соответствующая форма уравнений (1)-(6) приводит к подобию результатов, полученных при одинаковых направлениях ветра. Действительно, произведем замены размерных членов уравнений (1)-(6) на безразмерные, используя в качестве размерных параметров высоту области H и абсолютную величину скорости ветра
на верхней границе области Ur = ^/u^ (H) + u2r2 (H). Безразмерные характеристики будут
иметь следующий вид: x = x /H , y = y /H , z = z /H , Ui = ui / Ur, к = к / Ur2, e = eH / Ц3,
V = vT/(UrH), p = p/U2, vT =vT/(UrH), n =vc/(UrH), v=v/(UrH),
D0 = DJ [UrH ), c = c / Cs. Легко показатц что система (1)-(6) в безразмерн^1х переменных будет идентична исходной системе уравнений.
Граничные условия для скорости на верхней границе в случае безразмерных переменных будут иметь вид: U11 = ur1 (H) / Ur = cos (a), U21 = ur2 (H) / Ur = sin (a),
U3|Z=1 = 0, где a - угол, характеризующий направление ветра. Следовательно, безразмерные компоненты скорости ui , полученные при различных значениях скорости ветра Ur и при одинаковом направлении ветра a , будут одинаковы, а размерные значения компонент скорости при фиксированном a будут пропорциональны Ur: ui = ui (a)-Ur.
Граничные условия для безразмерной концентрации на стене здания будут определяться безразмерным потоком вещества F = (sqrt((Д • u) + (Д • u) ) + Дй). Поскольку
безразмерные компоненты скорости одинаковы при фиксированном направлении ветра на верхней границе области, то безразмерное поле концентрации, полученное в результате системы уравнений, будет зависеть только от направления ветра на верхней границе области.
Отмеченные свойства подобия решений системы (1)-(6) были подтверждены данными расчетов. На рис. 2 показаны поля скорости, вычисленные при a = 45°, Ur = 4,5 м/с и Ur = 1,0 м/с, а также поле безразмерной скорости, вычисленное при Ur = 4,5 м/с. Как
видно из данных, представленных на рисунке, размерное поле скорости на рис. 2б совпадает с безразмерным полем скорости на рис. 2в, что подтверждает подобие поля скорости по отношению к значению скорости на верхней границе области при заданном направлении ветра.
Эти результаты позволили ограничиться проведением вычислений при фиксированной скорости ветра на верхней границе области. Направление ветра на верхней границе
200-i
\\\\\\ \ WWW \ \ \ \ \ \
\ W w w w w w w w \ A\\\\\\\\w\
\ \ W W W ч» Ч Ч X 4 W \ \ 4 \ » -
www W w \\{± www \ \ W \ V
No.102
No.103
ww\ l www www www ; w\ww - / ww\wx^ \ \
WWWWV^A^' W W W W44»4» 4 w W W W W \ \ V
W W W W 4 4\ N W w"
N. N. \
-200 -150 -100
100 150 200
200-,
ч « •« •« ■» ■« •« »
■« ч ч ч ч ч
ЧЧ ЧЧЧЧЧЧЧЧЧЧЧЧ 4
150-
6).
ч 4 ч ч ч 4
ч ч 4 4 4 4 4. Ч Ч 4 4
100-
ч ч ч ч ч ч ч ч ч ч ч
4 4 4 4 4 4 4 4 4 4 ч
4 4 Ч Ч Ч 4 4 4 Ч 4 1« 4
ч ч ч ч ч ч
4 4 4 4 4 i,
ч ч v ч \ ,
ч< Ч \ \ V |
ii \ v v v i
ч ч ч ч ч s.
4
No.102
ч ч ч ч ч ч ч ч ч ч
No.103
4 4 4 4 4 Ч Ч Ч 4
Ч Ч Ч 4 4
4 4 4 4
« « Ч ч » 4 •« Ч Ч 4 «
ЧЧЧЧЧЧЧЧЧЧ
-2004
4 4 ч 4 Ч ч ч ч ч 4 ч ч ч 4 « 4 4 4 ЧЧ ЧЧЧЧЧЧЧЧЧ ч
4 4 « ч 4 4 4 4 4_«
-200 -150 -100 -50
0 X, M
100 150 200
150-
100-
50-
-100-
-150-
-200-
ч ч ч ч ч ч ч ч ч ч ч ч ч 4 ч ч ч ч ч ч
ч ч ч ч ч ч ч ч ч ч Ч Ч ч Л ч ч ч ч ч ч
ч в V ч ч ч ч ч ч ч Ч ч Ч * ч ч ч ч ч ч
ч ч ч ч ч ч ч ч ч ч "ч Ч Ч ч ч ч ч ч ч
ч ч ч ч ч ч ч ч Ч -ч Ч ^ -ч ч ч ч ч ч
ч ч ч ч ч * * ч ч ч
* I * * * * *
ч ч ч ч ч ч ч ч
ч ч « ч ч ч ч ч . . 1 No.102 t ? ч ч ч ч ч
ч ч « \ * t * J, -у v .
ч ч , No.103 - Ч - - ч
ч ч ч ч \ \ t ч * •
ч ч » * ч ч ч / 4 *
ч ч ч ч ч ч » ч \ ч ч •
ч ч ч ч ч ч ч ч ч ч Ч ч
ч ч ч ч ч ч ч ч ч Ч ч ч ч « ^ ■ч
ч ч ч ч ч ч ч ч ч ч ч ч ч ч ч ч * * •
ч ч ч ч ч ч ч ч ч ч ч ч ч * ч ч ч ч •ч
ч ч ч ч ч ч ч ч ч ч ч ч ч ч ч ч ч ч ч
-200 -150 -100
-50
0 х,м
50
100 150 200
Рис. 2. Рассчитанное поле скорости ветра вокруг зданий для следующих случаев: а - иг = 4,5 м/с, размерная скорость; б -
иг = 1,0 м/с, размерная скорость; в -
безразмерная скорость; показаны данные в каждом десятом узле расчетной сетки
области изменялось с шагом 45° от 0 до 315°. Полученные поля концентрации для различных направлений ветра осреднялись в зависимости от повторяемости направления ветра. Согласно данным Днепродзержинской метеостанции, различные направления ветра характеризуются следующими значениями повторяемости: северный ветер - 15%, северовосточный - 13%, восточный - 11%, юго-восточный - 14%, южный - 12%, юго-западный - 10%, западный - 10%, северозападный - 15%.
4. Результаты расчетов
В расчетах предполагалось, что все окна на северной стене здания открыты, и общая площадь окон Swindows равняется площади стены S ,, (a = S . d / S „ = 1). Значение a = 1 не
wall \ windows wall /
является реалистичным, поэтому результаты, представленные ниже, пересчитаны для другой относительной площади окон (a = 1/ 4) путем уменьшения концентраций, полученных при a = 1, в 4 раза.
На рис. 3 приведено поле приземной среднегодовой концентрации загрязнения. На рис. 4 приведено поле среднегодовой концентрации на стене здания №102. Из сопоставления данных на рис. 2 и 3 видно, что приземное поле концентрации ассиметрично, тогда как поле концентрации на стене здания №102 характеризуется большой степенью симметрии. В среднем приземная безразмерная концентрация загрязнения у стены здания № 102, создаваемая эмиссией из здания № 103, равна 0,2. Среднегодовая приземная безразмерная концентрация между зданиями № 102 и № 103 (исключая полосу шириной 3 м, прилегающую к зданию №103) составляет 0,25. Среднегодовая концентрация на расстоянии 25 м от восточной стены здания № 103 достигает значения 0,25, а на расстоянии 20 м от здания № 103 эта величина достигает значения 0,5. Такая разница в значениях концентрации по разные стороны здания № 103 объясняется наличием соседнего блокирующего здания с западной стороны, тогда как с восточной стороны здания № 103 другие строения отсутствуют.
X, м
Рис. 3. Изолинии среднегодовых значений безразмерной концентрации С / С3 на стене здания № 102
40
-60
-аи --1-1-г--1-г-т-1--т-т-г—
-100 -80 -60 -40 -20 0 20 40 60 80 100
Рис. 4. Изолинии среднегодовых приземных значений безразмерной концентрации С / С3; шаг изолиний равен 0,1; максимальное значение равно 1,9
Таким образом, приведенные на рисунках результаты позволяют сделать вывод о том, что в случае разрушения окон в здании №103 концентрации загрязнения в окрестности здания №102 и в других местах территории, прилегающей к зданию №103, заметно возрастут и могут представлять опасность для находящихся там людей.
Представим основные характеристики быстродействия модели в различных вариантах расчетов. Расчеты проводились на кластере Института проблем математических машин и систем НАН Украины. Кластер состоит из 4 узлов, на каждом из которых установлено два 4-ядерных процессора Intel(R) Xeon(R) CPU E5405, 2.00Ghz и 16 Гб оперативной памяти. В параллельном режиме расчетов на 8 ядрах один расчет, в котором осуществлялось численное решение системы уравнений (1)-(6), составлял в среднем 94151 с, что приблизительно составляет одни сутки. При этом время расчетов варьировалось в зависимости от направления ветра от 32155 с до 138398 с, то есть больше, чем в 4 раза. При расчете на одном ядре время расчетов составляло в среднем 1058070 с, что говорит о чрезвычайно высокой эффективности распараллеливания.
5. Выводы
В работе представлены результаты применения программного комплекса OpenFoam для оценки среднегодовых характеристик загрязнения в окрестности зданий бывшего ураново-
го производства «Приднепровский химический завод» (ПХЗ). Была использована упрощенная модель на основе стационарных уравнений несжимаемой жидкости, которая позволяет описывать поля ветра вокруг зданий при нейтральной стратификации. Использован также «конструктор моделей» OpenFoam, с помощью которого существующие модели были адаптированы для описания выветривания загрязнения из производственных зданий путем введения дополнительных членов в исходные уравнения Алгоритмы, реализованные в OpenFoam, позволяют решать исходные уравнения с применением методов параллельных вычислений на произвольных вычислительных сетках. В данной работе коэффициент ускорения за счет применения распараллеливания был практически равен количеству расчетных ядер.
Благодаря высокой эффективности используемых численных методов решения уравнений гидродинамики, в данной работе получены оценки среднегодовых характеристик загрязнения вблизи производственных зданий путем осреднения результатов расчетов при различных направлениях ветра с весовыми коэффициентами, определяемыми повторяемостью соответствующих направлений ветра. Показано, что, благодаря пропорциональной зависимости интенсивности источника от скорости ветра вблизи зданий, а также за счет подобия поля скорости по отношению к невозмущенной скорости ветра при фиксированном невозмущенном направлении ветра, концентрация вещества в области зависит только от направления невозмущенного ветра. В работе отмечено, что на территориях, прилегающих к производственному зданию № 103 ПХЗ, на расстояниях 10-20 м, за счет выветривания радиоактивности через разбитые окна могут достигаться высокие уровни приземной концентрации загрязнения - в среднем около одной четверти от концентрации внутри здания.
СПИСОК ЛИТЕРАТУРЫ
1. Ковалец И.В. Численная гидродинамическая модель атмосферной дисперсии загрязнений вокруг зданий / И.В. Ковалец // Зб. наук. праць 1нституту проблем моделювання в енергетищ iм. Г.С. Пухова. - 2011. - № 57. - С. 3 - 10.
2. Open Foam - The Open Source CFD ToolBox User Guide, Version 2.2.2 [Електронний ресурс] // OpenFoam Foundation. - 2013. - 28th September. - Режим доступу: http://openfoam.org.
3. Lavrova T. Radioecological assessment and remediation planning at the former uranium milling facilities at the Pridnieprovsky Chemical Plant in Ukraine / T. Lavrova, O. Voitsekhovych // Journal of Environmental Radioactivity. - 2013. - Vol. 115. - P. 118 - 123.
4. Численное моделирование воздушного распространения радона вокруг урановых хвостохрани-лищ / И.В. Ковалец, М.И. Железняк, А.В. Халченков [и др.] // Электронное моделирование. - 2010. - T. 32, № 3. - C. 67 - 82.
Стаття над1йшла до редакцп 18.03.2014
