CC BY
37
УДК 004.9:504:519.6 А.В. ХАЛЧЕНКОВ*
ПОДГОТОВКА ВЫЧИСЛИТЕЛЬНОЙ КОНЕЧНО-ОБЪЕМНОЙ СЕТКИ ДЛЯ ТРЕХМЕРНОГО CFD-МОДЕЛИРОВАНИЯ ПРИЗЕМНОГО СЛОЯ АТМОСФЕРЫ В ОБЛАСТИ ГОРОДСКОЙ ЗАСТРОЙКИ С ИСПОЛЬЗОВАНИЕМ ПРОГРАММНОГО КОМПЛЕКСА OPENFOAM
1нститут проблем математичних машин та систем НАН Украши, Кив, Украша
Анотаця. Дослгджена задача тдготовки горизонтально-неструктурованог розрахунковог стки за допомогою програми GMSH для тривим1рного г1дродинам1чного (CFD) моделювання приземного шару атмосфери в окол1 м^ськог забудови з використанням програмного комплексу OpenFoam та з урахуванням геопросторових даних, а також тформацп про забудову м1сцевост1. Розроблеш ал-горитми використання втьно доступних Г1С-даних програми Google Earth i глобальних набор1в категорт землекористування та топографа.
Ключовi слова: Г1С, розрахункова гiдродинамiка, неструктурована стка, OpenFoam, GMSH, то-пографiя, мiська забудова.
Аннотация. Исследована задача подготовки горизонтально-неструктурированной вычислительной сетки с помощью программы GMSH для трехмерного гидродинамического (CFD) моделирования приземного слоя атмосферы в окрестности городской застройки с использованием программного комплекса OpenFoam и информации о застройке местности. Разработаны алгоритмы использования свободно доступных ГИС-данных программы Google Earth и глобальных наборов категорий землепользования и топографии.
Ключевые слова: ГИС, вычислительная гидродинамика, неструктурированная сетка, OpenFoam, GMSH, топография, городская застройка.
Abstract. The problem of preparing a horizontally unstructured computational grid with the help of the GMSH program for three-dimensional hydrodynamic (CFD) modeling of the atmospheric surface layer in urban zone using the OpenFoam software package and information on the building area was investigated. The algorithms of using freely available GIS-data of the Google Earth program and the global data sets of land using and topography were developed.
Keywords: GIS, computational fluid dynamics (CFD), unstructured grid, OpenFoam, GMSH, topography, urban zone.
1. Введение
Программный комплекс OpenFoam [1], кроме ряда готовых гидродинамических моделей, предоставляет пользователю доступ к исходному коду, что позволяет легко и удобно комбинировать и модифицировать гидродинамические модели. Эта особенность позволяет создавать гидродинамические модели практически произвольной сложности и выгодно отличает OpenFoam от проприетарного программного обеспечения (ПО) с закрытым кодом. В программном комплексе OpenFoam присутствуют все необходимые инструменты для успешного моделирования приземного слоя атмосферы. Присутствующие в программном комплексе OpenFoam элементы конструктора моделей позволяют создавать новые модели турбулентности, что позволяет описывать в конечно-объемном виде систему гидродинамических уравнений с использованием моделей турбулентности, специализированных для атмосферных задач и правильно воспроизводящих характерные асимптотики приземного слоя атмосферы. Наряду с указанными преимуществами, программный комплекс OpenFoam обладает существенным недостатком: имеющиеся в наличии инструменты построения вычислительной сетки не предназначены для решения атмосферных задач и не позволяют создавать вычислительные сетки с учетом особенностей подстилающей по-
© Халченков А.В., 2017
ISSN 1028-9763. Математичш машини i системи, 2017, № 2
верхности (топография, шероховатость). Отдельную сложность представляет задача явного разрешения городской застройки, расположенной в условиях сложной топографии.
В комплекте поставки программного комплекса OpenFoam присутствует утилита blockMesh для построения параметрических блочно-структурированных сеток, предназначенная, в первую очередь, для решения задач с несложной геометрией. Например, она успешно использовалась в работе [2], где моделировалось распространение загрязнений в окрестностях зданий. Возможностей утилиты blockMesh было достаточно для этого исследования, поскольку не учитывалась топография (предполагалось, что в близлежащих окрестностях застройки ровная подстилающая поверхность). Среди учебных примеров, входящих в комплект поставки OpenFoam, присутствует пример моделирования потока ветра над поверхностью со сложной топографией. В этом примере для моделирования сложной поверхности используется утилита snappyHexMesh, позволяющая генерировать гексаэдральную сетку на основе импортируемой геометрии в формате (STL) или (OBJ). Формально, использование этой утилиты позволяет решить проблему одновременного учета наличия застройки и сложной топографии, но требует наличия файла с соответствующей геометрией, содержащей одновременно и топографию, и городскую застройку. При этом создание сложной координатно-привязанной геометрии с учетом топографии - это отдельная задача, сложно реализуемая без привлечения стороннего, как правило, проприетарного ПО.
Целью работы является создание алгоритма действий, позволяющего построить трехмерную, согласованную, координатно-привязанную вычислительную сетку с учетом особенностей подстилающей поверхности (топография, шероховатость) при наличии явно разрешаемой городской застройки с использованием исключительно свободно распространяемого ПО.
2. Постановка задачи. Состав входных, выходных и промежуточных данных
Процедура численного CFD-моделирования состоит из следующих четырех основных этапов:
• построение геометрии - крупномасштабное трехмерное моделирование всей области расчета с назначением характерных свойств отдельным геометрическим объектам: точкам, линиям, граням, объемам;
• построение расчетной сетки - процедура триангуляции в соответствии с геометрией с учетом свойств геометрических объектов. То есть, заполнение расчетной области конечными объемами с учетом областей сгущения и назначение типов граничных условий для граней, ограничивающих расчетную область;
• непосредственно решение системы нестационарных гидродинамических уравнений;
• обработка полученных результатов и визуализация.
Данная работа ограничивается первыми двумя пунктами из приведенных основных этапов. Оба эти этапа частично можно выполнить с использованием программы GMSH [3], предназначенной для создания трехмерных вычислительных сеток. Преимуществом этой программы, распространяемой по GNU (GPL) лицензии, является наличие собственных инструментов для построения крупномасштабной геометрии, возможность назначения граничных условий для геометрических граней и наличие собственных инструментов для построения сетки конечных элементов в соответствии с геометрией. Дополнительным преимуществом использования этой программы является наличие в стандартном комплекте поставки OpenFoam утилиты gmshToFoam, которая позволяет преобразовать выходной формат программы GMSH непосредственно в сформированный набор данных для Open-Foam.
Наряду с указанными преимуществами, программа GMSH не предназначена для моделирования атмосферных процессов, в ней отсутствуют собственные инструменты, позволяющие выполнить привязку расчетной сетки к координатам и учесть особенности подстилающей поверхности. Возникает потребность в решении трех дополнительных подзадач:
• привязка трехмерной вычислительной сетки к географическим координатам;
• учет топографии;
• определение параметра шероховатости.
В работе [4] была исследована задача построения качественной координатно привязанной двумерной сетки на основе карт Google Earth [5], в частности, там был изложен метод экспорта контуров, построенных в программе Google Earth, в программу GMSH для последующего построения двумерной сетки.
С учетом вышесказанного, в данной работе входными данными для моделирования будем считать созданные в программе Google Earth наборы полигонов, определяющие контуры отдельных зданий, значения высот для каждого здания, а также свободно доступные глобальные наборы топографических данных и данных о категориях землепользования.
В качестве промежуточных данных будут использоваться внутренние форматы программы GMSH: формат GEO для построения и модификации геометрии и формат MSH, содержащий список граничных условий, список узлов и список элементов (поверхностных и объемных). Этот файл содержит данные в обычном текстовом виде, каждый блок информации (список) начинается и оканчивается соответствующим ключевым словом. Например, блок узлов начинается с ключевого слова $Nodes, далее следует строка с общим количеством узлов, далее каждая строка содержит номер соответствующего узла и его координаты, оканчивается блок узлов ключевым словом $EndNodes. Блок элементов ограничен ключевыми словами $Elements и $EndElements. Каждая строка этого блока содержит номер элемента, тип элемента (треугольный, квадратный, 6-узлов.призма и пр.), список характеристик, описывающих принадлежность к определенным граничным условиям или геометрическим объектам, список узлов вершин элемента. Простая и известная структура файла позволяет легко извлекать и модифицировать содержимое этого файла.
Выходные данные будут состоять из папки polyMesh, получаемой путем конвертации MSH файла, в которой хранится ряд файлов (points, faces, owner, neighbour, boundary), позволяющих полностью описать структуру расчетной сетки OpenFoam, готовой для выполнения гидромеханического моделирования.
3. Создание трехмерной координатно привязанной сетки
В качестве примера в этой работе используется территория Приднепровского химического завода, контуры основных объектов которого (здания, хвостохранилища) приведены на рис. 1. На основе этих контуров в соответствии с [4] построим двумерную сетку, состоящую из треугольных элементов (рис. 2). Использование треугольных элементов значительно упрощает последующий процесс учета топографии (при использовании четырехугольных элементов необходимо обеспечивать нахождение всех узлов элемента в одной плоскости, что весьма не просто при интерполяции на сложную поверхность, такую, как произвольная топография). Путем подбора параметров и редактирования файла с геометрией GEO необходимо достичь создания приемлемой двумерной сетки, на основе которой в дальнейшем будет строиться трехмерная вычислительная сетка.
Трехмерная сетка на основе созданной двумерной сетки строится путем вытягивания поверхностей (экструзии). Для этого в файл GEO необходимо добавить команду:
Extrude {0,0,H} {Surface{S1, S2,.., Sn}; Layers{sl};Recombine;},
где Н - расстояние, на которое вытягиваются поверхности (в нашем случае по оси Z); ^1, S2,.., Sn} - список, состоящий из номеров поверхностей, которые участвуют в вытягивании; sl - количество параллельных слоев, которые получатся при последующем генерировании сетки.
Рис. 1. Контуры, привязывающие географические объекты Приднепровского химического завода к координатам
Рис. 2. Результат построения неструктурированной двумерной сетки из треугольных элементов (вся территория промплощадки)
За начальный уровень удобно принять уровень самой высокой постройки, далее воспользоваться командой Extrude и вытянуть все поверхности вверх до необходимого уровня, то есть заполнить геометрией весь объем выше зданий. Далее, для заполнения пространства между зданиями, нужно последовательно выполнять команду Extrude, опускаясь от уровня самого высокого здания, до уровня здания пониже и так далее до самой земли, выбрасывая из списка экструзии поверхности, соответствующие зданиям. Преимуществом процедуры экструдирования является возможность строить трехмерную сетку, заведомо состоящую из указанного количества параллельных слоев, что удобно при моделировании и дальнейшей модификации вычислительной области.
В результате выполненной процедуры вытягивания у нас образуется ряд геометрических объемов, которые необходимо преобразовать в один внутренний физический объем. Выполняется это с помощью команды
Physical Volume("intemal")={V1, V2,.., Vm},
где {V1, V2,.., Vm} - список геометрических объемов, образовавшихся при вытягивании поверхностей, «internal» - название физического объема (необходимо задать, но в дальнейшем это название использоваться не будет).
Далее необходимо создать и назвать ряд физических поверхностей, которые будут соответствовать различным граничным условиям последующего моделирования. То есть, необходимо сгруппировать и назвать все внешние поверхности, которые ограничивают физический объем, по типам будущего граничного условия. Делается это с помощью ряда команд:
Physical Surface("bottom") = { Sb1, Sb2,...};
Physical Surface("inlet") = { Si1, Si2,,...};
Physical Surface("top")={ St1, St2,...}
и т.д.,
где «bottom», «inlet», «top» - примеры названий будущих граничных условий; Sb1, Sb2 -номера поверхностей, ограничивающих физический объем со стороны соответствующего граничного условия (для последующего удобства задания топографии удобно крышу каждого здания назвать индивидуально, что позволит легко идентифицировать узлы и элементы, из которых она состоит).
Далее необходимо открыть отредактированный файл GEO в программе GMSH, перейти во вкладку «Mesh», выполнить команду «3D» и получить трехмерную конечно-объемную сетку, состоящую из слоев прямых призм с треугольным основанием (рис. 3). Для дальнейшей работы полученная трехмерная сетка сохраняется в формате MSH.
а б
Рис. 3. Трехмерная вычислительная сетка: общий вид (а), детализированная окрестность хвостохранилищ с близлежащими зданиями (б)
4. Учет топографии
Чтобы учесть топографию, в первом приближении необходимо извлечь из MSH-файла блок узлов, описанный в разделе 2, и выполнить для каждого узла следующую процедуру:
• взять координаты X и ^ некоего ьго узла;
• определить уровень поверхности грунта над уровнем моря Ы в точке с соответствующими координатами X и путем интерполяции значений из глобального набора топографических данных (например, SRTM 30-sec [6]);
• прибавить уровень поверхности грунта к ^ координате узла и обновить координату Zi в файл MSH.
Указанный подход очень прост в реализации, но есть недостаток: геометрия верхних поверхностей зданий (крыш) и верхней границы расчетной области перестают быть плоскими (рис. 4 а) и повторяют геометрию грунта (рис. 4 б).
Во втором приближении можно выровнять верхнюю границу расчетной области. Для этого нужно выполнить описанную выше процедуру для уровней, не превышающих высоту максимально высокого здания, а далее прибавлять не /?,, а параметр Аг, определяемый с помощью выражения
Аг = Л,.(Я-71.)/Я, (1)
где Аг - необходимое изменение уровня Zi, соответствующее / -му узлу, Н - высота вычислительной области (без учета топографии). Такой подход также прост в реализации, но пригоден только в случае несложной топографии (рис. 4 в), поскольку крыши зданий повторяют рельеф подстилающей поверхности.
В случае сложной топографии необходимо выполнять следующую процедуру «выравнивания» крыш:
• считывается информация об элементах, описанная в разделе 2;
• поскольку в процессе задания граничных условий мы для каждой крыши здания создали свою физическую группу (назвали индивидуально), легко выбираются элементы, которые образуют крыши;
• имея набор элементов, которые образуют крышу, определяются узлы, входящие в каждую из крыш;
• для каждой крыши усредняется координата Z, что обеспечивает выравнивание каждой крыши (то есть крыши становятся плоскими и горизонтальными).
К сожалению, в случае больших перепадов высот грунта на горизонтальных расстояниях масштаба длины зданий подобная модификация сетки приводит к грубым артефактам в узлах, близких к пересечению крыш и стен зданий (например, окрестность левого верхнего угла здания на рис. 4 г).
Чтобы исправить это, нужно после предыдущей процедуры «выравнивания крыш» выполнить следующую дополнительную процедуру «сглаживания» сетки (рис. 5):
• выбираются все узлы слоя, в котором расположена самая низкая крыша;
• определяются узлы, расположенные ближе некоторого подобранного расстояния к крыше;
• вычисляется новая высота этих узлов путем интерполирования между узлами крыши и остальными узлами соответствующего слоя. В результате получаем гладкий слой узлов с ровными крышами;
Рис. 4. Стадии построения сетки с учетом криволинейной топографии в окрестности прямоугольного здания
• при необходимости, имея высоты узлов на двух слоях (грунт и уровень самых низких крыш), высота промежуточных узлов вычисляется путем деления высот на равные отрезки в зависимости от количества промежуточных слоев.
Далее вместо уровня грунта подставляется «сглаженный» слой (нижайших крыш), вместо слоя нижайших крыш берется слой соответствующей следующей по высоте крыши, и выполняется аналогичная процедура сглаживания. Указанная процедура «сглаживания» итерационно повторяется вплоть до слоя, соответствующего высоте самых высоких зданий. Высота узлов, расположенных в слоях выше крыш, определяется по формуле (1).
В результате описанной процедуры получается гладкая вычислительная сетка (рис. 4 д) с ровной и горизонтальной верхней границей, с ровными и горизонтальными слоями ячеек, примыкающих к крышам, и с нижней границей, имеющей сложную форму, воспроизводящую топографию подстилающей поверхности.
Экспорт вычислительной сетки в OpenFoam выполняется с помощью утилиты gmshToFoam, входящей в стандартный комплект поставки OpenFoam.
5. Определение параметра шероховатости
В среде OpenFoam реализован ряд пристеночных функций, используя которые можно получать правильные пристеночные профили скорости ветра без непосредственного разрешения турбулентного пристеночного слоя мелкой расчетной сеткой. Наиболее естественной пристеночной функцией при решении атмосферных задач для турбулентной вязкости есть функция nutkAtmRoughWallFunction, которая обеспечивает корректный логарифмический слой скорости и зависит от заданного значения параметра шероховатости.
Пространственное распределение параметра шероховатости можно получить с помощью глобальных наборов данных о категориях землепользования (Modis, GlobeLand30, и пр.), сопоставив каждую приведенную в них категорию землепользования с некой усредненной шероховатостью, характерной для этой категории (например, препроцессор MAKEGEO модели атмосферного переноса Calpuff [7]). Поскольку используемая в OpenFoam неструктурированная сетка состоит из элементов, существенно отличающихся размером, параметр шероховатости каждого элемента определяется таким образом:
• если в элемент попадают узлы регулярной сетки глобального набора данных, то в качестве шероховатости для элемента применяется среднее значение шероховатости в этих узлах (рис. 6 а);
• если не попадают, то в качестве шероховатости для элемента применяется результат интерполяции значений в ближайших четырех узлах глобального набора в центр масс элемента (рис. 6 б).
Рис. 5. Высота красных узлов, образующих крышу, определяется путем усреднения (выравнивается крыша), высота узлов с фиолетовой заливкой определяется путем интерполяции между красными и зелеными узлами
Рис. 6. Интерполяция значений параметра шероховатости из регулярной сетки глобального набора на неструктурированную вычислительную сетку (звездочками отмечены узлы, которые участвуют в определении шероховатости для выделенного элемента)
Чтобы подставить значения шероховатости в модель OpenFoam, необходимо найти координаты центров ячеек нижней грани. Проще всего это выполнить с помощью утилиты writeCellCentres (входящей в стандартный комплект поставки OpenFoam). Из полученных в результате ее работы файлов необходимо выбрать значения в блоке, соответствующие граничному условию для нижней границы. Далее область расчета экспертно делится на подобласти следующих категорий.
• Категория I - явно моделируемые области застройки, в которых обтекания зданий напрямую воспроизводятся на расчетной сетке и внутри которой размеры ячеек сравнимы с масштабом шероховатости из глобального набора (около z0~1 м) - для нее задается минимальное значение шероховатости.
• Категория II - области, достаточно отдаленные от областей 1-й категории (степень отдаленности определяется высотой застройки), в которых обтекание элементов шероховатости является подсеточным процессом (то есть размер сетки много больше параметра шероховатости z0); значение z0 задается путем интерполяции значения из глобального набора данных.
• Категория III - промежуточные области, для которых значения шероховатости можно задать путем интерполирования значений между ближайшими областями категорий I и II.
Определение границ области застройки и промежуточной области выполняется с использованием тех же инструментов, что использовались при определении границ крыш и области близлежащих узлов при учете топографии.
6. Демонстрационный гидродинамический расчет, выполненный на неструктурированной трехмерной вычислительной сетке с учетом топографии и категорий землепользования
С помощью приведенного алгоритма была создана трехмерная вычислительная сетка, приведенная на рис. 7. Моделировался один стационарный метеорологический сценарий, соответствующий нейтральной стратификации и невозмущенному логарифмическому слою на входной границе расчетной области. Для твердых поверхностей, ограничивающих расчетную область, использовались пристеночные функции. Для выходной границы расчетной области использовалось граничное условие, соответствующее нулевому градиенту для всех переменных. В качестве источника загрязнений использовался равномерный единичный поток загрязнения (радона) с поверхности трех хвостохранилищ [8].
Результаты выполненного расчета приведены на рис. 8.
а б
Рис. 7. Трехмерная вычислительная сетка с топографией и явным разрешением городской застройки: общий вид (а), окрестность застройки (б)
а б
Рис. 8. Результаты расчета: поле скоростей на высоте 2 метра от поверхности земли (а),
поле приземной концентрации (б)
7. Выводы
В данной статье описаны основные этапы подготовки трехмерной неструктурированной вычислительной сетки для последующего решения задач вычислительной гидродинамики приземного слоя атмосферы с использованием программного комплекса OpenFoam.
Для построения и координатной привязки вычислительной сетки предложено совместное использование программы GMSH и программы Google Earth. Разработан алгоритм, позволяющий привязать явно разрешенную городскую застройку к топографии подстилающей поверхности, сохраняя при этом форму и ортогональность верхних и боковых поверхностей зданий, но привязав при этом нижнюю кромку зданий к криволинейной топографии. Приведена пошаговая инструкция, позволяющая устранить возможные сеточные артефакты, неизбежно возникающие при совмещении геометрических объектов различной кривизны.
В качестве информации о подстилающей поверхности предложено использование свободно доступных глобальных наборов данных (SRTM 30 для топографии, Modis или GlobeLand 30 для категорий землепользования). Предложена методика практического назначения параметра шероховатости для трех подобластей области расчета в зависимости от соотношения между горизонтальным размером ячеек сетки и значением параметра шероховатости из глобального набора. Описаны основные сложности и пути их решения,
возникающие при построении сетки и при аппроксимации данных о подстилающей поверхности.
БЛАГОДАРНОСТИ
Данная работа выполнена при поддержке гранта Президента Украины докторам наук № Ф74.
СПИСОК ЛИТЕРАТУРЫ
1. OpenFOAM Programmer's Guide. Version 3.0.1. OpenFOAMFoundation, 2015 [Электронный ресурс]. - Режим доступа: http://foam.sourceforge.net/docs/Guides-a4/ProgrammersGuide.pdf.
2. Халченков А.В. Использование вычислительной гидродинамической модели OpenFoam для оценки влияния загрязненных зданий на прилегающие территории / А.В. Халченков, И.В. Ковалец // Математичш машини i системи. - 2014. - № 2. - С. 97 - 104.
3. Geuzaine C. GMSH a three-dimensional finite element mesh generator with built in pre and post processing facilities / C. Geuzaine, J.F. Remacle // International Journal for Numerical Methods in Engineering. - 2009. - N 79 (71). - P. 1309 - 1331.
4. Халченков А.В. Создание расчетной сетки для экологических моделей на основе свободно доступных ГИС-данных / А.В. Халченков // Математичш машини i системи. - 2016. - № 2. - С. 109 -116.
5. Справочный центр - Google Планета Земля [Электронный ресурс]. - Режим доступа: http://support.google.com/earth/?hl=ru#topic=4363013.
6. An assessment of the SRTM topographic products, Technical Report JPL D-31639 / E. Rodriguez, C.S. Morris, J.E. Belz [et al.] // Jet Propulsion Laboratory. - Pasadena, California, 2005. - 143 p.
7. Scire J.S. A user's guide for the CALPUFF dispersion model (Version 5) / Scire J.S., Strimaitis D.G., Yamartino R.J. - Earth Tech. Inc., Concord. - 2000. - 521 p. [Электронный ресурс]. - Режим доступа: http: //www .src.com/calpuff/calpuff1.htm.
8. Atmospheric dispersion of radon around the uranium mill tailings of the former Pridneprovsky Chemical Plant in Ukraine / I. Kovalets, C. Asker, A. Khalchenkov [et al.] // Journal of Environmental Radioactivity - 2017 - doi: 10.1016/j.jenvrad.2017.03.025.
Стаття над1йшла до редакцп 31.03.2017
