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

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

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

Алгоритмы построения сеток разрабатывались в работах российских ученых С. К. Годунова, Г. П. Прокопова, H.H. Яненко, А. Ф. Сидорова,.

A.A. Самарского, H.H. Калиткина, С. А. Иваненко, В. Д. Лисейкина, Л. М. Дегтярева, A.A. Чарахчьяна, О. В. Ушаковой, В. И. Мажукина,.

B.Ф.Тишкина и др., иностранными исследователями — A. Winslow, J.F.Thompson, Z.U.A. Warsi, C.W. Mastin, P.R. Eiseman, P.L.George, P. Knupp, G. Liao, и др.

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

Разбиение области на ячейки следует осуществлять таким образом, чтобы получить как можно точнее численное решение физической задачи. Классификация работ по различным способам построения сеток приведена, например, в [193] (Thompson, Warsi, 1982). Методы построения сеток могут быть разделены на три основных вида: 1) построение с помощью алгебраических преобразований, с использованием различных видов интерполяции-или специальных функций преобразований, 2) посредством решения дифференциальных уравнений (при этом дифференциальные уравнения могут быть различных типов: эллиптические, гиперболические, параболические, смешанного типа и др.), 3) посредством решения вариационных задач, основанные на минимизации функционалов.

Алгебраические методы просты, обеспечивают быстрое построение сетки и контроль густоты и наклона координатных линий с помощью коэффициентов в интерполяционных формулах. Однако, они обладают рядом существенных недостатков. В областях сложной формы координатные линии (поверхности) могут перехлестываться и выходить за границу области, что приводит к вырождению ячеек сетки. Они переносят особенности границы (например, изломы) вглубь области. Для устранения этих недостатков были развиты методы построения нерегулярных сеток, состоящих из треугольных ячеек в двумерном случае и тетраэдральных в трехмерном. Подробные описания некоторых алгебраических методов приведены в [193](Thompson, Warsi, 1982), [194](Thompson и др., 1985), [133](Eiseman, 1985), [157](Knupp, Steinberg, 1993), [169](Liseikin, 1999), [143](Handbook of Grid Generation, 1999), [138](Prey, George, 2000). Построение квазиравномерных сеток с помощью спецальных преобразований рассматривалось в [49](Калиткин и др., 2005). В [24, 25](Ваби-щевич, 1989,1991) сетка строилась посредством декомпозиции двумерной области на подобласти с использованием в приграничной подобласти гранично-адаптивной сетки.

При построении сеток в областях со сложной геометрией границ, как правило, используются методы, основанные на решении эллиптических и параболических дифференциальных уравнений [30, 32](Годунов, Прокопов и др., 1972,1976), [58](Прокопов, 1988), [44](Иваненко, 1997), [174](Nakamura, 1982), [56](Михалин, 1995), [194](Thompson и др., 1985), [157](Knupp, Steinberg, 1993), [74,169](Liseikin и др., 1999,2005), [143] (Handbook of Grid Generation, 1999), [201](Winslow, 1966)3~Досто-инством этих методов является гладкость координатных линий сетки. Принцип максимума, который выполняется для этих систем, обеспечивает корректность постановки задачи и невырожденность сетки для достаточно широкого класса двумерных областей с не слишком сложной формой границы (не сильно изогнутой и без острых углов, направленных внутрь области).

Гиперболические дифференциальные уравнения более просты при реализации, т.к. используют при построении сеток маршевые методы, позволяют получать ортогональные сетки [188] (Steger, Chaussee, 1980), [189](Tai и др., 1995), [143, Chapt. 5](Chan, 1999). Они удобны для построения сеток, когда внутренняя граница области (например, обтекаемый профиль) задана, а внешняя находится на «бесконечности». Но если вся граница области задана жестко, то эти методы не применимы из-за некорректной постановки задачи. К тому же гиперболические уравнения обладают свойством переносить особенности решения (например, изломы на границе) вдоль характеристик. В случае использования квазилинейных гиперболических уравнений негладкости и даже разрывы решения могут возникать и при гладких исходных данных [58] (Прокопов, 1988).

