авторефераты диссертаций БЕСПЛАТНАЯ  БИБЛИОТЕКА

АВТОРЕФЕРАТЫ КАНДИДАТСКИХ, ДОКТОРСКИХ ДИССЕРТАЦИЙ

<< ГЛАВНАЯ
АГРОИНЖЕНЕРИЯ
АСТРОНОМИЯ
БЕЗОПАСНОСТЬ
БИОЛОГИЯ
ЗЕМЛЯ
ИНФОРМАТИКА
ИСКУССТВОВЕДЕНИЕ
ИСТОРИЯ
КУЛЬТУРОЛОГИЯ
МАШИНОСТРОЕНИЕ
МЕДИЦИНА
МЕТАЛЛУРГИЯ
МЕХАНИКА
ПЕДАГОГИКА
ПОЛИТИКА
ПРИБОРОСТРОЕНИЕ
ПРОДОВОЛЬСТВИЕ
ПСИХОЛОГИЯ
РАДИОТЕХНИКА
СЕЛЬСКОЕ ХОЗЯЙСТВО
СОЦИОЛОГИЯ
СТРОИТЕЛЬСТВО
ТЕХНИЧЕСКИЕ НАУКИ
ТРАНСПОРТ
ФАРМАЦЕВТИКА
ФИЗИКА
ФИЗИОЛОГИЯ
ФИЛОЛОГИЯ
ФИЛОСОФИЯ
ХИМИЯ
ЭКОНОМИКА
ЭЛЕКТРОТЕХНИКА
ЭНЕРГЕТИКА
ЮРИСПРУДЕНЦИЯ
ЯЗЫКОЗНАНИЕ
РАЗНОЕ
КОНТАКТЫ

Численное решение задач гидроаэромеханики на графических процессорах

САНКТ-ПЕТЕРБУРГСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ

На правах рукописи

Карпенко Антон Геннадьевич

ЧИСЛЕННОЕ РЕШЕНИЕ ЗАДАЧ

ГИДРОАЭРОМЕХАНИКИ НА ГРАФИЧЕСКИХ

ПРОЦЕССОРАХ

Специальность 01.02.05 — Механика жидкости,

газа и плазмы

Автореферат

диссертации на соискание учёной степени

кандидата физико-математических наук

Санкт-Петербург — 2013

Работа выполнена на кафедре гидроаэромеханики математико-механического факультета Санкт-Петербургского государственного университета.

Научный руководитель: доктор физико-математических наук, профессор Матвеев Сергей Константинович

Официальные оппоненты: доктор технических наук, профессор Гурьев Юрий Владимирович (ВУНЦ ВМФ «Военно морская академия им. Н.Г. Кузнецова»), заведующий кафедрой «Механика и гидромеханика»

кандидат физико-математических наук, доцент Бестужева Алла Николаевна, доцент кафедры прикладной математики (Петербургский государственный университет путей сообщения)

Ведущая организация: Федеральное государственное унитарное предприятие Российский федеральный ядерный центр — Всероссийский научно-исследовательский институт экспериментальной физики ФГУП «РФЯЦ - ВНИИЭФ»

Защита состоится ““ 2013 г. в часов на заседании диссертационного совета Д 211.232.30 по защите диссертаций на соискание ученой степени канди дата наук, на соискание ученой степени доктора наук при Санкт-Петербургском государственном университете по адресу: 198504, Санкт-Петербург, Петродво рец, Университетский пр., 28, математико-механический факультет, ауд. 405.

С диссертацией можно ознакомиться в Научной библиотеке им. М.Горького Санкт-Петербургского государственного университета по адресу: 199034, Сантк Петербург, Университетская наб., 7/9.

Автореферат разослан “_“ 2013 года.

Ученый секретарь диссертационного совета Д 211.232.30, д.ф.-м.н., проф.

Кустова Е.В.

Общая характеристика работы

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

Проведение экспериментальных исследований часто оказывается невозможным.

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

Современными задачами вычислительной гидроаэромеханики является расчет турбулентных течений подходами метода крупных вихрей LES и прямо го численного моделирования DNS, а также расчет нестационарных течений в каналах со сложной геометрией и расчет множества вариантов для оптимизации конструкции. Для этих задач необходимы все более подробные сетки и время расчета сильно возрастает, что делает расчеты малоэффективными. В последнее время появляется необходимость в универсальных расчетных схемах, позволяю щих рассчитывать как сжимаемые течения, так и несжимаемые в рамках одного подхода.

Одним из способов ускорения вычислений является распараллеливание.

В традиционном подходе производится декомпозиция расчетной области и каж дая область рассчитывается на своем центральном процессорном устройстве (ЦПУ). Но для обеспечения согласованного решения необходима пересылка дан ных между процессорами. На пересылку данных тратится значительная часть времени. В последнее время появились графические процессорные устройства (ГПУ) ориентированные на научные вычисления. Они представляют собой век торные процессоры с большим количеством арифметико-логических устройств (АЛУ), в отличие от ЦПУ у которых всего несколько АЛУ. Для ГПУ необходи мо разрабатывать новые схемы расчета учитывающие особенности архитектуры.

