автореферат диссертации по информатике, вычислительной технике и управлению, 05.13.18, диссертация на тему:Метод численного исследования обтекания пространственных конфигураций путём решения уравнений Навье-Стокса на основе схем высокого порядка точности

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

004600159

Волков Андрей Викторович

МЕТОД ЧИСЛЕННОГО ИССЛЕДОВАНИЯ ОБТЕКАНИЯ ПРОСТРАНСТВЕННЫХ КОНФИГУРАЦИЙ ПУТЁМ РЕШЕНИЯ УРАВНЕНИЙ НАВЬЕ-СТОКСА НА ОСНОВЕ СХЕМ ВЫСОКОГО ПОРЯДКА ТОЧНОСТИ

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

АВТОРЕФЕРАТ Диссертации на соискание ученой степени доктора физико-математических наук

Москва 2010

1 АПР 2010

004600159

Работа выполнена в Федеральном государственном унитарном предприятии Центральном Аэрогидродинамическом институте им. проф. Н.Е. Жуковского (ФГУП ЦАГИ)

Официальные оппоненты: доктор физико-математических наук,

профессор Толстых Андрей Игоревич, доктор физико-математических наук, в.н.с. Меньшов Игорь Станиславович, доктор технических наук, профессор Крицкий Борис Сергеевич Научный консультант - доктор физико-математических наук,

Чернышев Сергей Леонидович Ведущая организация - ФГУП «Центральный институт авиационного

моторостроения им. П.И.Баранова» г.Москва

Защита диссертации состоится 3 июня в 15.00 на заседании диссертационного совета Д 215.001.01 в Военно-воздушной инженерной академии имени профессора Н.Е.Жуковского по адресу 125190, г.Москва, ул.Планетная, д.З.

С диссертацией можно ознакомиться в библиотеке Военно-воздушной инженерной академии имени профессора Н.Е.Жуковского

Автореферат разослан 1$.02> 2010 г.

Ученый секретарь Диссертационного совета Д 215.001.01

Кандидат физико-математических наук j ц ,{

x/yW-^f A.C. Ненашев

ОБЩАЯ ХАРАКТЕРИСТИКА РАБОТЫ

Актуальность проблемы.

Решение проблем повышения качества и сокращения сроков проектирования новых JIA всё в большей мере связывается с совершенствованием численных методов расчёта их аэродинамических характеристик. Модель течения, основанная на полных уравнениях Навье-Стокса или Рейнольдса, позволяет обеспечить новый уровень проектирования. В последние полтора десятилетия было разработано большое количество подходов к решению этих уравнений. Наиболее удачные численные схемы, реализованные в известных коммерческих вычислительных программах (FLUENT, CFX, STAR-CD, NUMECA), получили широкое распространение и относительно успешно используются при решении многих прикладных задач. В подавляющем большинстве современные программы базируются на методах конечных разностей (МКР) или конечного объёма (МКО). Для решения задач аэродинамического проектирования пространственных конфигураций JIA современные расчётные схемы, базирующиеся на МКР или МКО, требуют использования значительных компьютерных и человеческих ресурсов, что приводит к заметному удорожанию результатов вычислений или, в некоторых случаях, к невозможности проведения численного эксперимента. Точность повсеместно используемых численных схем, как правило, не превышает второй порядок. Известно, что использование численных схем высокого порядка точности (выше второго) заметно снижает потребные вычислительные ресурсы и значительно расширяет круг решаемых задач. Однако, практическая реализация таких высокоточных схем встречает ряд принципиальных трудностей, решению которых посвящена данная диссертация.

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

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

Таким образом, актуальность работы определяется развитием авиационной техники и необходимостью разработки методов, направленных на повышение качества и сокращения сроков проектирования ЛА. Этой цели служит развитие численных схем высокого порядка точности, которые позволяют, с одной стороны, существенно уменьшить необходимые вычислительные ресурсы (память, время) по сравнению с существующими подходами и, с другой стороны, создают возможность построения алгоритмов, основанных на оптимальном автоматическом построении расчётных сеток с анизотропной адаптацией к особенностям течения.

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

Цель диссертационной работы состоит в развитии нового метода численного решения уравнений Навье-Стокса и Рейнольдса на основе схемы высокого порядка точности с использованием неструктурированных анизотропных адаптивных сеток. Особое внимание уделяется практической реализации метода и его сравнению с уже используемыми подходами МКР и МКО.

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

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

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

Метод РМГ высокого порядка точности для решения пространственных задач аэродинамики впервые практически реализован на гексаэдральных неструктурированных сетках.

Выполнено непосредственное сравнение требуемых компьютерных ресурсов при выполнении одних и тех же расчётов методами МКО и РМГ.

Автор защищает следующие результаты:

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

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

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

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

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

Апробация работы. Разработанный метод прошёл тщательное тестирование путём сравнения результатов расчёта с имеющимися экспериментальными данными и/или аналитическими решениями.

Основные результаты диссертации содержатся в 12 статьях, опубликованных в ведущих рецензируемых научных журналах, представленных в перечне ВАК. Результаты докладывались на многочисленных международных и всероссийских научно-технических конференциях, в том числе: на Международном конгрессе по авиационным наукам 1САБ (1996), 3-ей Европейской конференции по механике жидкости ЕШОМЕСН (1997), 6-ой международной конференции по генерации сеток для вычислительной аэродинамики

(1998), 16-м Конгрессе Международной ассоциации математического и компьютерного моделирования IMACS (2000), конференциях Американского Института Авиации и Аэронавтики AIAA (2007, 2009), международной конференции по высокоскоростным течениям WEHSFF (2007), 5-ом Европейском конгрессе по численным методам в прикладной науке ECCOMAS (2008), франко-российском семинаре ONERA-ЦАГИ (2009), международной школе-семинаре «Модели и методы аэродинамики» (2004, 2006) и на 12-ти школах-семинарах ЦАГИ «Аэродинамика летательных аппаратов» (1998 - 2009).

В данную диссертацию включены результаты исследований, поддержанные РФФИ (Проекты №96-01-00209-л, № 96-01-00210-л, № 00-01-00070-а, № 03-01-00236-а, №06-01-00283-а, №09-01-00243-а).

Объем работы. Диссертация состоит из введения, четырёх глав, заключения, списка литературы, включающего 156 наименований, и приложений. Общий объем - 189 страниц.

КРАТКОЕ СОДЕРЖАНИЕ РАБОТЫ

Во введении описано современное состояние рассматриваемой проблемы. Отмечается, что современные подходы к решению задач обтекания конфигураций JIA требуют привлечения значительных вычислительных ресурсов, что препятствует реальному прогрессу в решении многих практически важных задач. Снижение потребных ресурсов является одной из актуальных проблем вычислительной аэродинамики. Пути решения этой проблемы связаны с использованием анизотропных адаптивных сеток и с повышением порядка точности численных схем. Сетки с изотропной адаптацией ячеек не дают существенной экономии узлов. Анизотропные сетки (ячейки в которых вытянуты в заданном направлении) существенно экономичнее, но аппроксимация современньми подходами уравнений законов сохранения на таких сетках весьма проблематична. Известные численные коды, используемые в настоящее время в процессе проектирования, основаны на методе конечного объёма (МКО) второго порядка точности. Хорошие результаты МКО достигаются лишь на изотропных равномерных сетках. Повышение порядка точности МКО приводит к необходимости расширения шаблона аппроксимации и к требованию повышения качества расчётных сеток. Таким образом, одно из решений проблемы эффективности численных схем связано с использованием схем высокого порядка точности, позволяющих выполнять расчёты на анизотропных адаптивных сетках. Схемы с высоким порядком дискретизации исследовались на компактном (Толстых 1990) и расширенном (Visbai & Gaitonde 2002) шаблонах с использованием стандартного метода конечных разностей. Для получения монотонного решения с

высоким порядком аппроксимации в рамках МКР и МКО были разработаны также методы ENO/WENO (Harten et al. 1987; Liu et ah 1994; Jiang & Shu 1996). Все эти подходы имеют солидный теоретически« фундамент, но, к сожалению, применяются в основном на структурированных и неструктурированных сетках с почти равномерными распределением размеров и форм соседних ячеек. Возможность применения таких методов на анизотропных неструктурированных сетках не столь очевидна. Кроме того, в случае расширения шаблона аппроксимации открытыми остаются вопросы распараллеливания алгоритма расчёта и его точности в случае использования многоблочных сеток.

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

Одним из наиболее перспективных подходов к высокоточной аппроксимации на базе МКЭ является метод Галеркина с разрывными базисными функциями (РМГ). В последние годы этот метод вызывает повышенный интерес многих исследователей вследствие его общности, гибкости и надежной теоретической обоснованности. Впервые метод был предложен Reed & Hill (1973) для решения уравнения, описывающего перенос нейтронов, а первый анализ был дан в работе Saint & Raviart (1974). Численное решение 2-D уравнений Эйлера и Навье-Стокса на треугольных неструктурированных сетках этим методом впервые представлено в работах Bassi&Rebay (1977). Наиболее полное теоретическое описание метода с решением 1-D и 2-D модельных задач приведено в работах Cockburn&Shu (1998, 2001). Подробное исследование РМГ можно также найти в работах Lyapunov & Wolkov (2000), Волкова и Ляпунова (2006, 2007, 2009) и в диссертации Ляпунова (2008).

Практическая реализация РМГ для решения нелинейных уравнений Эйлера, Навье-Стокса и Рейнольдса встречается с принципиальными трудностями, по причине которых этот метод ещё не получил широкого распространения. В настоящей диссертационной работе было выделено несколько основных проблем, решение которых и

позволило успешно применить РМГ для расчёта пространственных конфигураций:

1) большое количество арифметических операций на одну степень свободы:

В отличие от МКО, в РМГ необходимо точное интегрирование не только потоков через грани контрольного объёма, но и вычисление объёмных интегралов внутри ячеек. Традиционный подход к интегрированию базируется на использовании квадратурных точек Гаусса. Суммирование значений подынтегральных функций в этих точках с определёнными весами обеспечивает получение точного значения в случае полиномиального представления подынтегральной функции. Чем выше порядок полинома, тем большее число гауссовых точек необходимо использовать. Присутствие объёмных и поверхностных интегралов в рассматриваемом методе приводит к необходимости использования чрезмерно большого количества этих точек на гранях контрольного объёма и внутри него.

2) немонотонность схемы высокого порядка в областях разрывов решений;

Решение задач обтекания может характеризоваться наличием в поле течения областей с большими локальными возмущениями, в частности, скачков уплотнения. Численное решение таких задач с использованием схем повышенного порядка точности приводит к нефизичным осцилляциям решения. В МКО для устранения нефизичных осцилляции используют ограничители градиентов решения, определяемые по значениям в ячейке и её соседях. Такой способ не применим в РМГ для решения задач неявным методом.

3) постановка граничных условий высокого порядка точности;

Решение дифференциальных уравнений с высоким порядком

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

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

4) проблема решения больших нелинейных систем сеточных уравнений;

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

В настоящее время число 3-0 реализаций РМГ для нелинейных уравнений законов сохранения весьма ограничено. Одна из первых успешных реализаций РМГ на неструктурированных гексаэдральных сетках представлена в публикациях автора диссертации.

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

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

В §§1.3 и 1.4 первой Главы рассмотрена общая теория РМГ. Предполагается, что вся расчётная область разбита на контрольные объёмы (элементы, ячейки). Решение в каждой ячейке хранится в примитивных переменных О = (р, и, V, №, /;). Система уравнений

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

г^+у-(г(и(д))-Р1,(и(д),уи))-8 = о, (1)

здесь и =(р, рг/, ру, р\у, рЕ) - вектор консервативных переменных; р -плотность, и, - компоненты вектора скорости течения, р - давление, Е - удельная полная энергия; Р(и) - невязкий, а Р (I), VI))- вязкий потоки; в - источниковый член, возникающий при использовании модели

турбулентности, матрица Г =

есть якобиан преобразования от

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

линейная комбинация которых и определяет решение в ячейке:

к,

