автореферат диссертации по информатике, вычислительной технике и управлению, 05.13.18, диссертация на тему:Вариационные методы построения сеток и приложения в газовой динамике

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

Азаренок Борис Николаевич

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

05.13.18 - Математическое моделирование, численные методы и комплексы программ

АВТОРЕФЕРАТ

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

1 4 ВИЗ

Москва - 2008

003458967

Работа выполнена в Вычислительном Центре им. А.А.Дородницына Российской академии наук.

Официальные

оппоненты: доктор физико-математических наук, Геннадий Павлович Прокопов,

доктор физико-математических наук, профессор Владимир Дмитриевич Лисейкин,

доктор физико-математических наук, профессор Николай Георгиевич Бураго

Ведущая

организация: Институт математики и механики УрО РАН

Защита состоится " 19 " марта 2009 г. в 11 ч. 00 мин. на заседании Диссертационного совета Д002.058.01 в Институте математического моделирования РАН по адресу: 125047 г. Москва, Миусская пл., д. 4А.

С диссертацией можно ознакомиться в библиотеке ИММ РАН. Автореферат разослан " " декабря 2008 г.

Ученый секретарь Диссертационного совета доктор физико-математических наук

Н.В.Змитренко

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

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

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

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

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

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

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

Целью работы является:

— разработка метода построения гексаэдральных сеток в областях со сложной геометрией с возможностью управления формой ячеек для использования в реальных физических и инженерных приложениях;

— разработка метода построения подвижных адаптивных сеток, подстраивающихся к особенностям решения моделируемой задачи;

— разработка численного метода решения задач двумерного течения невязкого газа, в том числе при наличии химической реакции, на подвижных сетках и его применение;

— разработка метода консервативной интерполяции с одной гек-саэдральной сетки на другую.

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

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

'Иваненко С.А. Вариационные методы построения сеток//Ж. вычисл. ма-тем. и матем. физ. 2003. Т. 43. № 6. С. 830-844.

метрического тензора можно воспроизвести любое заданное невырожденное отображение, а, следовательно, и сетку. Свойство универсальности функционала позволяет получать ячейки сетки произвольной заданной формы. При построении сетки предложено вместо невырожденности гексаэдральной ячейки с линейчатыми гранями потребовать невырожденность двух 12-гранных ячеек с треугольными гранями, что сводится к требованию положительности алгебраических объемов 10 тетраэдров. С помощью вычислительного эксперимента показана очень низкая вероятность появления вырожденных ячеек гексаэдральной сетки при выполнении этих условий. В практических примерах построения сеток предлагаемым вариационным методом выполнение этого условия обеспечивало невырожденность гексаэдральных сеток. Построена конечномерная функция, аппроксимирующая функционал и имеющая бесконечный барьер на границе множества невырожденных 12-гранных ячеек, минимизация которой значительно более экономична и эффективна по сравнению с предложенной ранее С.А.Иваненко процедурой аппроксимации функционала на 24 тетраэдрах2. Предложено необходимое условие невырожденности гексаэдральной ячейки, используемое при проверке сетки на невырожденность совместно с достаточными условиями невырожденности О.В.Ушаковой3. Предложен алгоритм перераспределения узлов сетки по граничным поверхностям и ребрам области. С использованием свойств универсальности и инвариантности функционала предложен алгоритм ортогонализации координатных линий и сгущения координатных поверхностей сетки к границе области, а также гладкого сопряжения приграничных слоев ячеек сетки к ячейкам, расположенным внутри области. Созданный комплекс программ обеспечивает построение сеток с управлением формы гексаэдральных ячеек в областях со сложной

2Иваненко С.А. Адаптивно-гармонические сетки. М. ВЦ РАН. 1997. 181 С.

3УшаковаО.В. Условия невырожденности трехмерных ячеек. Формула для объема ячеек//Ж. вычисл. матем. и матем. физ. 2001. Т. 41. 6. С. 881-894.

геометрией.

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

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

4ГодуновС.К. Разностный метод численного расчета разрывных решений уравнений гидродинамики//Матем. сб. 1959. Т. 47. Вып. 3. С. 271-306.

5ГодуновС.К., Забродин A.B., Иванов М.Я., КрайкоА.М., Прокопов Г.П. Численное решение многомерных задач газовой динамики. М.: Наука. 1976.

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

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

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

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

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

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

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

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

Апробация работы. Результаты диссертации были доложены на: VIII,IX Всероссийских совещаниях "Проблемы построения сеток для решения задач математической физики", 2000, 2002г.; XII,XV,XVI Всероссийских конференциях "Теоретические основы и конструирование численных алгоритмов решения задач математической физики" памяти К.И.Бабенко, 1998, 2004, 2006г.; 7th Russian-Japanese Intern. Sympos. on CFD, Moscow Lomonosov Univ., 2000; Intern. Confer. "OFEA'2001. Optimization of finite-element approximations, splines and wavelets", St.-Petersburg, 2001; Confer, on Numerical Methods for Fluid Dynamics, University of Oxford, UK, 1998; 8th International Symposium on CFD, Bremen, Germany, 1999; 2d Intern. Sympos. on Finite Volumes for Complex Applications - Problems and Perspectives, Duisburg, Germany, 1999; 7th,9th,10th Intern. Conference on Numerical Grid Generation, Whilster, Canada, 2000, San Jose, California, USA, 2005, Forth, Crete, Greece, 2007; 1st Intern. Conference on CFD, Kioto, Japan, 2000; 2nd Intern. Confer. Applied Mathematics for Industrial Flows, Ciocco, Italy, 2000; XI Международной конференции по