Отметим также методы, где сетка строится с помощью конформного отображения единичного квадрата на физическую область. Здесь решаются уравнения Лапласа с граничными условиями Неймана, см. [116](Chakravarty, Anderson, 1979), [172](Mastin, Thompson, 1984), [194](Thompson и др., 1985), [204] (Zhang и др., 2008). Но класс конформных отображений достаточно узок (согласно теореме Римана конформное отображение задается отображением трех точек на границе) и, как правило, эти сетки мало пригодны для расчетов.

Важными для приложений являются ортогональные сетки, когда координатные линии (поверхности) пересекаются под углом 90°, поскольку аппроксимация дифференциальных уравнений на таких сетках существенно проще. Методы построения двумерных ортогональных сеток рассматривались, например, в [66] (Сидоров, Шабашова, 1981), [57,60] (Прокопов, 1974,1998), [82](Ascoli и др., 1987), [182] (Ryskin, Leal, 1983), [194](Thompson и др., 1985), [152](Kang, Leal, 1992), [128](Duraiswami, Prosperetti, 1992), [143, Chapt.6,7](Khamayseh и др., Еда, 1999), [131](Еда, 1996), [76](Akcelik и др., 2001), [203] (Zhang и др., 2006). В [66](Сидоров, Шабашова, 1981), [60](Прокопов, 1998) было показано, что для ортогональных сеток соответствующая система дифференциальных уравнений имеет смешанный тип и ее применение для построения сеток может приводить к неустойчивым задачам. Поэтому на практике используется регуляризация уравнений, когда к ним добавляются эллиптические уравнения с малым весом [66](Сидоров, Шабашова, 1981), [60](Прокопов, 1998).

С помощью бигармонических уравнений двумерные сетки строились J в [104,185](Bell и др., 1982) и трехмерные сетки в [83](Altas и др., 2002).

Вариационные методы используются при построении сеток, удовлетворяющих ряду требований. Среди них невырожденность, гладкость, квазиравномерность, квазиортогональность и др. Первый вариационный подход в одномерном случае построения квазиравномерных сеток рассматривался в [65](Сидоров, 1966), а двумерном случае конструирование сеток с помощью конформных отображений проводилось в [29] (Годунов, Прокопов, 1967) и с помощью обратного гамонического отображения в [201](Winslow, 1966).

В [60](Прокопов, 1988) рассматривались различные способы задания функционалов при построении квазиортогональных сеток, когда функционал ортогональности регуляризуется посредством добавки эллиптических функционалов. В [21, 67](Сидоров, Ушакова и др., 1997,2003) использовалась сумма функционала равномерности, ортогональности и адаптивности с весовыми добавками для получения оптимальных сеток. В [34](Годунов и др., 1995) рассматривались квазиизометрические отображения, при которых отношение расстояния между любыми достаточно близкими точками к расстоянию между их образами ограничено сверху и снизу равномерно, а в [119](Chumakov, Chumakov, 1998,), [140](Godunov и др., 2007) были предложены алгоритмы построения квазиизометрических сеток. В [111](Brackbill, Saltzman, 1982) рассматривался функционал гладкости (для него следствием уравнений Эйлера есть обращенные уравнения Лапласа), а так же его комбинация с функционалами ортогональности и адаптивности. В [127](Dulikravich, Kennon, 1986) при конструировании двумерного целевого функционала использовалась локальная мера гладкости, определяемая как сумма квадратов разности площадей соседних по обоим направлениям ячеек, и мера локальной ортогональности. В [123](De Almeida, 1999) рассматривался функционал, взятый из теории упругости. Из распространенных функционалов для построения сеток следует также отметить функционалы длины и объема (площади в двумерном случае), см. монографии [115](Castillo и др., 1991), [157](Knupp, Steinberg, 1993), [143](Handbook of Grid Generation, 1999), [169](Liseikin, 1999), в которых имеются ссылки на статьи. Во многих работах рассматривается линейная комбинация этих функционалов вместе с функционалами гладкости и ортогональности. В [117](Chen, Jiang, 2008) сетка строилась с помощью минимизации целевой функции, включающей контроль за углами и длинами сторон ячеек сетки.