Необходимо разрабатывать подходы к описанию течения жидкости и газа имею щие параллелизм на уровне данных. Это позволит использовать весь потенциал ГПУ и существенно ускорить вычисления. Вновь разработанные схемы стано вятся конкурентоспособными со схемами разработанными для ЦПУ.

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

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

1. Разработать, реализовать и протестировать схемы расчета нестационарного одномерного течения газа на графических процессорах;

2. Разработать, реализовать и протестировать схемы расчета нестационарного вязкого сжимаемого течения газа на неструктурированных сетках в много мерном случае, используя ГПУ;

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

4. Продемонстрировать возможности разработанного подхода на примере ре шения практической задачи теплогидродинамической конвекции (ТГК) при транспортировке застывающих наливных грузов (ЗНГ).

Основные результаты:

1. Разработаны, реализованы и протестированы схемы расчета нестационар ного одмерного течения газа на графических процессорах. Показана воз можность и преимущество использования ГПУ. Разработан, реализован и протестирован программный комплекс расчета на ГПУ нестационарного вязкого сжимаемого течения газа на неструктурированных сетках в много мерном случае.

2. На основе разработанного комплекса решен ряд верификационных задач, для которых имеется аналитическое решение или подтвержденные экспе риментальные данные. Продемонстрировано преимущество и существен ное уменьшение времени расчета.

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

4. Разработанный комплекс использован при решении практических задач теплогидродинамической конвекции (ТГК) при транспортировки застыва ющих наливных грузов (ЗНГ).

Научная новизна:

1. Разработаны модификации схем расчета нестационарных одномерных те чений газа на графических процессорных устройствах ГПУ.

2. Разработаны схемы расчета нестационарного вязкого сжимаемого газа на неструктурированных сетках в многомерном случай с использованием ГПУ.

3. Разработан универсальный подход расчета течений газа при малых числах Маха и несжимаемых течений, а также течений со свободной конвекцией.

Подход основан на предобусловливании системы уравнений сжимаемого газа.

4. Впервые аналитически выведены матрицы диагонализации якобиана пото ка предобусловленной системы в многомерном случае.

Практическая значимость Разработка универсального программного комплекса расчета задач гидроаэромеханики на ГПУ позволит повысить быст родействие и обеспечит оперативность вычислений. Результаты работы и разра ботанные схемы расчета были внедрены в пакет программ трехмерного имита ционного моделирования на супер-ЭВМ ЛОГОС, являющегося отечественным аналогом ANSYS. Пакет программ инженерного анализа ЛОГОС предназначен для моделирования процессов тепломассопереноса (аэродинамика, газодинами ка, гидродинамика, теплопроводность) и разрабатывается Российским федераль ным ядерным центром — Всероссийский научно-исследовательский институт экспериментальной физики ФГУП «РФЯЦ-ВНИИЭФ» г. Саров.

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

Апробация работы. Основные результаты работы докладывались на:

• XIII Международный семинар "Супервычисления и математическое моде лирование" (Саров, ФГУП «РФЯЦ-ВНИИЭФ», 3-7 октября, 2011);

• XVIII Школа-семинар молодых ученых и специалистов под руководством академика РАН А.И. Леонтьева "Проблемы газодинамики и тепломассооб мена в новых энергетических технологиях" (Звенигород, 23–27 мая, 2011);

• XXII Всероссийский семинар с международным участием по струйным, отрывным и нестационарным течениям (Санкт-Петрбург, 22-25 июня, 2010);

• XVII Школа-семинар молодых ученых и специалистов под руководством академика РАН А.И.Леонтьева "Проблемы газодинамики и теплообмена в аэрокосмических технологиях" (Жуковский, 25–29 мая, 2009);

• Международная конференция по механике и баллистике «VI Окуневские чтения» (Санкт-Петербург, 23-27 июня, 2008);

• XVI Школа-семинар молодых ученых и специалистов под руководством академика РАН А.И.Леонтьева "Проблемы газодинамики и теплообмена в энергетических установках" (Санкт-Петербург, 21-25 мая, 2007);

• III Всероссийская научная конференция "Проектирование научных и ин женерных приложений в среде MATLAB" (Санкт-Петербург, 23-26 октября, 2007).

Публикации. Основные результаты по теме диссертации изложены в печатных изданиях, 3 из которых изданы в журналах, входящих в перечень ре цензируемых научных журналов и изданий, 7 — в трудах российских, междуна родных конференций и семинаров. В работе [1] автору принадлежит разработ ка и реализация метода постановки дозвуковых граничных условий для одно мерной системы уравнений газовой динамики, а также расчет тестовых задач.