(}((, X) - £ и. (г) сру (х) = и; (/) Ф, (х), (2)

где и .(0 - вектор коэффициентов разложения, которые должны быть

определены в процессе решения, Щ - количество базисных функций в ячейке. Величина К)■ связана с максимальной степенью базисного полинома К в пространственном случае через следующее соотношение: Кг (К) = + [)(К + 2){К + 3). При пространственной реализации РМГ

система базисных функций определена в физическом пространстве и имеет вид:

. . (х-хп)а (у-у,,)13 (г-гпУ %{Х) ---—> = 0<а+р+у<Х,

где х0, >'о, 2п - некая внутренняя точка элемента, а величины Ъх, Аг, /?-определяют размер ячейки вдоль соответствующих осей.

Система сеточных уравнений для коэффициентов иД/) из (2)

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

^ = Г11М"1 Ф|. (Р - ^) п + £ УФ/ (Р - ) сЮ + £ ф,- в сЮ

. (3)

¿1

Здесь М - матрица интегралов от произведений различных комбинаций базисных функций, - контрольный объём, £ = ЭГ2 - граница контрольного объёма. Уравнение (3) состоит из объёмных и

поверхностных интегралов по границе ячейки. Значения зависимых переменных терпят разрыв на границе элементов, и ключевую роль здесь играют правила вычисления переменных и потоков на границах. Как и в МКО, в РМГ величина невязкого потока на грани между двумя ячейками определяется в результате решения задачи Римана о распаде произвольного разрыва. Проанализированы три приближённых подхода решения задачи Римана. Это линеаризованные методики Poy(Roe 1981), Ошера (Engquist & Osher 1981) и упрощённый вариант - схема Лакса-Фридрихса (Lax 1954). Каждый из методов имеет свои преимущества и недостатки, а выбор конкретного подхода должен осуществляться в зависимости от решаемой задачи.

Правила аппроксимации градиентов примитивных переменных и вязких потоков в настоящей работе были реализованы в соответствии с рекомендациями работы Cockburn & Shu 1998.

В §1.5 обсуждаются особенности применения модели турбулентности Спаларта - Алмараса (Spalart& Allmaras 1992) для решения задач турбулентного обтекания. Выбор модели турбулентности был обусловлен кругом рассматриваемых задач - задач внешней аэродинамики при больших числах Рейнольдса. Отметим, что с точки зрения исследования эффективности использования схем высокого порядка точности, этот выбор не принципиален, и подобные исследования могли бы быть выполнены с любой другой моделью турбулентности. Добавление дополнительного уравнения, ответственного за модель турбулентности, значительно усложняет процесс получения общего решения, и эта сложность тем заметнее, чем выше предполагаемый порядок точности расчётной схемы. В диссертационной работе описаны некоторые технические детали применения этой модели, позволившие добиться надежной сходимости алгоритма поиска решения.

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

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

Интегрирование внутри элементарных фигур в канонических координатах осуществляется с использованием квадратурных формул Гаусса, которые сводят интегрирование к сумме значений подынтегральной функции в предписанных (гауссовых) точках с определёнными весами. Каждая квадратурная формула для определённой элементарной фигуры является уникальной и может быть найдена в специализированных справочниках (см., например, Strang & Fix, 1973), или в Интернете (например, www.cs.laileuven.ac.be/~nines/research/ecf).

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

JJj£r*-'dQ-c$r\di:=:0, r = {x,y,z}, к = 1,2,3...,

a s

где пх, пу, п2 - компоненты нормали к границе элемента. Данное

выражение является следствием теоремы о дивергенции.

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

заданными в соответствии с квадратурными формулами поверхностного интегрирования. Оптимизационный функционал строился как сумма квадратов разностей аналитического и численного значений интеграла от полинома заданной степени. Найденные решения обращали значение этого функционала в ноль с точностью до МО"22 и менее. Выполненные вычисления позволили сократить количество используемых внутренних точек и уменьшить общее расчётное время на 15 - 20 %.

Другой возможный подход к сокращению количества арифметических операций основан на аналитическом вычислении необходимых интегралов. В случае линейной зависимости потоков от зависимых переменных все интегралы сеточных уравнений (3) являются интегралами от полиномиальных функций, вычисления которых могут быть сделаны до начала основных итераций. Линейная связь потоков и зависимых переменных возникает, например, в задачах акустики, где решаются линеаризированные уравнения Эйлера: Г (11) - ^ фу. Такое

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

даУф,Р(и)сЮ = ]]{Уф,Г Ф;сЮ = Г, Д^ф, <р.<Ю,

11 Я 12

с^ф,Р(и)пс11 = с||ф,Г; <р. пс11 = с||ф(ф) п с1£.

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

м,

функциям: = ^Г,-Ф_,-• Здесь общее количество элементов

м

разложения М/ может быть отлично от количества базисных функций в элементе. Коэффициенты разложения могут быть найдены по правилам нахождения коэффициентов в ряде Тейлора (с последующим переразложением по базисным функциям элемента). Как правило, для обеспечения необходимой точности число элементов разложения М/ должно превышать число базисных функций К/в элементе. Однако число М/ можно уменьшить до К) посредством решения следующей линейной системы уравнений относительно новых коэффициентов разложения Г': /к, ^ /л/, \

2 V / и

Ф^О, ¡ = 1,..., К,

В реализованном РМГ на неструктурированных гексаэдральных сетках применение безквадратурного подхода оказалось эффективным лишь на грубых многосеточных уровнях, где основная цель решателя

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

§1.8 посвящен методам учёта кривизны обтекаемой границы. В настоящее время при решении задач вычислительной аэродинамики широкое распространение получил метод конечного объёма второго порядка точности. Кусочно-линейное представление обтекаемой границы не нарушает эту точность. Достижение порядка точности выше второго при решении уравнений Навье-Стокса и Рейнольдса требует более детального описания границ обтекаемых тел. Порядок полинома базисной функции должен согласовываться с порядком кривой, описывающей обтекаемую границу. Традиционный метод учёта кривизны обтекаемой поверхности базируется на изопараметрической трансформации граничных элементов в элементы простых форм (треугольник, квадрат, куб). Количество дополнительных узлов, лежащих на поверхности, определяет порядок изопараметрического преобразования. В случае «плоской» реализации РМГ преобразования второго или третьего порядков строились на базе одной или, соответственно, двух дополнительных точек на граничных рёбрах. При «пространственной» реализации метода четырёхугольная сторона гексаэдрального граничного элемента дополнялась центральным узлом, что позволило построить трикубичное преобразование гексаэдра в куб. Термин «трикубичное» означает, что вдоль основных координатных осей преобразование описывается полиномом третьей степени. Отметим, что учёт кривизны границы обтекаемого тела потребовал от сеточного генератора представления дополнительной информации о поверхности тела в виде возможности вычисления проекции центральной точки четырёхугольного граничного элемента на реальную обтекаемую поверхность.

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

В §2.1 отмечается, что использование схем порядка точности выше первого может приводить к нефизичным осцилляциям решения вблизи скачков уплотнений и в областях неразрешённых особенностей течений. Термин «неразрешенные особенности течения» означает область расчёта,

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

Нефизичные осцилляции приводят к невозможности получения решения на мелких и адаптированных сетках. Для их устранения используют либо добавление искусственной вязкости, либо применение монотонных схем первого порядка точности. Необходимость использования схем более высокого порядка точности в гладких областях течения обусловило создание так называемых «гибридных» схем или схем переменного порядка точности. Первая гибридная схема была предложена в работе Федоренко (1962), в которой также был описан алгоритм переключения между монотонной схемой и схемой высокого порядка. Другие гибридные схемы с плавным переходом от низкого порядка к высокому были предложены в работе Гольдин, Калиткин, Шишова (1965). Впоследствии, в рамках конечно-разностной идеологии построения численных схем, было разработано большое количество различных гибридных методов (обзор можно найти в книге Куликовский, Погорелов, Семенов, 2001).

В МКО для устранения нефизичных осцилляции в решении используют ограничители градиентов решения (лимитеры): U(t,x) = U0 +l-grad(U)(x-x0), предложенные Колганом (1972), а

позднее Van Leer (1979). Для предотвращения осцилляций решения коэффициент X уменьшается вплоть до нулевого значения. При этом восполнение решения в ячейке становится кусочно-постоянным, а точность схемы падает до первого порядка. Другой подход реализован в схемах ENO/WENO, в которых шаблон аппроксимации в процедуре реконструкции выбирается на основе локальной гладкости численного решения.

В §2.2 обсуждаются проблемы применения традиционных ограничителей решения в методе конечного элемента. Рассматривается

Эи Эи

одномерная схема РМГ для задачи переноса — + — = 0. В качестве

а! дх

набора базисных функций используются ортогональные полиномы Лежандра, определённые в промежутке от -1 до 1: <Pi=l, Ф2 =х,

Ф3 = -ЦЗдг2 -lj. Реконструкция решения в каждом элементе

представляется выражением: U (/, .г) = ф, (jc) + и2ф2 (х) + н3ф3 (х).

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

реконструкцию решения в виде коэффициента X перед базисными функциями высокого порядка: U(/, х) = г/,ф| (л-) +А-(г/2Ф2 (х)+г/3ф3 (х)).

При этом происходит ограничение вклада в решение составляющих высоких порядков, вплоть до использования кусочно-постоянного восполнения решения (при л=0). Впервые в РМГ такой подход был применен в работе Cockburn, Hou, Shu (1990) при использовании явной схемы Рунге-Кутта интегрирования по времени. Главный недостаток такого подхода состоит в невозможности его применения в неявных схемах. Введение лимитеров, принимающих нулевые значения, приводит к вырождению матрицы Якоби для системы сеточных уравнений относительно неизвестных значений ii\, г/2 и аь Действительно, система аппроксимирующих уравнений в ячейке представляется в виде системы нелинейных уравнений:

^,(г/,Дг/2Дг/з) = 0

=> i F2 (г/|, 'kih, Хи^) = 0. Fi{ul ,hi2 ,Хи3) = 0

Здесь Fi=F2=F3:=0 - сеточные уравнения для определения г/ь и2, щ в рассматриваемой ячейке, выписанные в соответствии с алгоритмом РМГ.

ЭF;

Очевидно, что при Х=0 Якобиан системы этих уравнений

|ФГ|+|И,

Э iij

вырождается. Таким образом, устранение осцилляций решения возможно на основе использования других известных подходов, а именно, добавление искусственной вязкости и использование гибридной схемы. Отметим, что метод добавления искусственной вязкости в РМГ опробован рядом исследователей Persson & Peraire (2006); Barter & Darmofal (2007), в то время как разработка и исследование гибридной схемы впервые представлены в диссертации автора.

Описанная в §2.3 гибридная схема в дальнейшем называется РМГ -X схемой. Ниже используется следующее буквенное обозначение с зависимостью от двух аргументов: РМГ (К, X), где К - максимальная степень базисных полиномов. Основная идея схемы состоит в использовании специфических базисных функций, которые получаются путем комбинации с коэффициентом X базисных функций высоких гармоник и кусочно-постоянных базисных функций с коэффициентом (1-Я,). Условно схема может быть представлена следующим образом:

РМГ (К, >.)= X РМГ (АО + (1 - X) РМГ (0). Здесь буквенная последовательность РМГ (К) обозначает классическую схему Галёркина с разрывными базисными функциями, описанную в §1.3. Особенность новой схемы заключается в возможности плавного

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

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

Для реализации РМГ (К, X) схемы каждый сеточный элемент разбивается на множество подэлементов: Е = е|ие2иеД!() и с каждым

подэлементом ассоциируется только одна узловая точка. Все базисные функции определяются в параметрическом пространстве для стандартного треугольника с вершинами Р^.О), Рг(0, 1) и Рз(0,0). Рассматриваются базисные функции «шапочного» типа, когда функция принимает значение единицы в некой узловой точке и ноль - во всех остальных.

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

для РМГ (К, X) схемы: а) - К= 1; б) - К=2

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

элемент (Ns=3) (рис.1,а) и в данном элементе вводятся три линейные базисные функции:

у2=у, V|/3 = 1-jc-JV .

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

= -х + 2х2, у2 =~У + 2У2, V|/3 = 1-Зх-Зу + 2х2 +2у2 + 4ху,

ц/4 =4х-4х2 ~4ху, v|f5 =4ху, \|/6 =4у-4у2 -4ху.

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

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

Итоговые базисные функции ф схемы РМГ (К, X) определяются как взвешенная сумма кусочно-линейных (или квадратичных) и кусочно-постоянных базисных функций: ф ,= X у,- + (1 - Я) ^ . При этом решение в каждом элементе Е представляется как линейная комбинация базисных функций: U(х) = фу (х). Как было описано в §1.3, система сеточных

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

if Ф^жЦ ».('-*)•■<=-

N

I

-Г Уф; (F-Fv)d£2- f ф ,■ SdQ

J£ic, jQe,

= 0.

Свойства РМГ (К, X) схемы рассмотрены в §2.4 диссертации. Здесь на примерах решения задач конвекции и конвекции-диффузии, имеющих гладкое решение, выполнены оценки порядка точности предложенной гибридной схемы. Показано, что сконструирована схема, в которой плавное изменение параметра X приводит не к дискретному, а к плавному

первого порядка

изменению порядка ошибки от монотонной схемы точности при Л=0 до схемы с порядком К+1 при Я=1.

В § 2.5 описаны и проанализированы некоторые известные сенсоры осцилляции решения: СоскЬигп&БЬи (1998), Кп\'ос1опо\'а Н а1. (2004), РегБвоп & Регапе (2006), и предложен альтернативный сенсор областей сингулярности решения.

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

2

течения: Я = уЕ + (у-\)

и ■+• V"

. Коэффициенты а, первого разложения

II(х) = а,ф,- (х) находятся в соответствии со стандартной процедурой

МКЭ: а

Н ф сЮ. Используя значения энтальпии Н' в Аф-

фиксированных точках элемента (его вершины, середины рёбер, центральная точка) можно получить другое разложение Я(х) = Ь/ф;(х) путём решения следующей системы линейных уравнений:

ь/ф /=Я/, /'=1,..., К/.

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

1, если я2

У

Х =

г

1 -+- -

1к

(

I

-Ь,

[0,

Здесь 50 и

к -

Ч

которые, как

,3 = ^,0

если — к - - Л'о + л" если > ^о + к некоторые настраиваемые параметры, показала практика, достаточно консервативны по отношению к выбору решаемой задачи и типу расчетной сетки, и во всех выполненных расчётах использовались следующие параметры: 50=-1.5, к= 1.

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

настроечных параметров, позволяет отличить гладкое решение на адаптивной сетке от разрывного, полученного на равномерной сетке. Таким образом, настоящий алгоритм может быть применён в качестве индикатора зон скачкообразного изменения решения при расчёте РМГ {К, А.) схемой на анизотропных адаптивных сетках.

В §§ 2.6-2.8 приведены примеры использования РМГ {К, X) схемы на последовательности анизотропных адаптивных сеток при решении задач конвекции с разрывом, расчёте обтекания профиля транс- и сверхзвуковым потоком невязкого сжимаемого газа (решение уравнений Эйлера) и расчёте обтекания профиля трансзвуковым турбулентным потоком (решение уравнений Рейнольдса).

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

В §2.8.А представлена серия расчетов для известного тестового случая «Case-6» Cook at. al. (1979) обтекания профиля RAE 2822 при числах М=0.725 и Re=6.5'106. Экспериментальный угол атаки составляет 2.92°, расчетный угол был выбран из условия равенства экспериментального и расчетного коэффициентов подъемной силы Су=0.738. Расчеты производились на последовательности адаптивных сеток с использованием как кусочно-линейной РМГ (1Д), так и кусочно-квадратичной РМГ (2, Я,) реконструкции решения в ячейке. Финальные расчетные сетки, полученные в процессе расчета с РМГ(1Д) и РМГ (2, X), представлены на рис. 2.

Со схемой второго порядка точности (А=1) было выполнено 10 итераций адаптации, в то время как схема третьего порядка точности (Ä=2) приблизилась к такому же количеству степеней свободы за 8 итераций (см. рис. 3). Результаты, полученные с использованием обеих схем, практически совпадают. Результат сравнения расчетного распределения давления с экспериментальными данными представлен на рис. 3.

а)

Рис. 2. Финальная расчетная сетка, полученная в процессе расчета: а) - РМГ (1, А.); б) - РМГ (2, X).

Су

Рнс. 3. Зависимость коэффициента

подъемной силы от количества степеней свободы расчетных схем РМГ(1, X), РМГ (2, X). Представлен расчет расчет турбулентного обтекания профиля ЯЛЕ 2822 при М=0.725,11е=6.510й