Одна из важнейших задач при построении сетки заключается в обеспечении ее невырожденности. Это обусловлено необходимостью корректной аппроксимации дифференциальных уравнений, описывающих физическое явление. Невырожденная сетка не содержит самопересекающихся ячеек, координатных линий или поверхностей, слипшихся ячеек или узлов, и в большинстве случаев ячеек с разной ориентацией ребер и граней [198](Ushakova, 2004). В вариационном методе невырожденная сетка строится с помощью гомеоморфного отображения параметрической области V из пространства переменных. ¦(в двумерном случае это квадрат или прямоугольник, а в трехмерном — куб или прямоугольный параллелепипед) с заданной квадратной (кубической) сеткой на физическую область из пространства переменных х=(а-1,.хп). Если отображение х (£): V—+Q, (Р, П — замыкания областей V, О) сохраняет ориентацию, т. е. оно гладкое, класса С1, и ^ якобиана отображения J=detx'(?) везде сохраняет знак, то оно может использоваться также и для построения криволинейной системы координат. Соответственно изначально ставилась задача поиска таких функционалов (или соответствующих им уравнений Эйлера), чтобы функции, доставляющие им минимум (являющиеся решением уравнений Эйлера)^обеспечивали гладкое гомеоморфное отображение параметрической области на физическую.

Для двумерного случая согласно теореме Радо (сформулирована в [180] (Rado, 1926) и доказана в [155](Kneser, 1926)) гармоническое отображение односвязной ограниченной области f2i на односвязную ограниченную выпуклую область Cl2 является диффеоморфизмом при условии заданного гомеоморфизма границы области dfti на dfh ¦ Поскольку в общем случае физическая область Г2 невыпуклая, то на практике рассматривают гармоническое отображение £(х): где V — параметрический квадрат (прямоугольник) с заданной квадратной сеткой. Для этого отображения условия теоремы Радо выполнены. Для построения сетки в физической области ?1 проводится замена переменных, уравнения Лапласа обращаются и решается задача Дирихле для квазилинейной системы дифференциальных уравнений с целью нахождения обратного гармонического отображения х (?):7 В [201](Winslow, 1966) обращенные уравнения Лапласа применялись для построения неструктурированных треугольных сеток, а в [192](Thompson и др., 1974) для получения структурированных четырехугольных сеток. В литературе использование этой системы уравнений для построения сеток принято называть методом Winslow, а задаваемое ими отображение обратным гармоническим, см. [194](Thompson и др., 1985), [157](Knupp, Steinberg, 1993), ^¡-(Иваненко, 1997), [169](Liseikin, 1999). Этот подход обеспечивает получение невырожденных сеток для довольно широкого класса областей'. В силу известных свойств уравнений Лапласа координатные линии криволинейной сетки получаются гладкими, а сама сетка квазиравномерной.

Однако, практика построения сеток показала, что для областей с гладкими, но сильно изогнутыми границами, см. [158](Knupp, Luczak), [16](Азаренок, 2009), использование обращенных уравнений Лапласа не обеспечивает получение невырожденных структурированных сеток при разумно приемлемом для практических вычислений количестве ячеек сетки, а качество невырожденных сеток является неудовлетворительным. Причина вырождения сеток состоит в ошибках аппроксимации квазилинейных дифференциальных уравнений, которые существенно возрастают в случае сильно изогнутых границ области [16](Азаренок, 2009). Если граница содержит направленные внутрь области изломы (острые углы), то при любой степени измельчения сетки четырехугольные ячейки вырождаются в окрестности направленных внутрь углов [157](Knupp, Steinberg, 1993), [44](Иваненко, 1997). В этом случае причина вырождения состоит в неустойчивости непрерывного отображения в окрестности особенностей границы области.