Вычислительной Механике и Современным Прикладным Программным Системам (ВМСППС'2001), Москва, 2001г.; 9th Intern. Confer, on Hyperbolic Problems, Theory, Numerics, Applications, California Institute of Technology, Pasadena, California, USA, 2002; Intern. Confer, on Scientific Computing and Partial Differential Equations, Hong Kong, 2002; Workshop "Grid Generation: Theory and Applications'^.: ВЦ РАН, 2002г., X Всероссийском семинаре "Современные проблемы численного моделирования", Новороссийск, 2003; III Intern. Workshop on Scientif. Comput. and Applications, City Univ. of Hong Kong, 2003; VII,VIII Международных конференциях "Забабахинские Научные Чтения", РФЯЦ-ВНИИТФ, Снежинск, 2003 и 2005; Всероссийских конференциях "Прикладная геометрия, построение расчетных сеток и высокопроизводительные вычисления", ВЦ РАН, Москва, 2004 и 2006; 4-й Международной школе-семинаре "Внутрикамерные процессы, горение и газовая динамика дисперсных систем", Балтийский гос. университет Военмех, С.-Петербург, 2004г.; V,VII Международных конференциях по неравновесным процессам в соплах и струях (NPNJ-2004 и NPNJ-2008), Самара, 2004, и Алушта, 2008; XIV,XV Международных конференциях по Вычислительной Механике и Современным Прикладным Программным Системам (ВМСППС-2005 и ВМСППС-2007), Алушта, 2005г. и 2007г.; International conference "Numerical geometry, grid generation and scientific computing", ВЦ РАН, Москва, 2008; на семинарах ВЦ РАН им. А.А.Дородницына, ИПМ РАН им. М.В.Келдыша, Институте Математики им. С.Л.Соболева СО РАН, Институте Математики и Механики УрО РАН, Институте Математического Моделирования РАН, Институте Вычислительной Математики РАН, Karlsruhe University, Germany; Hong Kong Baptist University; Hong Kong University of Science and Technology.

Работа над диссертацией проводилась в рамках проектов РФФИ: "Конструирование алгоритмов построения адаптивных сеток на основе теории гармонических отображений"(1999—2001г., код проекта 99-01-00264), "Разработка алгоритмов построения

многомерных сеток и их приложения в задачах математической физики"(2002-2004г., код проекта 02-01-00236); в рамках проекта Отделения Математических Наук РАН "Современные вычислительные и информационные технологии решения больших задач" (2005-2007г.).

Результаты диссертации использовались в совместных с зарубежными учеными исследованиях в Department of Mathematics of Hong Kong Baptist University (Hong Kong Research Grant Council, Project code HKBU 2045/02P and HKBU 201/03P), International Research Team on Complex System, Chinese Academy of Sciences.

Публикации. Основные результаты опубликованы в 14 рецензируемых журнальных статьях, рекомендованных ВАК, 2 рецензируемых журнальных статьях, монографии (2 работы), 4 препринтах ВЦ РАН, 9 статьях в трудах всероссийских и зарубежных конференций, 23 публикациях тезисов докладов на всероссийских и зарубежных конференциях.

Структура и объем работы. Диссертация состоит из введения, 4 глав, заключения и списка цитируемой литературы. Диссертация содержит 258 страниц, в общей сложности 100 рисунков и 8 таблиц. Список цитируемой литературы содержит 170 наименований.

Содержание работы

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

Алгоритмы построения сеток разрабатывались в работах российских ученых С.К. Годунова, Г.П. Прокопова, H.H.Яненко, А.Ф.Сидорова, A.A.Самарского, В.Д.Лисейкина, С.А.Иваненко, Л.М.Дегтярева, A.A.Чарахчьяна, О.В.Ушаковой,

B.И. Мажукина, H.A. Дарьина, В.Ф.Тишкина и др., иностранными исследователями - A.Winslow, J.F.Thompson, Z.U.A.Warsi,

C.W. Mastin, P.R. Eiseman, P. Knupp, G.Liao, и др.

Классификация работ по различным способам построения сеток приведена, например, в Thompson et al6. Методы построения сеток могут быть разделены на три основных вида: 1) алгебраические, с использованием различных видов интерполяции или специальных функций преобразований, 2) посредством решения дифференциальных уравнений (при этом дифференциальные уравнения могут быть различных типов: эллиптические, гиперболические, параболические, смешанного типа и др.), 3) вариационные, основанные на минимизации функционалов.

В диссертации рассматриваются структурированные (по другому регулярные) сетки. Для них упорядочивание узлов задается простым образом, посредством матрицы с двумя индексами в двумерном случае и матрицы с тремя индексами в трехмерном.

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

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

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

6Thompson J.F., WarsiZ.U.A. Boundary-fitted coordinate systems for numerical solution of partial differential equations//J. Сотр. Phys. 1982. V. 47. № 2. P. 1-108.

случае с работы А.Ф.Сидорова7 и в двумерном случае с работ С.К.Годунова, Г.П.Прокопова8 и A. Winslow9.

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

квадрат или прямоугольник, а в трехмерном - куб или прямоугольный параллелепипед) с заданной квадратной (кубической) сеткой на физическую область il из пространства переменных х=(ж1,..., хп). Если отображение (V, iî - замыкания

областей V, Q) сохраняет ориентацию, т.е. оно гладкое, и якобиан отображения J=detx'(£) везде сохраняет знак, то оно может использоваться также и для построения криволинейной системы координат. Если дополнительно обратное отображение £(х) : Q.—>V является гладким, то оно называется диффеоморфизмом. Соответственно изначально ставилась задача поиска таких функционалов (или дифференциальных уравнений), чтобы функции доставляющие им минимум (являющиеся решением дифференциальных уравнений) обеспечивали гладкое гомеоморфное отображение параметрической области V на физическую О..

Для двумерного случая согласно теореме Радо (сформулированной Radö10 и доказанной Kneser11) гармоническое отображение односвязной ограниченной области Q. на односвязную ограниченную выпуклую область V является диффеоморфизмом при условии заданного гомеоморфизма дО. на дР.

Поскольку в общем случае физическая область Q, невыпуклая,

7Сидоров А.Ф. Об одном алгоритме расчета оптимальных разностных се-ток//Тр. Матем. ин-та АН СССР. М. 1966. Т. 74. С. 147-151.

8ГодуновС.К., Прокопов Г.П. О расчетах конформных отображений и построении разностных сеток//Ж. вычисл. матем. и матем. физ. 1967. Т. 7. № 5. С. 1031-1059.

9Winslow A.M. Numerical solution of the quasi-linear Poisson equation in a nonuniform triangle mesh//J. Comput. Phys. 1966. V. 1. P. 149-172.

10Radö T. Aufgabe 41, Jahresber//Deutsche Math.-Verein. 1926. V. 35. P. 49.

uKneserH., Lösung der Aufgabe 41, Jahresber//Deutsche Math.-Verein. 1926. V. 35. P. 123-124.

то рассматривают гармоническое отображение £(х): ii—где V - параметрический квадрат (прямоугольник) с заданной квадратной сеткой. Для этого отображения условия теоремы Радо выполнены. Для построения сетки в физической области проводится замена переменных, уравнения Лапласа обращаются и решается краевая задача для нелинейной системы дифференциальных уравнений с соответствующими граничными условиями для нахождения обратного отображения х(£): "Р—>0. В литературе использование этой системы уравнений для построения сеток принято называть методом Winslow9, а задаваемое ими отображение гармоническим. Этот подход обеспечивает получение невырожденных сеток для довольно широкого класса областей. В силу известных свойств уравнений Лапласа координатные линии криволинейной сетки получаются гладкими, а сама сетка квазиравномерной.

Однако, практика построения сеток показала, что для областей с гладкими, но изогнутыми границами12, или направленными внутрь изломами границ2, использование обращенных уравнений Лапласа не обеспечивает получение невырожденных структурированных сеток при разумно приемлемом для практических вычислений количестве ячеек сетки, а качество невырожденных сеток является неудовлетворительным. Для построения невырожденных сеток в двумерных областях, состоящих из четырехугольных ячеек, был предложен вариационный барьерный метод13. Другой важной задачей является дополнительный контроль за координатными линиями сетки, иными словами, за формой ячеек. Для этого вводилась замена координат в параметрической области Р8'14. Замену координат в параметрической области V иногда

12KnuppP., LuczakR. Truncation error in grid generation: a case study, Numerical Methods for Partial Differential Equations. 1995. V. 11. P. 561-571.

13Иваненко С.А., Чарахчьян A.A. Криволинейные сетки из выпуклых четы-рехугольников//Ж. вычисл. матем. и матем. физ. 1988. Т 28. X« 4. С. 503-514.

14ГодуновС.К., Прокопов Г.П. Об использовании подвижных сеток в газодинамических расчетах//Ж. вычисл. матем. и матем. физ. 1972. Т. 12 JV« 2. С. 429-440.

удобно представлять в виде использования криволинейной сетки в канонической области С (вместо квадратной в?)и отображать ее на область П15 [2,4,5,46].