Емельянов В.Н. в этой работе выступал в роли научного консультанта и разра ботчика идеи постановки характеристических граничных условий. В работе [2] автору принадлежит обзор архитектуры графических процессоров, использова ния их для расчетов и технологии программирования. Соавторам принадлежит реализация схемы решения уравнения Лапласа, расчет течения в каверне. В ра боте [3] автору принадлежит разработка и реализация метода контрольного объ ема на графических процессорах, решение одномерной задачи Сода и задачи об ударной трубе в трехмерной постановке. Смирнов П.Г. и Тетерина И.В. выпол нили расчет течения в плоском канале. Волкову К.Н. принадлежит подробное описание алгоритма расчетов и общей архитектуры ГПУ. В работах [4], [5] ав тором представлен обзор возможности и перспективы использования ГПУ для расчета течений жидкости и газа. В работе [6] Емельянову В.Н. принадлежит общая постановка задач. В работах [10], [12] Азарову М.А. принадлежит раз работка математической модели течения в лазерном резонаторе, а автору расчет протекающих газодинамических процессов.

Объем и структура работы. Диссертация состоит из введения, четырех глав, заключения и приложения. Полный объем диссертации 178 страниц текста с 84 рисунками и 3 таблицами. Список литературы содержит 98 наименований.

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

Первая глава посвящена описанию архитектуры графических процес сорных устройств ГПУ, их отличию от традиционный центральных процессор ных устройств ЦПУ. Описаны преимущества и особенности использования ГПУ для вычислений. Дан обзор инструментов и технологий программирования ГПУ.

Произведен обзор литературы по использованию ГПУ для расчетов течений жидкости и газа.

Вторая глава посвящена использованию графических процессоров для решения задач гидроаэромеханики. На первом этапе производится численное ре шение системы уравнений нестационарного одномерного течения газа на ГПУ.

Рассматривается численное решение краевой задачи для системы уравнений од номерной нестационарной газовой динамики:

u U F (U ) = 0, где U = u, F (U ) = u2 + p +, t x (E + p)u E U (x, 0) = U (0) (x), (1) U (0, t) = Ul (t), U (L, t) = Ur (t).

Рассматривается эволюция начального распределения функции консервативных переменных U (0) (x) в области x [0, L] во времени. Вначале и в конце рассмат риваемой области задаются значения функций от времени. Такой постановкой можно описывать течение в ударной трубе после разрыва мембраны. Поэтому часто эту задачу называют задачей об ударной трубе.

Для дискретизации уравнения (1) используется метод контрольного объ ема. С точки зрения дискретизации по времени разностные схемы метода кон трольного объема могут быть явными или не явными. Явные схемы оказываются лучше согласованы с конечной скоростью распространения возмущений, харак терной для гиперболических уравнений, ограничивая их перенос одним шагом сетки за один шаг по времени. Наоборот, для дискретизации параболических уравнений, которые характеризуются мгновенной скоростью распространения возмущений, рационально использовать неявные схемы, снимающие жесткие ограничения на шаг интегрирования по времени.

Ограничение на шаг по времени вытекает из условия, что волны, образо вавшиеся после распада разрыва на грани контрольного объема, не достигли его центра, а при более слабом ограничении - другой грани. Схема Эйлера перво го порядка обеспечивает положительность решения если CCF L = 2|max|t 1, r где r = 2 xi - расстояние между центром ячейки и ближайшей гранью. Это условие должно быть выполнено для всех ячеек, поэтому за шаг интегрирования берется наименьший.

Различные версии метода контрольного объема различаются способом вычисления потоков на грани Fi+ 2. В работе Годунова предложено вычислять поток путем решения задачи о распаде произвольного разрыва на грани между соседними ячейками, так как решение является автомодельным. Для расчета по тока на грани требуется не вся информация из этого решения, а только значения на грани в следующий момент времени после распада разрыва. Таким образом, для расчета потоков на грани необходимо решать задачу о распаде разрыва на каждой грани, по полученным значениям параметров на грани рассчитывать зна чения потоков. Процедура расчета задачи о распаде разрыва является основной и наиболее затратной. Такой алгоритм ресурсоемкий, но он надежный и позволяет считать сильные разрывы и использовать уравнение состояния несовершенного газа. Эта схема получила название ”схема Годунова“ и имеет первый порядок аппроксимации по пространству.

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

1 Fi+ 2 = (F (UR ) + F (UL)) Ai+ 2 (UR UL), 1 2 где 2 Ai+ 2 (UR UL ) можно рассматривать как некую диссипативную добавку в схему с центральными разностями.

Для тестирования алгоритмов расчета используется набор стандартных верификационных задач. Эти задачи либо имеют точное аналитическое решение, либо достоверное численное решение. Реальные течения газа характеризуются широким диапазоном изменения параметров давления, плотности и скорости.

Поэтому для тестирования алгоритмов расчета используют понятие ”вычисли тельный газ“. Начальные условия для тестовых задач задаются в физических переменных и сведены в табл. 1. Во всех тестовых задачах сетка имеет внутренних ячеек, а граничные условия в фиктивных ячейках такие же как и в прилегающих внутренних. Контрольные времена выбраны таким образом, что возникающие при решении волны не успевают выйти за границу расчетной об ласти, что упрощает постановку условий на границах расчетной области в чис ленных реализациях задачи.