Для построения невырожденных структурированных сеток в двумерных областях в [42](Иваненко, Чарахчьян, 1988) был предложен вариационный барьерный метод. В нем аппроксимация функционала гладкости проводится таким образом, что при вырождении одного из четырех треугольников, на которые разбивается четырехугольная ячейка двумя диагоналями, площадь треугольника и, следовательно, якобиан линейного отображения, стоящий в знаменателе разностного функционала, стремятся к нулю, а сам функционал к бесконечности. Это является препятствием на пути вырождения четырехугольных ячеек в процессе минимизации дискретного функционала. Эта идея использовалась в [160](Knupp и др., 2002) при построении лагранжевых сеток с регулированием формы ячеек.

Другой важной задачей является дополнительный контроль за координатными линиями сетки, иными словами, за формой ячеек. Для этого в [29, 30](Годунов, Прокопов, 1967,1972) вводилась замена координат в параметрической области. Замену координат в параметрической области V иногда удобно представлять в виде использования криволинейной сетки в канонической области С (вместо квадратной в V) и отображение ее на О [194](Thompson и др., 1985). Эллиптические уравнения со второй пораматеризацией из [194], задаваемой локальным отображением, использовались, например, в [145](Hansen и др., 2004) для выглаживания сеток и в [153](Kaul, 2003), [199] (Villamizar и др., 2007) для управления сеточными линиями около границы области. Для управления сеточными линиями в [140](Годунов и др., 2007) использовалась суперпозиция двух отображений, квазиизометрического и конформного, а в [186](Spekreijse, 1999) композиция алгебраического и обратного гармонического отображений.

В [72] (Чарахчьян, 1999) с целью усиления контроля за формой ячеек в одном координатном направлении в качестве целевой разностной функции использовалась сумма разностного функционала гладкости из.

42](Иваненко, Чарахчьян, 1988) и функция, зависящая от расстановки узлов по выделенному направлению. Использование управляющего метрического тензора вместе со свойством барьерности дискретного функционала позволило воспроизводить произвольно заданные невырожденные сетки в двумерных областях [46] (Иваненко, 2000).

Трехмерный случай оказался значительно сложнее двумерного. Непростым является вопрос о нахождении условий, при которых исследуемое гладкое отображение является гомеоморфизмом. Например, для гармонического отображения теорема Радо не обобщается на трехмерный случай, см. [166](Liao, Liu, 1995), [162](Laugesen, 1996), [120](Clement и др., 1996). До настоящего времени не получено (и, по-видимому, не существует) условия, которое является одновременно необходимым и достаточным, обеспечивающего невырожденность гексаэдральной ячейки (называемой также линейчатой), грани которой есть линейчатые поверхности второго порядка [69, 198](Ушакова, 2001,2004). Тем не менее^ важность проблемы построения структурированных трехмерных сеток в реальных сложных областях обусловила существенный прогресс вариационных методов, см. [143](Handbook of Grid Generation, 1999), [64](Сахабутдинов и др., 1989), [195](Thompson, 1987), [161](Knupp, 2003), [21, 70](Ушакова, 2003,2007), [9,15, 96])(Азаренок, 2006,2008). В [195] (Thompson, 1987) использовалась система эллиптических уравнений для получения блочных структурированных сеток с ортогонализа-цией около границы. В [161](Кпирр, 2003) рассматривался дискретный функционал, характеризующий качество метрики сетки, а в [21,70](Ушакова, 2003,2007) функционал, полученный суммированием функционалов равномерности и ортогональности с весовыми коэффициентами. В.

146](Hansen и др., 2005) для выглаживания неструктурированных гек.