Трехмерный случай оказался значительно сложнее двумерного. Непростым является вопрос о нахождении условий, при которых исследуемое гладкое отображение является гомеоморфизмом. Например, для гармонических отображений теорема Радб не обобщается на трехмерный случай16,17.

Поэтому для построения гомеоморфного отображения параметрической области V на физическую Q в диссертации используется следующий подход. Глобальное гомеоморфное отображение х(£): V—^il ищется в виде склейки гладких гомеоморфных отображений Xi(£) каждой кубической ячейки из области V в гек-саздральную ячейку из области il. Для этого используются топологические теоремы о достаточных условиях гомеоморфности гладких, непрерывных и кусочно-гладких отображений. Гладкий гомеоморфизм х(£): обеспечивается благодаря барьерно-

му свойству дискретного функционала.

Глава 1 диссертации посвящена описанию вариационного метода построения гексаэдральных сеток. Основные результаты опубликованы в [5] (J.Comput.Phys., 2006), [26] (Препринт ВЦ РАН, 2006), [1] (Математ. Моделирование, 2008).

В разделе 1.1 дается определение невырожденной сетки, конструируемой с помощью отображения х(£) параметрической области V на физическую ii.

Определение 1. Сетка называется невырожденной, если она не содержит самопересекающихся ячеек, координатных линий или поверхностей, слипшихся ячеек или узлов, и якобиан отображения J— detx'(£)>0 почти всюду в V.

"Thompson J.F., WarsiZ.U.A., Mastin^W. Numerical Grid Generation. North-Holland, N.Y. etc. 1985.

16LiuH., LiaoG. A note on harmonic maps//Appl. Math. Lett. 1996. V. 9. № 4. P. 95-97.

17Laugesen R.S. Injectivity can fail for higher-dimensional harmonic extensi-ons//Complex Variables. 1996. V. 28. P. 357-369.

14

Требование 7>0 необходимо для построения невырожденной криволинейной системы координат в П. Однако, бывают исключения, когда <7 равен нулю на подмножестве V меньшей размерности (например, в изолированных точках, кривых или поверхностях). Гомеоморфизм отображения при этом сохраняется и полученные при этом отображении ячейки можно также называть невырожденными, если они допустимы для расчета основной дифференциальной задачи. Следовательно, задача построения сетки может быть решена заданием гомеоморфного отображения х(£): Р—>$7 при заданном отображении границ областей дР—>дП.

В разделе 1.2 приведена вариационная постановка задачи построения сетки и дается вывод функционала Т> в п-мерном случае. Для дополнительного управления формой ячеек сетки строится отображение канонической области С, в которой задана криволинейная сетка, на область Г2. Задание криволинейной сетки в С осуществляется с помощью отображения V, с заданной кубической сеткой, на область С.

Таким образом, рассматривается гомеоморфное гладкое, класса С1, отображение х(Х): Еп—канонической области С из пространства переменных Х=(Х\..., Хп) на физическую область из пространства переменных х=(х1,..., хп) с матрицей Якоби аг~дхг/дХЗ, г, .7=1,2,.. .,п. Используются еще два гомео-морфных гладких, класса С1, отображения х(£),Х(£): Мп—Л" параметрической области V из пространства переменных £=(£*,..., £п) на области и и С соответственно с матрицами Якоби Ьг~дх1/д^ и с^=дХ1 /д^3. Области С,Р,Г1 предполагаются одно-связными и ограниченными. На рис. 1 показан случай п=3. Если в области С определена криволинейная невырожденная сетка изс, то ее образом в £7 будет невырожденная сетка и>п. Для вывода функционала отображение х(Х) задается через композицию двух отображений: обратного £(Х) и прямого х(£).

Каждое отображение порождает метрику. Соответствующие метрические тензоры суть Ь=ота, д=ЬтЬ, 0=стс, где ат,Ьт,ст-транспонированные матрицы. Поскольку отображение Х(£) за-

сеткой, С - каноническая область с криволинейной сеткой, О, -физическая область.

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

Используется следующая вариационная постановка задачи построения сетки.

Постановка 1. Ставится задача построения гомеоморфного гладкого, С1, отображения х(Х): К"—>КП канонической области С на физическую область при условии, что отображение границ областей задано. Осуществляющие отображение функции х(Х) находятся с помощью минимизации функционала, который есть интеграл от функции, зависящей от инвариантов тензора С'1д.

Для конструирования функционала используется нормированное отношение инвариантов I\=tv(G~1g) к In= det(G_1 д) в соответствующих степенях, для того чтобы получить безразмерную величину, которое интегрируется по n-мерному единичному кубу

о " о

(1)

Функционал (1) рассматривался С.А.Иваненко1 с целью построения сеток. Следует отметить, что идея использовать отношение инвариантов метрики в подынтегральном выражении принадлежит В.Д.Лисейкину18'19. Функционал (1) инвариантен относительно ортогональных преобразований и растяжения координат х и X. Если рассматривать только отображение х(£)20, инварианты мониторной метрики д: I\ — trg, 7n=det<? и подынтегральное выражение в функционале вида Ii/(In)l/n, то получим меру отклонения сетки от конформной19.

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

Постановка 2. Ставится задача построения гомеоморфного гладкого, С1, отображения x(£):R"—>-М" параметрической области V на физическую область при заданном отображении границ областей. Осуществляющие отображение функции х(£) находятся в результате минимизации функционала (1), где Gij(£) - элементы симметричной положительно определенной матрицы, заданной в каждой точке области V.

Здесь также, как и в постановке 1, вводится вторая параметризация, т.е. зависимость Х(£), но для этого нет необходимости рас-

18ЛисейкинВ.Д. О построении структурированных сеток на n-мерных по-верхностях//Ж. вычисл. матем. и матем. физ. 1991. Т. 31. № 11. С. 1670-1683.

19Liseikin V.D. Grid Generation Methods. Springer-Verlag. New York. 1999.

20Для этого в (1) следует использовать евклидову метрику G, т.е. положить Gi,=l и Gi3= 0 при i^j.

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

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

Остается открытым вопрос об условиях, при которых в трехмерном случае гладкое отображение х(£): Р—>íí, задаваемое посредством минимизации функционала V, является гомеоморфизмом. Между тем, известны достаточные условия того, что исследуемое отображение является глобальным гомеоморфизмом, используя лишь информацию локального характера об отображении и его свойствах на границах образа и прообраза21. Эти условия могут быть использованы для обеспечения невырожденности сеток, поскольку в рассматриваемом в диссертации подходе сетка строится посредством реализации гомеоморфного кусочно-гладкого отображения, т.е. склеенного из гладких отображений каждой ячейки сетки из области Р в соответствующую ячейку сетки из области Q. Достаточное условие невырожденности разностной сетки формулируется в разделе 1.4 и состоит в следующем. При заданном отображении границы области дР на dfí,

21 Бобылев H.A., ИваненкоС.А., КазунинА.В. О кусочно-гладких гомео-морфных отображениях ограниченных областей и их приложениях к теории сеток//Ж. вычисл. матем. и матем. физ. 2003. Т. 43. № 6. С. 808-817.