На рис. 1 представлено распределение плотности для тестовых задач.

Сплошной линией показано точное аналитическое решение, красной линей от мечено решение с первым порядком аппроксимации по пространству, синей со вторым.

Таблица 1: Тестовые задачи.

Тест L uL pL R u R pR tf in 1 Задача Сода 1.0 0.75 1.0 0.125 0.0 0.1 0. 2 Две волны разрежения 1.0 -2.0 0.4 1.0 2.0 0.4 0. 3 Сильная ударная волна 1.0 0.0 1000.0 1.0 0.0 0.01 0. 4 Движущийся контактный разрыв 1.4 0.1 1.0 1.0 0.1 1.0 2. Density, t=0.20062 dh=0.01 Density, t=0.12007 dh=0. 1 0.8 0. 0.6 0. 0.4 0. 0.2 0. 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 L L Density, t=0.011893 dh=0.01 Density, t=2.0012 dh=0. 1. 1. 1. 1. 1. 1. 1. 1. 1. 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 L L Рис. 1: Распределение плотности.

Расчеты производились на одном из модулей Tesla S1070 с частотой ядер 1.44 ГГц (количество ядер 256) и на одном ядре ЦП AMD Phenom 2 с часто той ядра 3 ГГц. Для сравнения скорости вычисления расчеты проводились на сетках с различным числом ячеек. При переходе от сетки к сетке число ячеек увеличивается на порядок (1024 ячеек для сетки 1, 30720 ячеек для сетки 2 и 307200 ячеек для сетки 3). Сетка наилучшей разрешающей способности содер жит около 3 миллионов ячеек (сетка 4). Время, необходимое для расчета одного шага (производилось осреднение по 10 шагам), а также ускорение S решения задачи Сода приводятся в табл. 2 (время приводится в миллисекундах). Вариант 1 соответствует расчету по схеме Годунова, вариант 2 — расчету по схеме Рое.

Таблица 2: Задача Сода. Время и ускорение вычислений.

№ Сетка 1 Сетка 2 Сетка 3 Сетка ЦПУ ГПУ ЦПУ ГПУ ЦПУ ГПУ ЦПУ ГПУ S S S S 1 1.63 0.13 12.43 47.70 0.20 245.25 460.64 0.92 502.50 4627.61 8.06 574. 2 0.14 0.07 1.87 5.51 0.17 33.17 43.58 0.57 76.00 436.09 5.22 83. На следующем этапе производится расчтет трехмерных сжимаемых тече ний газа на неструктурированных сетках с использование ГПУ. Подход, описан ный для одномерного случая, можно обобщить на трехмерный случай. Для глад ких решений можно повысить порядок аппроксимации схемы по пространству.

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

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

За центр ячейки принимается ее геометрический центр тяжести. За центр грани выбирается ее геометрический центр тяжести.

Для верификации разработанного программного комплекса был выпол нен расчет ряда тестовых задач, которые имеют либо аналитическое решение, либо достоверные экспериментальные данные. Это позволяет говорить о приме нимости подхода к решению задач течений жидкости и газа для рассматривае мых в тестовых задачах режимах.

Взаимодействие прямой ударной волны с клином. В рассматривае мой задаче прямая ударная волна с числом Маха ударной волны MS = 1.7 и интенсивностью I = 3.2 распространяется вдоль оси OX рис. 2 слева. Ударная волна набегает на клин, расположенный под углом = 25. В общем случае, в зависимости от MS и существует четыре возможных типа взаимодействия.

При заданных параметрах происходит нерегулярное отражение с образованием прямого скачка уплотнения.

Рис. 2: Слева схема расчетной области. Справа поле градиента плотности вблизи передней кромки клина.

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

Расчет трансзвукового обтекания профиля NACA0012 под углом ата ки. Выполняется расчет обтекания плоского крыла трансзвуковым потоком воз духа. Схема расчетной области показана на рис. 3 слева внизу. На внешней гра нице расчетной области задаются граничные условия типа Римана, число Маха потока M = 0.75, поток натекает на профиль под углом атаки = 4.

Рис. 3: Слева вверху прямоугольная расчетная сетка. Справа вверху треугольная расчетная сетка. Слева внизу схема расчетной области. Справа внизу поле чисел Маха.

Расчеты производились на различных сетках. Справа рис. 3 вверху сетка полученная триангуляцией области, слева вверху прямоугольная сетка О-типа. В данном течении локально образуется сверхзвуковая зона над профилем. На рис.

3 справа внизу представлено поле чисел Маха.

Обтекание бесконечно тонкой пластинки вязким потоком. Рассчиты валось обтекание бесконечно тонкой пластинки длиной L = 1м вязким потоком воздуха. Скорость набегающего потока равна V = 10 м/с. Для разрешения по граничного слоя к пластинке выполнено сгущение. Для того чтобы сетка у краев пластинки не была сильно вытянута производилось сгущение к краям.