4 1 саэдральных сеток используются уравнения Бельтрами, являющий: ся V уравнениями Эйлера функционала энергии для гармонических отображений между Римановыми многообразиями (теория гармонических отображений представлена, например, в [132], Eells, Lemaire, 1988, [184],.

Schoen, Yau1, 1978). В [165](Lee, Soni, 2004) для этих целей использовались уравнения с управляющими коэффициентами, задаваемыми с помощью компонент метрического тензора дополнительного отображения из [194](Thompson и др., 1985). В [9,15,96](Азаренок, 2006,2008) был разработан вариационный барьерный метод построения гексаэдральных сеток с использованием функционала из [47] (Иваненко, 2003).

Как правило, при моделировании физического процесса существенное и резкое изменение физических параметров происходит на небольших участках рассматриваемой области. В этих зонах необходимо сильно измельчать сетку, для того чтобы получить численное решение с требуемой точностью. С другой стороны расчет на очень мелкой равномерной (или квазиравномерной) сетке для всей физической области привел бы к неоправданно большим затратам ресурсов ЭВМ, времени счета и оперативной памяти. Поэтому актуальным и важным разделом сеточных методов является построение адаптивных сеток, сгущающихся в зонах больших градиентов решения физической задачи. Вг настоящей работе рассматриваются только сетки с подвижными узлами. Такие методы позволяют сохранять регулярную структуру сетки, что, в свою очередь, существенно облегчает задачу аппроксимации дифференциальных уравнений, описывающих физический процесс. Описание некоторых методов построения адаптивных подвижных сеток и их приложения можно найти в монографиях [194](Thompson и др., 1985), [101](Baines, 1994), [44] (Иваненко, 1997), [143] (Handbook of Grid Generation, 1999), [169](Liseikin и др. 1999, 2005), [27](Гильманов, 2000), [49](Калиткин и др., 2005), [77](Advances in Grid Generation, 2007), обзорах [193] (Thompson, Warsi, 1982), [134](Eiseman, 1987), [144](Hawken и др., 1991), [54](Liseikin, 1996), [45](Иваненко, Прокопов, 1997), [191](Tang, 2005) и др.

В одномерном случае подход построения подвижных адаптивных сеток, основанный на оценке для ошибки (т.е. разности приближенного и точного решений) развивался, например, в [17] (Бахвалов, 1969),.

41](Емельянов, 1994), [68](Тихонов, Горбунов, 1964), [108](DeBoor, 1973), [183](Russel, Christiansen, 1978), [168](Lipnikov, Shashkov, 2006).

Был предложен принцип равномерного распределения узлов, основанный на создании такого распределения узлов, чтобы некоторая положительная весовая функция, умноженная на шаг сетки, была равна коэффициенту (в частном случае константе), в свою очередь зависящему от решения физической задачи, см. [200](White, 1979), [130](Dwyer и др., 1980), [III](Brackbill, Saltzman, 1982), [134](Eiseman, 1987), [101](Baines, 1994).

Методы построения сеток, основанные на оценке ошибки аппроксимации, рассматривались, например, в [173] (Miller, Miller, 1981), [38](Дегтярев, Иванова, 1993), [49](Калиткин и др., 2005).

В [36,37,55](Самарский, Мажукин, и др. 1988, 1989, 1993), [51] (Ма-жукин, Королева, 2007) был предложен и использовался метод динамической адаптации сеток для нестационарных задач. Этот подход основан на преобразовании координат, в результате которого решается система дифференциальных уравнений, где неизвестными являются как искомые сеточные функции, так и координаты узлов подвижной сетки.

В двумерном случае развивался подход связанный с получением априорных оценок. Для линейных эллиптических задач это было сделано в работе [100](Babuska, Rheinbold, 1978). Но его применение для нелинейных задач затруднено тем что ошибка решения не ограничивается ошибкой невязки уравнения. Другой подход основан на оценках ошибки интерполяции [176](Oden и др., 1986). В одномерном случае находится ошибка, которая есть разность между линейной и квадратичной аппроксимациями решения на интервале сетки, при том что решение задано в узлах сетки. Шаг «оптимальной» сетки находится из условия постоянства ошибки интерполяции решения на каждом интервале. Это требует вычисления вторых производных от искомого решения задачи и такая процедура была предложена в [205](Zienkiewicz, Morgan, 1983). Данный метод обобщается на многомерный случай, см., например, [78](Ait-Ali.

Yahia и др., 1996), [143, Chapt. 35](Handbook of Grid Generation, 1999).

