Диссертация: Программа моделирования воздействия среды на движущееся твердое тело

Для получения полностью дискретизированных уравнений осталось дискретизировать уравнения (2.23), (2.24) по элементам разностной сетки. Для этого воспользуемся формулами (2.4)_(2.13):

где:

Уравнение Пуассона(2.24) в дискретизированной форме принимает следующий вид:

Для успешного решения уравнения Пуассона во всех ячейках требуется знать следующие граничные значения:

Они задаются следующим образом:

Уравнение Пуассона (2.29) должно быть модифицировано в соответствии с рассчитанными выше граничными значениями (2.31). После этого оно принимает следующий вид:

где:

Решение системы линейных уравнений

Формула (2.32) представляет собой систему из уравнений с неизвестными. Метод Гаусса является одним из самых распространенных подходов к решению подобных систем, однако он нечасто используется при решении задач вычислительной гидродинамики. Сейчас в численных методах чаще всего применяется итеративный подход. В частности, в разрабатываемом приложении используется метод релаксации.

Пусть дана система из линейных уравнений:

где:

, , ;

Матрица может быть резделена на главную диагональ и верхнюю и нижнюю треугольные матрицы и

, , ;

Таким образом, система линейных уравнений может быть записана следующим образом:

где - фактор релаксации, константа, которая выбирается в зависимости от коэффициентов матрицы . Для задач вычислительной гидродинамики рекомендуется использовать значение [14].

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

где:

- правая часть уравнения Пуассона (2.32) для ячейки

При каждой итерации высчитываетсяневязка (residual) вычислений по следующей формуле:

где выражается следующей формулой:

Существует два условия прекращения работы алгоритма:

Превышение значения максимального количества итераций ;

Достижение невязки вычислений заданной степени точности ().

В разрабатываемом приложении установлены следующие значения количества итераций и точности:

Общая схема алгоритма

Общая схема алгоритма численного решения уравнений Навье-Стокса представлена в Приложение А. Схема алгоритма численного решения уравнений Навье-Стокса

Переход от двух измерений к трем измерениям

В трехмерном случае уравнения Навье-Стокса (1.5), (1.6) записываются следующим образом:

Уравнения движения:

Уравнение неразрывности:

где - Z-компонента вектора скорости.

Дискретизация уравнение Навье-Стокса теперь происходит в трехмерном прямоугольном параллелепипеде, ориентированном параллельно осям , и

Трехмерное пространство разбивается на маленьких прямоугольных параллелепипедов со сторонами и . Расположение переменных в разностной трехмерной сети изображены на Рис. 8.

Рис. 8. Расположение переменных в трехмерной ячейке

Граничные условия практически не подвергаются измерению при переходе к трем измерениям. Дискретизация членов уравнения Навье-Стокса производится аналогично двумерному случаю (2.4)_(2.13).

После дискретизации по времени получаются следующие уравнения движения:

где - Z-компонента вектора скорости в момент времени , а членаналогичен членам (2.27) и (2.28).

При объединении с уравнением неразрывности (2.43) получается трехмерное уравнение Пуассонадля давления в интервал времени:

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

Программная Реализация

Выбранные технологии

Для решения проблемы разработки под различные операционные системы Java был выбран в качестве языка программирования. Для создания пользовательского интерфейса, который на каждой операционной системе выглядел бы натурально, была выбрана библиотека пользовательского интерфейса Swing с бесплатными дополнениями внешнего вида «QuaquaLookandFeel»[15] для операционной системы Mac OS X и «JGoodiesLooks»[16] для операционных систем MicrosoftWindows. Также для создания меню были использованы бесплатные библиотеки «L2DProdCommonComponents»[17]и «JGoodiesForms»[18]. Для контроля всех зависимостейи сборки проекта под операционные системы MacOSXи MicrosoftWindow был использован фреймворкMaven[19].

В качестве инструмента визуализации используется графическая библиотека OpenGL версии 2.0. Для того, чтобы уделить большего внимания выполнению вычислительного модуля, было принято решение использовать готовый графический движок, совместимый с OpenGL и Java. В качестве такого движка был выбран JMonkeyEngine3[20].

В качествеIDEдля разработки приложения было использовано средство разработки IntelliJIdeaот компании JetBrains[21].Git[22]был использован в качестве системы контроля версий.

Алгоритм создания стартовых значений