На рис. 4 слева представлено поле скорости вблизи кончика пластины.

Видно, что образуется пограничный слой и его толщина увеличивается вдоль пластины. Несмотря на то, что пластинка бесконечно тонкая, вблизи кромки об разуется зона торможения и повышенного давления. Это происходит из за того, что поток тормозится об образовавшийся пограничный слой. Постановка этой задачи несколько отличается от классической задачи Блазиуса. В задаче Блази уса решаются уравнения пограничного слоя на бесконечной пластине. Несмотря на это, на рис. 4 справа приведено сравнение полученного решения с решени ем задачи Блазиуса. График приведен в безразмерных координатах, сплошной линией отмечено аналитическое решение, а красной и зеленой линией обозна чено решение в двух срезах на некотором расстоянии от кромки. Видно, что решение также является автомодельным и количественно совпадает с решение задачи Блазиуса.

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Vx/Vx besk Рис. 4: Слева поле модуля скорости. Справа распределение безразмерной скорости в сечениях x = 0.01 м и x = 0.1 м.

Произведено сравнение времени расчета на ГПУ и центральном процес соре. Расчеты производились на сетках различного размера от 10000 ячеек до 10 млн., с учетом и без учета вязкости, с различными порядками аппроксима ции по времени и по пространству. Расчеты производились на одном из модулей Tesla S1070 с частотой ядер 1.44 ГГЦ (количество ядер 256) и на одном ядре ЦП AMD Phenom 2 с частотой ядра 3 ГГц. На данном оборудовании удалось полу чить ускорение в 42 раза, но для реальных вязких задач необходимо производить расчет со вторым порядком аппроксимации по пространству. В этом случае уда лось достигнуть ускорения лишь в 22 раза. Графические процессоры постоянно развиваются, увеличивая производительность. Этот подход и этот алгоритм без доработок можно использовать на более современном оборудовании и получить более высокие результаты.

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

Разработан и реализован параллельный алгоритм расчета вязкого сжима емого газа на неструктурированных сетках. Показана возможность применения ГПУ для задач вычислительной газодинамики при сложной геометрии расчетной области. Для верификации алгоритма решен ряд тестовых задач, которые отра жают основные исследуемые режимы течений. Получено существенное ускоре ние времени расчета по сравнению с временем расчета на ЦПУ. Полученные результаты позволяют применять явные схемы для расчета реальных задач тех ники со сложной геометрией на неструктурированных сетках.

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

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

Стандартные методы расчета сжимаемых течений, описанные в предыду щих главах, оказываются неэффективными при расчете течений с малыми чис лами Маха. Это происходит из за жесткости системы уравнений для сжимаемого газа. Один из методов решения этой проблемы состоит в изменении основной системы уравнений таким образом, чтобы устранить жесткость системы.

Первой работой описывающей подход видоизменения основной системы уравнений или предобусловливание была работа Chorin. Применение сжима емых алгоритмов для расчета системы уравнений несжимаемаемой жидкости выполнено путем добавления в уравнение неразрывности производной от дав p ления по времени умноженной на коэффициент. Для стационарных течений t этот член стремится к нулю, поэтому решение стремится к решению исходной системы уравнений для несжимаемой жидкости. Добавление производной по времени делает уравнения глобально гиперболическими (если рассматривать не вязкий случай) и позволяет решать их маршевыми методами.

Этот подход можно обобщить если рассматривать систему уравнений в матричной форме. Рассмотрим систему одномерных уравнений и домножим ее на некоторую матрицу P.

U Q F (2) + = 0, P Q t x U где якобиан перехода от консервативных переменных U = V U (, u, v, w, E) к физическим Q = (p, u, v, w, T ). Матрица = P V называется матрицей предобусловливания и выбирается таким образом, чтобы устранить жесткость системы при малых числах Маха. На первый взгляд изменение системы уравнения приведет к решению отличному от решения исходной системы, но в случае расчета стационарных течений первый член стремится к нулю. В случае расчета нестационарных течений используется ме тод двойных шагов. Все существующие методы отличаются выбором матрицы предобусловливания. В данной работе матрица предобусловливания имеет вид: 0 0 0 T v 0 T vx x = vy 0 0.

T vy vz 0 0 T vz H 1 vx vy vz T H + Cp Параметр водится следующим образом:

1 T =.

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

Для совершенного газа регулирующий параметр Ur можно ввести следу ющим образом: c, если |v| c Ur = |v|, если c |v| c c, если |v| c Если локальная скорость потока меньше скорости звука то включается меха низм предобусловливания и параметр Ur задается равным локальной скорости потока. Если же скорость потока выше скорости звука, то предобусловливание системы уравнений не нужно и параметр задается Ur равным скорости звука, при этом система уравнений переходит в исходную систему для сжимаемых по токов и собственные числа якобиана становятся равны собственным числам для непредобусловленной системы уравнений.