если граница каждого элемента сетки д~Рг (Р* - куб в области V), гомеоморфно отображается на границу элемента сЮ» (0,г -гексаэдр в области и) и для всех элементов якобиан ,/г отображения х$(£): Рг—■положителен, тогда глобальное кусочно-гладкое отображение х(£) -.V—является гомеоморфизмом. Следовательно, необходимо иметь условия невырожденности гексаэд-ральной ячейки.

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

Г=(1 - - £>2 + €2Г3] + (1 - + (1 - *2)п]} +

+№[(1 - ^2)Г6 + £2ГТ] + (1 - ^)К2Г8 + (1 - С2)Г5]}, (2)

где Гг=(хг, у%) ~ координаты вершин ячейки (см. рис. 2а). Невырожденность ячейки означает положительность якобиана отображения (2)

«7=г£1-(г?2хг£з) > 0 (3)

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

вого (б) и второго (в) типа с теми же вершинами.

19

Для рассматриваемого в диссертации алгоритма построения сеток предлагаются следующие условия невырожденности линейчатой ячейки. Линейчатая ячейка заменяется на два 12-гранника первого и второго типа с треугольными гранями (см. рис. 2). Каждый 12-гранник состоит из 5 тетраэдров. Единичный параметрический куб разбивается на 10 базисных тетраэдров подобно разбиению двух 12-гранников (см. рис.26,в): 8 угловых тетраэдров при вершинах куба и 2 внутренних. Трилинейное отображение (2) заменяется на набор линейных отображений гл(£) 10 базисных тетраэдров, на которые разбивается параметрический куб, в соответствующие тетраэдры, из которых составлены два 12-гранника. Условие невырожденности ячейки (3) заменяется на условие невырожденности 10 линейных отображений г'1(£), т.е. невырожденности двух 12-гранников. Это равносильно условию положительности алгебраических объемов 10 тетраэдров, на которые разбиваются 12-гранники. Несмотря на то, что для отдельно взятой ячейки множество невырожденных линейчатых ячеек не совпадает с множеством невырожденных 12-гранников, практика построения сеток рассматриваемым методом показала, что это условие обеспечивает невырожденность гексаэдральной сетки.

Приводится ряд необходимых условий невырожденности гексаэдральной ячейки. Часть из них рассматривалась О.В.Ушаковой3. Предложено необходимое условие НУ4, которое состоит в проверке знака якобиана отображения (2) на отрезках, соединяющих соответствующие точки противоположных граней линейчатой ячейки, разбитых на 10x10 четырехугольников, при том что каждая грань параметрического куба разбивается на 10x10 квадратов. Поскольку вдоль каждого отрезка 3 есть квадратичная функция одного параметра С, то зная значения 1 на концах и посередине отрезка, легко определить принимает ли он отрицательные значения на всем отрезке. Это необходимое условие вместе достаточным условием 2 из работы О.В.Ушаковой3 применяется для проверки невырожденности построенной гексаэдральной сетки.

В разделе 1.6 дискретизация функционала V в каждой гекса-

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

В разделе 1.7 описывается квазиньютоновская процедура минимизации дискретного функционала T>h. В матрице вторых производных удерживаются только диагональные элементы, что значительно сокращает время вычислений. Необходимые для этого расчетные формулы приводятся в разделе 1.8. Минимизация этой функции Vh является значительно более экономичной и эффективной процедурой по сравнению с предложенной ранее С.А.Иваненко процедурой, использующей аппроксимацию функционала на другом множестве 24 базисных тетраэдров2.

В разделе 1.9 показывается, что Vh имеет бесконечный барьер на границе множества невырожденных сеток, состоящих из 12-гранных ячеек первого или второго типа. Смысл барьерного свойства дискретного функционала Vh состоит в следующем. Пусть при перемещении узлов сетки у некоторой 12-гранной ячейки один из 5 тетраэдров, на которые она разбивается, близок к вырождению, т.е. его объем стремится к нулю, оставаясь положительным. В знаменателе функционала V присутствует якобиан отображения J. Для линейного отображения базисного тетраэдра в рассматриваемый тетраэдр из Q, якобиан отображения J равен шести объемам тетраэдра. Поэтому знаменатель одного из слагаемых в дискретном функционале Vh, отвечающего этому тетраэдру, стремится к нулю, а само слагаемое стремится к +оо и, следовательно, Vh тоже. Это является бесконечным барьером на пути вырождения тетраэдра и 12-гранной ячейки при минимизации При этом числитель этого слагаемого не может стремиться к нулю (что сняло бы бесконечный барьер у Vh), поскольку тогда все узлы сетки стянутся в одну точку, а это противоречит условию расстановки узлов по границе области дО,.

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

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

В разделе 1.11 описывается способ сгущения и ортогонолиза-ции сетки к границе ¿Ю с помощью задания управляющей метрики С. Для этого используется свойство универсальности функционала £>, т.е. способность воспроизводить заданную сетку. Маршевым методом задается несколько слоев ячеек сетки около 30, а затем за несколько переходных слоев ячеек сетка трансформируется к квазиравномерной. Важным вопросом является получение первоначальной невырожденной сетки, которая в дальнейшем сглаживается посредством минимизации V11. В областях со сложной геометрией построение начальной невырожденной сетки представляет собой отдельную задачу. Для того чтобы попасть в допустимое множество невырожденных сеток, используется ре-гуляризованный функционал, представленный в разделе 1.12. Этот функционал является аналогом двумерного регуляризован-ного функционала22. Примеры построения сеток приводятся в разделе 1.13. На рис. 3 представлена сетка в межлопаточном канале турбины. Сетка строится со сгущением координатных поверхностей к лопаткам и ортогонализацией координатных линий в окрестности поверхности лопаток.

В разделе 1.14 обсуждается вопрос о том, что функционал V не является единственно возможным. Любая функция <р(Е) (Е -подынтегральное выражение в V), монотонно возрастающая при Е>1, также обладает свойством воспроизведения произвольной заданной невырожденной сетки при соответствующем назначении

22ГаранжаВ.А., КапоринИ.Е. Регуляризация барьерного вариационного метода построения разностных сеток//Ж. вычисл. матем. и матем. физ. 1999. Т. 39. № 9. С. 1489-1503.

D

Рис. 3. Сетка 81x41x31 в межлопаточном канале турбины.

управляющей метрики.

Глава 2 посвящена описанию вариационного метода построения подвижных адаптивных гексаэдральных сеток и некоторых вопросов построения одномерных и двумерных адаптивных сеток. Основные результаты опубликованы в [2, 9, 13] (ЖВМиМФ, 2000,2003,2008), [11] (Int. J. Numer. Meth. Fluids, 2002), [10] (SIAM J. Numer. Anal.,2002), [4] (J. Comput. Phys., 2007), [23] (Препринт ВЦ РАН,2007), [18] (Advances in Grid Generation, Chapt.4,2007).

В разделе 2.1 приводится постановка задачи и вариационный функционал Т> в n-мерном случае, используемый для построения адаптивных подвижных сеток. Пусть задано многообразие М. в пространстве Mn+m переменных х=(ж1,... ,хп, 71,..., /т)=(ж\..., хп+тп). Здесь п - размерность евклидова пространства Мп переменных х=(ж1,..., хп), m - число компонент мониторной вектор-функции f=(/1,... ,/т), по которой проводится адаптация сетки. Каждая компонента мониторной функции зависит от х: fp=fp(x1,... ,хп).