Одной из сложностей при разработкепрограммы было задание начальных данных для алгоритма численного решения уравнений Навье-Стокса. Требовалось преобразовать трехмерную сцену, созданную с помощью jMonkeyEngine 3 в трехмерный массив ячеек, с которым работает алгоритм. Прямое преобразование (без отображения сцены) не представляется возможным, т.к. хранение трехмерной сцены jMonkeyEngine 3 осуществляется с помощью дерева объектов[23]. Было разработано следующее решение:

Пусть даны координаты двух точек

Эти точки задают прямоугольный параллелепипед, определяющий пространство эксперимента. Также известна сторона ячейки трехмерной разностной сетки - .

Сначала определяются переменные - размер разностной сеткиcучетом дополнительных двух рядов по каждой из координат для описания граничных условий:

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

Из центра ячейки пускается луч в направлении, параллельном положительному направлении оси . Возможности движка jMonkeyEngine3 позволяют определять координаты столкновения луча с трехмерным объектом. На основании этих данных формируются зоны жидкости и зоны препятствий.

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

Твердое тело - граничное условие свободного прилипания(2.15)

Граница пространства эксперимента - граничное условие свободного прилипания(2.16)

Место присоединения трубы - граничное условие свободного потока(2.17)

Структура файлов и каталогов

При создании проекта с именем в выбранной директории создается следующая иерархия каталоговРис. 9:

Рис. 9. Структура файлов и каталогов проекта SimulationTest

В файле, где - название проекта, хранятся общие настройки проекта, включающие в себя:

Название, дату создания проекта

Координаты двух точек, ограничивающих пространство эксперимента (3.1)

Информация об исследуемой сред

Длительность эксперимента

Фиксированный временной шаг , если отключено автоматическое определение временного шага

Начальная скорость потока из трубы

Сечение пространства эксперимента, если моделирование производится в двух измерениях

В каталогеassetsхранятся файлы, относящиеся к отображению 3D-сцены на движке jMonkeyEngine 3:

Файл solid.j3oопределяет геометрию и расположение в пространстве твердого тела

Файл flowSource.j3oопределяет геометрию и расположение в пространстве источника потока.

Данные файлы генерируются стандартным механизмом сохранения сцен движка jMonkeyEngine3. Для каждого класса объектов (твердого тела и источника потока) потребовалось реализовать интерфейса jme3.export.Savable[23].

По умолчанию модель твердого тела - сфера с радиусом 1 метр. Когда пользователь загружает свою модель твердого тела, то она также копируется в эту директорию. Трехмерная модель тела задается в формате Wavefront.OBJ [24] - одном из самых распространенных форматов представления моделей, который на данный момент поддерживается большим количеством 3D-редакторов [25]

В каталогеsimulationхранятся файлы, содержащие информацию о разностной сетке.Название файла задано следующей формулой:

где:

- название проекта;

- дата начала моделирования

- шаг по времени

Файл создается путем сериализации объекта классов model.solver.CellArray2Dи model.solver.CellArray3D. Для каждого класса был реализован интерфейс java.io.Serializable[26].

Для каждой ячейки хранятся следующие параметры:

Давление

Компоненты вектора скорости

Тип ячейки. Всего существует три типа ячеек:

Ячейка жидкости

Ячейка - препятствие

Ячейка, находящаяся в месте присоединения трубы

В каталогеviewsхранятся изображения, визуализирующие результаты эксперимента. Их можно просмотреть как напрямую (стандартным программным обеспечением для просмотра фотографий), так и с помощью экрана просмотра визуализации результатов эксперимента. Формат названия файла для иллюстраций аналогичен формату названия файлов разностной сети (3.3)

Визуализация результатов эксперимента

Визуализация результатов эксперимента представляет собой набор двумерных изображений. При моделировании двумерного случая создается одно изображение разностной сетки для каждого шага по времени. При моделировании трехмерного случая возможно построение проекций на координатные плоскости . В таком случае для каждого временного интервала создается изображений проекции на плоскость , изображений проекции на плоскость и изображений проекции на плоскость . Препятствие представлено черным цветом, место присоединения трубы - синим, а цвет ячейки с жидкостью представлен линейным градиентом от зеленого к красному (Рис. 10). Зеленый цвет соответствует давлению , а красный - давлению . Оба значения давления задаются пользователем.

Рис. 10. Картина распределения давления.

Предположим, выбрана проекция разностной сетки на плоскость в сечении . Каждой ячейке разностной сетки соответствует квадрат 3х3 пикселя. Для визуализации скорости выбираются такие ячейки , для которых верно следующее условие:

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