Метод контрольного объема в многомерном случае с использованием матрицы предобусловливания можно записать так:

Ni dQi I V + (Fij + Fij )Sij = 0.

dt Vi j V Вязкие потоки Fij рассчитываются явно, таким же образом как и в случае без использования предобусловливания. Второй порядок по пространству достига ется интерполированием на грань. Схема расчета потока на грани контрольного объема запишется следующим образом:

1 1 Fi+ 2 = (FR + FL ) M | |M Q, 2 где в диссипативной добавке производится диагонализация якобиана потока пре Fx добусловленной системы уравнений A = 1. В работе впервые выведены Q матрицы диагонализации M и M для многомерного случая.

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

Часто бывает интересно стационарное решение задачи течения в канале.

Решение такой задачи может быть получено методом установления. Суть мето да состоит в том, что решение стационарной задачи отыскивается как предел к которому стремится решение нестационарных уравнений при t. Инте грирование системы уравнений по времени ведется до тех пор, пока не будет достигнут стационарный режим. О достижении стационарного режимов можно судить по тому, как отличается решение на слое n от решения на слое n + 1.

Это можно делать, оценивая L2 норму приращения сеточного решения. Реше ние можно считать сошедшимся, если L2, где некоторая заданная малая величина.

Будем рассматривать течение в канале круглого сопла, профиль которого 1/ задается выражением y = 1+x. Течение в канале сопла происходит без от рыва и входная граница находится слева, а выходная справа. На входе задаются параметры соответствующие параметрам в камере, на выходе задается давление окружающей среды. Для течений в канале удобно задавать перепад давления между входом и выходом.

При перепаде давления p = 1600 Па, число Маха в критическом сече нии приближенно будет равно M = 0.3. При таких скоростях потока, подходы с использованием предобусловливания и без предобусловливания дают одинако вую скорость сходимости. Это говорит о том, что для чисел Маха M 0.3 нет необходимости использовать предобусловливание.

При уменьшении перепада давления до p = 175 Па, число Маха в кри тическом сечении приближенно будет равно M = 0.1. Решение полученное по схеме Роэ при таких числах Маха начинает отличатся от решения полученно го с использованием предобусловливания. Это связано с тем что в схеме Роэ используются абсолютные значения давления и жесткостью системы уравнений при малых числах Маха. На рис. 5 представлены графики сходимости реше ния при использовании 20 ячеек. Линия 1 соответствует невязке скорости, а линия 2 невязке давления. Видно, что решение без использования предобуслов ливания не сходится, тогда как с использованием сходится. Стоит отметить, что распределения параметров в данных случаях отличаются. Это говорит о плохой сходимости алгоритма без предобусловливания при числах Маха M 0.1.

Рис. 5: Сходимость решения. Слева без использования предобусловливания, справа с предобусловливанием Сходимость метода без предобусловливания при числе Маха M 0.3 за метно ухудшается, а с алгоритмом предобусловливаниея скорость сходимости не изменяется. Для описания свободно-конвективных течений используют систему уравнений Навье-Стокса с добавлением источникого члена потенциальных сил в уравнение импульсов. Таким течениям характерны малые скорости потока и поэтому для их расчета можно применять метод предобусловливания системы уравнений сжимаемой жидкости.

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

При описании свободно-конвективных течений используется приближение Бу синесска. Основная идея которого состоит в том, что зависимость плотности от температуры учитывается только при вычислении массовых сил. Расчетная область является двумерной и замкнутой, течение в области является нестаци онарным, но в окрестности пластины, после некоторого промежутка времени, параметры слабо меняются и можно считать течение установившимся. В началь ный момент времени воздух покоится имеет равномерную температуру T = К, затем пластина приобретает температуру T = 300 К и из за разности плотно сти нагретого и холодного газа возникает конвективное течение. На пластинке выставлялись условия прилипания и температура стенки.

0. 0. 0. 0. 0. 0. 0. 0. 0. 0 0.5 1 1.5 2 2.5 Рис. 6: Профиль безразмерной температуры в сечении x = 0.03 м.

Эта задача может быть решена аналитически в приближении погранич ного слоя. В этом подходе вводят безразмерные температуру, скорость и координату. Стоит отметить, что данная постановка задачи отличается от по становки в приближении пограничного слоя. В данном случае решаются полные уравнения Навье-Стокса. Поэтому профили скорости у кромки пластинки будут отличаться от автомодельного решения. На рис. 6 представлено сравнение про филей безразмерной температуры в сечении y = 0.03 м. полученного численно с автомодельным решением. Из графиков видно, что численное решение хорошо совпадает с автомодельным решением.