Рассматривается гомеоморфное гладкое, класса С1, отображение х(Х) : Mn—>Rn+rn канонической области С из пространства переменных Х=(Х1,..., Хп) на многообразие М с матрицей Якоби alj=dxl /дХ\ i=1,..., n+m, j=1,..., п. Используется также вспомогательное аналогичное отображение х(£): En—>En+m параметрической области V из пространства переменных £=(£*,... , £п) (как правило, это куб или прямоугольный параллелепипед) на многообразие М. с матрицей Якоби Ьг~дх%/д%>, г=1,..., n-fm, j=l, ...,п, и вспомогательное отображение X(£):Rn—>Rn области V на С с матрицей Якоби сг~дХг/д^, г, j=l,..., п. Если в области С определена криволинейная невырожденная сетка, то ее образом в М. будет также невырожденная сетка, а проекция последней на физическую область fi в пространстве Rn переменных х является адаптивной сеткой.

На рис.4 показан пример для случая п=2, m= 1, т.е. скалярной мониторной функции /. Здесь параметрический квадрат V в плос-

кости отображается на каноническую область С в плоскости X1, X2, и каждая из этих областей отображается на многообразие М, которое суть поверхность в трехмерном пространстве х1, х2, /. Заданная криволинейная сетка в области С с помощью отображения х(Х) трансформируется в сетку в М, а проекция последней на плоскость ж1, х2 является адаптивной сеткой в физической области

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

ся задача построения гомеоморфного гладкого, С1, отображения х(£): Е"—параметрической области V на многообразие М при заданном гомеоморфизме границ дТ-^дМ. Осуществляющие отображение функции х(£) находятся в результате минимизации функционала (1), где - элементы симметричной положи-

тельно определенной матрицы, заданной в каждой точке V.

Идея конструирования адаптивных сеток с помощью записи функционалов на мониторном многообразии была предложена В.Д.Лисейкиным18. Если в (1) использовать в качестве С евклидову метрику, то получим функционал конформности на монитор-ной поверхности, предложенный В.Д.Лисейкиным19, но с другими степенями (подынтегральное выражение суть 1\/{1п)1/п).

В разделе 2.2 приведены одномерный и двумерный функционалы, получаемые непосредственно из функционала общего вида (1), если элементы управляющей метрики йц—бц (5^ - символ Кронекера). Такого вида функционалы были предложены В.Д.Лисейкиным18. Функционал с векторной мониторной функцией использовался при моделировании газодинамических течений на адаптивных сетках в [11,13]. В разделе 2.3 проводится аппроксимация одномерного и двумерного функционалов. Трехточечная модель адаптации представлена в разделе 2.4. Для нелинейного уравнения переноса показано, что при соблюдении некоторых условий значения сеточной мониторной функции в ячейках подвижной сетки остаются неизменными. Тогда задачу построения адаптивной сетки можно рассматривать отдельно от основной задачи, как если бы использовалась аналитическая мониторная функция. Это позволяет исследовать некоторые важные свойства дискретного функционала Т>к, которые сохраняются и в общем случае, например, при расчете газодинамических течений с ударными волнами.

В разделе 2.5 рассматриваются свойства функционала Т>н в одномерном и частном двумерном случае с использованием трехточечной модели адаптации. Показано, что при наличии разрывов у мониторной функции необходимо "замораживать" производные

от этой функции для сохранения бесконечного барьера у одномерного функционала Барьерное свойство позволяет сильно сгущать сетку в окрестности разрыва мониторной функции, что в задачах моделирования газодинамических течений с разрывным решением должно приводить к существенному повышению точности расчетов. В двумерном случае замораживание производных от мониторной функции обуславливает барьерное свойство итерационной процедуры минимизации функционала Т>н, также препятствующее вырождению ячеек сетки. Этот прием использовался С.А.Иваненко2, а теоретический анализ был проведен в [9,10]. Проведенный анализ свойств дискретных функционалов в одномерном и двумерном случаях показывает, что они являются несогласованными между собой, т.е. осуществляемое с их помощью сгущение сетки к разрыву мониторной функции происходит по разному внутри области О и на ее границе дП. Следовательно, возникает задача согласованной расстановки узлов сетки в П и на Она решается в разделе 2.7.

В разделе 2.6 на одномерном примере показано, что при использовании уравнений Эйлера для построения сеток бесконечный барьер на границе множества невырожденных сеток отсутствует, что приводит к вырождению ячеек во время адаптации. В разделе 2.7 описаны 5 способов перестроения узлов на границе области <90 при адаптации. Для согласованной расстановки узлов сетки внутри П и на ее границе дП решается задача условной минимизации дискретного функционала Т>,г при наличии ограничений типа равенств, задающих <9П. Примеры использования различных методов расстановки граничных узлов для аналитических мониторных функций приведены в разделе 2.8. Они иллюстрируют несогласованное и согласованное сгущение узлов внутри области и вдоль ее границы дО. при использовании различных способов расстановки точек на д(2.

В разделе 2.9 дан вывод трехмерного функционала V, записанного на многобразии М в пространстве К3+т. Он имеет следующий вид [2,23]:

раметрического куба V в физическую область Г2. Якобиан этого отображения суть

где х=(х1,х2,х3).

В разделе 2.10 описана дискретизация функционала (4), осуществляемая усреднением аппроксимаций функционала на 10 базисных тетраэдрах для каждой линейчатой ячейки (см. рис.2). Там же приведена процедура минимизации дискретного функционала XРасчетные формулы даны в разделе 2.11. В разделе 2.12 рассмотрен алгоритм расстановки граничных узлов, условная и параметрическая минимизация В разделе 2.13 приведены примеры построения адаптивных сеток. На рис. 5 представлена сетка 81x61x41 в межлопаточном канале турбины. Заданием управляющей метрики С обеспечивается сгущение координатных поверхностей сетки к лопаткам турбины. Проводится адаптация сетки по аналитически заданной мониторной функции /, определяющей форму слоя высоких градиентов.

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

\Zdet5 = х^г(х^хх^з)

Рис. 5. Адаптивная сетка (а); фрагмент граничной поверхности к=41 (б).

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

В главе 3 представлен метод расчета двумерных течений невязкого газа, включая случаи течений с выделением химической энергии, на подвижных сетках. Основные результаты опубликованы в [54] (Препринт ВЦ РАН,1997), [6, 8, 9, 13] (ЖВМ и МФ, 2000,2003,2005), [12] (Comput. Methods in Applied Mech. and Engin.,2000), [11] (Intern. J. for Numer. Meth. in Fluids,2002), [15] (Comm. Math. Sei.,2003), [7] (J. Comput. Phys.,2005), [16] (Comput. Fluid Dynamics J.,2001). Метод включает в себя элементы схемы С.К.Годунова4'5: аппроксимацию уравнений газовой динамики, записанных в виде интегральных законов сохранения, и решение задачи о распаде разрыва для определения потоков через границы ячейки. При этом зона горения газа разрешается сгущением узлов подвижной адаптивной сетки.

В разделе 3.1 приводится система уравнений, описывающая одномерное течение невязкого газа при наличии химической реакции. Обсуждается вопрос постановки граничных условий для различных режимов горения газа. Численная схема расчета на подвижной сетке рассматривается в разделе 3.2. В ней аппроксимируются уравнения движения газа, записанные в виде интегральных законов сохранения. Задача о распаде разрыва на подвижной сетке для течений газа с детонационными волнами рассмотрена в разделе 3.3. Здесь используется метод решения нелинейной системы уравнений для негорючего газа5 в момент времени in+1/2, и решается задача о распаде разрыва для уравнения химической кинетики. В разделе 3.4 приводится условие устойчивости на шаг по времени при расчете на подвижной сетке. Система уравнений для двумерного течения газа при наличии химической реакции выписана в разделе 3.5 в форме законов сохранения в интегральной форме. Разностная схема решения уравнений на подвижной сетке рассматривается в разделе 3.6.

В разделе 3.7 приведены результаты расчетов течения газа без химической кинетики: задачи о распаде разрыва, сверхзвукового течения в плоском канале, обтекания крылового профиля, нестационарного течения в плоском канале со ступенькой, двумерной задачи о взрыве. На примере расчета задачи о сверхзвуковом течении газа в плоском канале показано, что использование адаптивных сеток позволяет сэкономить память компьютера в 25 раз и уменьшает время расчета в 50-60 раз по сравнению с расчетом на квазиравномерной фиксированной сетке. В разделе 3.8 приведены результаты расчетов течения газа с химической кинетикой: одномерной детонации в режиме Чепмена-Жуге, одномерного течения неустойчивой пересжатой волны, неустойчивой пересжатой волны в плоском канале (см. адаптивную сетку на рис.6).

В главе 4 описан алгоритм консервативной интерполяции с одной гексаэдральной сетки на другую. Основные результаты опубликованы в [3] (Матем. Моделирование,2008), [27] (Препринт ВЦ РАН,2006), [18] (Advances in Grid Generation, Chapt. 12, 2007).

В разделе 4.1 дается постановка задачи интерполяции для кусочно-постоянного распределения консервативной функции в ячейках сетки. Она сводится к определению объема фигуры пересечения ячеек старой сетки ш0 и новой и)п. Для определенности рассматривается интерполяция массы т. Гексаэдральная ячейка с линейчатыми гранями заменяется на две 12-гранные с треугольными гранями (см. рис.2). С помощью поверхностного интеграла вычисляется объем 12-гранников и показывается, что полусумма объемов двух 12-гранников равна объему гексаэдральной ячейки. Это равенство между объемами линейчатой ячейки и десяти тетраэдров, составляющих два 12-гранника, было получено О.В.Ушаковой3 непосредственным вычислением объемного интеграла. Построение фигуры пересечения двух линейчатых ячеек заменяется построением фигуры пересечения четырех пар 12-гранников. В каждой паре участвуют по одному 12-граннику от новой и старой гексаэдральных ячеек. Объем фигуры пересечения линейчатых ячеек берется равным среднему от найденных 4