Рис. 4. Распределение давления: - экспериментальные результаты а=2.92°, Сс=0.738, Сх=0.0127; - -расчётРМГ(1, X) а=2.485°, С>=0.738,Сл=0.0151.

В §2.8.Б рассмотрено трансзвуковое отрывное обтекание профиля NACA64A010 при М=0.8, Re=2-10\ а=6.2° (Johnson & Bachalo 1980).

В §2.9 приведён пример использования РМГ (К, X) схемы на неструктурированной гексаэдральной сетке для расчета крыла в трансзвуковом турбулентном потоке.

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

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

В §3.1 представлен явный алгоритм поиска решения, традиционно используемый в РМГ (см. Cockburn&Shu 1998). В частности, рассмотрено семейство так называемых оптимальных TVD схем Рунге-Кутта (Shu & Osher 1988), то есть схем с наибольшим допустимым числом Куранта, обеспечивающим при этом сохранение по времени полной вариации решения. В качестве примера использования явной TVD схемы интегрирования уравнений РМГ аппроксимации приводится решение задачи о распространении цилиндрической акустической волны малой амплитуды. Оценивается поле давления к моменту времени, когда акустическая волна подошла к границам расчётной области. Решение этой задачи требует одновременного использования схемы высокого порядка, как по времени, так и по пространству. Задача имеет аналитическое решение, что даёт возможность оценить реальный порядок точности разработанной численной схемы. Для реконструкции решения в ячейке использованы полиномиальные базисные функции четвёртого порядка (К-А). Интегрирование по времени выполнено шести шаговой TVD схемой Рунге-Кутта пятого порядка точности. Результаты, полученные на вложенной последовательности сеток, демонстрируют сходимость расчётного поля давления к аналитическому полю с порядком точности, превышающим пятый.

В §3.2 рассмотрен неявный алгоритм поиска решений, допускающий использование бесконечно больших чисел Куранта. Отмечается, что неявный метод решения систем сеточных уравнений весьма эффективен, однако требует грандиозных затрат оперативной памяти компьютера, требующейся для хранения матрицы невязки Якоби. В диссертационной работе описываются некоторые подходы, позволяющие снизить использование требуемых вычислительных ресурсов. Неявные методы решения систем сеточных уравнений, полученных в результате аппроксимации методами МКР и МКО, получили широкое распространение (например, Venkatakrishnan 1988; Егоров и Зайцев 1991). Эффективность неявного метода в рамках РМГ была апробирована на решении задачи расчёта многоэлементного профиля с элементами взлётно-посадочной механизации крыла на последовательности адаптивных сеток (Волков и Ляпунов 2006, 2007).

§3.3 посвящен описанию наиболее эффективного метода решения -многосеточного метода. В настоящее время в промышленных программах решения уравнений Навье-Стокса и Рейнольдса с целью ускорения итерационного процесса сходимости к стационарному решению используют многосеточный метод, предложенный Федоренко в 1961. Суть метода состоит в ускорении передачи информации между далеко расположенными друг от друга ячейками сетки. Это осуществляется через укрупнение сеточных ячеек путём их агломераций. В РМГ такой подход малоэффективен, так как изначально предполагается использование грубых сеток с большим количеством переменных в каждой ячейке. В последнее время философия построения многосеточного метода была пересмотрена и расширена (R6nquist& Patera 1987). В работах Helenbrook et al. (2003), Жуков и др. (2004) метод был адаптирован к использованию в конечно-элементном методе аппроксимации течений. По аналогии с решением задачи на последовательности измельчённых сеток, в РМГ решается задача на фиксированной сетке, но с использованием различного набора базисных функций. В настоящей работе полиномиальный многосеточный метод впервые применён к решению пространственных уравнений Навье-Стокса в рамках РМГ. Предложены новые операторы перевода решения между различными уровнями многосеточного метода (операторы интерполяции и сборки), а также локально-неявный сглаживатель ошибки решения.

Предложен следующий алгоритм решателя, подробно описанный в §3.3.1. Определяется последовательность уровней многосеточного метода l = \,...,¿, такая, что уровень £+1 в сравнении с уровнем I имеет либо большее количество ячеек, либо большее количество базисных функций (степеней свободы). При традиционном сеточном измельчении (h-измельчение) используется последовательность вложенных гексаэдральных сеток, генерируемых алгоритмом агломерации соседних ячеек. При измельчении по степеням свободы (р-измельчение) используются наборы из базисных функций, описанных в таблице 1.

Таблица 1. Набор базисных функций в элементе

А=() К= 1 К= 2 К= 3

Базис набора Базис наборов Базис наборов

К= 0 К= 0, 1 К= 0, 1,2

дополняется дополняется дополняется

функциями: функциями: функциями:

1 x,y,z 2 Ъ х-,у г, ху, XZ, yz, 2 2 2 xy,xz,yx, y2z, zzx, z2y, xyz

Далее описание полиномиального алгоритма полностью совпадает с описанием традиционного многосеточного алгоритма по схеме полной аппроксимации (FAS), представленной в работе Brandt (1980).

На самом верхнем (мелком) многосеточном уровне L решается

Да,