Показана необходимость изменения схем расчета сжимаемых течений для течений с малыми числами Маха. Разработан и реализован параллельный алгоритм расчета низкоскоростных и несжимаемых течений на неструктуриро ванных сетках с использованием ГПУ. Для этого подхода аналитически выве дены матрицы диагонализации якобиана предобусловленной системы для трех мерного случая. Преимущества подхода предобусловливания показаны на при мере решения квазиодномерного течения в канале сопла. Применение ГПУ и предобусловливание дает новое направление развития явных схем для расчета низкоскоростных и несжимаемых течений. Использование ГПУ позволяет при менять явные схемы и предобусловливание сжимаемой системы уравнений для расчета течений с малыми числами Маха, несжимаемых течений и течений со свободной конвекцией. Использование ГПУ делает конкурентными такие под ходы в сравнении с алгоритмами метода коррекции давления, SIMPLE, PISO и др., специально разрабатываемых для несжимаемых течений. Для верификации алгоритма предобусловливания решен ряд тестовых задач.

В четвертой главе произведен расчет практической задачи о конвекции мазута в железнодорожной цистерне.

Жидкие нефтепродукты, в том числе и мазут перевозятся в цистернах по железной дороге в холодное время года. При охлаждении у мазутов и масел происходит резкий рост вязкости. Такие вязкие нефтепродукты при перевозках охлаждаются и загустевают настолько, что их выгрузка без разогрева становится невозможной. При температурах, близких к температурам застывания, большин ство застывающих наливных грузов (ЗНГ) становятся структурированными рео логическими средами. Вязкость таких сред зависит от напряжения сдвига. При высоких же температурах, когда вязкость меньше 0.3 · max можно применять модель Ньютоновской жидкости.

На станции загрузки жидкие топлива заливаются горячими с темпера турой около Tж = 70 C. Температура окружающей среды в зимнее время в некоторых регионах может опускаться до Tв = 40 C. При этом, время ожида ния цистерн к отправке достигает несколько десятков часов. При таких условиях в цистерне из за разности температур развивается естественная или теплогид родинамическая конвекция (ТГК). При этом имеет место существенный отток тепла из цистерны. В работе Моисеева В.И. предлагается подавлять ТГК около стенки цистерны встречным восходящим потоком жидкости. Этот поток жид кости предлагается организовывать с помощью тепловых аккумуляторов (ТА).

Тепловые аккумуляторы представляют собой тела с большой теплоемкостью.

Тепловые аккумуляторы предлагается установить по контуру в цистерну.

В начале загрузки температура нефтепродуктов равна Tж = 70 C и ТА нагре ваются до той же температуры. В процессе охлаждения цистерны температура нефтепродуктов уменьшается, но температура ТА остается почти такой же и снижается намного медленнее. Тем самым, около ТА образуется ТКГ направ ленная вверх, против течения около стенки цистерны рис. 7. Этим предлагается снижать интенсивность ТГК около стенки цистерны и уменьшать теплоотдачу.

Рис. 7: Схема расчетной области для цистерн без тепловых аккумуляторов (слева) и с ТА (справа).

Будем рассматривать течения мазута марки М40. Это один из самых рас пространенных видов топлива. Зависимость вязкости нефтепродуктов от темпе ратуры хорошо описываются экспоненциальной зависимостью:

(T ) = AeBT, где для мазута марки М40 параметр A = 4.2127 · 104, а параметр B = 0.0325.

Так как цистерна достаточно длинная, будем рассматривать двумерную задачу.

Расчетная сетка является неструктурированной и к стенке цистерны и ТА вы полнено сгущение ячеек для разрешения пограничного слоя. Очевидно, задача имеет плоскость симметрии, поэтому рассматривается только половина обла сти, а на границе выставляются условия симметрии. Так как время охлаждение достаточно длительное, будем рассматривать задачу в квазистационарной поста новке. Будем рассчитывать ТГК в некоторый момент времени, когда жидкость уже остыла до температуры Tж = 50 C, а тепловые аккумуляторы до температу ры TТА = 55 C. Поэтому, в начале расчета в обоих случаях жидкость покоится и имеет температуру Tж = 50 C. На стенке цистерны выставляется температура равная Tст = 0 C, а на стенке тепловых аккумуляторов TТА = 55 C. Произ водится расчет установившегося течения при таких параметрах. Это течение соответствует некоторому моменту времени охлаждения. В результате расчета выполняется сравнение интенсивности течения около холодной стенки цистер ны в обоих случаях.

На рис. 8 представлено поле скорости для варианта без установки теп ловых аккумуляторов и с установкой. Видно, что восходящий поток от ТА на правлен в противоположную сторону от течения около стенки. Встречный поток замедляет течение около стенки, тем самым снижая теплоотдачу.

Рис. 8: Векторное поле скорости для цистерн без тепловых аккумуляторов (слева) и с ТА (справа) В заключении приведены основные результаты работы, которые заклю чаются в следующем:

1. Разработаны, реализованы и протестированы схемы расчета нестационар ного одномерного течения газа на графических процессорах. Показана воз можность и преимущество использования ГПУ. Разработан, реализован и протестирован программный комплекс расчета на ГПУ нестационарного вязкого сжимаемого течения газа на неструктурированных сетках в много мерном случае.