Принцип равномерного распределения узлов сетки вдоль каждой из координат для уравнений Пуассона из [192](Thompson и др., 1974) с весовыми коэффициентами был предложен в двумерном случае для построения адаптивных сеток в [79, 80](Anderson, 1987, 1990). На трехмерные адаптивные сетки эта идея была обобщена в [81](Anderson, Munipalli, 1996) с использованием функции Грина.

В [75](Яненко и др., 1977) для построения адаптивных сеток в задачах газовой динамики был предложен функционал, сочетающий меру близости сетки к лагранжевой, меру деформации и меру концентрации сетки. В [129](Dvinsky, 1991) было предложено использовать теорию гармонических отображений многообразий для построения двумерных адаптивных сеток. Эта идея с использованием различных функций при конструировании мониторной метрики применялась в [190] (Tang, Tang, 2003). В [112](Brackbbill, 1993) построение адаптивных сеток выводилось из вариационной формулировки диффузионного метода [202](Winslow, 1981) и функционала с управлением по направлениям. Принцип равномерного распределения ошибки применялся, например, в [l02](Baines, 1998), [113](Russel и др., 1999) и др. В [20](Дегтярев и др., 2001) использовался функционал в виде комбинации функционалов ортогональности и адаптации сетки. В [114](Russel и др., 2002), [103] (Baines и др., 2004) рассматривался закон сохранения площади области для получения скорости движения узлов сетки. В [23](Бураго, 2004) для адаптации сетки к решению было предложено использовать уравнения нелинейной термоупругости. В [53,169](Liseikin, 1991,1999)^было предложено для построения адаптивных сеток функционал^ зависящий от отношения инвариантов метрического тензора, записывать на графике поверхности мониторной функции или мониторном многообразии. В [43,44](Иваненко, 1993,1997) рассматривался алгоритм минимизации дискретного аналога для этого функционала в двумерном и трехмерном случаях. В [46,47](Иваненко, 2000,2003,2004) было предложено в функционале, записанном на мониторном многообразии из [53](Liseikin, 1991), использовать дополнительное отображение для управления формой ячеек. В [5, 89](Азаре-нок, 2002,2003) рассматривался способ движения узлов по граничному контуру физической области, использующий процедуру условной минимизации дискретного функционала из [53](Liseikin, 1991) в двумерном случае. В [197] (Tu, Thompson, 1991) использовались уравнения Пуассона из [194](Thompson и др., 1985) с управляющей функцией для построения блочно-структурированных гексаэдральных адаптивных сеток. В [142](Hagmeijer, 1994) было предложено адаптировать сетку в параметрической области. Деформационный метод для расчета двумерных сеток был предложен в [106](Bochev и др., 1996), а его модификация с использованием линий уровня была рассматривалась и применялась для двумерных и трехмерных сеток в [167](Liao и др., 2000). В [175](Nakahashi, Deiwert, 1986) трехмерные адаптивные сетки строились на основе принципа равномерного распределения, когда наряду с контролем положения узлов в пространстве посредством весовых функций использовался контроль за скошенностью сетки с помощью механической модели, в которой узлы сетки связаны пружинами, работающими на растяжение и изгиб. В [12,14,99](Азаренок, 2007,2008) был разработан метод построения трехмерных адаптивных сеток с использованием функционала из [47,149](Иваненко, 2003, Ivanenko 2004).

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

— л «18.