Рис. 6. Адаптивная сетка для расчета неустойчивой пересжатой волны в плоском канале в момент ¿=60.44 (а) и ¿=61.59 (б).

объемов фигур пересечения 12-гранников.

В разделе 4.2 дается общее описание метода. Алгоритм интерполяции состоит из 4 этапов. В разделе 4.3 описывается этап I, построение £оп, линии пересечения поверхностей 12-гранных ячеек С1п новой сетки шп и ячеек Í)G старой сетки ш0. На этом этапе находится ломанная линия С0п, каждое звено которой суть линия пересечения двух треугольных граней ячеек Г1п и íl0. Предполагается, что линия Соп является контуром, без точек ветвления. Может существовать несколько таких линий Сгоп.

В разделе 4.4 описывается этап II. Здесь находится fion, фигура пересечения ячеек й0 и Í2n, которая является симплексом. Этот этап осуществляется одновременно с этапом I. На каждой из рассматриваемых треугольных граней ячеек ft0 (или Qn) строится многоугольник, высекаемый линией Lm и лежащий внутри ячейки Г\п (или íl0). Заключенная внутри него многоугольная плоская поверхность есть часть границы искомой фигуры flon. Поскольку может существовать несколько линий Сгт, то и многоугольников на рассматриваемой грани может быть больше одного (но не более двух). Объединение всех многоугольных поверхностей образует границу Г20П. Отметим, что область fion может состоять из двух симплексов. В разделе 4.5 описывается этап III, вычисление объема и массы фигуры пересечения fion.

В разделе 4.6 описывается этап IV, алгоритм перебора ячеек сетки ш0. Необходимо найти окружение рассматриваемой новой ячейки fín, которое состоит из старых ячеек fio, 9=1,2,.. .,gmai, имеющих в пересечении с ячейкой Í2n непустое множество. Для текущей ячейки fin используется оптимальный алгоритм перебора ячеек сетки ш0 с тем, чтобы проверять только ячейки окружения

ni

В разделе 4.7 рассмотрен вопрос построения слоев фиктивных ячеек. Ячейки и грани старой и>0 и новой шп сеток не должны совпадать. Иначе задача определения линии пересечения совпадающих граней ячеек является неоднозначно определенной, что при вычислениях приводит к делению на ноль. Поэтому невозможно

интерполировать величины с сетки на саму себя. Чтобы избежать ситуаций прерывания следует слегка сдвинуть одну сетку (например, новую ип) относительно другой. Поскольку при сдвиге граничные ячейки сетки шп выйдут за расчетную область, вводится дополнительный внешний слой фиктивных ячеек для старой сетки ш0 по всей границе области с экстраполяцией в них значений плотности23. Такую расширенную сетку будем обозначать ш1хЬ. Теперь ячейки сетки шп будут целиком лежать внутри области, занимаемой ячейками и потери массы после интерполяции не происходит.

В разделе 4.8 проводится анализ ошибки интерполяции. Она складывается из ошибок двух типов. Первая связана с ошибками округлений при вычислениях. При увеличении разрядной сетки компьютера эта ошибка будет стремиться к 0. Вторая, систематическая ошибка, связана с тем, что сетка и>п сдвинута относительно ша. Чтобы понять причину возникновения систематической ошибки рассматривается пример интерполяции с прямоугольной равномерной сетки и>о на эту же сетку сип, сдвинутую относительно старой. Проводится оценка систематической ошибки.

В разделе 4.9 приводятся численные примеры интерполяции. На рис. 7а представлена старая сетка ш0 с числом узлов 29x29x49 в прямоугольной области 0<х,у<1, 0<2<4 (сетка была построена О.В.Ушаковой) с заданным распределением плотности в виде

где 1, pt=100, а=4. На рис. 76 показана новая прямоугольная сетка шп : 31x31x51, на которую осуществляется интерполяция массы.

Проекция координатной поверхности j=1 на плоскость х, z показана на рис.8а. На рис. 86,в представлены изолинии плотности в слое 1 сеток ш0 и ton соответственно. Относительное изменение

23Dukowicz J.K., Padial N.T. REMAP3D: A conservative three-dimensional remapping code. Los Alamos report. 1991.