2. На основе разработанного комплекса решен ряд верификационных задач, для которых имеется аналитическое решение или подтвержденные экспе риментальные данные. Продемонстрировано преимущество и существен ное уменьшение времени расчета задач гидроаэромеханики на графических процессорах.

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

4. Разработанный комплекс использован при решении практических задач теплогидродинамической конвекции (ТГК) при транспортировке застыва ющих наливных грузов (ЗНГ).

В Приложении А подробно описан метод решения задачи о распаде про извольного разрыва методом последовательных приближений.

В Приложении Б приведены аналитические выводы собственных чисел и матриц диагонализации якобиана потока предобусловленной системы уравне ний. Также приведен вид матрицы обратной к матрице предобусловливания в многомерном случае.

Публикации автора по теме диссертации 1. Емельянов В.Н., Карпенко А.Г. Метод постановки дозвуковых граничных условий для системы уравнений газовой динамики // Вестник тихоокеан ского государственного университета. 2013. N 2 (29). С. 29-38. (Из перечня ВАК) 2. Волков К.Н., Емельянов В.Н., Карпенко А.Г., Курова И.В., Серов А.Е., Смир нов П.Г. Численное решение задач гидродинамики на графических процес сорах общего назначения // Вычислительные методы и программирование.

2013. Т. 14. N 1. С. 82-90. (Из перечня ВАК) 3. Волков К.Н., Емельянов В.Н., Карпенко А.Г., Смирнов П.Г., Тетерина И.В.

Реализация метода конечных объемов и расчет течений вязкого сжимаемого газа на графических процессорах // Вычислительные методы и программи рование. 2013. Т. 14. N 1. С. 183-194. (Из перечня ВАК) 4. Emelyanov V.N., Karpenko A.G., Volkov K.N. Development of advanced CFD tools and their application to simulation of internal turbulent ows // Proceedings of the 5th European Conference for Aeronautics and Space Sciences (EUCASS 2013), Munich, Germany. P. 15.

5. Волков К.Н., Емельянов В.Н., Курова И.В., Карпенко А.Г. Современные вычислительные технологии в разработке и оптимизации технических устройств различного назначения // Труды XIII международного семи нара «Супервычисления и математическое моделирование». Издательско полиграфический комплекс ФГУП «РФЯЦ - ВНИИЭФ». г. Саров. 2011. С.

161-171.

6. Емельянов В.Н., Карпенко А.Г. Реализация алгоритма решения трехмерных уравнений нестационарной газовой динамики на графических процессо рах // Тезисы XVIII Школы-семинара молодых ученых и специалистов под руководством академика РАН А.И. Леонтьева: Проблемы газодинамики и тепломассообмена в новых энергетических технологиях. М: Издательский дом МЭИ. 2011. С. 56-58.

7. Жигалко Е.Ф., Карпенко А.Г., Цветков А. И. Деструкция вихря // Труды семинара «Компьютерные методы в механике сплошной среды». Издатель ство С.Петерб. ун-та. 2011. С. 3-6.

8. Карпенко А.Г. Решение задачи о распаде произвольного разрыва используя векторные алгоритмы на графических процессорах // Тезисы XXII Всерос сийский семинар с международным участием по струйным, отрывным и нестационарным течениям. Балт. гос. техн. ун-т. СПб. 2010. С. 84-86.

9. Карпенко А.Г. Разработка алгоритма расчёта течений сжимаемого газа на неструктурированных сетках // Труды XVII Школы-Семинара молодых уче ных и специалистов под руководством академика РАН А.И.Леонтьева: Про блемы газодинамики и теплообмена в аэрокосмических технологиях. М:

Издательский дом МЭИ. 2009. Т. 1. С. 353-356.

10. Азаров М.А., Карпенко А.Г. Гидродинамика и теплообмен в импульсно периодических лазерных системах // Сборник трудов Международной кон ференции по механике и баллистике «VI Окуневские чтения». Балт. гос.

техн. ун-т. 2008. Т. 3. С. 79-84.

11. Карпенко А.Г. Математическое моделирование нестационарных газоди намических процессов и теплопереноса в рабочих полостях импульсно периодических систем // Труды XVI Школы-семинара молодых ученых и специалистов под руководством академика РАН А.И.Леонтьева: Проблемы газодинамики и теплообмена в энергетических установках. М: Издатель ский дом МЭИ. 2007. Т.2. С. 29-32.

12. Азаров М.А., Емельянов В.Н., Карпенко А.Г. Математическое моделирова ние нестационарных газодинамических процессов и теплопереноса в рабо чих полостях импульсно-периодических систем // Труды III научной кон ференции "Проектирование Инженерных и научных приложений в среде Matlab". 2007. С. 5-24.



 




 
2013 www.netess.ru - «Бесплатная библиотека авторефератов кандидатских и докторских диссертаций»

Материалы этого сайта размещены для ознакомления, все права принадлежат их авторам.
Если Вы не согласны с тем, что Ваш материал размещён на этом сайте, пожалуйста, напишите нам, мы в течении 1-2 рабочих дней удалим его.