Введение
Авторам данной статьи довелось выполнить довольно редкий по своему характеру проект. Требовалось разработать коммерческую программу моделирования процессов движения воздушной среды в чистом помещении. Чистое помещение — это производственное помещение, отвечающее определенным требованиям по чистоте воздуха, температуре и скорости его движения. Основной показатель чистоты — это класс чистоты, который определяется ГОСТом по концентрации частиц в воздухе. Потоки воздуха в чистом помещении направляются так, чтобы обеспечить эффективное удаление пыли и аэрозолей из помещения. Требования могут также ограничивать градиенты температуры в пространстве и во времени. Программа Cleanroom предназначена для использования в качестве инструмента проектировщика чистых помещений. С ее помощью проектировщик должен выполнять размещение оборудования и элементов вентиляции, а по результатам моделирования процессов в воздушной среде определять степень соответствия варианта размещения предъявляемым требованиям по чистоте.
Моделирование концентрации частиц и вычисление класса чистоты здесь мы рассматривать не будем. В разработанной программе класс чистоты вычисляется для помещения в целом в соответствии с методом Густавссона, изложенным в книге www.asincom-group.ru/book.htm. В его реализации не было принципиальных затруднений. Очень полезной была бы реализация моделирования перемещения пыли в пространстве помещения, но эта задача не уложилась в бюджет проекта.
Движение воздуха описывается системой уравнений Навье-Стокса, которой посвящено большое количество литературных источников, и есть простое изложение основ даже здесь habrahabr.ru/post/171327/. Уравнения связывают между собой скорость движения воздуха, его температуру и давление. В чистом помещении параметры приточной вентиляции, вытяжной вентиляции и нагрева воздуха оборудованием составляют граничные условия для протекающих аэродинамических процессов.
Уравнения Навье-Стокса представляют собой систему дифференциальных уравнений в частных производных, которая, за исключением нескольких простых случаев, не может быть решена аналитически. Поэтому, как правило, эта система решается численно методом конечных элементов. Для решения подобных задач существуют универсальные приложения, относящиеся к классу CFD (Computational Fluid Dynamics). Это чрезвычайно сложные и дорогие программные продукты (Ansys, SolidWorks, Flow Vision и др.), которые позволяют исследовать течение жидкостей и газов в разнообразных условиях. Универсальность этих приложений приводит к тому, что пользователь должен обладать знаниями физики и математики на очень высоком уровне, а освоение этих программ требует многих месяцев обучения и практической работы.
В нашем случае требовалось разработать специализированную программу, которая обеспечивает выполнение следующих функций:
- Ввод геометрических размеров помещения и технологических компонентов. Ввод температуры всех поверхностей и отдельных регионов поверхности.
- Ввод геометрических размеров регионов на поверхности, соответствующих элементам приточной и вытяжной вентиляции, задание параметров движения воздуха через границу.
- Вычисление векторного поля скорости воздуха и скалярного поля температуры. Графическое отображение этих полей.
В данной программе система уравнений Навье-Стокса использовалась в приближении Буссинеска, т.е. с упрощением зависимости плотности воздуха от его температуры. Технологические элементы размещения и помещение представляются либо параллелепипедами, либо комбинациями параллелепипедов. На грани элементов размещения накладываются прямоугольные регионы со своими параметрами. Грани элементов размещения и наложенные на них регионы задают граничные условия: температуру и скорость движения воздуха. Решение системы уравнений методом конечных элементов дает искомый результат: скалярное поле температуры и векторное поле скорости воздуха внутри помещения. Таким образом, разрабатываемая программа является специализированным решателем одной и той же системы уравнений Навье-Стокса, для которой изменяются только граничные условия при различных конфигурациях помещения и размещенного оборудования. На рис. 1 показано ее окно с результатами моделирования. Графически моделируемое помещение и поля скорости и температуры отображаются в виде сечений в трех плоскостях.
Рис. 1. Окно приложения с результатами моделирования
Данное приложение, не смотря на его очевидную наукоемкость, было разработано за относительно небольшой срок. Это оказалось возможным благодаря использованию сторонних библиотек с открытым кодом Freefem++ и NetGen. Цель данной статьи видится в том, что наш опыт позволит другим разработчикам обойти некоторые «подводные камни» в разработках подобных приложений с привлечением указанных библиотек.
Реализация метода конечных элементов
Поиск подходящей библиотеки для решения уравнений Навье-Стокса выполнялся как среди специализированных программных средств для решения задач гидродинамики, так и среди универсальных средств решения систем дифференциальных уравнений в частных производных. При этом пришлось учитывать следующую особенность выполняемого проекта. Разрабатываемая программа Cleanroom — это коммерческий продукт, который должен функционировать под управлением Windows. Возможность приобретения платных библиотек была исключена техническим заданием. Поэтому сторонние библиотеки должны были быть бесплатными и иметь лицензию, допускающую их использование в коммерческом продукте.
Последнее ограничение исключало возможность использования сторонних программных продуктов с лицензией GNU GPL. Поэтому пришлось отказаться, например, от известной и широко используемой библиотеки OpenFoam (официальный сайт www.openfoam.org/ и имеется кратное описание здесь habrahabr.ru/company/intel/blog/170675/ ). Другие библиотеки также либо имеют не подходящую лицензию, либо проекты заморожены и их использование связано с большим риском. Единственной библиотекой, которая подходила по всем параметрам, оказалась Freefem++.
Freefem++ (официальный сайт www.Freefem.org/ff++/) — это программа численного решения дифференциальных уравнений в частных производных. Она основывается на методе конечных элементов. Это приложение, которое принимает на входе файл с текстом программы (скриптом), описывающей решение задачи на специальном C-подобном языке. Входной язык Freefem++ включает в себя основные алгоритмические конструкции: объявление переменных, операторы цикла, оператор ветвления, операторы ввода-вывода в стиле потокового ввода-вывода C++, объявление функций, элементарные типы данных, массивы, файлы и т.п. Для решения систем дифференциальных уравнений в языке предусмотрены специальные типы данных и функции. Такими типами являются триангуляционные сетки, пространство конечных элементов на базе этих сеток. Имеются специальные способы записи производных от функций и интегралов. Предусмотрена специальная конструкция solve, с помощью которой записывается вариационная формулировка задачи и решается система уравнений.
Для решения системы уравнений Навье-Стокса задаются две вариационные формулировки, которые используются последовательно. При помощи первой решается система уравнений для скорости воздуха и давления при заданных значениях температуры для заданного момента времени. Вторая вариационная формулировка позволяет пересчитать поле температуры по вычисленному полю скорости. Этот процесс циклически повторяется для последующих дискретов времени. В состав первой вариационной формулировки включаются граничные условия для скорости воздуха. Для непроницаемых стенок скорость задается равной нулю (эффект прилипания), а для регионов приточной и вытяжной вентиляции указывается направление и скорость движения воздуха. Вторая вариационная формулировка включает граничные условия, заданные в виде значений температуры стенок.
Геометрические элементы, которым надо приписать свойства, представляются пользователю как грани пространственных фигур и регионы поверхности на них (см. рис. 2). Поэтому для Freefem++ необходимо сформировать сетку конечных элементов так, чтобы треугольники на границе сетки полностью укладывались в регионы поверхностей, как показано на рис. 3. В типе данных mesh3 для трехмерной сетки конечных элементов для граничных треугольников проставляются метки (идентификаторы). Эти метки используются в вариационной формулировке для задания соответствующей функции на границе.
Рис. 2. Регионы поверхности, на которых должны быть заданы граничные условия
Рис. 3. Граничные треугольники пространства конечных элементов
Генерация сетки конечных элементов
К сожалению, Freefem++ не имеет собственного генератора 3D сеток. Для этой цели авторами Freefem++ предлагается использовать сторонний модуль TetGen (http://tetgen.berlios.de). Модуль TetGen подключается в скрипте стандартным для Freefem++ способом (оператор load «tetgen»). В языке имеется несколько предопределенных функций генерации 3D сеток, которые фактический реализуются именно этим модулем. Однако модуль имеет своеобразную лицензию, которая требует согласования с автором возможности использования TetGen в коммерческих приложениях. Поэтому от этого генератора 3D-сеток пришлось отказаться и искать ему замену.
В языке Freefem++ имеется оператор импорта 3D-сетки из текстового файла формата mesh. Поэтому может быть использован практически любой сторонний генератор. Однако в нашей разработке сторонний генератор сеток конечных элементов должен быть бесплатным продуктом с лицензией, позволяющей использовать его в коммерческом приложении. В конечном счете был выбран NetGen (http://sourceforge.net/projects/netgen-mesher/).
NetGen предназначен для генерации и отображения 3D-сеток и состоит из нескольких библиотек. Основная библиотека nglib реализует непосредственно генерацию сеток. Имеется два способа задания пространства, которое должно быть заполнено тетраэдальными элементами:
- Задание пространства при помощи операторов конструктивной блочной геометрии (Constructive Solid Geometry, CSG). В этой технологии искомое пространство описывается совокупностью примитивных объектов, над которыми задаются такие операции как объединение, пересечение, разность.
- Описание поверхности, являющейся границей искомого пространства, в файле формата STL. В этом формате информация о поверхности представляется в виде координат треугольных граней поверхности и их нормалей.
Не смотря на привлекательность использования конструктивной блочной геометрии, от нее пришлось отказаться. Согласно ТЗ, поверхность помещения и любого элемента размещения может содержать регионы с собственными параметрами моделирования (температурой и скоростью воздуха). В реализации конструктивной блочной геометрии имеется возможность задать только свойства поверхности примитивной пространственной фигуры в целом. А найти разумный способ задания свойств отдельных регионов не удалось.
Для решения поставленной задачи оказался пригодным только вариант с формированием конечных элементов по геометрическому описанию в формате STL. Для того, чтобы сформировать этот файл, необходимо предварительно выполнить триангуляцию поверхностей, составляющих границу пространства. Причем триангуляция должна быть выполнена с учетом наложенных на поверхность регионов.
Генерация описания поверхности в формате STL
Построение триангуляции поверхностей было решено выполнять при помощи Freefem++, которая, в отличие от реализации генерации 3D-сеток, не использует для этих целей дополнительные сторонние модули. Так как в разрабатываемом приложении все элементы размещения и помещение представляются параллелепипедами, то для получения триангуляции их поверхностей надо было выполнить две операции:
- Выполнить 2D триангуляцию каждой грани с учетом регионов, наложенных на грань.
- Вершинам полученных треугольников приписать пространственные 3D координаты.
Исходные данные, необходимые для выполнения 2D триангуляции во Freefem++, представляются объектами специального типа — border. В таком объекте указывается параметрическая запись одного фрагмента границы. Для формирования триангуляционной сети предназначена функция buildmesh, которой на вход передаются все границы разделяемой поверхности. При этом для каждого фрагмента границы указывается количество отрезков треугольников, которые будут опираться на границу. В данном случае это количество должно быть указано равным единице, т.е. каждый фрагмент границы целиком будет одним ребром треугольника из триангуляции.
Исходные данные, вводимые пользователем, описывают накладываемые регионы поверхностей (рис. 4-а). В результате такого наложения на одной грани элемента размещения может образоваться множество границ произвольной сложности. Для выполнения процедуры триангуляции необходимо преобразовать геометрические данные о регионах в систему направленных отрезков как показано на рис. 4-б. Соответствующий алгоритм реализован авторами в Cleanroom. Полученная система границ описывается в скрипте, который передается на вход Freefem++. Данный скрипт выводит результаты триангуляции в файл (пример триангуляции на рис. 4-в). После этого разработанная программа, считав триангуляцию из файла, каждому треугольнику приписывает идентификатор, соответствующий набору физических свойств того региона поверхности, которому принадлежит треугольник. Очевидно, что триангуляция всех граней всех элементов размещения может быть выполнена одним скриптом за один вызов Freefem++.
Рис. 4. Преобразование регионов поверхности в систему границ и триангуляцию поверхности
Общий алгоритм использования Freefem++ и NetGen
Сведем воедино последовательность обращения к Freefem++ и NetGen в процедуре моделирования. Напомним, что исходные данные для моделирования — объектная модель элементов размещения с физическими свойствами регионов поверхностей. Выполняются следующие действия:
- Формирование системы границ регионов в виде направленных отрезков на каждой грани всех элементов размещения (см. рис. 4).
- Формирование файла со скриптом для Freefem++ с описанием границ поверхностей, подлежащих триангуляции.
- Вызов Freefem++. На его входе — подготовленный файл скрипта. Freefem++ выполняет триангуляцию заданных поверхностей. На выходе — файл 2D-триангуляции для граней элементов размещения.
- Преобразование 2D-треугольников в 3D-треугольники. Формирование файла в формате STL, в котором указываются координаты вершин треугольников и нормали к ним.
- Выполнение последовательности вызовов функций NetGen. На входе — STL-файл. NetGen генерирует 3D-сетку конечных элементов. На выходе — mesh-файл с описанием тетраэдров сетки и граничных треугольников.
- Подготовка файла со скриптом на входном языке Freefem+ для решения системы дифференциальных уравнений. В скрипте импортируется сетка конечных элементов из mesh-файла. В скрипт записываются граничные условия для меток граничных треугольников.
- Запуск Freefem++. На входе — mesh-файл и файл со скриптом. Freefem++ решает систему уравнений Навье-Стокса. На выходе — файл со значениями компонент скорости воздуха в узлах сетки и файл со значениями температуры в узлах сетки.
- Чтение mesh-файла с сеткой конечных элементов, файла со значениями скорости и файла со значениями температуры. Отображение искомых полей в графическом виде.
Важная особенность использования NetGen для генерации сетки конечных элементов заключается в следующем. В NetGen предусмотрена возможность указания набора предопределенных отрезков (ребер — в терминологии NetGen) на поверхности пространства, которые сохраняются при формировании сетки конечных элементов. В нашем случае обязательно надо задать в качестве ребер границы прямоугольных регионов, заданных пользователем. Тогда NetGen генерирует сетку конечных элементов так, что граничные треугольники плотно заполняют регионы поверхности, т.е. каждый граничный треугольник сгенерированной сетки полностью попадает в один прямоугольный регион. Благодаря этому, для граничного треугольника определены заданные физические свойства (скорость воздуха и температура), которые были указаны для региона. Если же ребра не задавать, то NetGen исходные треугольники из STL-описания не сохранит и полученные в результате граничные треугольники сети конечных элементов будут произвольным образом покрывать поверхность пространства, пересекая границы регионов с заданными свойствами. А это означает, что треугольнику нельзя будет приписать физические свойства для граничного условия.
При реализации передачи сетки конечных элементов, сгенерированной при помощи NetGen, на вход Freefem++ выявилась недокументированная особенность Freefem++, которая заключается в следующем. В mesh-файле с сеткой конечных элементов, который импортируется во Freefem++, тетраэдальные элементы описываются в виде строк с номерами вершин. Порядок следования вершин в документации на Freefem++ не оговаривается. Однако оказалось, что если вершины тетраэдра перечислены в том порядке, в котором их выдает NetGen, Freefem++ завершает работу из-за исключительной ситуации. Анализ исходного кода Freefem++ показал, что вершины каждого тетраэдра должны быть заданы в таком порядке, чтобы в результате вычисления объема тетраэдра во Freefem++ не получалось отрицательное число. Поэтому вершины тетраэдров из сетки, получаемой с помощью NetGen, необходимо переставлять соответствующим образом. К сожалению, выяснение этого факта отняло довольно много времени при разработке.
При решении тестовых задач выяснилось, что различные методы решения системы алгебраических уравнений, в которую преобразуется система дифференциальных уравнений в методе конечных элементов, имеют недокументированные ограничения на максимальное количество конечных элементов. Только метод, обозначаемый во входном языке Freefem++ как CG (Conjugate Gradient), позволил решать реальные задачи без возникновения исключительных ситуаций.
Место Freefem++ и NetGen в структуре приложения
В Cleanroom процедуры формирования конечных элементов и решения системы уравнений являются наиболее ресурсоемкими частями, и именно к ним предъявляются требования по быстродействию. Требования по быстродействию ко всем остальным функциям несопоставимо ниже. Поэтому для сокращения времени разработки было принято решение разрабатывать Cleanroom на C#. Для бесшовного использования NetGen, т.е. включения библиотеки на C++ в программу на C#, была разработана тонкая прослойка на C++/CLI. В ней осуществляется вызов функций из NetGen и преобразования, необходимые для адаптации данных к слою бизнес-логики. При этом в сборке на C++/CLI используется слой данных, разработанный на C#. Такое соединение в одном месте двух разных средств разработки (C++ и C#) оказалось чрезвычайно удобным. Так как Freefem++ является загрузочным модулем, то он вызывается как отдельное приложение, и связь с ним осуществляется через файлы. Разработанная структура программы показана на рис. 5.
Рис. 5. Структура приложения Cleanroom
Вспомогательные инструментальные средства
Для отладки приложения использовались различные инструменты, которые позволяют визуально отображать STL-файлы, mesh-файлы и получаемые поля процессов. Для просмотра поверхности, описываемой в формате STL, использовался GLC-player (http://www.glc-player.net/). Он достаточно устойчив к ошибкам в STL-файле и иногда позволяет по дефектам изображения поверхности эти ошибки выявить.
В состав проекта Freefem++ входит программа ffmedit, которая является адаптацией программы Medit (http://www.ann.jussieu.fr/frey/software.html) для Freefem++. Эта программа визуализирует mesh-файл. С ее помощью можно детально рассмотреть пространство конечных элементов. В ней можно даже посмотреть поверхности равных значений функции, вычисленной на этом пространстве. Для этого необходимо сформировать файл значений поля в вершинах сетки элементов в соответствии с определенными правилами.
В составе NetGen имеется утилита, использующая эту же библиотеку, которая позволяет просматривать и поверхности заданные в STL, и файлы с сеткой конечных элементов. В этой программе можно сгенерировать 3D-сетку элементов по STL описанию, выполнить повторное разбиение, просмотреть отдельные тетраэдры сетки.
Заключение
Использование Freefem++ и NetGen позволило в короткие сроки разработать весьма сложную по своей математической функциональности программу. При этом лицензии на Freefem++ и NetGen не ограничивают лицензию разрабатываемого продукта. Имеются некоторые проблемы использования Freefem++ и NetGen, но они решаемы. Поэтому Freefem++ и NetGen целесообразно использовать в разработках программных продуктов, аналогичных Cleanroom.
Автор: InrecoLAN
