Содержание к диссертации
Введение
1 О моделировании горения 10
1.1 Обзор литературы по проблеме исследования 10
1.2 Методы решения нестационарных уравнений газовой динамики . 24
1.3 Метод крупных частиц 26
1.4 Распараллеливание и переносимость программного кода . 29
1.4.1 Использование многопроцессорных систем 29
1.4.2 Переносимость программного кода 33
2 Математическое моделирование многомерных процессов в химически реагирующих проницаемых средах 47
2.1 Предположения и допущения 47
2.2 Основные динамические уравнения 49
2.2.1 Уравнения движения газа в отсутствии химически реагирующего вещества 49
2.2.2 Уравнения движения газа в присутствии химически реагирующего вещества 51
2.2.3 Начальные и граничные условия 5?
2.3 Полная математическая модель 54
3 Модификация метода крупных частиц 57
3.1 Выбор сеточной области 57
3.1.1 Неравномерные сетки 57
3.1.2 Равномерные сетки 59
3.2 Построение разностной схемы 60
4 Анализ математической модели и ее применение 67
4.1 Апробация математической модели многомерных процессов в химически реагирующих проницаемых средах 67
4.2 Исследование устойчивости по входным данным на примере решения задачи о горении пороха в замкнутом объеме
4.3 Исследование сходимости решение при уменьшении временного шага и измельчении сетки 79
1 4.4 Параллелизм в задачах численного моделирования 81
4.4.1 Параллелизм типа "коллективного,, решения 83
4.4.2 Геометрический параллелизм 84
4.4.3 Масштабируемость модели 8*5
4.5 Применение математической модели к различным типам задач . 91
4.5.1 Влияние высокого давления на скорость горения 91
4.5.2 Прохождение ударных волн через проницаемую преграду . 95
4.5.3 Моделирование работы донного газогенератора 97
Заключение 100
Литература
- Методы решения нестационарных уравнений газовой динамики
- Уравнения движения газа в отсутствии химически реагирующего вещества
- Неравномерные сетки
- Исследование устойчивости по входным данным на примере решения задачи о горении пороха в замкнутом объеме
Введение к работе
Актуальность работы. Научное исследование горения началось в XVIII в. вместе со стремительным развитием химии. На первоначальном этапе горение определялось как соединение с кислородом горючих веществ (в первую очереди содержащих водород и углерод).
Выяснение химической сущности горения подготовило базу для развития энергетики и термодинамики, поскольку горение - основной поставщик газов высокой температуры и энергии.
Этап изучения горения и взрывов, начавшийся в конце XIX в. и продолжающийся до настоящего времени, был связан с появлением двигателей вну* треннего сгорания, развитием внутренней баллистики артиллерийских орудий и взрывного дела, это связано с широким внедрением в технику реактивных двигателей. Это во многом стимулировало быстрое развитие науки о горении.
Большой вклад в развитие теории горения внесли В.А. Михельсон, P. Vieille, М. Berthelot, Е. Jouguet, J. Taffanel, P. Daniell, D.L. Chapman.
На современном этапе исследований процесса горения в центре внимание стоит вопрос о скорости химического превращения.
Известные отечественные ученые (ДА. Франк-Каменецкий, Я.Б. Зельдович, Н.Н. Семенов, К.К. Андреев, А.Ф. Беляев, ЮА. Победоносцев, П.Ф. Похил, А.Г. Мержанов, Б.В. Новожилов) и зарубежные ученые (Th. von Karman, М. Summerfield, А.К. Oppenheim, G.H. Markstein, FA. Williams и др.) обогатили науку о горении и ее приложения.
Теория горения, как часть математической физики, включает и использует достижения многих родственных наук: теорию тепло- и массообмена, газодинамику реагирующих потоков, химическую кинетику, турбулентное движение газа и др.
Известные модели горения твердого топлива обнаруживают существенные недостатки: только одно- или двумерная постановка задачи, зависимость модели от химического состава используемого топлива, от геометрии заряда, использование эмпирических коэффициентов, для получения которых необходимо проведение эксперимента, громоздкая численная реализация. Кроме того, в имеющихся моделях не учитывалась проницаемость твердого топлива, то есть влияние окружающей среды на процесс горения.
Поэтому построение математической модели многомерных процессов ь химически реагирующих проницаемых средах является актуальной задачей.
Цель работы. Целью диссертационной работы является построение математической модели многомерных процессов в химически реагирующих проницаемых средах и создание комплекса программ для моделирования этих процессов.
В соответствии с целью работы должны быть решены следующие задачи; построение математической модели и исследование границ применимости разработанной математической модели; модификация метода крупных частиц для численной реализации построенной модели; разработка алгоритма решения задач для двумерного и трехмерного случаев; распараллеливание расчетной процедуры (для двумерного и трехмерного случаев); кроссплатформенная (переносимая) реализация алгоритма решения на языке программирования C++; .* решение прикладных задач путем численного моделирования на основе модели.
Научная новизна работы заключается: в построении новой математической модели горения твердого топлива (для двумерного и трехмерного случаев), на основе которой возможно числен »' ное моделирование задач о взаимодействии многомерных процессов в химически реагирующих проницаемых средах в широком диапазоне начальных и граничных условий; в адаптации метода крупных частиц для моделирования процесса горения; в решении новых прикладных задач с помощью разработанной модели.
Практическая значимость. Разработанная математическая модель процесса горения твердого топлива позволяет решать задачи о взаимодействии многомерных процессах в химически реагирующих проницаемых средах. Кроме того, численное моделирование процесса горения сокращает материальные и временные затраты на натурные эксперименты.
Распараллеливание алгоритма решения математической модели позволя" ет применять ее на многопроцессорных вычислительных комплексах, что в свою очередь, снимает ограничение по размеру моделируемой задачи. Крос-сплатформенность (переносимость) программного кода расчетной процедуры обеспечивает независимость реализации модели от операционной системы компьютера или многопроцессорного вычислительного комплекса различной архи, тектуры.
Разработанная модель может быть применена в различных отраслях производства (системы вооружений, внутренняя и промежуточная баллистика и т.д.)
Достоверность результатов. Достоверность полученных результатов вытекает из корректной постановки задачи и обоснованности применяемого ме* тода крупных частиц; обеспечивается проведением расчетов на ЭВМ с контролируемой точностью; подтверждается совпадением полученных решений с известными результатами для частных случаев.
Апробация работы. Основные результаты работы докладывались на конференциях: XII Международная конференция по вычислительной механике и современным прикладным программным системам (30 июня - 5 июля 2003 г., Владимир, Россия);
21 международный баллистический симпозиум (19-23 апреля 2004г., Аделаида, Австралия);
Международная конференция "Четвертые Окуневские чтения"(22-25 июн^ 2004 г., Санкт-Петербург, Россия);
Четвертая международная школа семинар "Внутрикамерные процессы, горение и газовая динамика дисперсных систем"(28 июня - 03 июля 2004 г., Санкт-Петербург, Россия).
Использование результатов. Разработанная модель внедрена в программный комплекс GasDynamicsTool, использующийся в научных и про- мышленных организациях, конструкторских бюро и высших учебных заведе-ниях. Среди них:
НИИ механики Московского государственного университета;
Институт автоматизации проектирования РАН;
Московский физико-технический институт; Sandia National Laboratory (USA). *
Диссертационная работа состоит из введения, 4 глав, заключения и списка литературы. Работа содержит 114 страниц, в том числе 26 рисунков и 6 таблиц. Список литературы включает 104 наименования.
Во введении обоснована актуальность темы диссертационной работы, сформулирована цель работы, показана новизна и практическая ценность по-лученных результатов, представлены положения, выносимые на защиту.
Первая глава диссертации посвящена анализу состояния проблемы исследования. Приведен обзор существующих моделей горения и проанализированы достоинства и недостатки этих моделей. Рассмотрены вопросы о распараллеливании и переносимости программного кода при организации расчетной процедуры. * Во второй главе осуществлено построение базовой математической модели многомерных процессов в химически реагирующих средах.
В третьей главе осуществлено разрешение математической модели с помощью метода крупных частиц. Кроме того, в главе рассмотрены различные методы построения неравномерных сеток. Показаны достоинства и недостатки применяемых в работе равномерных сеток по сравнению с неравномерными с точки зрения распараллеливания алгоритма расчета.
В четвертой главе произведена проверка математической модели на адекватность с использованием известных теоретических результатов. Исследована сходимость решения при измельчении сетки и уменьшения величины временного шага. Произведена проверка модели на масштабируемость. Приводятся примеры применения модели к различным классам задач.
В заключении приведены основные теоретические и практические результаты, полученные в диссертации.
Методы решения нестационарных уравнений газовой динамики
Среди существующих численных подходов для решений нестационарных уравнений газовой динамики широкое распространение получили конечно, разностные методы.
Конечно-разностные методы решения уравнений газовой динамики, подобно методам исследования сплошной среды, условно можно разделить на два класса: эйлеровы и лагранжевы. При этом конечно-разностная сетка будет оставаться неподвижной в эйлеровом случае и двигаться со средой в лагран-жевом. Приложение каждого из них к задачам газовой динамики имеет ка положительные, так и отрицательные моменты. Достоинства эйлеровых методов состоят в том, что можно изучать движение среды при значительных искажениях элементарных объемов жидкости. Однако возникают трудности с численной реализацией граничных условий на поверхностях, не совпадающих с поверхностями эйлеровой системы координат. Особенно сложно удовлетворить граничным условиям в задачах взаимодействия жидкости с деформируемы ; телом, когда движущаяся контактная поверхность произвольным образом пересекает неподвижные ячейки эйлеровой сетки. Напротив, в лагранжевом подходе постановка граничных условий на контактных поверхностях осуществляется достаточно просто при соответствующем выборе начальной системы координат: граница контакта жидкости с деформируемым телом должна состоять из лагранжевых координатных линий конечно-разностной сетки. Но лагранжеу. метод неприменим в задачах, где среда сильно искажается. В этом случае ла-гранжевы ячейки, двигаясь вместе со средой, также будут вытягиваться, выворачиваться, перехлестываться. Движение жидкости приводит к значительному искажению сетки, следствием чего становится нарушение порядка аппроксимации производных конечными разностями и невозможность дальнейшего расче-та.
Таким образом, при рассмотрении задач, где искажения элементарного объема жидкости незначительны, целесообразно использовать лагранжев способ описания. Конечно-разностные методы интегрирования уравнений газовой динамики в переменных Лагранжа подробно обсуждаются в работах [66,74,83,94,96]. При значительных искажениях элементарного объема жидко-сти необходимо применять другие подходы: методы подвижных координат [47], совместный эйлерово-лагранжев метод [60], смешанный метод [75], метод полос [73], модифицированный метод маркеров и ячеек [104], методы перестройки, основанные на законах сохранения массы, импульса и энергии [95], произвольный лагранжевы-эйлерово метод [80,93].
В последнее время интенсивно развиваются так называемые смешанные лагранжево-эйлеровы методы, описание которых и многочисленные примеры можно найти в работах [14,17,18,20,22,23,34]. Они объединяют положительные качества эйлерова и лагранжева подходов: наряду с тем, что допускаются сильные искажения среды, можно следить за контактными подвижными поверхностями. Поэтому именно смешанные лагранжево-эйлеровы численные методы находят все большее приложение в задачах газовой динамики.
Ф. Харлоу [79] и др. предложили в 1955 г. оригинальный численный метод частиц в ячейках для расчета нестационарных задач гидродинамики. Можно отметить наиболее существенные моменты этого подхода.
- Дискретное представление среды. В методе частиц в ячейках сплош ная среда трактуется дискретным образом и представляется как некая совокупность (ансамбль) частиц фиксированной массы, которые движутся из ячейки в ячейку в неподвижной сетке координат. При этом масса, импульс и энергия каждой частицы вычитаются из соответствующих величин прежней ячейки (откуда частица ушла) и прибавляются к значениям в новой ячейке, куда частица переместилась. Плотность в каждой ячейке здесь определяется как частное от деления общей массы всех частиц этой ячейки на ее объем (площадь - в двумерном случае), и, таким образом, закон сохранения массы всегда удовлетворяется при этом автоматически с определенной точностью.
- Эйлерово-лагранжевы представление. По существу, указанный под ход использует совместное эйлерово-лагранжево представление. Область решения здесь разбивается неподвижной, фиксированной по пространству (эйлеровой) расчетной сеткой, а сплошная среда трактуется дискретной моделью - рассматривается совокупность частиц фиксированной массы (лагранжева сетка частиц), которые движутся через эйлерову сетку ячеек. Частицы служат для определения параметров самой жидкости, в то время как эйлерова сетка используется для определения параметров поля.
Уравнения движения газа в отсутствии химически реагирующего вещества
В математической модели взаимодействие газа, выделяющегося при горении, с окружающей средой определяется следующими параметрами: давлением (Р, [Па]), плотностью (р, [кг/м3]), вектором скорости (V, [м/с]), показателем адиабаты (7), коволюмом (а, [м3/кг]), удельной теплоемкостью при постоянном объеме (су, [Дж/кг - К]), внутренней энергией (є, [Дж]).
Для описания поведения сплошной среды используются интегральные зэ? коны сохранения массы, импульса и энергии [44]: где Q - некоторый объем пространства, S - площадь поверхности, ограничи — вающей Г2, р - плотность газа, V - скорость движения газа, Р - давление, є -внутренняя энергия. Уравнение состояния имеет вид [81]: -4) где 7 - показатель адиабаты газа, а - коволюм газа м3/кг.
Каждая компонента газа характеризуется показателем адиабаты, теплоемкостью при постоянном объеме (су), коволюмом. Вдоль линии тока параметры, характеризующие внутренние свойства газа, не изменяются, т.е. при лагранжевом рассмотрении вдоль линии тока имеем
или, переходя к привычным эйлеровых координатам, получим: д У - - — + div V) — ydivV = О С/6 или + Vgradj = 0. (2.5)
Если умножить (2.5) на р и сложить с уравнением неразрывности: i-div(pV) = 0У умноженным на 7 5 то получим уравнение, эквивалентное по форме уравнению неразрывности !t + divfrpV) = 0. (2.6) Перепишем уравнение (2.6) в интегральной форме: jt j 1PdQ = -jlpV- dS. (2.7) n s
Рассуждая также, как и в случае с показателем адиабаты у, получаем ещі два дополнительных уравнения для определения коволюма а и теплоемкости при постоянном объёме:
Уравнения (2.7, 2.8, 2.3) представляет собой замкнутую систему из восьми уравнений для определения девяти переменных: трех компонент вектора ско-рости V, давления Р, плотности р, внутренней энергии є, показателя адиабаты 7, коволюма а, удельной теплоемкости при постоянной объеме су. Эта система описывает поведение газа, образующегося при горении твердого топлива.
Данная система уравнений описывает поведение сплошной среды в объеме Q в отсутствии процесса горения.
При горении твердого топлива в объеме W, величина Q является функцией от времени. Структурная концентрация c(t) будет изменяться в соответствии с законом (2.2). Объем, занимаемый твердым топливом (WTB) будет уменьшаться со временем, пока не превратится в 0. Обозначим изменяющийся во времени объем, доступный газу через Q = Q(t). Следовательно, площадь поверхности, замыкающей объем Q(t), будет тоже функцией от времени. Обозначим ее как S(t). S(t) = [ОД]2/3 (2.9)
В определяющих уравнениях (2.3) из-за горения твердого топлива появятся дополнительные члены. Запишем уравнения (2.3) для изменяющегося во времени объёма: - [ pdQ = - I pVdS + Mi, Q(t) S(t) (2.10) dt n(t) - J pedtl I pVdtl = - і (pVV + P)dS, S(t) {peV + PV)dS + Ei, Q(t) S(t) где Mi, Ei - приход массы и энергии в изменяющийся объем в единицу времени за счет химически реагирующего твердого топлива. Приход массы в единицу времени есть произведение плотности твердого топлива на изменение объема 0,(t), происходящее из-за горения топлива. dQ(t) (2.11) dt ТВ Mi = p где /9ТВ - плотность твердого топлива. Приход энергии в объём Q(t) осуществляется за счет химической реакции горения твердого топлива: Ei = Mi Q, ще Q - тепловыделение, [Дж/кг]. Последнее выражение с учетом (2.10) перепишется в следующем виде: dQ(t) (2.12) dt El = Ртв Q Модифицируем уравнения (2.7, 2.8), описывающие поведение показателя адиабаты, коволюма, удельной теплоемкости при постоянном объеме для переменного во времени объема:
Неравномерные сетки
В отличие от неравномерных сеток, разбиение расчетной области на равномерные ячейки обладает целым рядом преимуществ. Разностная схема для равномерной сетки записывается более простыми формулами, то есть алгоритм расчета становится более прозрачным для поиска ошибок при его программировании. Использование равномерной сетки дает выигрыш по скорости при расчете на одинаковом количестве ячеек по сравнению с неравномерной сет кой.
Расчет на равномерной сетке занимает меньше оперативной памяти, так как в этом случае не приходится держать в памяти координаты узлов сетки. Координата любого узла регулярной сетки вычисляется явно.
При разбиении расчетной области на равные части для распараллеливания процесса расчета выгоднее использовать регулярную сетку. В этом случае в каждой подобласти будет примерно равное количество ячеек, что обеспечит равномерную загрузку процессоров.
При использовании неравномерной сетки количество ячеек в подобластях может существенно различаться. Это различие появляется со временем при использовании адаптивных сеток. В некоторых подобластях происходит сгущение сетки, в других — разрежение. Использование адаптивных сеток приводи к появлению проблемы сшивки на границах подобластей при перестройке расчетной сетки. Все это негативно отражается на скорости расчета.
Понятия сетки и сеточного узла являются основными при построении большой группы приближенных методов решения прикладных задач газовой динамики, называемых сеточными методами. В таких методах непрерывное пространственное распределение искомых величин и описание их непрерывного изменения во времени представляют совокупностью их значений в узлах временной сетки. При этом производные искомых функций, входящие в дифференциальные уравнения краевые условия, приближенно заменяют (аппрок-симируют) в каждом узле конечными разностями. В итоге исходную математическую формулировку задачи сводят к системе уравнений (в общем случае нелинейных) относительно неизвестных узловых значений. Такие уравнения называют разностными, а их систему вместе с правилами их построения называются разностной схемой.
Описанный подход применяется к одному из наиболее применяемых ва - риантов метода сеток — методу конечных разностей приближенного решения прикладных задач.
Приведенная в параграфе 2.3 математическая модель будет разрешаться численно с использованием метода крупных частиц [9], так как нахождение аналитического решения системы нелинейных интегральных уравнений является очень сложной (если не невозможной) задачей.
Рассмотрим некоторую область пространства W, разбитую на множество равномерных кубических ячеек с ребром h. Величины с целыми индексами относятся ко всей ячейке. Величины с дробными индексами і + 1/2, j + 1/2 и т.д. относятся к граням ячейки i,j,k. Верхний индекс п - номер временного шага.
В начальный момент времени все величины ячейки заполняются исходя из начальных условий для газа (2.14) и твердого топлива.
В начале расчета необходимо определить некоторые величины. Структурная концентрация определяется из закона горения (2.2), записанного в конечно-разностном виде: откуда Объем ячейки, доступный газу, и площадь грани вычисляются по формулам: Щі,к = (1 - clj,k)h » s?j,k = ( )
Введем в рассмотрение характеристику ячейки, называемую проницаемостью и будем вычислять ее по следующей формуле: \,j,k = №i,j,k)
Метод крупных частиц является развитием метода "частиц в ячейках", предложенного Ф. Харлоу в 1955 году. Основная идея метода крупных частиц, предполагает расщепление физического процесса на подпроцессы. Вначале изучается состояние подсистем, находящихся в ячейках, в предположении их замо-роженности или неподвижности, затем рассматривается смещение всех частиц пропорционально их скорости и шагу по времени без изменения внутреннего состояния подсистем с последующим пересчетом расчетной сетки в исходное состояние. В отличие от метода "частиц в ячейках", метод крупных частиц рассматривает не дискретную модель частиц, а непрерывные потоки массы через границы эйлеровых ячеек. Эти потоки определяются из закона сохранения массы, записанного в разностной форме для каждой ячейки (крупной частицы), совпадающей в данный момент с эйлеровой ячейкой.
Исследование устойчивости по входным данным на примере решения задачи о горении пороха в замкнутом объеме
Отметим, что метод геометрического параллелизма является методом статической балансировки загрузки, заранее определяя обрабатываемую каждым процессором часть сетки. Статическая балансировка эффективна пр ; условии, что априорной информации достаточно для предварительного распределения общей вычислительной нагрузки поровну между процессорными узлами. Метод коллективного решения является методом динамической балансировки загрузки. В его рамках перед началом вычислений не известно, какие именно узлы сетки будут обработаны тем или иным процессором. Процессоры получают задания динамически, по мере выполнения уже поступивших, чтг обеспечивает равномерную загрузку процессорных узлов при наличии большого набора независимых заданий.
Параллелизм типа „коллективного решения" удобен при проведении вычислений, распадающихся на большое количество однотипных задач, каждая из которых решается независимо от остальных. Передачи данных между такими задачами нет, а значит, полностью отсутствует необходимость их взаимной синхронизации.
Выделяется один управляющий процессор, все остальные используются в качестве обрабатывающих узлов, то есть счетных. Каждый обрабатывающий процессор выполняет элементарные задания - решение системы ОДУ для очередного узла сетки с соответствующими локальными параметрами. На управляющий процессор возлагаются функции распределения элементарных заданий по обрабатывающим процессорам и сбора полученных результатов.
В начале очередного шага каждый из процессоров ожидает новую порцию данных, обрабатывает ее, возвращает результат и снова переходит к ожиданию очередного задания, пока не получит в ответ вместо задания сообщение, что все узлы сетки уже обработаны.
Отсутствие необходимости выполнения синхронизации между элементарными задачами позволяет передавать разным процессорам разное количество расчетных узлов по мере окончания обработки данных. Тем самым решается проблема равномерной загрузки процессоров, даже если время решения системы уравнений для разных узлов сетки или производительность процессоров сильно отличается.
В случае неоднородной вычислительной нагрузки при расчете разных узлов пространственной сетки использование принципа „коллективного решения" потенциально позволяет значительно снизить время простоя и повысить эффективность распараллеливания по сравнению с рассмотренным далее мето дом геометрического параллелизма. Преимущества этого метода в полной мере могут быть реализованы, если данные для обработки изначально сосредоточены на одном из процессоров, который в этом случае может выполнять управляющие функции. При исходных данных, изначально распределенных между процессорами случайным образом, использование метода требует предварительного сбора данных, соответствующих всем расчетным точкам, на одном из процессоров. Необходимость предварительного копирования данных со всех процессоров на один и последующего возврата результатов с этого одного процессора процессорам-„владельцам" точек в значительной мере снижают эффективность данного метода, что делает его мало пригодным для практического применения при решении многих задач численного моделирования.
Исходную задачу мы можем разбить на группу областей, независимых друг от друга на каждом расчетном шаге и пересекающихся только по границе разбиения. Т.е. мы рассчитываем (п -+ 1)-м временной слой в каждой области, затем согласуем границы и переходим к расчету следующего слоя.
Однако, при таком подходе, когда мы делим расчетную область на непересекающиеся подобласти, у нас возникают проблемы с пересчетом значений на границах между данными областями, поэтому предлагается следующий достаточно логичный шаг, - делить исходную область на взаимно перекрывающиеся подобласти (4.15).