рь, если z<z(, — 0.43

Pb + (Pt ~ Pb){ 1 ~ exp[-a(z - zb + 0.43)2]}, иначе,

Рис. 7. Криволинейная сетка ша:29х29x49 (а), прямоугольная равномерная сетка шп: 31x31x51 (б).

массы и объема после интерполяции составило

£m = = 2.464-ю-4, sv = — 1-518 • 1СГ6, т0 Vo

где m0, V0 - масса и объем области. Заметим, что ev=em когда р= 1.

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

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

35

■0.5

(а) ' (б) (в)

Рис. 8. Проекция координатной поверхности ^'=1 сетки ш0 на плоскость х, г (а); изолинии р в слое ]=\ сетки и0 (б) и сетки шп (в).

-л о I -I

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

36

областях сложной формы.

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

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

сеток позволяет сэкономить время расчета в 50-60 раз и оперативную память компьютера в 25 раз по сравнению с расчетом на фиксированных сетках.

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

Публикации по теме диссертации

Статьи в рецензируемых журналах, рекомендованных ВАК.

1. АзаренокБ.Н. Вариационный метод построения гексаэдраль-ных сеток с управляющей метрикой//Матем. моделирование. 2008. Т.20. № 9. С. 3-22.

2. АзаренокБ.Н. Вариационный метод построения пространственных адаптивных сеток//Ж. вычисл. матем. и матем. физ. 2008. Т.48. № 5. С. 100-119.

3. АзаренокБ.Н. Об одном методе консервативной интерполяции на гексаздральных сетках//Матем. моделирование. 2008. Т.20. № 2. С. 59-75.

4. Azarenok B.N. A method of constructing adaptive hexahedral moving grids//J. Сотр. Phys. V. 226. Issue 1. 2007. 2007.

P. 1102-1121.

5. AzarenokB.N. A variational hexahedral grid generator with control metric//J. Сотр. Phys. V. 218. Issue 2. 2006. P. 720-747.

6. Азаренок Б.Н. Об одной схеме расчета детонационных волн на подвижных сетках//Ж. вычисл. матем. и матем. физ. 2005. Т.45. № 12. С. 2260-2282.

7. AzarenokB.N., TangT. Second-order Godunov-type scheme for reactive flow calculations on moving meshes//J. Сотр. Phys. V. 206. Issue 1. 2005. P. 48-80.

8. Азаренок Б.Н. Расчет задачи о взрыве на подвижной адаптивной сетке//Ж. вычисл. матем. и матем. физ. 2003. Т. 43. N® 6. С. 856-865.

9. Азаренок Б.Н. О применении вариационного барьерного метода в гиперболических задачах газовой динамики//Ж. вычисл. матем. и матем. физ. 2003. Т. 43. № 7. С. 1072-1096.

10. AzarenokB.N. Variational barrier method of adaptive grid generation in hyperbolic problems of gas dynamics//SIAM J. Numer. Anal. 2002. V. 40. № 2. P. 651-682.

11. IvanenkoS.A., AzarenokB.N. Application of moving adaptive grids for numerical solution of nonstationary problems in gas dynamics//Int. J. Numer. Meth. Fluids. 2002. V. 39. P. 1-22.

12. AzarenokB.N. Realization of a second-order Godunov's method //Comput. Methods in Applied Mechanics and Engineering. 2000. V. 189. № 3. P. 1031-1052.

13. Азаренок Б.Н., Иваненко С.А. О применении адаптивных сеток для численного решения нестационарных задач газовой динамики//Ж. вычисл. матем. и матем. физ. 2000. Т. 40. № 9. С. 1386-1407.

14. Азаренок Б.Н. О решении стационарных задач магнитной гид-родинамики//Матем. моделирование. 1999. Т. 11. № 11.

С. 91-108.

Статьи в рецензируемых журналах.

15. AzarenokB.N., IvanenkoS.A., TangT. Adaptive mesh redistribution method based on Godunov's scheme//Comm. Math. Sci. 2003. V. 1. № 1. P. 152-179.

16. AzarenokB.N., IvanenkoS.A. Application of moving adaptive grids for simulation of supersonic gas flow//Comput. Fluid Dynamics Journ. Japan. 2001. V. 10. № 3. P. 400-404.

Главы в коллективных монографиях.

17. AzarenokB.N., IvanenkoS.A. Grid optimization and adaptation. In: Advances in Grid Generation. (Ushakova O.V. ed.). Nova Science Publishers. New York. 2007. Chapt. 4. P. 85-125.

18. AzarenokB.N. Conservative remapping on hexahedral meshes. In: Advances in Grid Generation. (Ushakova O.V. ed.). Nova Science. New York. 2007. Chapt. 12. P. 337-379.

Прочие публикации.

19. АзаренокБ.Н. К вопросу о построении пространственных адаптивных сеток//Ргос. of International conference "Numerical geometry, grid generation and scientific computing", A.A. Do-rodnicyn Computing Center RAS, Moscow, June 10-13, 2008, P. 67-74.

20. АзаренокБ.Н. Об одном вариационном методе построения структурированных сеток в двумерных областях//Тезисы. VII Международная Конференция по Неравновесным Процессам в Соплах и Струях (NPNJ 2008). 25-31 мая 2008 г. Алушта. М.: Вузовская книга. 2008. С. 36-37.

21. AzarenokB.N. A variational approach to hexahedral mesh generation. In: SoniB., HauserJ., EismanP., and Thompson J., eds. Proceedings of the 10th International Conference on Numerical Grid Generation, Forth, Crete, Greece, 16-20

40

September, 2007. Birmingham. Alabama: International Society of Grid Generation. P. 3-12.

22. АзаренокБ.Н. О построении адаптивных подвижных гекса-эдральных сеток//Тезисы. Материалы XV Международной Конференции по Вычислительной Механике и Современным Прикладным Программным Системам (ВМСППС-2007). 2531 мая 2007 г. Алушта. М.: Вузовская книга. 2007. С. 38-39.

23. АзаренокБ.Н. О построении подвижных адаптивных пространственных сеток. М.: ВЦ РАН. 2007. 50 С.

24. АзаренокБ.Н. О построении пространственных сеток//Тези-сы доклада на XVI Всероссийской конференции "Теоретические основы и конструирование численных алгоритмов для решения задач математической физики с приложением к многопроцессорным системам" посвященная памяти К.И.Бабенко. 4-10 сентября 2006. Дюрсо. С. 3.

25. АзаренокБ.Н. К вопросу о построении гексаэдральных се-ток//Труды Всероссийской Конференции "Численная геометрия, построение расчетных сеток и высокопроизводительные вычисления". М.: ВЦ РАН. 4-7 июля 2006. С. 100-107.

26. АзаренокБ.Н. Об одном вариационном методе построения пространственных сеток. М.: ВЦ РАН. 2006. 51 С.

27. АзаренокБ.Н. Алгоритм консервативной интерполяции на гексаэдральных сетках. М.: ВЦ РАН. 2006. 58 С.

28. AzarenokB.N. First-Order Algorithm of Conservative Interpolation on Hexahedral Meshes//In: PapadopolousP., SoniB., HauserJ., EismanP., ThompsonJ. (eds.) Proceedings of the 9th International Conference on Numerical Grid Generation, San Jose, California, 12-15 June 2005. Birmingham, Alabama: International Society of Grid Generation, pp. 3-12.

29. АзаренокБ.Н. Об одном алгоритме построения регулярных гексаэдральных сеток//Тезисы. Материалы XIV Международной Конференции по Вычислительной Механике и Современным Прикладным Программным Системам (ВМСППС-2005). 25-31 мая 2005. Алушта. М.: Вузовская книга. 2005. С. 31-32.

30. АзаренокБ.Н. Вариационный метод конструирования регулярных гексаэдральных сеток с управлением формы ячеек/ /Тезисы. VIII Международная Конференция "Забабахин-ские Научные Чтения". 5-10 сентября 2005. Снежинск.

С. 216-217.

31. АзаренокБ.Н. О расчете течений газа с детонационными волнами на подвижных адаптивных сетках//Труды Всероссийской конференции "Прикладная геометрия, построение расчетных сеток и высокопроизводительные вычисления", ВЦ РАН. Москва. 28 июня - 1июля 2004. Т. 1. С. 87-96.

32. АзаренокБ.Н. О расчете детонационных волн на адаптивных сетках//Тезисы. 4-я Международная Школа-Семинар "Внут-рикамерные процессы, горение и газовая динамика дисперсных систем", Балтийский гос. университет Военмех. С.-Петербург. 27 июня - 3 июля 2004. Т. 2. С. 150-152.

33. АзаренокБ.Н. Моделирование распространения детонационных волн с использованием адаптивных подвижных сеток/ /Тезисы. Материалы V Международной Конференции по Неравновесным Процессам в Соплах и Струях (NPNJ-2004). 5-10 июля 2004. Самара.' С. 17-18.

34. АзаренокБ.Н. Консервативная интерполяция гидродинамических параметров при расчетах на гексаэдральных сетках/ /Тезисы. Материалы V Международной Конференции по Неравновесным Процессам в Соплах и Струях (NPNJ-2004). 5-10 июля 2004. Самара. С. 18-19.

35. АзаренокБ.Н. О консервативной интерполяции газодинамических параметров на гексаэдральных сетках//Тезисы доклада на XV Всероссийской Конференции "Теоретические основы и конструирование численных алгоритмов для решения задач математической физики с приложением к многопроцессорным системам "посвященная памяти К.И.Бабенко. 8-11 сентября 2004. Дюрсо. С. 4-5.

36. АзаренокБ.Н. Расчет течений с детонационными волнами на подвижных адаптивных сетках//Тезисы. VII Международная Конференция "Забабахинские Научные Чтения". 8-12 сентября 2003. Снежинск. С. 115-116.

37. АзаренокБ.Н. Алгоритм консервативной интерполяции газодинамических полей на гексаэдральных сетках//Тезисы. VII Международная Конференция "Забабахинские Научные Чтения", 8-12 сентября 2003. Снежинск. С. 114-115.

38. AzarenokB.N. Moving adaptive meshes and their application in hyperbolic problems of gas dynamics//Abstract of III Intern. Workshop on Scientif. Comput. and Applications. City University of Hong Kong. Jan. 6-9. 2003.

39. AzarenokB.N. Application of Moving Adaptive Meshes in Hyperbolic Problems of Gas Dynamics//Proceedings of the workshop "Grid Generation: Theory and Applications". Moscow, June 24-28, 2002. Computing Center Russian Academy of Sciences, P. 135-144.

40. AzarenokB.N., IvanenkoS.A. Moving Adaptive Meshes and Godunov's Scheme//Abstracts of 9th International Conference on Hyperbolic Problems, Theory, Numerics, Applications. California Institute of Technology. Pasadena. California. March 25-29. 2002. P. 153-154.

41. AzarenokB.N., TangT. Adaptive mesh redistribution method in hyperbolic problems of gas dynamics//International Conference

43

on Scientific Computing and Partial Differential Equations. Book of Abstracts. Hong Kong Baptist University. December 12-15. 2002. P. 9-10.

42. АзаренокБ.Н. К расчету двумерных течений газа на адаптивных подвижных сетках//Тезисы доклада на IX Всероссийском совещании по проблемам построения сеток для решения задач математической физики. 8-11 сентября 2002. Дюрсо. С. 4-5.

43. AzarenokB.N. Adaptive moving grids in problems of gas dynamics, Grid Generation: New Trends and Applications in Real-World Simulations, S.A.Ivanenko et al. (Eds.)//Proceedings of the minisimposium in Intern. Confer. "Optimization of finite-element approximations, splines and wavelets", St.-Petersburg, June 25-29, 2001. P. 30-44.

44. AzarenokB.N., IvanenkoS.A. Moving grid technology for shock waves simulation//Abstracts of 7th Russian-Japanese Intern. Sympos. on CFD, July 31-Aug. 6, 2000, Russ. Acad, of Science, Mosc. Lomonosov Universty.

45. AzarenokB.N. Adaptive Moving Grids in Supersonic Flow Simulation/Numerical Grid Generation in Computational Field Simulations, Proceedings of the 7th International Conference, September 25-28, 2000, Whistler, British Columbia. Edited by B.K. Soni, J. Haeuser, J.F. Thompson, P. Eiseman. P. 629-638.

46. AzarenokB.N. Application of adaptive moving grids for Simula tion of supersonic gas flow in channel//Abstracts of the 2nd Intern. Confer. Applied Mathematics for Industrial Flows, 12-14 Oct., 2000, Ciocco, Italy. P. 25.

47. AzarenokB.N., IvanenkoS.A. Hypersonic flow simulation based on moving grid technology//Abstracts of the 1st Intern. Conference on CFD, 10-14 July, 2000. Kioto. Japan.

48. AzarenokB.N. Second-order Godunov's method for supersonic problems//Abstracts of the 1st Intern. Conference on CFD, 10-14 July, 2000. Kioto. Japan.

49. АзаренокБ.Н. О вариационном барьерном методе построения сеток при решении гиперболических задач газовой динамика/Тезисы докл. на VIII Всероссийском совещании по проблемам построения сеток для решения задач математической физики. Пущино. 21-28 июля. 2000.

50. AzarenokB.N., IvanenkoS.A. Moving grid technology for finite volume methods in gas dynamics//In: Finite Volumes for Complex Applications II - Problems and Perspectives, R.Vilsmaeier, F.Benkhaldoum and D.Hanel (Eds.), Hermes, 1999, 795-802.

51. AzarenokB.N., IvanenkoS.A. Moving Grid Technology for Shock Waves Simulation/ 8th International Symposium on Computational Fluid Dynamics, Bremen, Germany, 5-10 September, 1999.

52. АзаренокБ.Н. О схеме С.К.Годунова второго порядка аппроксимации/ /Тезисы доклада на XII Всероссийской "Теоретические основы и конструирование численных алгоритмов для решения задач математической физики" памяти К.И.Бабенко, 9-13 сент. 1998, Новороссийск, Дюрсо.

53. AzarenokB.N. Realization of a Second-Order Godunov's Method//Conference on Numerical Methods for Fluid Dynamics, 31 March - 3 April 1998, University of Oxford, UK.

54. АзаренокБ.Н. Об одной реализации схемы С.К. Годунова высокого порядка аппроксимации. М.: ВЦ РАН. 1997. 22 С.

Азаренок Борис Николаевич Вариационные методы построения сеток и приложения в газовой динамике

Подписано в печать 10.11.2008

Формат бумаги 60x84 1/16 Уч.-изд. л. 2. Усл.-печ. л. 2,75 Тираж 100 экз. Заказ 33

Отпечатано на ротапринтах в Вычислительном центре им. A.A. Дородницына Российской академии наук 119991, Москва, ул. Вавилова, 40