Численное моделирование распространения детонационных волн в газе имеет свои трудности, поскольку ширина зоны реакции значительно меньше характерного размера счетной области. Для разрешения зоны горения потребовалось бы использовать квазиравномерную разностную сетку с десятками и даже сотнями тысяч узлов в направлении движения волны. В [73](Черкашин, 1974) при расчете без выделения фронта детонации на грубых сетках энерговыделение, имитирующее процесс детонации, производилось в точках области по критерию равенства внутренней энергии продуктов взрыва энергии в точке Жуге. При этом правильная скорость «размазанной» детонационной волны автоматически не формируется. При выделении фронта волны ее скорость определялась точно, но при этом необходимо было проводить расчеты в области с изменяющимися границами, т. е. на подвижных сетках. Здесь использовалась схема С. К. Годунова первого порядка аппроксимации, см. [32](Годунов и др., 1976). В [118](Chorin, 1976) выделение химической энергии учитывалось при решении задачи о распаде разрыва. В [39] (Дерюгин и др., 1990) моделирование одномерного движения волны проводилось на лагранжевой сетке методом С. К. Годунова. Одним из распространенных алгоритмов расчета является метод дробных шагов, когда на первом шаге вычисления проводятся на основе уравнений движения идеального газа, а на втором решается обыкновенное дифференциальное уравнение для уравнения химической кинетики, см., например, [39](Дерюгин и др., 1990), [109](Bourlioux и др., 1990), [122](Colella и др., 1986), [105](Ben-Artzi, 1989), [147](Hwang и др., 2000), [178](Pember R.B., 1993), [179](Quirk, 1994), [196](Ton, 1996). Во всех указанных работах расчет проводился на основе схем типа Годунова, т. е. с использованием решения задачи о распаде разрыва. Некоторые другие алгоритмы приведены в монографии [177](Oran, Boris, 1987).

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

При расчете реальных задач, в частности гидродинамических течений, часто возникает потребность в некоторый момент времени перейти от вычислений на одной сетке к расчету на другой. Для этого разрабатываются специальные алгоритмы консервативной интерполяции. В двумерном случае такие работы осуществлялись, например, в [124](Вико-шсг, 1984), [181](Ramshaw, 1985). В трехмерном случае для интерполяции с одной гексаэдральной сетки на другую в [126](Вико" шсг, РасИа1, 1991) был предложен метод второго порядка точности. Интерполяционная задача сводилась к построению фигуры пересечения гекса-эдральных ячеек старой и новой сеток. Линии пересечения линейчатых поверхностей в пространстве находились при решении дифференциальных уравнений. Эти уравнения могут иметь особенности в областях, где грани ячеек старой и новой сеток пересекаются и почти параллельны, в местах соприкосновения граней и т. д. В [141](Сгапс1у, 1999) гексаэдраль-ная ячейка заменялась 24-гранником с треугольными гранями, имеющим такой же объем. 24-гранник делился на 48 тетраэдров и определялись объемы фигур пересечения каждого из 48 тетраэдров новой ячейки с каждым из 48 тетраэдров старой ячейки. Этот интерполяционный метод имеет первый порядок точности. В диссертации рассматривается алгоритм консервативной интерполяции при замене гексаэдральной ячейки на два 12-гранника с треугольными гранями. [10,13, 94, 98](Азаренок, 2005;2008).

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

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

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

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

3) Предложен новый способ аппроксимации функционала на 10 тетраэдрах для каждой гексаэдральной ячейки. Предложены новые приближенные условия невырожденности гексаэдральной ячейки, которые в практических расчетах обеспечивают невырожденность гексаэдраль-ных ячеек сетки.

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

5) Предложен метод расстановки узлов на границе области, заключающийся в условной минимизации функционала.

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

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

8) На основе предложенной трехточечной модели адаптации проведен теоретический анализ свойств дискретных функционалов в одномерном, двумерном и трехмерном случае при наличии разрывов у мониторной функции. Анализ показал необходимость замораживания производных от мониторной функции при минимизации функционала.

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

10) Разработан численный метод расчета двумерных нестационарных газодинамических течений с выделением химической энергии на подвижных сетках.

11) Разработан новый алгоритм консервативной интерполяции на гексаэдральных сетках.

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

ЗАКЛЮЧЕНИЕ

.

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

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

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