система уравнений Навье-Стокса или Рейнольдса: М—-+ R, (q^J = 0.

А t

Здесь q - искомый вектор решения, R(q) - вектор-функция невязки решаемых уравнений или результат пространственной аппроксимации законов сохранения без учёта производной решения по времени, М -блочно-диагональная матрица, состоящая из интегралов от произведений различных базисных функций элемента. Решение данной системы уравнений может быть выполнено явным или неявным способами, описанными выше. Эти способы - итерационные, и применение одной или нескольких итераций уменьшает лишь высокочастотную ошибку решения, но не позволяет решить задачу полностью. Поэтому данную операцию принято называть сглаживанием, а соответствующий алгоритм - сглаживателем. Применение сглаживателя позволяет получить коррекцию решения на верхнем слое L: q£ = q; + AqL.

На более грубом (расположенном ниже) многосеточном уровне I = L -1 к решаемой системе уравнений добавляется правая часть:

M1T + R'(4'>F" ГДе F' =R'(I^l4(+,]) + i(it,[F£+1-R(+1(q(+1)].

Здесь [•], I'+1 [■] - операторы переноса решения и невязки решения с

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

Я*"""' = q, + qCulT, qCOT={lL[q(4-ir'[qt]]j-

Здесь фигурные скобки означают осреднение значения в рассматриваемой ячейке и её соседях.

Как было указано выше, многосеточные уровни I и 1-\ могут различаться либо разным количеством ячеек в сетке при одинаковом количестве базисных функций в реконструированном решении, либо разным порядком полиномиальной реконструкции, но на фиксированной сетке. В зависимости от этих обстоятельств варьируется и вид операторов переноса решения и невязки решения, подробно описанных в §§ 3.3.1 и 3.3.2.

В случае перевода решения между уровнями с разным количеством сеточных ячеек используются традиционные операторы интерполяции и сборки. Для перевода решения между уровнями с разным количеством базисных функций в реконструируемом решении были разработаны новые операторы интерполяции и сборки 1'+|:

С к ]=[ЧрГ Ммч< > С к+, ]=|Х,ЧГ ,

1г+.[»г+,] = мчр[ми,]"

Здесь М - массовые матрицы, каждый член которых определяется в соответствии со следующими правилами:

= /ф^'ф^п , м;-; = , м^ = |ф1'-,ф;сю , м^ = {Ф;Ф;-мп.

Ключевым алгоритмом любого многосеточного метода является так называемый сглаживатель - алгоритм подавления низкочастотной составляющей ошибки решения.

Сглаживатель на основе явного метода является наиболее экономичным алгоритмом. Скорость сходимости итерационного процесса может быть увеличена за счёт применения неявного метода на грубых многосеточных уровнях. Использование неявного сглаживателя на самом мелком уровне даёт весьма заметное ускорение сходимости, однако приводит к значительному перерасходу оперативной памяти. С целью экономии памяти в настоящей работе (§3.3.4) было предложено на верхнем многосеточном уровне использовать упрощенный неявный метод. Его суть состоит в том, что в процессе поиска решения в матрице Якоби учитываются только диагональные блоки, отражающие взаимное влияние переменных внутри только одного элемента. При этом не учитывается влияние соседних элементов. Для конечно-элементного метода высокого порядка точности такой подход наиболее эффективен, так как решение формируется во многом за счет большого количества переменных внутри ячеек. Поэтому в настоящей работе этому подходу было дано название «локально-неявный» метод, которое, по всей видимости, точнее передает суть этого численного алгоритма.

В § 3.3.5 приведён пример оценки эффективности полиномиального многосеточного метода в задаче об обтекании сферы невязким потоком при числе М=0.15. Семейство неструктурированных гексаэдрапьных расчётных сеток было построено с использованием программы НЕХРЯЕВБ (компании КиМЕСА) (см. рис. 5). Алгоритм агломерации позволил построить систему из трёх вложенных сеток, решения на которых получалось с использованием наименьшего, кусочно-постоянного базиса. Далее уровни многосеточного алгоритма измельчались путём обогащения базиса, как это было описано ранее. Таким образом, расчёты РМГ с К=3 (кусочно-кубический базис)

выполнялись с использованием шестиуровневого алгоритма. Типичная история сходимости (зависимость L2 нормы невязки решения от времени CPU) представлена на рис. 6 для случая кусочно-квадратичной аппроксимации решения. Рисунок 6 демонстрирует эффект ускорения решения при использовании многосеточного алгоритма. Отметим, что без применения многосеточного алгоритма при кусочно-кубичном базисе получение решения за разумное время становится невозможным.

Рис. 5. Фрагмент гексаэдральной сетки около сферы

ихКПЧ ч. и РМГ К=2)

- 5-т* шагов ясхамаРу ««-Кутта

: ЪЛИНОМИ 1ЬНЫЙ MHO метод Время ■ 1 . . эоеточиый 3PU(cat)

Рис. 6. История сходимости. Эффект ускорения решения при использовании многосеточного алгоритма.

В Главе 4 проведены численные исследования точности новой схемы аппроксимации на последовательности вложенных сеток. Дано сравнение результатов расчёта и потребных вычислительных ресурсов по предложенной схеме высокого порядка точности с теоретическими и (или) экспериментальными результатами и результатами, полученными с использованием стандартной схемы МКО второго порядка точности.

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

Известно, что наличие скачка уплотнения ограничивает порядок точности численной схемы до первого (см., например, Иванов & Крайко 1978). При этом возникает очевидный вопрос о целесообразности применения схем высокого порядка для решения задач с разрывом. Однако возможность адаптивного изменения размеров ячеек в поле течения реабилитирует необходимость развития высокоточных подходов.

Использование подобных схем возможно в области течения не примыкающей к скачку уплотнения (в «гладкой» области). Очевидно, что в уравнениях газовой динамики ошибка, зарождающаяся на скачке уплотнения, добавляется к ошибке в гладкой области течения:

|и(х)-и;(фо(ь:)+о(ь:)'.

Здесь и(х) - аналитическое решение, и^ (х) - численное решение, полученное на сетке п из вложенной последовательности сеток, Ь™ -размер ячеек в области скачка уплотнения, Ь" - размер ячеек в гладкой области течения, р - порядок точности численной схемы. Используя измельчённую сетку в окрестности скачка уплотнения, величину ошибки можно сделать сравнимой с ошибкой в гладкой области течения, где используется схема с высоким порядком точности, но на сетке с крупными ячейками. Такой подход приводит к существенной экономии количества расчётных узлов и, как следствие, к уменьшению компьютерных ресурсов.

В § 4.2 приводится оценка фактического порядка точности РМГ в задаче о ламинарном обтекании кругового цилиндра при М=0.15 и 11е=40. Рассмотрены две вложенных последовательности из 4 сеток. Первая последовательность состояла из сеток с ячейками треугольной формы. Каждая последующая сетка получена из предыдущей посредством удвоения количества узлов в обоих направлениях. Самая крупная сетка содержала 10 узлов на поверхности цилиндра и 5 узлов в направлении от обтекаемой поверхности к внешней границе течения, радиус которой во всех случаях составлял 300 диаметров цилиндра. Другая последовательность сеток состояла из ячеек сложной формы, которые образовывались путем искусственной агломерации некоторых соседних треугольных элементов. Фрагменты таких сеток приведены на рис.7.

Предельное значение сопротивления, экстраполированное на бесконечно мелкую сетку, составило 1.5041881. Оценки порядков ошибки выполнены при использовании базисных функций третьего РМГ(АГ=3) и четвёртого РМГ(А=4) порядков. Величины сопротивления Сх и порядка ошибки р на последовательности из сеток треугольных ячеек и ячеек сложной формы отражены в таблице 2. Из представленных результатов видно, что, во-первых, несмотря на большие значения ошибок, полученных на сетках с агломерированными ячейками, порядок ошибки близок к порядку, полученному на треугольных сетках. Кроме того, для всех случаев полученные значения близки к теоретическому результату: К+1.

около кругового цилиндра, а) сетка №-1: 10x5, б) сетка №-4: 80x40

Таблица 2. Порядки ошибок в задаче о ламинарном обтекании цилиндра

Треугольные ячейки Ячейки сложной формы

К=3 К=4 А=3 К=4

Сх Р Сх Р Сх Р Сх Р

10x5 1.5344122 1.4905387 2.1969095 1.9978443

20x10 1.5059944 4.06 1.5045982 5.06 1.5371651 3.75 1.5152094 4.68

40x20 1.5042885 4.17 1.5041962 5.66 1.5053013 4.69 1.5045512 4.73

80x40 1.5041953 3.80 1.5041883 5.40 1.5042349 4.4 1.5041995 4.80

Аналогичные результаты получены и в задаче о турбулентном обтекании пластины при числах М=0.1 и 11е=Ы06. Так же как и в предыдущем случае, порядок точности оценивался по величине сопротивления. С целью исключения сингулярной особенности течения в передней части пластины сопротивление вычислялось на серединном участке пластины. При использовании кусочно-линейной и квадратичной реконструкций решения в ячейках в данном численном эксперименте также удалось получить оптимальный порядок точности 0[Ь(Л+|)].

В §4.5 описана серия расчетов обтекания сферы дозвуковым потоком газа при М=0.15, которая была выполнена на последовательности вложенных сеток с использованием кусочно-линейного, квадратичного и кубичного полиномиального базисов. Все расчёты с высоким порядком точности (К > 2) выполнялись с учётом кривизны граней гексаэдров, непосредственно примыкающих к сфере. На рис. 8 представлена зависимость величины коэффициента сопротивления от числа переменных задачи. Известно, что в невязком течении величина сопротивления должна быть равна нулю, и, таким образом, величина отличия расчётного коэффициента сопротивления от нуля косвенно характеризует качество расчётной схемы. Анализ результатов показал,

что скорость стремления к нулю коэффициента сопротивления по мере сгущения сеток выше, чем 2К+1, что, возможно, говорит о достижении суперсходимости (см. Barth 2004) конечно-элементного метода. На самой густой из использованных сеток, содержащей 8 148 ячеек, при применении кусочно-кубичного базиса для реконструкции решения коэффициент сопротивления четверти сферы составил величину 3.2x10"6. Получить такое малое значение коэффициента сопротивления методом конечного объёма чрезвычайно затруднительно. Например, использование кода FLUENT, даже на сетке в 240 ООО узлов, не позволило получить сопротивление меньше, чем 2.1-10"4.

|Сх|

' л РМГ(К=1 рмг (К=2) )

N. РМГ (f =3)

J _. I

5-104 105 .. 2-10*

N

Рис. 8. Зависимость величины коэффициента сопротивления сферы от числа переменных в РМГ для кусочно-линейного, квадратичного и кубичного базисов

В §4.4 представлено сравнение расчётов ламинарного обтекания пластины (М=0.35, 11е=76000) методами РМГ и МКО в условиях эквивалентного количества степеней свободы. Расчёты выполнялись на гексаэдральных сетках разной густоты, имеющих только одну ячейку в г-направлении. Некоторые параметры сеток, представленные в таблице 3, показывают: (1) - общее количество ячеек в сетке, (2) - расстояние от поверхности пластины до первого слоя ячеек, (3) - коэффициент увеличения расстояния между последовательными рядами сеточных точек и общее количество узлов (4) - в х и (5) - в у-направлениях.

Таблица 3. Параметры расчётных

Сетка (1) (2) (3) (4) (5)

№1 6739 1-Ю'5 1.02 25 128

№2 281 4-10"5 2.50 5 32

сеток

Отметим, что количество элементов в сетке №1 более чем в 20 раз превышает количество элементов в сетке №2. Такой выбор обеспечивает возможность проводить сравнения численных схем в условиях близкого количества степеней свободы. На рис. 9 представлены расчётные профили скорости и даны сравнения с аналитическим решением Блазиуса. Расчётные профили изображены комбинацией линий и маркеров. Каждая линия построена через множество точек в точном соответствии с формулой для реконструкции решения в ячейке. Маркеры же на каждой линии расположены приблизительно в серединах ячеек. Таким образом, количество маркеров на линии отражает густоту сетки.

N =6739 N =281x20=5620 02

Рис. 9. Сравнение и и V компонент профилей скорости с аналитическим решением Блазиуса (зелёная линия). Частота маркеров отражает густоту сетки

Продольная ¿/-компонента скорости, показанная на рис. 9.а), хорошо согласуется с аналитическим решением для обеих схем. Поперечная F-компонента оказывается более чувствительной к качеству сетки (рис. 9.6). Результаты, полученные с РМГ К=Ъ на сетке №2, лучше согласуются с аналитическим решением, чем результаты, полученные традиционным МКО на самой мелкой сетке №1, несмотря на эквивалентное количество задействованных степеней свободы. Отметим, что осциллирующее поведение расчётной кривой, полученной МКО, является результатом использования недостаточно мелкой сетки.

Сравнение расчётов трёхмерного ламинарного течения в изогнутой трубке методами РМГ и МКО в условиях эквивалентного количества степеней свободы представлено в §4.6. Пространственное ламинарное течение воды внутри изогнутой под прямым углом трубки постоянного сечения было экспериментально исследовано Enayet et al. (1982) с помощью лазеро-доплеровского измерителя скорости. Результаты расчёта и эксперимента сравнивались при числе Рейнольдса Re=500, соответствующему ламинарному течению. Скорость потока на

входе в трубу - 10 м/с. Установившееся течение в трубке определяется балансом вязких и невязких сил и характеризуется наличием пары вихрей, закрученных в противоположные стороны, образующихся вниз по потоку после прохождения трубочного изгиба. Имеются экспериментально измеренные профили скорости в пяти сечениях вдоль трубки: сечение 1 - 0.58 диаметра трубки вверх по потоку от изгиба, сечение 2 - 30° по изгибу, сечение 3 - 60° по изгибу, сечение 4 - 75° по изгибу, сечение 5 - один диаметр вниз по потоку от изгиба. Для численных исследований построено три неструктурированных гексаэдральных сетки. Грубая сетка №1 содержит 2 500 ячеек. При использовании РМГ К=3 на такой сетке задействованное количество степеней свободы равно 50 000. С целью сравнений результатов РМГ и МКО построена также сетка №2 с количеством ячеек 62 000 и сетка №3, содержащая 140 000 ячеек. Некоторые фрагменты сеток представлены на рис. 10.

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

На рис. 11 профили скорости вдоль сечений трубки, полученные с использованием схемы РМГ К=3 на сетке №1, сравниваются с аналогичными результатами, полученными МКО на сетке №3. Отметим, что результаты, полученные МКО на сетке №2 с числом элементов, обеспечивающим равное количество степеней свободы с РМГ, весьма далеки от экспериментальных данных. На рис. 11 экспериментальные результаты представлены кружками. Линия с квадратными маркерами демонстрирует результаты РМГ(А-З), в то время как пунктирная линия с треугольными маркерами показывает результаты МКО на сетке №3.

о о о Эксперимент л.....л-л Сетка 139,902; МКО ■—■—■ Сетка 2,502; РМГ К=3

Рис. 11. Профиль скорости в 4-х сечениях трубки. Сопоставление эксперимента и численных результатов, полученных с РМГ К=3 на сетке №1 и МКО на сетке №3

Густота маркеров отражает густоту сетки, так как на одну ячейку приходится один маркер. Представленное сравнение демонстрирует эквивалентную точность расчётов по МКО и РМГ, а также их хорошее согласование с экспериментальными данными. Рассчитанный профиль скорости в сечении 1, как правило, во всех случаях находится в удовлетворительном согласовании с экспериментальными данными и поэтому здесь не представлен.

Таблица 4 представляет сравнение вычислительных затрат. Использование схемы РМГ 4-го порядка точности на сетке №1 требует почти в 3 раза большего времени по сравнению со схемой МКО на сетке №2, обеспечивающей использование эквивалентного количества степеней свободы. Однако, принимая во внимание большую точность результатов РМГ схемы, целесообразно сравнить общие времена расчётных схем, требуемые для получения результатов с одинаковой точностью, которая, как видно из рис.11, обеспечивается МКО на сетке №3.

Отметим, что в расчёте на 50 итераций многосеточного алгоритма обе схемы затрачивают примерно одинаковое время. Расчёт новым методом на грубой сетке позволяет использовать почти в два раза больший шаг по времени, а, следовательно, до сходимости одинакового уровня требуется в два раза меньшее количество итераций. Таким образом, использование РМГ схемы для решения этой задачи позволило сократить не только объём оперативной памяти, но и общее расчётное время.

Таблица 4. Сравнение вычислительных затрат различных схем аппроксимации

Сет- Кол-во Время на Общее

ка Кол-во степ. Отнош. 50 Мв кол-во Отнош

№ ячеек своб. памяти итераций (секунды) итерации врем.

РМГ 1 2502 50040 0.36 1185 4375 0.47

к=з

МКО 2 62689 62689 0.45 496 7880 0.86

МКО 3 139902 139902 1.00 1115 9155 1.00

В §4.7 оцениваются возможности схемы РМГ при решении задачи о распространении сферической акустической волны. Известно, что решение задач аэроакустики требует применения численных схем высокого порядка точности, обладающих минимальными дисперсионными и диффузионными свойствами. Исследуемая схема, как оказалось, обладает такими свойствами и может быть использована в практических приложениях. Свойства схемы были изучены на примере классической задачи о распространении акустического импульса, имеющего сравнительно малое начальное возмущение. Задача рассматривается в кубической расчетной области. Решение сеточных уравнений осуществляется методом Рунге-Кутта 5-го порядка точности с глобальным шагом по времени. Представленные результаты соответствуют времени, когда передний фронт возмущения еще не дошел до границ расчетной области. На рис. 12 приводится сравнение результатов для РМГ К=\, 2, 3 схем и традиционной схемы МКО 2-го попядка точности (ЖГМЕСА), полученных на последовательности равномерных структурированных сеток. Количество ячеек в сетке определяется как произведение числа разбиений расчётной области по каждому из направлений Л^хЛг,хЛг-. На представленных рисунках в подписи к соответствующей кривой отражается не только размерность использованной сетки, но и количество степеней свободы, задействованных при выполнении соответствующего расчёта. В 4-х символьном обозначении: NЛxN,xN.^x.m число т указывает количество использованных переменных в ячейке. Видно, что расчёты,

101301.5

Р

101301.0

101300.5

101300.0

101299.5

101299.0

ооо Точное решение

— МКО; 51=86x86x86x1=636056 ■——■ РМГ К=1; 31=54x54x54x4=629856 • - • РМГК=2; =40x40x40x10=640 ООО . • • РМГК=3; 51=32x32x32x20=655360

Рис. 12. Сравнение точного и численного решений для акустического импульса. Эквивалентное число степеней свободы

со

•

\А

•

i_i i i i i i i i i i. I i i i i

представленные на рис. 12, выполнялись на индивидуальных сетках, подобранных таким образом, чтобы задействованное количество степеней свободы было эквивалентным. Линии снабжены маркерами, из расчёта один маркер на ячейку. Поэтому плотность маркеров также отражает густоту сеток.

Сравнение компьютерных ресурсов для этого случая приведено в таблице 5. Схема 4-го порядка точности обеспечила хороший результат в 5.6 раза быстрее, использовав при этом в 1.4 раза меньшую память. Отметим, что получить такую же точность результатов со схемой РМГ А=1 и, тем более, с использованием традиционной конечно-объемной схемы не представилось возможным из-за чрезмерно большой потребной компьютерной памяти.

Таблица 5. Сравнение компьютерных ресурсов в случае близкой

точности результатов для схем РМГ К=2 и К~3

Сетка Кол-во степеней свободы Исп. память (КЬ) Отн. памяти Время CPU Отн. врем.

РМГ К= 2 60x60x60 2 160 000 743 1.4 13ч 48мин 5.6

РМГ А'=3 32x32x32 655 360 528 1.0 2ч 28мин 1.0

В §4.8 выполнено сравнение схем МКО второго порядка точности и РМГ (К= 1) в задаче о турбулентном обтекании изолированного крыла дозвуковым потоком вязкого газа (М=0.4, а=0.6°, Re=7.17xl06). В качестве тестового случая был выбран известный эксперимент Müller at al. (1996) с крылом LANN. Расчёты выполнены на неструктурированных гексаэдральных сетках, обеспечивающих использование эквивалентного количества переменных. Грубая сетка с числом элементов 190 213 использована для РМГ, в то время как расчёты МКО выполнялись на сетке с числом элементов 625 076. Сравнения результатов расчёта распределённых характеристик обтекания демонстрируют, что оба подхода обеспечивают примерно одинаковые результаты, несмотря на использование заметно более грубой сетки в РМГ.

Сравнение схем МКО и РМГ при решении задачи об обтекании конфигурации «крыло + фюзеляж» представлено в §4.9. Самолётная конфигурация «крыло+фюзеляж» DLR-F4, выбранная в качестве тестовой модели конференцией AIAA CFD Drag Prediction Workshop (DPW), была рассчитана РМГ(А=1) на сетке, содержащей 230 тысяч узлов. Объём использованной сетки был ограничен возможностями персонального компьютера с 2-мя гигабайтами оперативной памяти и 32-битной операционной системой. Однако новый метод позволил получить результаты, удовлетворительно согласующиеся с экспериментальными данными и с расчётами, выполненными по известным коммерческим кодам МКО на более густых расчётных сетках. Отметим, что по классификации DPW сетка с объёмом ячеек 1.5-н 2 млн. считается грубой, но даже для расчёта на таких сетках требуется вычислительная система с большей оперативной памятью и 64 битной операционной системой, позволяющей работать с памятью, превышающей ~ 3 Гб.

Расчётная сетка №1 для конфигурации DLR-F4, сгенерированная кодом HEXPRESS™ (пакета NUMECA), содержит 229 739 ячеек и 261 752 узлов. Фрагмент поверхностной сетки и общий вид изображены на рис. 13.

Рис. 13. Фрагмент поверхностной сетки: а) конфигурация «крыло+фюзеляж» включая плоскость симметрии, б) внешняя граница расчётной области

Средняя аэродинамическая хорда крыла рассчитываемой конфигурации составляет 0.1412 м, длина фюзеляжа - 1.92 м, а полуразмах - 0.5856 м. Внешняя граница расчётной области представляет собой полуцилиндр радиусом 5 м и высотой 20 м. Ячейки, прилегающие к обтекаемой границе, имеют дополнительное разбиение, состоящее из 4-х слоев. Расстояние до первого слоя узлов составляет 2-10"5м, а толщина последующих слоев увеличивается с коэффициентом геометрической прогрессии 1.4.

Расчёты выполнены при числах М=0.75 и Re=3-106. Начальный коэффициент турбулентной вязкости на бесконечности в модели турбулентности Спаларта-Алмараса был равен 5. Расчёты по МКО были выполнены кодом Fine/Hexa™ (пакета NUMECA) на сетке №2, содержащей более 1 млн. ячеек. Количество слоев в приграничных ячейках использованной сетки было равным 12, а коэффициент геометрической прогрессии 1.22. На рис. 14 и в таблице 6 представлены сравнения результатов расчётов выполненных РМГ на сетке №1 и МКО на сетке №2. Видно, что результаты по схеме РМГ находятся в лучшем согласовании с экспериментом DLR.

Сопоставление расчётных (РМГ Ä=l) и экспериментальных интегральных характеристик обтекания по углам атаки представлено на рис. 15. Сплошная красная линия отражает результаты расчёта, а маркерами изображены экспериментальные значения (Redeker, 1994), полученные в различных АДТ (DLR, ONERA, DRA). Пунктирные линии ограничивают область результатов получаемых различными современными МКО кодами, представленными в материалах конференции DPW (Lee-Rausch et al. 2003).

Отклонение расчётной поляры от экспериментальной составляет величину 2-10"3 (20 каунтов). Такое различие может быть объяснено несовершенством используемой модели турбулентности и недостаточным количеством задействованных степеней свободы задачи. Тем не менее, это отклонение минимально среди других результатов, полученных МКО на сетках с большим количеством ячеек и с использованием больших компьютерных ресурсов.

Таблица 6. Сравнение интегральных характеристик расчёта при а=0

Су Сх

МКО 0.5728 0.0325

РМГ 0.4773 О.ОЗОЗ

Эксперимент DLR 0.4812 0.0278

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

Рис. 15. Интегральные характеристики обтекания конфигурации 01^-Р4

Рис. 16. Течение в области стыка крыла с фюзеляжем

ОСНОВНЫЕ РЕЗУЛЬТАТЫ РАБОТЫ И ВЫВОДЫ

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

2. Для эффективного решения пространственных уравнений конечно-элементной аппроксимации уравнений Эйлера и Навье-Стокса впервые реализован полиномиальный многосеточный подход с локально-

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

3. Предложено два основных пути сокращения арифметических операций РМГ в расчёте на одну степень свободы:

о использование аналитического способа интегрирования потоков на

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

расчёта поверхностных и объёмных интегралов.

4. Впервые метод РМГ высокого порядка точности реализован для решения пространственных задач аэродинамики и акустики с использованием неструктурированных гексаэдральных сеток и приведено его широкое сравнение с МКО.

5. Проведены численные исследования точности новой схемы аппроксимации на последовательности вложенных сеток. В процессе решения уравнений Эйлера и Навье-Стокса на неструктурированных сетках продемонстрирована оптимальная сходимость результатов к точным значениям с порядком К+\, где К - максимальный порядок базисного полинома.

6. Рассмотрены разнообразные тестовые случаи (вязкое дозвуковое обтекание цилиндра, сферы, пластины; пространственное течение вязкой жидкости в изогнутой трубке; распространение акустических волн в невязком газе; турбулентное трансзвуковое обтекание профиля и крыла), демонстрирующие хорошее согласование расчётов по предложенной схеме высокого порядка точности с теоретическими и (или) экспериментальными результатами даже при использовании экстремально грубых сеток. При сравнении аналогичных по точности решений, полученных методами РМГ высокого порядка и стандартной схемой МКО второго порядка, оказалось, что наблюдается устойчивый выигрыш схемы РМГ в использовании вычислительных ресурсов (память, время).

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

8. Показано, что разработанный метод позволяет проводить анализ особенностей местной аэродинамики ЛА, например, области стыка крыла с фюзеляжем, что позволяет использовать настоящий подход на этапе

предварительного проектирования и открывает перспективы улучшения качества проектирования новых ЛА.

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

ОСНОВНЫЕ РЕЗУЛЬТАТЫ ОПУБЛИКОВАНЫ В РАБОТАХ (в научных журналах, представленных в перечне ВАК)

1. А.В. Волков. Особенности применения метода Галеркина к решению пространственных уравнений Навье-Стокса на неструктурированных гексаэдральных сетках // Ученые записки ЦАГИ. т. XL №6,2009.

2. А.В. Волков, С.В. Ляпунов. Монотонизация метода конечного элемента в задачах газовой динамики // Ученые записки ЦАГИ. т. XL №4,2009.

3. В.Д. Боксер, А.В. Волков, А.В. Петров. Применение тангенциального выдува струй для снижения сопротивления сверхкритических профилей при больших дозвуковых скоростях // Ученые записки ЦАГИ. т. XL №1,2009.

4. N.B. Petrovskaya, A.V. Wolkov, S.V. Lyapunov. Modification of Basis Functions in High Order Discontinuous Galerkin Schemes for Advection Equation // Applied Mathematical Modelling. Volume 32, Issue 5 pp.826-835, May 2008.

5. В.Д. Боксер, А.В. Волков, А.В. Петров, Г.Г. Судаков. Можно ли улучшить аэродинамику сверхкритического крыла? Общероссийский научно-технический журнал «Полет», стр. 70-76,2008.

6. N.B. Petrovskaya, A.V, Wolkov. The Issues of Solution Approximation in Higher Order Schemes on Distorted Grids // International Journal of Computational Methods. (IJCM) v. 4, No. 2 Page: 367 - 382, June 2007.

7. А.В. Волков, С.В. Ляпунов. Применение конечно-элементного метода Галеркина с разрывными базисными функциями к решению уравнений Рейнольдса на неструктурированных адаптивных сетках // Ученые записки ЦАГИ. т. XXXVIII, № 3.4, Стр. 22-30,2007.

8. А.В. Волков, С.В. Ляпунов. Исследование эффективности

использования численных схем высокого порядка точности для решения уравнений Навье-Стокса и Рейнольдса на неструктурированных адаптивных сетках // ЖВМиМФ, Том 46, №10, стр. 1894-1907,2006.

9. А.В.Волков, С.В.Ляпунов, А.Н.Храбров. Влияние установившегося

вращения на аэродинамические характеристики профиля при наличии отрыва потока // Ученые записки ЦАГИ. Том XXXIV, №3-4,2003.

10. А.В.Волков, С.В.Ляпунов. Метод расчета вязкого отрывного обтекания систем крыловых профилей // Ученые записки ЦАГИ, том XXIX №3-4, 1998.

11. A.V.Wolkov, S.V.Lyapunov. Numerical Prediction of Transonic Viscous Separated Flow Past an Airfoil // Theoretical and Computational Fluid Dynamics, vol.6, №1,1994.

12. А.В.Волков, С.В.Ляпунов. Метод расчета трансзвукового вязкого обтекания профиля с учетом изменения энтропии на скачках уплотнения // Ученые записки ЦАГИ, T.XXIV, №1, 1993.

Представлены в Трудах конференций:

13. Ch.Hirsch, A.Wolkov, B.Leonard. Discontinuous Galerkin Method on Unstructured Hexahedral Grids // AIAA paper 2009-01771.

14. А.В.Волков. Анализ компьютерных ресурсов необходимых при численном решении 3-D уравнений Навье-Стокса методом DG. Материалы ХХ-ой школы-семинара "Аэродинамика летательных аппаратов", ЦАГИ, 2009 г.

15. Charles Hirsch, Andrey Wolkov, Benoit Leonard. Discontinuous Galerkin Method on Unstructured Hexahedral Grids for 3D Euler and Navier-Stokes Equations. 5th. European Congress on Computational Methods in Applied Sciences and Engineering //ECCOMAS 200, held at the Lido Island in Venice (Italy) on 30 June - 4 July, 2008.

16. А.В. Волков. Учёт кривизны обтекаемой поверхности в методе высокого порядка точности. Материалы XIX-ой школы-семинара "Аэродинамика летательных аппаратов", ЦАГИ, 2008 г.

17. A.Wolkov, Ch.Hirsch, B.Leonard. Discontinuous Galerkin Method on Unstructured Hexahedral Grids for 3D Euler and Navier-Stokes Equations //AIAA paper 2007-4078.

18. A.B. Волков. Применение р-многосеточного метода к ускорению решения конечно-элементного метода аппроксимации при решении пространственных уравнений Навье-Стокса. Материалы XVIII-ой школы-семинара "Аэродинамика летательных аппаратов", ЦАГИ, 2007 г.

19. V.Vlasenko, A.Wolkov, Ch. Hirsch, Computationally effective Discontinuous Galerkin scheme for Linearized Euler Equations // WEHSFFC 2007. Moscow, November 19-22,2007.

20. A.B. Волков. Применение метода Галеркина с разрывными базисными функциями к решению пространственных уравнений Эйлера и Навье-Стокса. Труды VI международной школы семинар ММА «Модели и Методы Аэродинамики», Евпатория 2006.

21. А.В.Волков, С.В.Ляпунов. Численный метод высокого порядка точности для решения системы уравнений Навье-Стокса и Рейнольдса. ММА («Модели и Методы Аэродинамики»), Материалы Третьей Международной школы-семинара, Евпатория 2003.

22. S.V.Lyapunov, A.V.Wolkov. Application of Discontinuous Galerkin finite element method to the solution of partial differential equations. Part II. System of nonlinear equations: Euler equations // 16th IMACS World Congress. Lausanne-Switzerland August 21-25,2000.

23. A.V.Wolkov, S.V.Lyapunov. Application of Discontinuous Galerkin finite element method to the solution of partial differential equations. Part I. 2D scalar conservation laws // 16th IMACS World Congress. Lausanne-Switzerland August 21-25,2000.

24. V.S.Sakovich, A.M.Sorokin, A.V.Wolkov, S.V.Lyapunov. Anisotropic unstructured grid generation for 3D flow simulation problems. 6th Int. Conf. On numerical grid generation in computational fluid simulation, 1998.

25. A.V.Wolkov, S.V.Lyapunov. A Separated Flow Calculation About Airfoils and High-Lift Systems on Basis of Viscous-Inviscid Interaction Methods. EUROMECH, 3rd European Fluid Mechanics Conference, 1997, Gottingen, Germany

26. S.V.Lyapunov, A.V.Wolkov. Application of Viscous-Inviscid Interaction Methods for a Separated Flow Calculation About Airfoils and High-Lift Systems. ICAS Proceedings, ICAS-96-1.10.2,1996.

Издательский отдел ЦАГИ Заказ № 5364, Тираж 100 экз.

Используемые аббревиатуры и обозначения.

ВВЕДЕНИЕ.

1. МЕТОД ГАЛЁРКИНА С РАЗРЫВНЫМИ БАЗИСНЫМИ ФУНКЦИЯМИ И ОСОБЕННОСТИ ЕГО РЕАЛИЗАЦИИ ПРИМЕНИТЕЛЬНО К РЕШЕНИЮ УРАВНЕНИЙ НАВЬЕ-СТОКСА.

1.1 Методы аппроксимации уравнений Навье-Стокса.

1.2 Сравнение схем конечного объёма и конечного элемента на деформированных сетках.

1.3 Особенности дискретизации уравнений Навье-Стокса методом Галеркина с разрывными базисными функциями.

1.4 Методы приближённого решения задачи Римана в РМГ.

1.5 Особенности применения модели турбулентности Спаларта-Алмараса

1.6 Квадратурный и аналитический способы интегрирования потоков.

1.7 Базисные функции РМГ.

1.8 Методы учёта кривизны обтекаемой границы.

Совершенствование аэродинамических форм летательных аппаратов (ЛА) и их элементов достигается за счёт проведения большого объёма дорогостоящих трубных и лётных экспериментов. Высокая стоимость такого рода исследований делает целесообразным использование численных методов при проектировании ЛА. На этапе зарождения численных алгоритмов стоимость расчёта аэродинамических характеристик всегда оставалась пренебрежимо малой по сравнению со стоимостью эксперимента. Однако достоверность расчётных результатов была недостаточной из-за использования сравнительно простых моделей течения. Первые численные алгоритмы основывались на уравнениях идеального несжимаемого газа и ограничивались рассмотрением плоских (Серебрийский 1944, Павловец 1971) или упрощённых пространственных течений (Белоцерковский 1965). Совершенствование вычислительной техники позволило перейти к расчёту невязких сжимаемых течений путём решения уравнения для потенциала скорости и уравнений Эйлера (вагаЬеШап 1971, Колган 1972, Минайлос 1976, Ляпунов 1980, Вышинский 1983). Развитие асимптотических методов исследований задачи взаимодействия внешнего потенциального течения с пограничным слоем (Нейланд 1969; Сычев В., Рубан, Сычев Вик., Королев 1987) позволило перейти к более сложным моделям течений. При практическом проектировании элементов летательного аппарата широкое применение получили зональные подходы, базирующиеся на решении задачи о вязко-невязком взаимодействии (Брутян 1978, ЬеВаНеиег 1981, \Уо11соу & Ьуарипоу 1994, Волков и Ляпунов 1998). Эти методы позволили выполнять расчеты не только безотрывного обтекания, но и рассчитывать режимы с локальными зонами отрыва пограничного слоя. Широкую известность в промышленности получила вычислительная программа ВЬ\¥Б Карася и Ковалёва (1989), которая даёт возможность получать достоверные распределённые и интегральные аэродинамические характеристики компоновок современных JIA с учётом влияния вязкости на до- и трансзвуковых режимах обтекания.

Развитие вычислительной техники и дальнейший прогресс вычислительных методов позволил перейти к рассмотрению модели течения вязкого газа на основе уравнений Навье-Стокса. В последние полтора десятилетия было разработано большое количество подходов к решению этих уравнений (Anderson et al.l 996, Jespersen et al. 1997, Gerhold et al. 1997, Mavriplis & Venkatakrishnan 1997, ICrist et al. 1998, Kroll et al. 1998, Босняков, Власенко, Матяш, Михайлов 2007). Наиболее удачные численные схемы, реализованные в известных коммерческих вычислительных программах (FLUENT, CFX, STAR-CD, NUMECA), получили широкое распространение и относительно успешно используются при решении многих прикладных задач. Отметим, что в подавляющем большинстве численные коды базируются на методах конечных разностей (МКР) или конечного объёма (МКО). Точность повсеместно используемых численных схем, как правило, не превышает второй порядок, что в настоящее время является фактически стандартом.

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

Высокая стоимость вычислительных работ при прямом численном решении уравнений Навье-Стокса обусловлена необходимостью привлечения больших компьютерных ресурсов (оперативная память, время расчёта) и высокопрофессиональных вычислителей. Этот факт убедительно и продемонстрирован в трудах известной международной конференции AIAA CFD Drag Prediction Workshop (DPW) (Lee-Rausch et al. 2003, Rumsey et al. 2004, Morrison & Hemsch 2006), посвященной оценке величины аэродинамического сопротивления тестовой конфигурации «крыло+фюзеляж». Выбранная модель DLR F—4 была исследована в различных аэродинамических трубах. Оценка точности экспериментальных измерений на режиме М=0.85, Су= 0.5, Re=5T06 составляет величину 4-Ю-4. При такой погрешности в определении сопротивления погрешность величины аэродинамического качества компоновки составляет 0.2 единицы, что приемлемо при современной технике проектирования JIA. Очевидно, что не следует ожидать от расчётов идеальной сходимости с экспериментом, хотя бы по причине отсутствия совершенной модели турбулентности и отсутствия учёта многих других факторов, оказывающих влияние на формирование реального течения. Однако требуется получить точное численное решение в рамках имеющейся модели течения, описываемой уравнениями Навье-Стокса и конкретной моделью турбулентности. Теоретически такое решение может существовать на сетках с очень мелкими ячейками, дальнейшее измельчение которых не приводит к заметным изменениям решения. В рассматриваемом тестовом случае достаточное приближение к точному численному решению, как оказалось, достигается лишь на структурированных сетках с объёмом ячеек более 10 млн. Отметим, что построение сеток такого типа представляет собой трудоёмкий процесс, требующий высокой квалификации пользователя. Неструктурированные сетки могут быть сгенерированы значительно быстрее и с минимальным участием пользователя. Однако, добиться при этом приемлемого качества величины сопротивления удаётся лишь на сетках с объёмом более 30 млн. узлов.

Практический расчёт обтекания JIA в крейсерской конфигурации не является задачей, требующей решения именно уравнений Навье-Стокса и привлечения значительных вычислительных ресурсов. В настоящее время с расчётом подобной конфигурации надёжно справляется, например, вычислительная программа ВЬ\УТ, основанная на более простой модели течения. Очевидно, что уравнения Навье-Стокса необходимо применять для расчетов геометрически сложных конфигураций, течений со значительными отрывными зонами, например, для задачи расчёта взлётно-посадочных аэродинамических характеристик ЛА при существенно отклонённой механизации крыла или, например, для расчёта обтекания маневренных самолётов при больших углах атаки. Построение структурированных сеток для таких конфигураций является технически сложной проблемой. Простые оценки показывают, что необходимый объём неструктурированных сеток в этом случае может превысить сто миллионов узлов. Вычислительные ресурсы, требуемые для решения таких задач, пока доступны лишь отдельным крупным научным центрам. Поэтому снижение потребных вычислительных ресурсов является одной из актуальных проблем, решение которой позволит обеспечить реальный прогресс во многих практически важных приложениях.

Другая проблема современных методов решения уравнений Навье-Стокса, связана с необходимостью минимизации влияния человеческого фактора, неизбежно присутствующего при проведении вычислений. Часто результаты численных решений одних и тех же задач у разных пользователей заметно отличаются друг от друга. В основном различие связано с построением расчётных сеток. Современные пакеты генерации сеток предоставляют пользователю широкий набор возможностей, позволяя измельчать или укрупнять ячейки в заданных областях расчётного поля. Каждый пользователь делает это в соответствии со своими представлениями о физике течения, что в конечном итоге и обуславливает разницу результатов. Отметим, что такое «ручное» построение сеток может быть причиной «необъяснимых» существенных различий в результатах полученных одним и тем же пользователем на незначительно модифицированных геометриях. Решение проблемы сопряжено с задачей определения оптимального положения сеточных узлов или, другими словами, с необходимостью использования алгоритмов адаптации сеток и численных схем к особенностям течения. Идеи адаптации известны давно (см., например, 8акоу1сИ е1 а1. 1998), однако, в реальных расчётах они практически не используются по двум основным причинам. Во-первых, технология адаптации предполагает использование неструктурированных сеток, а как было описано выше, точность расчётов на неструктурированных сетках заметно хуже. Во-вторых, реальная адаптация приводит к необходимости расчётов на неравномерных сетках с сильно деформированными ячейками. Современные же схемы, базирующиеся на МКР и МКО, ориентированы на использование равномерных сеток с ячейками правильных форм.

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

Рис. 1. Фрагменты сетки с плавной вариацией форм и размеров ячеек пригодной для расчета МКО. Количество ячеек -110 ООО

Данная сетка содержит -110 тыс. ячеек и вполне годится для получения достоверных результатов при помощи МКО. Использование для решения этой задачи адаптивных сеток с сильной анизотропией форм ячеек позволяет существенно сэкономить их необходимое количество. Пример такой сетки изображён на рис. 2.

Рис. 2. Фрагменты адаптивной сетки с сильной анизотропией форм ячеек.

Количество ячеек ~ 28 ООО

Такое же покрытие поля течения обеспечивается значительно меньшим количеством ячеек ~ 28 тыс. Данную сетку характеризует высокая неравномерность форм и размеров ячеек, соседство крупных элементов с тонкими вытянутыми («щелочными») элементами. Использование МКО для решения задачи на такой сетке невозможно из-за потери точности аппроксимации и сложностей поиска решения получающейся системы сеточных уравнений. Таким образом, разработка схем расчёта на подобных высокоанизотропных сетках позволит, с одной стороны, снизить объёмы расчётных сеток, и, с другой стороны, устранить влияние человеческого фактора на получение окончательных результатов путём использования технологии автоматического построения сеток.

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

Точность

Рис. 3. Зависимость точности численной схемы от числа неизвестных дискретной задачи

Таким образом, одно из решений проблемы эффективности численных схем связано с использованием схем высокого порядка точности, позволяющих выполнять расчёты на анизотропных адаптивных сетках. Разработка высокоточных схем ускорит решение не только задач внешней и внутренней аэродинамики, но также позволит решать многие другие актуальные задачи. Среди них задачи акустики, моделирование турбулентности (LES/DES) и другие.

Схемы с высоким порядком дискретизации исследовались на компактном шаблоне в работах Толстых (1990), Lele (1992), Tarn & Webb (1993), Tolstykh & Lipavskii (1998), а на расширенном шаблоне со стандартным методом конечных разностей в работе Visbai & Gaitonde (2002). Для получения монотонного решения с высоким порядком аппроксимации были разработаны также методы ENO (Harten et al. 1987), WENO (Liu et al. 1994, Jiang & Shu 1996). Все эти методы имеют солидный теоретический фундамент, но, к сожалению, применяются в основном на структурированных и неструктурированных сетках с почти равномерным распределением размеров и форм соседних ячеек. Возможность применения таких методов на анизотропных неструктурированных сетках не столь очевидна. Кроме того, в случае расширения шаблона аппроксимации открытыми остаются вопросы распараллеливания алгоритма расчёта и его точности в случае использования многоблочных сеток.

В работе Ladeinde et al. (2006) впервые расчетная схема высокого порядка точности использована при решении задач внешней аэродинамики для тел сложных геометрических форм. Здесь рассматривалась возможность создания надежного, высокоточного алгоритма решения уравнений Рейнольдса на базе ENO/WENO схем с использованием блочно-структурированных сеток. В качестве примера рассмотрен расчет компоновки современного самолета.

В настоящее время имеется широкий выбор методов, альтернативных МКР или МКО. В работах Ekaterinaris (2000) и Wang (2007) дан обзор различных способов аппроксимации высокого порядка точности на неструктурированных сетках.

Одним из наиболее перспективных подходов к высокоточной аппроксимации как на структурированных, так и неструктурированных сетках является метод Галеркина с разрывными базисными функциями (РМГ). В англоязычной литературе метод имеет название Discontinuous Galerkin Method (DGM). В последние годы этот метод вызывает повышенный интерес многих исследователей вследствие его общности, гибкости и надежной теоретической обоснованности. Впервые метод был предложен Reed & Hill (1973) для решения уравнения описывающего перенос нейтронов, а первый анализ был дан в работе Saint & Raviart (1974). Численное решение 2-D уравнений Эйлера и Навье-Стокса на треугольных неструктурированных сетках этим методом впервые представлено в работах Bassi & Rebay (1997). Наиболее полное теоретическое описание метода с решением 1-D и 2—D модельных задач приведено в работах Cockburn & Shu (1998, 2001). Подробное исследование

РМГ можно также найти в работах Lyapunov & Wollcov (2000), Волкова и Ляпунова (2006, 2007, 2009) и в диссертации Ляпунова (2008).

Первые 3-D реализации метода относятся к решению задач аэроакустики, где решались линеаризованные уравнения Эйлера (Atkins & Lockard 1999, Remalci et al. 2002, Reymen et al. 2005, Peyret 2007). Практическая реализация РМГ для решения нелинейных уравнений Эйлера и Навье-Стокса встречается с принципиальными трудностями, по причине которых этот метод ещё не получил широкого распространения. В настоящее время можно выделить несколько основных проблем, которые необходимо решить для успешного применения РМГ.

1) Большое количество арифметических операций на одну степень свободы. В отличие от МКО в РМГ необходимо точное интегрирование не только потоков через грани контрольного объёма, но и вычисление объёмных интегралов внутри ячеек. Традиционный подход к интегрированию базируется на использовании квадратурных точек Гаусса. Суммирование значений подынтегральных функций в этих точках с определёнными весами обеспечивает получение точного значения в случае полиномиального представления подынтегральной функции. Чем выше порядок полинома, тем большее число гауссовых точек необходимо использовать. Присутствие объёмных и поверхностных интегралов в рассматриваемом методе приводит к необходимости использования чрезмерно большого количества этих точек и на гранях контрольного объёма и внутри него.

2) Немонотонность схемы высокого порядка в областях разрывов решений. Решение задач обтекания может характеризоваться наличием в поле течения областей с большими локальными возмущениями, в частности, скачков уплотнения. Численное решение таких задач с использованием схем повышенного порядка точности приводит к нефизичным осцилляциям решения. В МКО для устранения нефизичных осцилляций используют ограничители градиентов решения, определяемые по значениям в ячейке и её соседях. Такой способ неприменим в РМГ для решения стационарных задач неявным методом.

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

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

4) Проблема решения больших нелинейных систем сеточных уравнений.

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

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

В настоящее время число 3-0 реализаций РМГ для нелинейных уравнений законов сохранения весьма ограничено. В открытой печати можно найти информацию лишь о нескольких примерах успешной реализации метода. К примеру, применение РМГ на тетраэдальных сетках описано в работе Luo et al. (2006). Первая успешная реализация РМГ на неструктурированных гексаэдральных сетках представлена в настоящей работе и в публикациях: Wolkov, Hirsch, Leonard (2007); Hirsch, Wolkov, Leonard (2009); Волков (2009); Волков (2010).

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

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

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

Практическая значимость работы заключается в разработке эффективных методов расчёта уравнений Эйлера и Навье-Стокса. Дан набор рекомендаций по применению предложенного подхода. Созданные в процессе работы вычислительные программы используются при проведении расчётных исследований, как отдельных элементов ЛА, так и пространственной конфигурации в целом, в процессе выполнения инициативных научных исследований, а также при выполнении НИР ЦАГИ по контрактам с Роспромом. Разработанные методы открывают перспективы эффективного расчёта обтекания сложных конфигураций ЛА.

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

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

В качестве эффективного решателя системы сеточных уравнений впервые предложен полиномиальный многосеточный подход с локально-неявным сглаживателем.

Метод РМГ высокого порядка точности для решения пространственных задач аэродинамики впервые практически реализован на гексаэдральных неструктурированных сетках.

Выполнено непосредственное сравнение требуемых компьютерных ресурсов при выполнении одних и тех же расчётов методами МКО и РМГ.

Автор защищает следующие результаты:

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

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

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

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

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

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

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

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

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

Основные результаты диссертации опубликованы в работах Lyapunov & Wolkov (1996); Sakovich, Sorokin, Wolkov, Lyapunov (1998); Wolkov & Lyapunov (2000); Lyapunov & Wolkov (2000); Волков & Ляпунов (2004, 2006, 2007, 2009); Petrovskaya, Wolkov, Lyapunov (2006); Wolkov, Hirsch, Leonard (2007); Vlasenko, Wolkov, Hirsch (2007); Petrovskaya & Wolkov (2007, 2008); Hirsch, Wolkov, Leonard (2009); Волков (2009, 2010); Петровская & Волков (2010).

Результаты работы обсуждались на Международном конгрессе по авиационным наукам ICAS (1996), 3-ей Европейской конференции по механике жидкости EUROMECH (1997), 6-ой международной конференции по генерации сеток для вычислительной аэродинамики (1998), 16-м Конгрессе Международной ассоциации математического и компьютерного моделирования IMACS (2000), конференциях Американского Института Авиации и Аэронавтики AIAA (2007, 2009), международной конференции по высокоскоростным течениям WEHSFF (2007), 5-ом Европейском конгрессе по численным методам в прикладной науке ECCOMAS (2008), франко-российском семинаре ONERA-ЦАГИ (2009), международной школе-семинаре «Модели и методы аэродинамики» (2004, 2006), и на 12-ти школах-семинарах ЦАГИ «Аэродинамика летательных аппаратов» (1998 - 2009).

В данную диссертацию ' включены результаты исследований, поддержанные РФФИ (Проекты №96-01-00209-л, № 96-01-00210-л, № 00-01-00070-а, № 03-01-00236-а, №06-01-00283-а, №09-01-00243-а).

Автор выражает глубокую благодарность д.ф.м.н. С.В.Ляпунову, в тесной работе с которым были выполнены основные исследования, изложенные в диссертации. Автор считает своим долгом выразить благодарность профессору Г.А.Павловцу за постоянную поддержку и внимание к данной работе. Автор выражает признательность президенту компании NUMECA профессору Ч.Хиршу (Ch. Hirsch) за предоставление промышленной программы метода конечного объема и помощь в проведении настоящих исследований.

ЗАКЛЮЧЕНИЕ И ОСНОВНЫЕ РЕЗУЛЬТАТЫ

Решение проблем повышения качества и сокращения сроков проектирования новых ЛА всё в большей мере связывается с совершенствованием численных методов расчёта их аэродинамических характеристик. Модель течения, основанная на полных уравнениях Навье-Стокса, позволяет обеспечить новый уровень проектирования, однако современные подходы, основанные на методе конечного объёма (МКО) второго порядка точности, требуют привлечения очень больших компьютерных ресурсов, недоступных в настоящее время большинству пользователей.

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

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

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

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

4. Впервые метод РМГ высокого порядка точности реализован для решения пространственных задач аэродинамики и акустики с использованием неструктурированных гексаэдральных сеток и приведено его широкое сравнение с МКО.

5. Проведены численные исследования точности новой схемы аппроксимации на последовательности вложенных сеток. В процессе решения уравнений Эйлера и Навье-Стокса на неструктурированных сетках продемонстрирована оптимальная сходимость результатов к точным значениям с порядком К+1, где К максимальный порядок базисного полинома.

6. Рассмотрены разнообразные тестовые случаи (вязкое дозвуковое обтекание цилиндра, сферы, пластины; пространственное течение вязкой жидкости в изогнутой трубке; распространение акустических волн в невязком газе; турбулентное трансзвуковое обтекание профиля и крыла), демонстрирующие хорошее согласование расчётов по предложенной схеме высокого порядка точности с теоретическими и (или) экспериментальными результатами даже при использовании экстремально грубых сеток. При сравнении аналогичных по точности решений, полученных методами РМГ высокого порядка и стандартной схемой МКО второго порядка, оказалось, что наблюдается устойчивый выигрыш схемы РМГ в использовании вычислительных ресурсов (память, время).

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

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

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

1. Белоцерковский С.М. (1965) Тонкая несущая поверхность в дозвуковом потоке газа. М.: Наука.

2. Босняков С.М. (2007) Концепция программного продукта EWT ЦАГИ и основные этапы ее развития // Труды ЦАГИ, выпуск 2671, 2007, стр.3-19.

3. БрутянМ.А. (1978) Исследование отрывного обтекания симметричного профиля в несжимаемой жидкости // Ученые записки ЦАГИ, 6.

4. Власенко В.В. (2007) О математическом подходе и принципах построения численных методологий для пакета прикладных программ EWT ЦАГИ. стр.20-85 // Труды ЦАГИ, выпуск 2671, стр.3-19.

5. Волков A.B. (2006) Применение метода Галеркина с разрывными базисными функциями к решению пространственных уравнений Эйлера и Навье-Стокса // Труды VI международной школы семинар «Модели и Методы Аэродинамики». Евпатория 2006.

6. Волков A.B. (2009) Особенности применения метода Галеркина к решению пространственных уравнений Навье-Стокса на неструктурированных гексаэдральных сетках // Ученые записки ЦАГИ, т. XL, .№6.

7. Волков A.B. (2010) Применение многосеточного подхода к решению 3-D уравнений Навье-Стокса на гексаэдральных сетках методом Галеркина с разрывными базисными функциями // ЖВМиМФ, т.50, №3, с. 517 531.

8. Волков A.B. (2010) Методы решения сеточных уравнений конечно-элементной аппроксимации пространственных течений. Ученые записки ЦАГИ, т. XLI, .№3.

9. Волков A.B., Ляпунов C.B. (1998) Метод расчета вязкого отрывного обтекания систем крыловых профилей // Ученые записки ЦАГИ, том XXIX №3-4.

10. Волков A.B., Ляпунов C.B. (2004) Метод и результаты расчета аэродинамических характеристик крылового профиля при наличии льда разнообразных форм на передней кромке. Техника воздушного флота. Том LXXVII, № 1 (666).

11. Волков A.B., Ляпунов C.B. (2006) Исследование эффективности использования численных схем высокого порядка точности для решения уравнений Навье-Стокса и Рейнольдса на неструктурированных адаптивных сетках // ЖВМиМФ, Том 46, №10, стр. 1894-1907.

12. Волков A.B., Ляпунов C.B. (2007) Применение конечно-элементного метода Галеркина с разрывными базисными функциями к решению уравнений Рейнольдса на неструктурированных адаптивных сетках // Ученые записки ЦАГИ. том XXXVIII, № 3-4, стр. 22-30.

13. Волков A.B., Ляпунов C.B. (2009) Монотонизация метода конечного элемента в задачах газовой динамики // Ученые записки ЦАГИ. том XL, №4.

14. Вышинский B.B. (1983) К расчету пространственного околозвукового обтекания удлиненных тел // ЖВМиМФ, т.23, №1, с.160-169.

15. Годунов С.К. (1957) Разностный метод расчёта ударных волн. Успехи математических наук. Том. 12, № 1 (73), 176 177.

16. Годунов С.К., Рябенький B.C. (1973) Разностные схемы. (Введение в теорию). М.: Наука.

17. ГольдинВ.Я., Калиткин H.H., ШишоваТ.В. (1965) Нелинейные разностные схемы для гиперболических уравнений// ЖВМиМФ 5, №5, 938944.

18. Егоров И.В., Зайцев O.JI. (1991) Об одном подходе к численному решению двумерных уравнений Навье-Стокса методом сквозного счёта // ЖВМиМФ, т.31, №2.

19. Жуков В.Т., Феодоритова О.Б., Янг Д.П. (2004) Итерационные алгоритмы для схем конечных элементов высокого порядка// Матем. моделирование, том 16, номер 7, с. 117-128.

20. Иванов М.Я., Крайко А.Н. (1978) Об аппроксимации разрывных решений при использовании разностных схем сквозного счёта// ЖВМиМФ, №3, стр. 780-783.

21. Карась О.В., Ковалев В.Е. (1989) Применение обратного метода расчета трехмерного пограничного слоя к задаче обтекания крыла с учетом влияния вязкости // Уч. зап. ЦАГИ. Т. 20. № 5. С. 1-11.

22. Колган В.П. (1972) Применение принципа минимальных значений производных к постороению конечно-разностных схем для расчёта разрывных решений газовой динамики// Учёные записки ЦАГИ, 3, №6, 6872.

23. Корн Г., Корн Т. (1973) Справочник по математике для научных работников и инженеров / М.: Наука.

24. Куликовский А.Г., Погорелов Н.В., Семенов А.Ю. (2001) Математические вопросы численного решения гиперболических систем уравнений / М.: ФИЗМАТЛИТ.

25. Ландау Л.Д., Лифшиц Е.М. (1986) Гидродинамика. Теоретическая физика 6 /М.: Наука.

26. Ляпунов C.B. (1980) Программа расчёта обтекания профиля трансзвуковым потоком идеального газа // Труды ЦАГИ. Вып. 2064. М.

27. Ляпунов C.B. (2008) Разработка и исследование численных схем высокого порядка точности для решения уравнений газовой динамики на неструктурированных сетках / Диссертация на соискание ученой степени доктора физико-математических наук.

28. Меньшов И.С. (1991) Обобщенная задача о распаде произвольного разрыва // Прикладная математика и механика, т. 55, №1, 86-94.

29. Минайлос А.Н. (1976) Сверхзвуковое течение у тонкого трапециевидного крыла // Ученые записки ЦАГИ. т.7, №4.

30. Михайлов C.B. (2007) Объектно-ориентированный подход к созданию эффективных программ, реализующих параллельные алгоритмы расчета // Труды ЦАГИ, выпуск 2671, стр.86-108.

31. НейландВ.Я. (1969) К теории отрыва ламинарного пограничного слоя в сверхзвуковом потоке // Изв. АН СССР, МЖГ, №4, с.53-57.

32. ПавловецГ.А. (1971) Методы расчета обтекания сечений крыла идеальным несжимаемым потоком. Труды ЦАГИ, вып. 1344.

33. Петровская Н.Б., Волков A.B. (2010) Влияние геометрии сетки на точность реконструкции решения в конечно-объёмных и конечно-элементных схемах высокого порядка // Математическое моделирование, 22, №3, стр. 145 160.

34. Пинчуков В.И., Шу Ч.В. (2000) Численные методы высоких порядков для задач аэрогидродинамики / Новосибирск: Издательство Сибирского отделения РАН.

35. Родионов А.В. (1987) Повышение порядка аппроксимации схемы С.К. Годунова//ЖВМиМФ, 27, №12, 1853-1860.

36. Самарский А.А. (1989) Теория разностных схем. М.: Наука.

37. Серебрийский Я.М. (1944) Обтекание крыловых профилей произвольной формы// Труды ЦАГИ, вып. 553.

38. Сычев В.В., Рубан А.И., Сычев В.В., Королев Г.В. (1987) Асимптотическая теория отрывных течений/ Под ред. В.В. Сычева. М.: Наука.

39. Толстых А.И. (1990) Компактные разностные схемы и их применение в задачах аэрогидродинамики. М. : Наука.

40. Толстых А.И. (2000) О построении схем заданного порядка с линейными комбинациями операторов// ЖВМиМФ, 40, №8, 1206-1220.

41. Федоренко Р.П. (1961) Релаксационный метод решения разностных эллиптических уравнений// ЖВМ и МФ, т.1, №5, 922-927.

42. Федоренко Р.П. (1962) Применение разностных схем высокой точности для численного решения гиперболических уравнений// ЖВМиМФ 2, №6, 1122-1128.

43. Anderson W.K., RauschR.D., and Bonhaus D.L. (1996) Implicit Multigrid Algorithms for Incompressible Turbulent Flows on Unstructured Grids // Journal of Computational Physics, v. 128, No. 2, pp. 391-408.

44. Atkins H.L. and Lockard D.P. (1999) A high-order method using unstructured grids for the aeroacoustic analysis of realistic aircraft configurations// AIAA 991945.

45. Atkins H.L., Shu C.W. (1996) Quadrature-Free Implementation of Discontinuous Galerkin Method for Hyperbolic Equations // AIAA 96-1683.

46. Barter G.E., Darmofal D.L. (2007) Shock capturing with high-order, PDE-based artificial viscosity // AIAA-3823.

47. Barth T.J. (1991) A three-dimensional upwind Euler solver of unstructured meshes //AIAA-91-1548.

48. Barth T.J. (1992) Aspects of unstructured grids and finite-volume solvers for the Euler and Navier-Stokes equations. Special course on unstructured grid methods for advection dominated flows // AGARD-R-787, pp.6-1 6-61.

49. Barth T.J. (1993) Recent developments in high order k-exact reconstruction on unstructured meshes // AIAA Paper-0668.

50. Barth T.J. (2004) A posteriori error estimation and mesh adaptivity for finite volume and finite element methods. Springer series lecture notes in computational science and engineering, v. 41.

51. Barth T.J., Frederickson P.O. (1990) Higher order solution of the Euler equations on unstructured grids using quadratic reconstruction // AIAA-90-0013.

52. Barth T.J., Jespersen D.C. (1989) The design and application of upwind schemes on unstructured meshes// AIAA Paper-89-0366.

53. Barth T.J., Larson M.G. (2002) A-posteriori error estimation for higher order Godunov finite volume methods on unstructured meshes // NASA Technical Report NAS-02-001.

54. Bassi F. and RebayS. (1997) High-Order Accurate Discontinuous Finite Element Solution of the 2D Euler Equations // J. Comput. Phys., v. 138, pp. 251285.

55. Bassi F., RebayS. (1997) A High-Order Accurate Discontinuous Finite Element Method for the Numerical Solution of the Compressible Navier-Stokes Equations// Journal of Computational Physics v. 131, pp. 267-279.

56. Brandt A. (1980) Multilevel adaptive computations in fluid dynamics // AIAA Journal, V.18, №-10, October.

57. Chapman A., SaadY., WigtonL. (2000) High-order ILU preconditioners for CFD problems // Intematioanl Jouranal for numerical methods in fluids, v.33, 767788.

58. Cockburn B., Hou S., Shu C. (1990) TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws IV: the multidimensional case // Math. Comp., 54.

59. Cockburn B., ShuC-W. (1998) The Local Discontinuous Galerkin Finite Element Method for Convection-Diffusion Systems // SIAM J. Numer. Anal., v.175, pp. 2440-2463.

60. Cockburn B., Shu C-W. (2001) Runge-Kutta Discontinuous Galerkin Methods for Convection-Dominated Problems // Journal of Scientific Computing, v. 16, No. 3, September.

61. Cook P.H., McDonald M.A., & FirminM.C.P. (1979) Aerofoil RAE 2822 -pressure distribution, and boundary layer and wake measurements // AGARD-AR-138.

62. Ekaterinaris J.A. (2000) Implicit High-Order-Accurate-in-Space Algorithms for the Navier-Stokes Equations // AIAA Journal Vol. 38, No. 9, September 2000.

63. EnayetM., Gibson M., Taylor A., YianneskisM. (1982) Laser Doppler Measurements of Laminar and Turbulent Flow in a Pipe Bend // NASA Contract Report CR-3551.

64. EngquistB., OsherS. (1981) One-sided difference approximations for nonlinear conservation laws // Math. Comput. v.36, No. 154, 321-351.

65. Fidkowski K.J., Oliver T.A., Lu J., Darmofal D.L. (2005) p-Multigrid solution of high-order discontinuous Galerkin discretizations of the compressible Navier-Stokes equations//Journal of Computational Physics v.207 pp. 92-113.

66. FriedrichO. (1998) Weighted essentially non-oscillatory schemes for the interpolation of mean values on unstructured grids // Journal of Computational Physics, v. 144, pp. 194-212.

67. Garabedian P.R., Corn D.G. (1971) Analysis of transonic airfoils // Comm.Pure Appl.Math., vol.24, No 6.

68. Gerhold T., FriedrichO., Evans J., Galle M. (1997) Calculation of Complex Three-Dimensional Configurations Employing the DLR TAU-Code // AIAA Paper 97-0167.

69. Giles M., Pierce N.A. (1999) Improved lift and drag estimates using adjoint Euler equations // AIAA-99-3293.

70. Harten A. (1984) On a Class of High Resolution Total Variation-Stable Finite-Difference Schemes // SIAM J. Numer. Anal., v. 21.

71. Harten A., Engquist B, Osher S., Chakravarthy S. (1987) Uniformly high order essentially non-oscillatory schemes, III // Journal of Computational Physics, v.71, pp. 231-303.

72. Harten A., EngquistB., OsherS. and Chakravarthy S. (1987) Uniformly high order essentially non-oscillatory schemes, III // J. Comp. Phys., v. 71, pp. 231— 303.

73. Helenbrook B.T., Mavriplis D., Atkins H.L. (2003) Analysis of p-multigrid for continuous and discontinuous finite element discretizations // AIAA Paper 3989.

74. Hirsch Ch. (2007) Numerical Computation of Internal and External Flows / Elsevier.

75. Hirsch Ch., Wolkov A., B.Leonard B. (2008) Discontinuous Galerkin Method on Unstructured Hexahedral Grids for 3D Euler and Navier-Stokes Equations // 5th. European Congress on Computational Methods in Applied Sciences and Engineering (ECCOMAS 2008).

76. Houston P., Senior B., Slili E. (2002) hp-Discontinuous Galerkin finite element methods for hyperbolic problems: error analysis and adaptivity. Int. J. Numer. Meth. Fluid, v.40, pp. 153-169.

77. Hughes T.J.R., Brooks A.A (1979) Multidimensional upwind scheme with no crosswind diffusion, in Finite element methods for convection dominated flows. New York: ASME.

78. Hughes T.J.R., Engel G, Mazzei L, Larson M.G. (2000) Te continuous Galerkin method is locally conservative // Journal of Computational Physics, v.163, pp.467-488.

79. Hughes T.J.R., Franca L.P., and HulbertG.M. (1989) A new finite element formulation for computational fluid dynamics: VII The Galerkin/least-squares method for advective-diffusive systems // Comput.Methods Appl. Mech. Engrg, v.73, pp. 173-189.

80. Jespersen D.C., Pulliam T.H., and Buning P.G. (1997) Recent Enhancements to OVERFLOW // Tech. Rep. AIAA-97-0644.

81. Jiang G. and ShuC.-W. (1996) Efficient implementation of weighted ENO schemes // J. Comp. Phys., v. 126, pp. 202-228.

82. Johnson D.A., Bachalo W.D. (1980) Transonic flow past a symmetrical airfoil inviscid and turbulent flow properties // AIAA J. V. 18, p. 16-24.

83. Klaij C.M., van der Vegt J.J.W. and vanderVenH. (2006) Space-time discontinuous Galerkin method for the compressible Navier-Stokes equations, J. Comput. Phys., v.217, Issue 2 , pp. 589-611.

84. Krist S.L., Biedron R.T., and Rumsey C.L. (1998) CFL3D User's Manual (Version 5.0) // NASA TM-1998-208444.

85. Krivodonova L, BergerM. (2006) High-order accurate implementation of solid wall boundary conditions in curved geometries // Journal of Computational Physics v.211, pp. 492-512.

86. Krivodonova L., XinJ., Remacle J.-F., Chevaugeon N., Flaherty J.E. (2004) Shock detection and limiting with Discontinuous Galerkin methods for hyperbolic conservation laws// Applied Numerical Mathematics, v.48.

87. KrollN, Rossow C.-C., Becker K., Thiele F. (1998) MEGAFLOW A Numerical Flow Simulation System// 21stICAS Congress, ICAS.

88. Ladeinde F., Alabi K., Safta C., Cai X. and Johnson F. (2006) The First HighOrder CFD Simulation of Aircraft: Challenges and Opportunities // 44th AIAA Aerospace Sciences Meeting and Exhibit, AIAA-Paper 2006-1526, Jan.

89. Larson M.G., BarthT.J. (1999) A-posteriori error estimation for adaptive discontinuous Galerkin approximations of hyperbolic systems. // NAS Technical Report NAS-99-010.

90. LaxP.D. (1954) Weak solutions of nonlinear hyperbolic equations and their numerical computation// Comm. Pure Appl. Math, v.7 No.l, 159-193.

91. Le Balleuer J.C. (1981) Strong Matching Method for Computing Transonic Viscous Flows Including Wakes and Separations // La Recherche Aerospatiale, No.3, pp 21-45.

92. Lee-Rausch E., M. Buning P. G., Morrison J. H., Park M. A., Rivers S. M., Rumsey C. L. (2003) CFD Sensitivity Analysis of a Drag Prediction Workshop Wing/Body Transport Configuration AIAA 2003-3400.

93. Lele S.K. (1992) Compact finite difference schemes with spectral-like resolution. // J. Comp. Phys., v. 103, pp. 16^42.

94. Léonard B., Patel A., DelanayeM. and HischCh. (2000) All hexahedra unstructured mesh adaptive solver for turbulent flow simulations// Proc. Of the 1 st International Conference on Computational Fluid Dynamica (ICCFD), Kyoto, Japan, July.

95. Lesaint P., Raviart P.A. (1974) On a finite element method to solve the neutron transport equation, in Mathematical Aspects of Finite Elements in Partial Differential Equations, edited by C.de Boor. New York: Academic Press.

96. LiuX.D., OsherS. and ChanT. (1994) Weighted essentially-non-oscillatory schemes // J. Comp. Phys., v. 115, pp. 200-212.

97. LuoH., BaumJ.D., LôhnerR. (2006) A fast, p-multigrid discontinuous Galerkin method for compressible flows at all speeds // AIAA 2006-110.

98. LuoH., BaumJ.D., LôhnerR. (2007) A Hermite WENO-based limiter for Discontinuous Galerkin method on unstructured grids // AIAA 2007-510.

99. Lyapunov S.V., WolkovA.V. (1996) Application of the Viscous-Invisid Interaction Model to Calculation of Two-Dimensionai Separated Flows // TsAGI Journal, v.2, 1, 1996.

100. Lyapunov S. V., WolkovA.V. (1996) Application of Viscous-Inviscid Interaction Methods for a Separated Flow Calculation About Airfoils and HighLift Systems // ICAS Proceedings, ICAS-96-1.10.2.

101. Lyapunov S.V., WolkovA.V. (2000) Application of Discontinuous Galerkin finite element method to the solution of partial differential equations. Part I. 2D scalar conservation laws. 16th IMACS World Congress. Lausanne-Switzerland. August 21-25.

102. Mavriplis D.J. and Venkatakrishnan V. (1995) Agglomeration Multigrid for Two Dimensional Viscous Flows// Computers and Fluids, v. 24, No. 5, pp. 553570.

103. Mavriplis D.J. and Venkatakrishnan V. (1997) A Unified Multigrid Solver for the Navier-Stokes Equations on Mixed Element Meshes // International Journal for Computational Fluid Dynamics. No. 8, pp. 247-263.

104. Morrison J., HemschM. (2006) "Statistical Analysis of CFD Solutions from the 3rd AIAA Drag Prediction Workshop", 3rd AIAA APA Drag Prediction Workshop.http://aaac.Iarc.nasa.gov/tsab/cfdlarc/aiaadpw/Workshop3/presentations/index.html

105. MiillerU.R., SchulzeB., HenkeH. (1996) Computation of transonic steady and unsteady flow about LANN wing. Validation of CFD codes and assessment of turbulence models, ECARP Report 58, pp. 479-500, Vieweg.

106. Nastase C.R. Mavriplis D.J. (2006) Discontinuous Galerkin Methods Using an hp-Multigrid Solver for Inviscid Compressible Flows on Three-dimensional Unstructured Meshes // AIAA 2006-107.

107. Nielsen E.J., Andersen W.K., Walters R.W., Keyes D. (1995) Application of Newton-Krylov methodology to a three-dimensional unstructured Euler code // AIAA-95-1733.

108. Persson P.-O., PeraireJ. (2006) Sub-cell shock capturing for discontinuous Galerkin method// AIAA 112.

109. Petrovskaya N.B., Wolkov A.V. (2007) The Issues of Solution Approximation in Higher Order Schemes on Distorted Grids. International Journal of Computational Methods// (IJCM) v. 4, No. 2, pp: 367 382.

110. Petrovskaya N.B., Wolkov A.V., Lyapunov S.V. (2006) Modification of Basis Functions in High Order Discontinuous Galerkin Schemes for Advection Equation. Preprint 2006/26. School of Mathematics, Edgbaston, Birmigham B215 2TT, U.K.

111. Petrovskaya N.B., Wolkov A.V., Lyapunov S.V. (2008) Modification of Basis Functions in High Order Discontinuous Galerkin Schemes for Advection Equation. Applied Mathematical Modelling, v. 32, Issue 5, May 2008, pp.826-835.

112. Peyret C. (2007) hp Discontinuous Galerkin Method for Computational Aeroacoustics // AIAA-2007-3475.

113. ReedW.H. and Hill T.R. (1973) Triangular mesh methods for the neutron transport equation// Technical Report LA-UR-73—479. Los Alamos Scientific Laboratory.

114. RedekerG. (1994) DLR-F4 wing body configuration. A selection of experimental test cases for the validation of CFD code // AGARD-AR-303, v.2.

115. RemakiM., Habashi W.G., Ait-Ali-Yahia D. and Jay A. (2002). A 3D discontinuous Galerkin method for multiple pure tone noise problems // AIAA 2002-0229.

116. Reymen Y., De Roeck W., Rubio G., Baelmans M. and Desmet W. (2005) A 3D discontinuous Galerkin method for aeroacoustic propagation. // Twelfth International Congress on Sound and Vibration. ICSV12 Lisbon.

117. Roe P.L. (1981) Approximate Riemann problem solvers, parameter vectors, and difference schemes//J.Comput.Physics. v.43, No-2, 357-372.

118. Ronquist E.V., Patera A.T. (1987) Spectral element multigrid. I. Formulation and numerical results// Journal of Scientific Computing, v.2. No.4.

119. Rumsey C.L., Rivers S. M., Morrison J.H. (2004) Study of CFD variation on transport configurations from the second drag-prediction workshop // AIAA 2004394.

120. SaadY. (2004) Iterative Methods for Sparse Linear Systems. SIAM edition. Имеется в открытом доступе: http://www-users.cs.umn.edu/~saad/books.html.

121. SaadY. and Schultz M. (1986) GMRES: A Generalized Minimal Residual Algorithm for Solving Nonsymmetric Linear Systems// SIAM J. Sci. Statist. Comput. 7, 856-869.

122. Saint P. Le, RaviartP. (1974) On a finite element method for solving the neutron transport equation I in: C. de Boor (Ed.), Mathematical Aspects of Finite Elements in Partial Differential Equations, Academic Press, New York, pp. 89145.

123. Sakovich V.S., SorokinA.M., WolkovA.V., Lyapunov S.V. (1998) Anisotropic unstructured grid generation for 3D flow simulation problems // 6th Int. Conf. on numerical grid generation in computational fluid simulation.

124. ShakibF., Hughes T.J.R., JohanZ., (1989) A multi-element group preconditioned GMRES algorithm for nonsymmetric problems arising in finite element analysis// Comput. Methods Appl. Mech. Engrg, 87, pp 415-456.

125. ShuC.W. (2001) Different Formulations of the Discontinuous Galerkin Method for the Viscous Term //Advances in Scientific Computing.

126. ShuC.-W. (2001) High Order Finite Difference and Finite Volume WENO Schemes and Discontinuous Galerkin Methods for CFD // NASA/CR-2001-210865, ICASE Report No. 2001-11.

127. Shu C.-W., OsherS. (1988) Efficient implementation of essentially non-oscillatory shock capturing schemes // Journal of Computational Physics, v.11, pp. 439-471.

128. Shu C.-W., OsherS. (1989) Efficient implementation of essentially non-oscillatory shock capturing schemes II // Journal of Computational Physics, v.83, pp. 32-78.

129. Spalart P.R., Allmaras S.R. (1992) A One-Equation Turbulent Model for Aerodynamic flows. AIAA-92-0439.

130. Strang G., Fix G. (1973) An analysis of the finite element method. / Prentice Hall, New Jercey.

131. TamC.K.W. and Webb J.C. (1993) Dispersion-relation-preserving finite difference schemes for computational acoustics // J. Comp. Phys. v. 107, pp. 262.

132. Tolstykh A.I., Lipavskii M.V. (1998) On performance of methods with third-and fifth-order compact upwind differencing// J.Comp.Phys., v. 140, №2, 205-232.

133. Van LeerB. (1979) Towards the ultimate conservative difference schemes V. A second order sequel to Godunov's method// J.Comp.Phys., v.32.

134. Van der Vegt J.J.W, Van der Ven H.(2002) Space-time discontinuous Galerkin finite element method with dynamic grid motion for inviscid compressible flows. I. General formulation // J Comput Phys 2002;v.l82, pp.546-585.

135. Van der Ven H, Van der Vegt J.J.W. (2002) Space-time discontinuous Galerkin finite element method with dynamic grid motion for inviscid compressible flows. II. Efficient flux quadrature // Comput Meth Appl Mech Engrg; v.191, pp. 4747-4780.

136. Vassberg J.C., BuningP.G., RumseyC.L. (2002) Drag prediction for the DLR-F4 wing/body using OVERFLOW and CFL3D on an overset mesh. AIAA-2002-0840.

137. Venditti D.A., Darmofal D.L. (2003) Anisotropic grid adaptation for functional outputs: application to two-dimensional viscous flows // Journal of Computational Physics.

138. Venkatakrishnan V. (1988) Newton solution of inviscid and viscous problems // AIAA-99-0413.

139. Venkatakrishnan V. (1995) Comvergence to steady state solutions of the Euler Equations on Unstructured Grids with limiters // Journal of Computational Physics, v. 118, 120-130.

140. Venkatakrishnan V., Mavriplis D.J. (1993) Implicit solvers for unstructured meshes // Journal of Computational Physics, v. 105, pp.83-91.

141. Venkatakrishnan V., Allmaras S.R., Kamenetskij D.S., Johnson F.T. (2003) Higher order schemes for the compressible Navier-Stolces equations // AIAA-3987.

142. Visbal M.R. and Gaitonde D.V. (2002) On the use of higher-order finite-difference schemes on curvilinear and deforming meshes // J. Comp. Phys., v. 181, pp. 155-185.

143. Vlasenko V., WolkovA., HirschCh., (2007) Computationally effective Discontinuous Galerkin scheme for Linearized Euler Equations // West-East High Speed Flow Field Conference 2007. Moscow, November 19-22.

144. Wang Li, Mavriplis D.J. (2009) Adjoint-based h-p adaptive Discontinuous Galerkin methods for the compressible Euler equations //AIAA 2009-952.

145. Wang Z.J. (2007) High-order methods for the Euler and Navier-Stokes equations on unstructured grids // Progress in Aerospace Sciences v.43, pp 1-41.

146. WigtonL., Young D. (1985) GMRES acceleration of computational fluid dynamics codes// AIAA-1494.

147. Wolkov A., Hirsch Ch., Leonard B. (2007) Discontinuous Galerkin Method on Unstructured Hexahedral Grids for 3D Euler and Navier-Stokes Equations / 18th AIAA Computational Fluid Dynamics Conference // AIAA Paper 2007- 4078.

148. Wolkov A.V., Lyapunov S.V. (1994) Numerical Prediction of Transonic Viscous Separated Flow Past an Airfoil // Theoretical and Computational Fluid Dynamics, v.6, №1, 49-63.

149. Wolkov A.V., Lyapunov S.V. (1997) A Separated Flow Calculation About Airfoils and High-Lift Systems on Basis of Viscous-Inviscid Interaction Methods // EUROMECH, 3rd European Fluid Mechanics Conference, Book of Abstracts, Gottingen, Germany.