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

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

Ляпунов Сергей Владимирович

РАЗРАБОТКА И ИССЛЕДОВАНИЕ ЧИСЛЕННЫХ СХЕМ ВЫСОКОГО ПОРЯДКА ТОЧНОСТИ ДЛЯ РЕШЕНИЯ УРАВНЕНИЙ ГАЗОВОЙ ДИНАМИКИ НА НЕСТРУКТУРИРОВАННЫХ СЕТКАХ

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

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

□ОЗ168878

Москва-2008

003168878

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

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

профессор Крайко Александр Николаевич,

Ведущая организация - Вычислительный центр РАН (ВЦ РАН, Москва)

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

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

Автореферат разослан 2008 г

Ученый секретарь

Диссертационного совета Д 215 001 01

л

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

доктор физико-математических наук, профессор Хлопков Юрий Иванович доктор технических наук,

доцент Босняков Сергей Михайлович

А С Ненашев

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

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

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

Перспективными подходами к построению численных методов, удовлетворяющих перечисленным требованиям, являются применение неструктурированных вычислительных сеток, численных схем высокого порядка точности, адаптация сетки и численной схемы к решению В настоящей диссертационной работе предложены и исследованы методы адаптации неструктурированных сеток и локальной адаптации порядка точности численного метода к решению на базе метода конечного элемента Галеркина с разрывными функциями Впервые этот метод был предложен в работе (Reed, Hill 1973) для решения уравнения переноса нейтронов Дальнейшие многочисленные исследования были посвящены анализу математических аспектов метода, таких как скорость сходимости и пр (LeSamt, Raviart, 1974, Johnson, Pitkaranta, 1986), а также развитию метода (Cockburn, Shu, 1989, Cockburn, Lin, Shu, 1989, Cockbum, Hou, Shu, 1990, Cockburn, Shu, 1998), (Bassi, Rebay 1997), (Warburton, Lomtev, Kirby, Karniadakis 1998, Lomtev, Karniadakis 1997)

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

• Данный метод легко адаптируется к неструктурированным сеткам и, следовательно, удобен для работы с областями сложной геометрии

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

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

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

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

Цель диссертационной работы состоит в теоретической и методической разработке методов ускорения расчета и повышения точности результатов численного решения уравнений газовой динамики - уравнений Эйлера (невязкий случай) и уравнений Навье-Стокса и Рейнольдса (вязкий случай) путем адаптации не только неструктурированной расчетной сетки, но и локального порядка точности численной схемы Данные подходы основаны на модификации метода Галеркина с разрывными функциями (DG - Discontinuous Galerkin в англоязычной литературе), который является комбинацией метода конечного элемента и метода конечного объема типа метода Годунова Особое внимание уделяется анализу порядка точности получаемых схем, выявлению их преимуществ по сравнению со стандартными схемами типа Годунова Подробно рассмотрены возможности, которые обеспечивает схема DG с точки зрения адаптации к решению на неструктурированных сетках Приведены результаты исследований на примерах одномерных и двумерных задач Выбор этих задач в значительной степени обусловлен интересом к анализу порядка ошибки численного решения, что требует знания точного решения

Общая методика выполнения исследований состоит в

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

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

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

- выявлении положительных особенностей метода с точки зрения адаптации схемы к решению путем адаптации вычислительной сетки (Ь-гейпешеп!) или порядка точности схемы (р-геГтетаи)

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

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

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

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

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

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

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

Основные результаты проведенного автором исследования содержатся в 31 публикациях, опубликованных в российских научных изданиях и за рубежом, а также докладывались на международных и российских научно-технических конференциях, в том числе, на 6-ой международной конференции по генерации сеток для вычислительной аэродинамики (1998), 16-м Конгрессе Международной ассоциации математического и компьютерного моделирования IMACS (2000), 3-еи Европейской конференции по механике жидкости EUROMECH (1997), 5-ой Российско-китайской конференции по аэродинамике и механике полета (1997), Международных Конгрессах по авиационным наукам ICAS (1990, 1992, 1996), школах-семинарах «Аэродинамика летательных аппаратов (1998, 1999,2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007), Международной конференции молодых ученых и специалистов «Современные проблемы аэрокосмической науки и техники (2000), франко-российском семинаре ONERA-ЦАГИ (2002) и др

В данную диссертацию включены исследования, поддержанные РФФИ (Проекты № 98-01-00032-а, 1998, №00-01-00070-а, 2000, №02-01-00124-а, 2002, №03-01-00236-а, 2003, №06-01-00283-а, 2006)

Объем работы Диссертация состоит из введения, пяти глав, заключения и списка литературы, включающего 93 наименования Общий объем -127 страниц

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

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

В Главе I рассмотрена общая теория метода Галеркина с разрывными функциями (Ой) для численного решения одномерных законов сохранения Метод представляет собой симбиоз метода конечного элемента и конечного объема Метод применен к решению модельных задач для одномерных уравнений конвекции, теплопроводности и уравнения Бюргерса с вязкостью Данные примеры являются простейшими моделями конвективных, вязких членов уравнений газовой динамики, разрывных решений Основное внимание уделено анализу порядка ошибки численного решения

В § 1 1 рассмотрена общая теория метода Галеркина с разрывными функциями применительно к решению одномерных законов сохранения Конкретно рассматривается следующая задача определения функции и с начальным условием и периодическими граничными условиями для одномерного уравнения конвекции-диффузии с произвольной функцией невязкого потока /(и)

д,п + дх (/{и)-/лдхи) = 0 в области 0= (0,1) х (0,Т), ,,(/ = 0) = н0, хе (0,1),

ч(х = 0,0 = и(х = 1,0, дхи(х - 0,0 = дх11(х = 1,0

Здесь ¡1 - постоянный коэффициент диффузии

Это уравнение второго порядка сводится к системе уравнений первого порядка Данная система уравнений приводится к полудискретной системе обыкновенных дифференциальных уравнений путем умножения на базисные функции в каждом элементе (отрезке сетки) и интегрирования по элементу по частям (метод Галеркина) Базисными функциями являются полиномы Лежандра внутри каждого элемента Полученная система обыкновенных дифференциальных уравнений по времени решалась с помощью метода Рунге-Кутта Рассмотрены различные способы

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

В §12 рассмотрено численное решение одномерного уравнения конвекции Э,и + Эхи = 0 с начальным условием в виде П-образного импульса Точное решение соответствует перемещению этого импульса с единичной скоростью без изменения его формы Результаты расчетов показывают, что диссипативность схемы уменьшается с увеличением порядка полиномов, а точность решения повышается, что демонстрирует преимущества применения численных схем высокого порядка точности

На примере задачи для того же одномерного уравнения конвекции с гладким начальным условием были проведены исследования порядка ошибки численного решения в С-норме (модуль максимальной разницы между численным и точным решениями) Порядком ошибки называется величина N в асимптотической оценке ошибки Егг=0{Иы), где Ь-шаг сетки Эти исследования показали, что порядок ошибки близок к к+1, те при кусочно-постоянной аппроксимации численного решения (к = 0) получаем схему первого порядка точности, при кусочно-линейной аппроксимации (¿=1) - второго порядка и тд Эти же результаты свидетельствуют о возможности получения схемы произвольного порядка точности в рамках единообразного подхода

В §13 приведены результаты аналогичных исследований для одномерного уравнения теплопроводности вида д,и = дх(дхи) с тем же выводом - порядок численной ошибки близок к к+\ Здесь же приведены особенности аппроксимации вязких членов в рамках схемы Галеркина с разрывными функциями

В §1 4 приведен анализ решения одномерного уравнения Бюргерса с

вязкостью Э,г/ + Э;

2 у

= цд2хи Это уравнение при нулевом коэффициенте

вязкости (I дает разрывное решение при гладких начальных условиях и может моделировать решения со скачками уплотнения в газовой динамике И в этом случае порядок ошибки равен к+1 вне области разрыва решения

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

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

В §2.1 приведена постановка модельной задачи решения линейного уравнения конвекции-диффузии вида дхи + f'[x)dvii = цд~и в единичном

квадрате с разрывным условием вблизи левой границы квадрата. На границах квадрата, где течение «втекало» в область решения (вектор (1, f'(x) направлен внутрь области), использовались граничные условия для решения, соответствующие аналитическому решению. Аналитическое решение приведено на рис.1, где красному цвету соответствует значение решения +1, а синему - значение -1. Решение содержит криволинейный слой смешения, форма которого определяется выбором функции f(x), с переходом решения от -1 к +1. Ширина слоя определяется величиной коэффициента вязкости (I. С точки зрения газовой динамики эта задача моделирует развитие вязкого следа за профилем. Численная схема Галеркина с разрывными функциями применялась на неструктурированной сетке из треугольников, которая, на начальном этапе представляла собой триангуляцию структурированной сетки, как показано на рис. 1.

Рис. 1. Точное решение модельной задачи для двумерного уравнения конвекции-диффузии и расчетная сетка

Общая теория построения дискретной схемы аппроксимации данного двумерного уравнения с помощью метода Галеркина с разрывными функциями приведена в §2.2. Результирующие уравнения представляют собой систему линейных уравнений с блочной разреженной, вообще говоря, неструктурированной матрицей. Эта система уравнений решалась с помощью обобщенного метода минимальных невязок GMR.ES с предобуславливателем в виде неполного Ш разложения, описанного в §2.3.

В §2 4 приведены результаты численных исследований На рис 2 представлены в логарифмическом масштабе зависимости ошибки численного решения в норме Ь2 (средняя интегральная ошибка) и в норме С (максимальная ошибка) уравнения конвекции-диффузии в зависимости от числа степеней свободы (неизвестных переменных) при расчетах на разных сетках Приведены кривые для различных порядков к полиномиальных базисных функций (отмеченные на рис как Ой(к)) Эти зависимости иллюстрируют одно из преимуществ схем высокого порядка точности, а именно, начиная с некоторого числа переменных, схема высокого порядка точности будет давать меньшую ошибку, чем схема низшего порядка Анализ порядка ошибки, который равен удвоенному тангенсу угла наклона приведенных зависимостей, показал, что он (порядок) близок к к+1

Рис 2 Зависимости ошибок численного решения двумерного уравнения конвекции-диффузии от числа степеней свободы при различных порядках к полиномов базисных функций

Рассматриваемый метод Галеркина с разрывными функциями обладает большой гибкостью с точки зрения «приспосабливания» схемы к особенностям решения или расчетной сетки В частности, в данном параграфе была рассмотрена возможность использования полиномов различных порядков в различных ячейках Изменение порядка схемы в различных областях может быть использовано, как средство адаптации схемы к решению Такой способ адаптации в англоязычной литературе называется р-геАпетеп1, в отличие от Ь-геПпстсгИ: - адаптации с помощью изменения размеров и других геометрических характеристик ячеек

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

Макс ошибка

Пример сетки с «плохими» ячейками приведен на рис. 3. Здесь узлы сетки в одном сечении по оси х (.г = 0.55) были смещены по определенному закону. В результате получается искаженная сетка с двумя столбцами «плохих» ячеек. Три других графика на рис. 3 представляют поля ошибки решений, полученных различными способами. Численное решение, полученное с помощью кусочно-линейной аппроксимации (00(1), верхний правый график), дает максимальную ошибку 0.52. Отметим, что максимальная ошибка, полученная с 00(1) на неискаженной сетке, равнялась 0.32. Это иллюстрирует влияние «плохих» ячеек, вытянутых в направлении градиента решения, на его точность. Решение, полученное с помощью кусочно-квадратичной аппроксимации (00(2), нижний правый график), уменьшает максимальную ошибку до 0.31. Применение 00(1) во всей вычислительной области, за исключением двух столбцов плохих ячеек, где использовался метод Ов(2) (квадратичные полиномы), как показано на нижнем левом графике, дало максимальную ошибку, очень близкую к случаю применения метода 00(2) во всей области течения, а именно, 0.32. Решение перед столбцами «плохих» ячеек очень близко к решению с методом Ов( 1) во всей области течения, в то время как за этими столбцами - близко к решению методом 00(2) во всей области течения. Таким образом, применение методов более высокого порядка в областях с «плохими» элементами и более низкого - в остальной области течения позволяет повысить порядок точности решения.

Ошибка для 06(1) и 00(2)

в "плохих" ячейках. Макс. оши6ка=0.32 Ошибка для 00(2). Макс.оши6ка=0.31

Рис. 3. Исследование влияния локального изменения порядка полиномов в случае сетки с «плохими» ячейками

Увеличение порядка полиномов может также применяться локально не только в «плохих» ячейках, но и в областях с большими градиентами решения, где желательна более высокая точность. Результаты такого рода исследований приведены на рис. 4. Здесь представлены поля величин ошибок для метода 00(1) (верхний левый рис.) и 00(2) (верхний правый рис.). Левый нижний рис. представляет поле ошибки решения, когда метод 0С( 1) применялся всюду, за исключением ячеек, которые лежат в полосе ширины 0.05 (приблизительно два ряда ячеек) около линии скачка решения. В последнем случае ошибка в норме и максимальная ошибка стали гораздо меньше, чем в случае применения метода 00(1) во всей области (см. таблицу на рис. 4) при небольшом увеличении числа неизвестных (10%). Результаты аналогичного подхода при ширине полосы с квадратичной аппроксимацией решения (00(2)), равной 0.1, очень близки к результатам при применении метода 00(2) во всей области (которые помечены в таблице "ширина полосы - <3"), в то время как число неизвестных гораздо меньше, чем в случае применения метода ОО (2) во всей области.

Ошибка Ов(1)

Ошиькэ БЩ2)

1Ш|№На полосы ЧЮЛО Н*ИЗЮСТ«,И С111Иб)3 Б юрт 1, МЭМ.0ШИб>2

0 9Е00 20.79е-3 1203е-2

005 10560 3.41 е-3 279е-2

0.1 11520 1.55е-3 1.27е-2

О 19200 155е-3 127е-2

Ошибка 06(1)+06(2) в полосе шириной 0.05

Рис.4. Влияние локального увеличения порядка полипомов в области с большими градиентами решения

Еще один пример возможности адаптации метода 00 к решению приведен ниже. Если решение медленно меняется в каком-либо направлении, возможно использование полипомов низкого порядка для аппроксимации

решения в этом направлении и полиномов более высокого порядка для аппроксимации решения в ортогональном направлении Это направление может меняться от элемента к элементу В приведенном ниже примере в каждом элементе вводилось направление градиента решения (Л-направление) и нормальное к нему (^-направление) Численное решение аппроксимировалось полиномами по координатам 5 и Л7, причем аппроксимация по координате 5 была кусочно-постоянной, а по координате N - кусочно-квадратичной Три базисные функции были равны у,=1, 1'2=М, у3=Л'~ Таким образом, число неизвестных в таком подходе равно числу неизвестных в стандартном методе БС(1) На рис 5 приведено сравнение норм ошибок, полученных на равномерных сетках с помощью метода ОС(1) (кусочно-линейная аппроксимация), 00(2) (кусочно-квадратичная аппроксимация) и БСт - аппроксимация с помощью модифицированного базиса, описанного выше Точность, полученная с помощью метода Бвт, значительно выше, чем по методу 00(1) при одинаковом числе неизвестных на данной сетке и близка к точности, получаемой по методу Ой(2) Такой подход к выбору набора базисных функций можно рассматривать, как способ адаптации схемы к решению

Рис 5 Сравнение норм ошибок методов 00(1), 00(2), и ООт с адаптированным к решению набором базисных функций

В Главе III рассмотрены вопросы применения метода Галеркина с разрывными базисными функциями (00) к решению уравнений Эйлера -уравнений динамики идеального газа

§3 1 содержит изложение особенностей применения метода Бв к решению уравнений Эйлера для плоских течений Уравнения Эйлера формулируются относительно консервативных переменных в следующем виде

Э,и +V Р(и)=0

Здесь консервативные переменные и и декартовы составляющие Г(и) и g(u) функции потока Р(и) задаются в виде

р ' Р» ру

Р" р и' + р риу 2

ру р !(У Р V +Р

Ре. риИ ру/г

0)

Ь = е + р/ р, ¿> = (у-1)р(е-(г/2+у2)/2),

где р - плотность, и и V - компоненты скорости, р - давление, е - полная энергия единицы массы, Л - полная энтальпия единицы массы, у -отношение удельных теплоемкостей

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

В отличие от традиционных методов численного решения задач газовой динамики с помощью схем конечного объема, метод ОС оказался очень чувствительным к геометрическому описанию криволинейных границ обтекаемого тела В §3 2 представлены результаты расчета обтекания кругового цилиндра потоком идеального газа при числе Маха 0 1 Как обычно принято в методе конечного объема, криволинейная граница обтекаемого тела представлена полигоном (многоугольником)

Типичные результаты расчета распределения давления на поверхности цилиндра на сетке 128*64 изображены на рис 6а Здесь штриховой линией представлено точное решение для несжимаемой жидкости Сплошная линия соответствует решению, полученному с помощью стандартной, классической схемы конечного объема второго порядка точности, обозначенному как «схема КО 2-го порядка» Сплошная линия с кружочками соответствует расчету по исследуемой Б0(1) схеме с аппроксимацией решения с помощью кусочно-линейных функций Решение методом 00(1) носит осциллирующий характер Отметим, что это решение соответствует невязке, уменьшенной до машинного нуля, те полностью сошедшееся решение В классической схеме решение плавное и без осцилляций и близко к точному решению, в то время как решение с использованием БО схемы сильно осциллирующее Отметим, что мы имеем

большие ошибки в градиентах решений в каждом элементе, в то время как, в среднем, решение довольно близко к точному

Решение осциллирует только вблизи поверхности цилиндра На рис 6а приведено решение по 00(1) схеме на расстоянии 11=1 165 от центра цилиндра (радиус самого цилиндра равен 1) Видно, что в этом случае осцилляций нет Таким образом, можно предположить, что причина плохого решения по 00 схеме лежит в неадекватном представлении границы обтекаемого тела

Рис 6 Распределение давления на поверхности цилиндра при М=0 1 Сетка 128*64

Возможным способом решения этой проблемы является использование вблизи криволинейной границы т н параметрических элементов с одним криволинейным ребром Теоретические аспекты применения таких элементов содержатся в §3 3, а результаты расчетов - в §3 4 Результат расчета распределения давления по поверхности цилиндра на той же сетке, что и на рис 6а, приведен на рис 66 Видно, что применение параметрических элементов обеспечивает высокое качество решения Схема 00(1) с кусочно-линейной аппроксимацией решения обеспечивает также меньшую диссипацию, чем схема конечного объема второго порядка на той же расчетной сетке Использование схемы КО приводит к несимметричному обтеканию, являющемуся результатом возникшего отрыва потока в задней части цилиндра, в то время как картина течения, полученная по 00 схеме, выглядит более симметричной, в соответствии с точным решением

Были также проведены исследования порядка точности различных схем путем анализа поведения коэффициента сопротивления на последовательности сеток с удваивающимся количеством узлов в каждом направлении Результаты представлены на рис 7, где показана зависимость коэффициента сопротивления от количества неизвестных N в

логарифмическом масштабе для а) схемы конечного объема первого порядка (кусочно-постоянная реконструкция, схема Годунова), Ь) схемы конечного объема второго порядка (кусочно-линейная реконструкция), с) ОО(О) схемы (кусочно-постоянные элементы), с1) 00(1) схемы (кусочно-линейные элементы), е) 00(2) схемы (кусочно-квадратичные элементы). Для всех 00 схем расчеты были выполнены с квадратическими параметрическими элементами на границе. Как и ожидалось, обе схемы первого порядка демонстрируют первый порядок ошибки в сопротивлении. Обе схемы второго порядка дают очень близкие значения сопротивлений при одном и том же количестве неизвестных, в то время как порядок величины ошибки асимптотически стремится к 3. Схема Ов(2) дает порядок ошибки, асимптотически стремящийся к 5 (на рис.7 приведены треугольники, углы наклона гипотенузы которых соответствуют указанному порядку точности схемы). Такие высокие порядки величин могут быть следствием высокой регулярности используемых ссток. Следует отметить высокую точность метода 00(2) - величина коэффициента сопротивления составляет 10°, что весьма трудно получить при использовании схем меньшего порядка точности при разумном числе неизвестных задачи.

Рис.7. Коэффициент сопротивления в зависимости от количества неизвестных N в случае обтекания цилиндра при М=0.1

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

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

Ниже приведены два примера решения уравнений Эйлера на адаптивных сетках. На рис.8 представлено решение задачи о сверхзвуковом обтекании биплана Буземана при числе Маха невозмущенного потока М = 2.766. Начальная сетка около данной конфигурации состояла из 1741 ячеек. На рис. представлена сетка с 9224 ячейками и решение после 7 итераций адаптации сетки. Данный пример иллюстрирует возможность получения решения по методу ОС на весьма нерегулярных сетках, а также существенное улучшение качества решения (разрешения скачков и волн разрежения) в результате адаптации сетки к решению.

О 0 5 1 0 0.5 1

X X

Рис. 8. Сетка и решение (поле числа Маха) после 7 итераций адаптации сетки к решению для биплана Буземана. Число узлов сетки 4695, 9224 ячеек

Еще один пример расчета на адаптивных сетках представляет собой расчет обтекания профиля КАСА0012 при М=0.85, а=1.25°. Использовался метод ОС(1) с линейной реконструкцией решения и было выполнено 17 итерации адаптации сетки. Изменение коэффициента подъемной силы Су на последних итерациях не превышает 0.7%, коэффициента сопротивления Сх - 2%. На рис. 9 приведены сетка, распределение давления и поле чисел Маха при расчете на сетке, полученной в результате адаптации с числом

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

Рис. 9. Сетка после 17-ти итераций адаптации и решение (коэффициент давления и поле числа Маха) для профиля ЫАСА0012, М=0.85, а=1.25°

В Главе IV рассматривается решение уравнений Навье-Стокса для плоских ламинарных течений с использованием конечно-элементного метода Галеркипа (Бв) с разрывными базисными функциями.

В §4.1 изложены особенности применения метода БО к решению уравнений Навье-Стокса, которые записываются для вектора консервативных переменных и в следующей форме:

Э,и + V • Р;(и) = V • ^ (и).

Консервативные переменные и и декартовы составляющие ^(и) и §|(и) определяются согласно (1). Декартовы компоненты Ци) и функции вязкого потока Гч(и) задаются в виде

Г, =И

О

2их + \(их + V,.)

и[2и + Х(их + V,, ,)]+ у(ух + и )+—<?,

Рг

2Уг +\(Их + У1)

иЬх+к +Х(их + V )\+—ех

Рг

где д - коэффициент динамической вязкости, Рг - число Прандтля, Я=-2/3.

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

В §4 2 приведены примеры решения уравнений Навье-Стокса с помощью метода DG с различным порядком точности схемы Рассмотрены задачи обтекания плоской пластины при числе Рейнольдса Re=106 при малых скоростях М=0 1, кругового цилиндра при Re=40 при М=0 1, а также два случая ламинарного обтекания профиля NACA0012 для режимов М=0 5, а=0, Re=5000 и М=0 8, а=10°, Rc=73 В качестве примера на рис 10 приведены результаты расчета кругового цилиндра при Re=40, М=0 1 на достаточно редкой сетке 20*10 узлов при различных порядках базисных полиномов к в виде картины линий тока Данный режим обтекания по экспериментальным данным и результатам расчетов других авторов является стационарным и характеризуется наличием развитой отрывной зоны за цилиндром Представленные на рис 10 данные демонстрируют симметризацию картины течения и улучшение описания отрывной зоны по мере увеличения порядка базисных функций от 0 до 4

В таблице 1 приведено общее количество переменных (степеней свободы) при различных сочетаниях сеток и величин порядка базисных полиномов к В таблице 2 приведено угловое значение положения точки отрыва для различных сеток и различных значений к Здесь же дано экспериментальное значение этого угла 53 5°, приведенное в работе (Coutanceau, Bouard, 1977), а также расчетное значение из работы (Hauke, Hughes, 1998) Анализ результатов, представленных на рис 10, а также в таблицах 1, 2 показывает, что решения близкого качества и степени согласования с экспериментом могут быть получены на мелкой сетке 80*40 при А=1 и на крупной сетке 20*10 при к=3 или 4 В последних случаях общее число переменных в 3-5 раз меньше, чем в первом, что свидетельствует о преимуществах применения схем высокого порядка точности

В §4 3 приведен пример решения уравнений Навье-Стокса с использованием метода DG с адаптацией сетки к решению Методика адаптации неструктурированной сетки, описанная в главе 3, была применена к решению задачи вязкого обтекания профиля NACA0012 при двух режимах дозвуковое обтекание при числе Маха М = 0 5, а=0, Re=5000 и сверхзвуковое обтекание при числе Маха М=2 0, а=10°, Re=106 Расчеты данных случаев обтекания приведены также в работе (Bassi, Rebay, 1997) Использовался метод DG(1) с линейной реконструкцией решения Сетка адаптировалась к числу Маха На рис И приведены сетка и поле чисел Маха для второго из указанных случаев после 14 итераций адаптации сетки к решению Наблюдается, в частности, хорошее разрешение головной ударной волны

К=2 к=3

Рис. 10. Круговой цилиндр. Сетка 20* 10. Г1е=40, М=0.1. Линии тока

Таблица 1. Общее число переменных

сетка к=0 к=1 к=2 к=Ъ к=4

20*10 400 1200 2400 4000 6000

40*20 1600 4800 9600 16000

80*40 6400 19200 38400

Таблица 2. Круговой цилиндр, Яе=40. Угловое положение точки отрыва

сетка к= 0 k= 1 к= 2 А=3 к=4

20*10 36.0° 48.1° 54.0° 53.2° 53.4°

40*20 45.0° 55.6° 54.0° 53.7°

80*40 49.2° 54.4° 53.8°

Эксперимент (Coutanceau, Bouard, 1977) 53.5°

Расчет (Hauke, Hughes, 1998) 54.0°

Рис. 11. Сетка после 14 итераций адаптации и решение (поле числа Маха) для профиля КАСА0012, М=2.0, сс=10\ Яе=106

В Главе V рассматривается применение метода Галеркигга с разрывными базисными функциями к расчету вязких турбулентных плоских течений, описываемых уравнениями Навье-Стокса, осредненными по Рейнольдсу (уравнения Рейнольдса) Для замыкания системы уравнений Рейнольдса используется модель турбулентности Спаларта-Алмараса (8ра1аП, АНтагав, 1992)

В §5 1 описаны рассматриваемые уравнения Рейнольдса, кратко описана модель турбулентности Спаларта-Альмараса и особенности применения метода ОС в данном случае Уравнения Рейнольдса и уравнение Спаларта-Алмараса могут быть записаны совместно в консервативной форме Эта система уравнений имеет следующий вид

Э,и +V Р,(и) = V Р»+8 (2)

Консервативные переменные и и декартовы составляющие (".(и) и g,(u) невязкого потока Р,(и) задаются в виде

и =

Здесь, как и ранее, р - плотность, и и V - компоненты скорости, р -давление, е и И - полная энергия и энтальпия единицы массы, V - рабочая переменная в уравнении Спаларта - Алмараса, определяющая величину турбулентной вязкости

Декартовы компоненты £(и) и gv(u) вязкого потока ^(и) задаются в

виде

р ри ру

ри ри2+р рЫУ

ру риу ■ &=■ ру2 + р

ре риИ рук

Р7. ри\> рУУ

О

Т

ху

<7 6х

11Т}Х+УТ>Г+()У

Эу ду

Источниковый член Б в уравнении (2) содержит дополнительные слагаемые только в уравнении для турбулентной вязкости, связанные с производством, диффузией и диссипацией турбулентной вязкости, в соответствии с работой (8ра1аП, АНшагав, 1992)

Здесь компоненты тензора вязких напряжений л>хл„ и вектора теплового потока (}х, Qy определяются обычным образом в соответствии с гипотезой Буссинеска

, х/ 4 Эи zöv

2dv

зау

/ и Эг< 3v

я>

4dv

з а v

2 Эм

Jäc

= tU+M,)^-Ex+M, Рг

Рг,

J_ Рг

J_ Рг,

_1_ Рг

Составляющие градиента удельной внутренней энергии Е равны _ de Эи ^Эу _ Эе ^Эм * Эх Эх Эх г Эу Эу Эу Параметры модели турбулентности определяются в соответствии с рекомендациями авторов модели, которые в деталях воспроизведены в диссертации

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

В §5 2 приведены некоторые примеры решения уравнений Рейнольдса помощью метода DG Рассмотрены задачи турбулентного обтекания плоской пластины при числах Рейнольдса Re=106, 107 и малых скоростях М=0 1, кругового цилиндра при Re=3 6 10б и малых скоростях М=0 1, а также турбулентного обтекания профиля RAE2822 при М=0 725, а=2 92°, Re=6 5 1 06 и М=0 725, а=3 19°, Re=6 5 10б, соответствующих безотрывному и отрывному режимам обтекания

Как указывалось выше, одним из потенциальных преимуществ схемы DG является возможность варьирования порядка точности схемы в различных областях течения (изменение порядка полиномов, аппроксимирующих решение), а также адаптации схемы к решению с помощью изменения указанного порядка (p-refinement) Повышение порядка точности схемы будет наиболее эффективным в областях, где решение является гладким, но быстро меняющимся, например, в пограничных слоях В качестве первого шага в этом направлении рассмотрено численное решение задачи о вязком турбулентном обтекании профиля RAE2822 при условиях М=0 725, а=2 92°, Re=6 5 10б (CASE-6)

На рис 12 представлены профили скорости в пограничном слое в трех сечениях на верхней поверхности профиля Представлено 3 решения Решение №1 получено на грубой сетке из 3660 ячеек по схеме DG(1) 2-го

порядка точности (10980 неизвестных задачи) Данное решение не обладает достаточной точностью Решение №2 получено на измельченной сетке из 6448 ячеек по той же схеме (19334 неизвестных) Результаты проведенных расчетов на различных сетках позволяют утверждать, что это решение является достаточно точным Решение №3 получено на той же сетке, что и решение №1, однако в области пограничного слоя использовалась схема повышенного 3-го порядка Число неизвестных в этом случае равнялось 13869, те меньше, чем в расчете№2, а формы профилей скорости было получены с высокой точностью и практически совпадают с профилями скорости, полученными в расчете №2 Этот пример является иллюстрацией эффективности использования схемы повышенного порядка в областях гладкого быстро меняющегося решения, например в пограничном слое

у/5

3660 ячеек 10980 неизв схема 2-го порядка

--- 6448 ячеек, 19334

y/g неизв, схема 2-го порядка

3660 ячеек, 13896 неизв схема смешанного порядка

1

u/Ue

В

u/Ue

u/Ue

Рис 12 Профили скорости в пограничном слое

В качестве примера расчета обтекания тела сложной геометрии ниже приведены результаты расчета обтекания трехэлементного профиля корпорации Mcdonnell Douglas (MCD) Данный профиль представляет собой сечение крыла с отклоненными на 30° предкрылком и закрылком в посадочной конфигурации Этот профиль был подробно экспериментально исследован в малотурбулентной аэродинамической трубе НИЦ им Лэнгли HACA, и экспериментальные результаты, использованные ниже, взяты из работ [Chin, и др, 1993, Spaid, Lynch, 1996] Общий вид профиля представлен на рис 16 Эксперимент был проведен при числе М=0 2 и числе Re=910б Расчеты были проведены по схеме DG(1) с линейной реконструкцией решения в ячейках, с моделью турбулентности Спаларта-

Альмараса и применением методики адаптации сетки к решению, описанной выше.

На рис.13 приведены исходная (на первой итерации) и финальная (на седьмой итерации) расчетные сетки. Видно, что сгущение и обогащение сетки происходят в областях пограничных слоев и вязких следов. Результаты расчета поля числа Маха и его фрагментов приведены на рис.14. Обтекание характеризуется наличием заметных отрывов потока на нижней поверхности предкрылка и в нише основного профиля. Наблюдается хорошее разрешение вязких следов, сходящих с различных элементов профиля, что свидетельствует о достаточной эффективности процедуры адаптации. Данная процедура позволяет также разрешать и другие особенности течения, такие, например, как скачки уплотнения. Это иллюстрируется данными, представленными на рис.15, где приведены фрагменты сетки и поля чисел Маха вблизи передней кромки предкрылка при а=24°, что соответствует критическому углу атаки (максимальной подъемной силе). Данный режим характерен наличием малой местной сверхзвуковой зоны на передней кромке предкрылка, с максимальным местным числом Маха М=1.3 и замыкающим скачком уплотнения. Эта тонкая особенность течения также хорошо разрешена в результате адаптации сетки к решению.

Рис.13. Расчетные сетки около профиля МСЭ на первой и последней итерации адаптации сетки к решению

Рис.14. Поле числа Маха при обтекании профиля МСЭ при М=0.2, а=18°,

Ие=9-10(1

Рис.15. Фрагменты адаптивной сетки и поля чисел Маха вблизи передней кромки предкрылка профиля MCD при М=0.2, а=24°, Re=9106

Экспериментальные исследования обтекания профиля MCD включали измерения профиля полного давления в сечении х/с=0.975 при трех значениях угла атаки (а=8°, 16°, 21°). На рис.16 результаты измерения

коэффициента полного давления С „ = —сравниваются с

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

С практической точки зрения представляющей интерес характеристикой профиля с взлетно-посадочной механизацией является зависимость подъемной силы от угла атаки и величина максимальной подъемной силы. На рис.17 приведено сравнение этих зависимостей, полученных экспериментально [van Dam, 2002] и в расчете. Наклоны этих зависимостей на линейном участке и величины критического угла атаки согласуются хорошо. Разница в величине максимальной подъемной силы составляет около 2%.

А

коэффициента полного давления на закрылке профиля МСО при х/с=0 975,

а=8°, 16°, 21°

Рис 17 Сравнение расчетной и экспериментальной зависимостей коэффициента подъемной силы от угла атаки трехэлементного профиля

МСЭ

Для иллюстрации возможности применимости рассматриваемого метода к расчету пространственных течений, в работе приведен пример расчета трехмерного обтекания крыла ЬА^ при М = 0.82, Ке = 7.3-10 , а = 0.6°. Течение считалось полностью турбулентным. Расчетная сетка была построена с помощью промышленного пакета прикладных программ и она состояла из 190213 гексаэдральпых элементов. Расчет проводился по схеме СС( 1) с кусочно-линейной реконструкцией решения в ячейках, что дает общее количество неизвестных задачи равное 760852. На рис.18 представлено распределение коэффициента давления по верхней поверхности крыла и сравнение с экспериментальными данными для сечения г = 0.475. Наблюдается хорошее согласование результатов расчета с экспериментом.

Рис.18. Сравнение распределения коэффициента давления для крыла в сечении г - 0.475 с экспериментальными данными

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

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

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

- удобство работы с неструктурированными сетками;

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

- адаптация путем локальной модификации набора базисных функций;

- автоматическая адаптация сетки к решению.

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

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

3 Метод описан, протестирован на многочисленных примерах и подготовлен к практическому применению

ОСНОВНЫЕ РЕЗУЛЬТАТЫ ОПУБЛИКОВАНЫ В РАБОТАХ

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

2 А В Волков, С В Ляпунов Применение конечноэлементного метода Галеркина с разрывными базисными функциями к решению уравнений Рейнольдса па неструктурированных сетках Ученые записки ЦАГИ,

т XXXVIII, №3-4, 2007

3 S V Lyapunov, А V Wolkov 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

4 A V Wolkov, S V Lyapunov 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

5 N В Petrovskaya, A V Wolkov, S V Lyapunov Modification of basis functions in high order discontinuous Galerkin schemes for advection equations University of Birmingham, School of mathematics, Preprint 2006/26,2006

6 N В Petrovskaya, A V Wolkov, S V Lyapunov Modification of basis functions in high order discontinuous Galerkin schemes for advection equations Applied Mathemetical Modelling, Elsevier, 2006

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

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 А В Волков, С В Ляпунов, А Н Храбров Влияние установившегося вращения на аэродинамические характеристики профиля при наличии отрыва потока Ученые записки ЦАГИ, т XXXIV, №3-4, 2003 S V Lyapunov, А V Wolkov A separated flow calculation about airfoils and high-lift systems on basis of viscous-inviscid interaction methods EUROMECH, 3rd European Fluid Mechanics Conf, 1997 S.V Lyapunov, A V Wolkov Calculation of separated flows about airfoils and high-lift systems using viscous-inviscid interaction approach Proc Of the 5lh russian-chmise simposium on aerodynamics and flight dynamics, 1997 S V Lyapunov, A.V Wolkov Application of Viscous-inviscid Interaction Methods for a Separated Flow Calculation About Airfoils and High-Lift Systems 1CAS Proceedings, 1СAS-96-1 10 2,1996 S V Lyapunov, A V Wolkov Application of the Viscous-Invisid Interaction Model to Calculation of Two-Dimensional Separated Flows TsAGl Journal, vol 2, 1, 1996

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

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

S V Lyapunov Accelerated Method of the Euler Equation Solution in Transonic Airfoil Flow Problem ICAS Proceedings, ICAS-92-4 2 3, 1992 С В Ляпунов Ускоренный метод решения уравнений Эйлера в задаче о трансзвуковом обтекании профиля Мат моделирование, N4,1991

5 V Lyapunov Convergence acceleration and wave drag determination m transonic airfoil calculations. ICAS Proc | ICAS-90-6 9 2, 1990

N A Vladimirova, V V Vyshinsky, S V Lyapunov, Ya M Serebnisky Convergence acceleration of computational methods for two- and three-dimensional flows about bodies Fluid Mechanics, v 2, 1986 H А Владимирова, В В Вышинский, С В Ляпунов, Я М Серебрийский

06 ускорении сходимости методов расчета плоского и пространственного трансзвукового обтекания тел в неограниченном потоке Ученые записки ЦАГИ, т 16, №4,1985

М А Брутян, С В Ляпунов О второй вариации функционала Рябушинского ДАН СССР, т 258, №4, 1981

22 В В Вышинский, С В Ляпунов Расчет околозвукового осесимметричного обтекания мотогондолы с учетом вязкости Ученые записки ЦАГИ, т 19, №3, 1988

23 А В Волков, В С Сакович, А М Сорокин, Ю Б Лифшиц, С В Ляпунов Применение неструктурированных сеток в вычислительной аэродинамике Материалы IX школы-семинара "Аэродинамика летательных аппаратов", ЦАГИ, 1998

24 А В Волков, С В Ляпунов Применение ориентированных схем второго порядка для решения двумерных уравнений Эйлера на неструктурированных сетках Материалы X школы-семинара "Аэродинамика летательных аппаратов", ЦАГИ, 1999

25 А В Волков, С В Ляпунов Исследование схемы аппроксимации высокого порядка уравнений газовой динамики. Материалы XI школы-семинара "Аэродинамика летательных аппаратов", ЦАГИ, 2000

26 А В Волков, С В Ляпунов Исследование схемы аппроксимации высокого порядка уравнений газовой динамики Современные проблемы аэрокосмической науки и техники Международная научно-техническая конференция молодых ученых и специалистов, 2000

27 А В Волков, С.ВЛяпунов Применение схемы высокого порядка точности для решения двумерных уравнений Эйлера и Навье-Стокса на неструктурированных сетках Материалы XII школы-семинара "Аэродинамика летательных аппаратов", ЦАГИ, 2001

28 А В Волков, С В Ляпунов Численный метод высокого порядка точности для решения системы уравнений Навье-Стокса и Рейнольдса Материалы XIII школы-семинара "Аэродинамика летательных аппаратов", ЦАГИ, 2002

29 А В Волков, С В Ляпунов Разработка новой монотонной схемы высокого порядка точности на компактном шаблоне для решения уравнений Эйлера и Навье-Стокса Материалы XVI школы-семинара "Аэродинамика летательных аппаратов", ЦАГИ, 2003

30 А В Волков, С В Ляпунов Применение метода конечных элементов к расчету вязких турбулентных течений на неструктурированных адаптивных сетках Материалы XV школы-семинара "Аэродинамика летательных аппаратов", ЦАГИ, 2004

31 А В Волков, С В Ляпунов Разработка новой монотонной схемы высокого порядка точности на компактном шаблоне для решения уравнений Эйлера и Навье-Стокса Материалы XVI школы-семинара "Аэродинамика летательных аппаратов", ЦАГИ, 2005

32 А В Волков, С В Ляпунов О решении акустической задачи методом Галеркина с разрывными базисными функциями Материалы XVII школы-семинара "Аэродинамика летательных аппаратов", ЦАГИ, 2006

Отпечатано в типографии ЦАГИ

Подписано в печать 03 04 2008 Тираж 100 экз Заказ № 5237

ВВЕДЕНИЕ

ПРИНЯТЫЕ ОБОЗНАЧЕНИЯ

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

1.1. Общая теория метода Галеркина с разрывными функциями для решения одномерных законов сохранения

1.2. Одномерное уравнение конвекции 21 1.2.1 Постановка задачи 21 1.2.2. Модельная задача и анализ порядка численной ошибки

1.3. Одномерное уравнение теплопроводности 25 1.3.1 Постановка задачи 25 1.3.2. Модельная задача и анализ порядка численной ошибки

1.4. Уравнение Бюргерса с вязкостью 28 1.4.1 Точное решение уравнения Бюргерса 28 1.4.2. Анализ порядка численной ошибки

Течения невязких и вязких газов описываются, соответственно, системами уравнений Эйлера и Навье-Стокса. Эти системы уравнений представляют собой законы сохранения массы, импульса и энергии. Аналитическое решение этих уравнений возможно лишь в небольшом числе случаев, как правило, далеких от практики. Единственным общим средством решения этих уравнений являются численные методы. Эта область аэродинамики и прикладной математики столь обширна, что выделяется в отдельную ветвь науки, которая называется "вычислительная аэродинамика" (в англоязычной литературе - CFD - Computational Fluid Dynamics).

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

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

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

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

Большая часть упомянутых проблем в настоящее время нашла то или иное решение. За последние 25-30 лет в области разработки и применения численных методов решения уравнений газовой динамики был достигнут колоссальный прогресс. Обзоры успехов и проблем в области вычислительной аэродинамики можно найти в многочисленных публикациях (см., например Sloof, Schmidt, 1994; LeVeque, 1992; Fletcher, 1988; Куликовский, Погорелов, Семенов , 2001; Роуч, 1980). Автор диссертации в течение ряда лет также работал в данной области и полученные им результаты опубликованы, в частности, в работах (Брутян, Ляпунов, 1981, Владимирова, Вышинский, Ляпунов, Серебрийский, 1985, Волков, Ляпунов, 1993, Волков, Ляпунов, Храбров, 2003, Ляпунов, 1991, 1993, Lyapunov, Wolkov, 1994, 1996а,б, 1997а,б, Lyapunov, 1990, 1992).

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

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

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

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

• При аппроксимации на пространственной сетке в дискретных выражениях явно или неявно присутствуют слагаемые, называемые "искусственная (схемная) вязкость", которые обеспечивают устойчивость схемы. В идеальном случае эти слагаемые должны быть минимально возможными для обеспечения устойчивости. Наиболее популярными способами внесения таких слагаемых являются явный подход для центральноразностных схем (Jameson, Schmidt, Türkei, 1981; Türkei, 1988) и неявный подход с использованием схем с разностями против потока (Годунов, 1959; Steger, Warming, 1981; Van Leer, 1982; Roe, 1981).

• Первые используемые численные схемы обеспечивали первый порядок аппроксимации уравнений на гладких решениях (напр. Годунов, 1959). Такие схемы обладали избыточными диссипативными свойствами и слишком сильно размазывали разрывы решения (скачки уплотнения). В дальнейшем наиболее популярными стали схемы, обеспечивающие второй порядок аппроксимации на гладких решениях (Lax, Wendroff, 1964; Колган, 1972, Родионов, 1987).

• Для численных схем второго порядка точности и выше, в отличие от схем первого порядка, характерны нефизичные осцилляции на разрывах. Для их устранения требуется в окрестности разрывов переходить к монотонным схемам первого порядка точности, т.е. схемы должны иметь переменный порядок точности. Схемы такого типа получили название "схемы высокого разрешения" (high resolution schemes, Harten, 1983; Yee, 1989). Обычно это обеспечивается путем ограничения производных при реконструкции решения (limiters).

Пространственная сетка, используемая для дискретизации задачи, представляет собой систему областей, называемых обычно ячейками, которые покрывают всю область и не пересекаются, Иногда понятие сетки формулируется другим образом (например, система точек — узлов сетки или сетка с перекрывающимися элементами), однако в данной работе будут рассматриваться сетки в указанной выше формулировке. Обычно ячейки сетки представляют собой многоугольники в плоском случае (далее именуемом 20) и многогранники в пространственном случае (30).

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

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

Структурированные и неструктурированные сетки имеют свои преимущества и недостатки. Неструктурированные сетки требуют хранения информации о связях (соседях ячеек, гранях ребрах и вершинах ячеек) и использования специальных процедур доступа к этой информации. С этой точки зрения они менее эффективны, чем структурированные сетки, Обычно считается, что «накладные расходы» при использовании неструктурированных сеток по памяти и быстродействию составляют 20-40%. С другой стороны, генерация структурированной сетки около тел сложной геометрии — весьма непростая задача. Обычным способом решения этой задачи является применение блочно-структурированных сеток, когда область вне тела разбивается на блоки простой топологии, внутри каждого блока генерируется структурированная сетка, и (желательно) узлы сеток внутри блоков «сшиваются» на границах блоков. Для тел сложной геометрии разбиение области на блоки в значительной мере является искусством и требует обычно довольно много времени. Генерация неструктурированных сеток в меньшей степени опирается на искусство вычислителя, и эта процедура может быть более формализована.

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

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

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

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

Метод конечного объема является в настоящее время наиболее распространенным. Он основан на интегральном представлении законов сохранения и, поэтому, применим к построению разрывных решений. В соответствии с этим подходом, неизвестными величинами являются средние значения физических переменных внутри ячеек. Решение в каждой ячейке реконструируется с использованием указанных неизвестных средних значений, и дискретные уравнения получаются путем записи законов сохранения для каждой ячейки для реконструированного решения. Поскольку на границах между ячейками реконструкция решения является разрывной, для описания потоков используются точные или приближенные методы решения задачи Римана о распаде произвольного разрыва. Простейшим случаем реализации такой схемы является кусочно-постоянная реконструкция, что обеспечивает первый порядок схемы по пространству (Годунов, 1957, 1959). Схема второго порядка получается при кусочно-линейной реконструкции (Колган, 1972, Родионов, 1987). Как было указано выше, схемы высокого порядка могут приводить к нефизичным осцилляциям численного решения. Т.н. монотонные схемы свободны от этого недостатка. В работе (Годунов, 1959) доказано, что не существует монотонной линейной схемы порядка точности выше первого. Монотонизация схем высокого порядка (или обеспечение свойства TVD — невозрастания полной вариации) осуществляется обычно путем введения ограничителей на градиенты реконструкции решения в ячейках. Впервые эта процедура при моделировании двумерных течений была предложена в работе (Колган, 1972). Схемы такого типа, называемые в зарубежной литературе MUSCL - (аббревиатура от слов «монотонные ориентированные схемы для законов сохранения») были подробно рассмотрены в работах (van Leer, 1973, 1974, 1977а, 19776, 1979). По существу этот подход представляет собой способ уменьшения порядка точности схемы до первого в областях разрывов решения. Еще более высокий порядок точности получается при реконструкции решения с помощью полиномов второго (Collela, Woodward, 1984) и более высоких порядков в т.н. схемах ENO (Harten et al, 1987, Harten, Osher, 1987).

Третий способ дискретизации законов сохранения возможен с помощью методов конечного элемента. В соответствии с ним, численное решение аппроксимируется линейной комбинацией непрерывных функций (или имеющих требуемый порядок гладкости, в соответствии с порядком рассматриваемых уравнений) с неизвестными коэффициентами, т.е. решение принадлежит линейному подпространству т.н. базисных функций. Обычно базисные функции имеют компактный носитель, т.е. они равны нулю вне ограниченной области, состоящей из одной или нескольких ячеек и представляющей собой собственно конечный элемент. В большинстве методов конечного элемента в качестве базисных функций используются полиномы от координат. Далее возможны различные способы получения уравнений для определения коэффициентов при базисных функциях. Одним из таких часто используемых методов является метод Галеркина, который также называют методом Бубнова-Галеркина (Михлин, Смолицкий, 1965). В соответствии с ним дискретные уравнения представляют собой условия ортогональности невязки базисным функциям в функциональном пространстве Ьг, записанные в т.н. слабой или галеркинской формулировке, получаемой интегрированием указанного условия ортогональности по частям. При этом требования к гладкости базисных функций являются минимальными. Если невязка аппроксимируется также с помощью базисных функций, то такой подход называют собственно методом Галеркина. Если для аппроксимации невязки используется другое подпространство функций, то говорят о методе Петрова-Галеркина.

Классический метод Галеркина достаточно хорошо математически обоснован и широко применяется для решения эллиптических и параболических задач. Применение его к гиперболическим задачам, к которым, в частности, относятся уравнения Эйлера динамики идеального газа, сталкивается с определенными трудностями, такими как неустойчивость итерационных процедур решения задач с начальными данными. Источник этих трудностей заключается в том, что дискретные уравнения, получаемые с помощью метода Галеркина, являются «центрированными», т.е. они не учитывают скорость и направление распространения физических возмущений в гиперболических задачах. Иными словами, в метод необходимо добавить диссипативные члены, аналогично тому, как это было описано выше для конечно-разностных и конечно-объемных методов. Имеются различные способы добавления диссипативных членов в конечноэлементные методы. Одним из таких методов является т.н. метод SUPG (аббревиатура от английского названия Streamline Upwind Petrov-Galerkin - метод Петрова-Галеркина с несимметрией вдоль линий тока). Этот метод также иногда называют SD — Streamline Diffusion - метод с диффузией вдоль линий тока. Этот метод был предложен в работе (Hughes, Brooks 1979) для решения линейных задач конвекции-диффузии и в дальнейшем был развит и применен для решения более сложных уравнений вплоть до уравнений Навье-Стокса. Обзор вопроса и обширный список ссылок можно найти в (Johnson, 1992). По существу, это классический метод Галеркина с непрерывными базисными функциями, в котором пробные функции - функции, на которые скалярно умножаются уравнения движения для получения дискретных уравнений — отличаются от базисных функций, используемых для аппроксимации невязки на величину аддитивных слагаемых, обеспечивающих необходимую диссипацию.

Альтернативный подход к адаптации методов конечного элемента к гиперболическим задачам заключается в использовании базисных функций, заданных внутри конечного элемента и терпящих разрыв на его границе (т.н. несогласованные элементы). Диссипация в этом методе вводится в точности в соответствии с рецептами схем типа Годунова - путем решения задачи Римана о распаде разрыва на границе элемента. По существу такой метод представляет собой симбиоз методов конечного элемента и конечного объема. Данный подход в англоязычной литературе получил название DG (Discontinuous Galerkin - метод Галеркина с разрывными функциями). В дальнейшем оба названия — метод DG и метод Галеркина с разрывными функциями будут использоваться как синонимы. Впервые этот метод был предложен в работе (Reed, Hill . 1973) для решения уравнения переноса нейтронов. Дальнейшие многочисленные исследования были посвящены анализу математических аспектов метода, таких как скорость сходимости и пр. (LeSaint, Raviart, 1974, Johnson, Pitkaranta, 1986), а также развитию метода (Cockburn, Shu, 1989, Cockburn, Lin, Shu, 1989, Cockburn, Hou, Shu, 1990, Cockburn, Shu, 1998) и применению его к решению более сложных уравнений, таких как уравнения Эйлера (Bassi, Rebay 1997а) и Навье-Стокса (Bassi, Rebay 1997b, Warburton, Lomtev, Kirby, Karniadakis 1998, Lomtev, Karniadakis 1997).

Преимущества метода Галеркина с разрывными функциями заключаются в следующем.

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

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

Остановимся подробнее на двух последних аспектах. Если ошибка численного решения равна О(Ик), где к — шаг сетки, то число к называется порядком ошибки. Очевидно, что, начиная с некоторого числа неизвестных дискретной задачи, схемы высокого порядка будут обеспечивать меньшую ошибку решения, по сравнению со схемами низкого порядка, что делает такие схемы весьма привлекательными. Исследованию схем высокого порядка точности, в частности, посвящены работы (Толстых, 1990, 2000, То1з1укЬ, Ырауэкп, 1998). В схемах конечного объема повышение порядка обеспечивается повышением порядка полиномиальной аппроксимации решения внутри конечного объема. Поскольку в схемах конечного объема неизвестными являются средние значения решения внутри конечных объемов, повышение порядка аппроксимации требует расширения шаблона (например, в двумерных задачах на треугольных конечных объемах линейная аппроксимация требует шаблона из трех элементов, в качестве которых можно использовать соседние элементы, а для квадратичной аппроксимации требуется шаблон из шести элементов). Выбор шаблона в схеме конечного объема высокого порядка не очевиден и не однозначен, особенно на неструктурированных сетках. В схеме БО коэффициенты полиномиальной аппроксимации решения внутри элемента являются неизвестными, и формулировка дискретной задачи определения этих коэффициентов не требует определения какого-либо шаблона. В записи дискретных уравнений в невязком случае участвуют только элементы, соседние с данным элементом, а в вязком случае — также и соседи второго порядка (соседи соседей). Подробнее этот вопрос рассмотрен в содержании диссертации на конкретных примерах.

Выше была указана важность адаптации сетки к решению. Это означает, что сетка должна содержать большее количество элементов меньшего размера в областях сильного изменения решения и меньшее количество узлов в областях, где решение меняется слабо. Важна и форма ячеек, в частности в пограничных слоях ячейки должны быть вытянуты' вдоль слоев. Операции адаптации сетки путем сгущения и разрежения ячеек в англоязычной литературе получили название Ь-гейпетеп1. Возможен и другой способ адаптации, когда при неизменной сетке в зависимости от характера решения меняется порядок точности численной схемы. Такой способ адаптации называется p-refinement. Наконец, возможно и сочетание этих способов (hp-refinement).

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

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

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

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

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

Содержание диссертации изложено в пяти главах.

В первой главе приведены результаты применения метода Галеркина с разрывными функциями к решению одномерных модельных линейных и нелинейных задач. Рассмотрены линейное уравнение конвекции (переноса), теплопроводности, нелинейное уравнение Бюргерса. Исследован порядок ошибки схемы БЭ в зависимости от порядка базисных полиномов.

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

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

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

Наконец, в пятой главе приведены примеры применения метода БО к решеншо уравнений Рейнольдса, описывающих турбулентное течение, с моделью турбулентности Спаларта-Алмараса (8ра1аг1:, АПтагаэ, 1992). Приведены примеры расчета турбулентного обтекания пластинки, кругового цилиндра и профиля при больших числах Рейнольдса.

Основные результаты диссертации опубликованы в работах (Sakovich, Sorokin, Wolkov, Lyapunov, 1998, Волков, Сорокин, Сакович, Лифшиц, Ляпунов, 1998, Волков, Ляпунов, 1999, 2000а,б, 2001, 2002, 2003, 2004, 2005, 2006а,б, Lyapunov, Wolkov, 2000, Petrovskaya, Wolkov, Lyapunov, 2006a,б, Wolkov, Lyapunov, 2000).

Результаты работы обсуждались на на 6-ой международной конференции по генерации сеток для вычислительной аэродинамики (1998), 16-м Конгрессе Международной ассоциации математического и компьютерного моделирования IMACS (2000), 3-ей Европейской конференции по механике жидкости EUROMECH (1997), 5-ой Российско-Китайской конференции по аэродинамике и механике полета (1997), Международных Конгрессах по авиационным наукам ICAS (1990, 1992, 1996), школах-семинарах «Аэродинамика летательных аппаратов (1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007), Международной конференции молодых ученых и специалистов «Современные проблемы аэрокосмической науки и техники (2000), франко-российском семинаре ONERA-ЦАГИ (2002) и др.

В данную диссертацию включены исследования, поддержанные РФФИ (Проекты № 98-01-00032-а, 1998, №00-01-00070-а, 2000, №02-01 -00124-а, 2002, №03-01-00236-а, 2003, №06-01-00283-а, 2006).

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

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

ПРИНЯТЫЕ ОБОЗНАЧЕНИЯ

А Матрица системы линейных уравнений а Скорость звука

С Симметричная блочно-диагональная матрица в уравнении (3.4)

С, cu, С\2, С21, С22 Диссипативная матрица и ее элементы сьь Cb2> Cvb Cwi, cW2, cw3, а, Полуэмпирические константы в модели турбулентности к, Г],щ Спаларта-Алмараса

Cf Коэффициент сопротивления трения

CFL Число Куранта

Ср Коэффициент давления

Сх Коэффициент аэродинамического сопротивления

Схр Коэффициент сопротивления давления

Су Коэффициент подъемной силы

D Слагаемое, связанное с диссипацией турбулентной вязкости в модели турбулентности Спаларта-Алмараса d Расстояние до поверхности обтекаемого тела в пограничном слое е Полная энергия единицы массы

Е Удельная внутренняя энергия e¡ Длина i-ой стороны элемента (треугольника)

F Вектор невязких потоков в уравнениях Эйлера f Вектор правых частей системы линейных уравнений f(u) Функция невязкого потока в одномерном уравнении конвекциидиффузии f,g; fj,gj Составляющие вектора невязкого потока в уравнениях Эйлера и

Навье-Стокса

F¡ Вектор невязких потоков в уравнениях Навье-Стокса

Fv Вектор вязких потоков в уравнениях Навье-Стокса fv,gv Составляющие вектора вязкого потока в уравнениях Навье

Стокса fvi, fV2, fw, S, г Вспомогательные функции в модели турбулентности Спаларта

Алмараса h Шаг сетки по пространству h Полная энтальпия единицы массы hu, hq Функции невязких потоков величин и и q в уравнении конвекции-диффузии huq, hqq Функции вязких потоков величин и и q в уравнении конвекциидиффузии

Ij Интервалы разбиения отрезка в одномерной задаче

J Матрица Якоби невязки системы дискретных уравнений

К Число базисных функций

L Невязка m Размерность подпространства Крылова

М Матрица предобуславливателя для системы линейных уравнений М Число Маха

N Число интервалов (элементов) сетки в одномерном уравнении п Вектор единичной нормали к ребру элемента

N Размерность системы линейных уравнений р давление

Р Слагаемое, связанное с производством турбулентной вязкости в модели турбулентности Спаларта-Алмараса Давление торможения Полином Лежандра степени / Число Прандтля Турбулентное число Прандтля

Коэффициенты разложения вспомогательной функции д по базисным функциям Компоненты вектора теплового потока Вспомогательная переменная в одномерном уравнении конвекции-диффузии (уравнение (1.2)) Модуль скорости

Вектор невязки системы линейных уравнений Число Рейнольдса

Вектор невязки системы дискретных уравнений Число Рейнольдса по координате х Площадь ячейки

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

Матрица собственных векторов якобиана вектора потоков в уравнениях Эйлера

Температура

Время

Собственные вектора якобиана вектора потоков в уравнениях Эйлера (столбцы матрицы Т)

Неизвестная функция в одномерном и двумерном уравнениях конвекции-диффузии

Вектор консервативных переменных в уравнениях Эйлера

Скорость, нормальная к разрыву (границе элемента)

Коэффициенты разложения решения и по базисным функциям

Скорость набегающего потока

Компоненты вектора скорости

Базисные и пробные функции в методе Галеркина

Вектора подпространства Крылова

Вектор характеристических переменных \у=Т"1и

Вектор неизвестных системы линейных уравнений

Декартовы ортогональные координаты

Вектор поправки к решению системы линейных уравнений

Вспомогательная переменная в уравнении Спаларта — Алмараса, определяющая величину турбулентной вязкости

Компоненты тензора вязких напряжений

Коэффициент диффузии в уравнении конвекции-диффузии, коэффициент (динамической) вязкости в уравнениях Навье

Стокса

Угол атаки

Плотность у отношение удельных теплоемкостей у Вспомогательная переменная в уравнении Спаларта - Алмараса.

V Коэффициент молекулярной кинематической вязкости

Параметрические координаты а,р,к Степень полинома

61 Толщина вытеснения пограничного слоя а,/, Р,/ Параметры схемы Рунге-Кутта

Лк Собственные значения матрицы Якоби вектора потоков в уравнениях Эйлера Шаг по времени

V, Коэффициент турбулентной кинематической вязкости

Коэффициент (динамической) турбулентной вязкости в уравнениях Навье-Стокса <|Хз Составляющие вектора единичной нормали к ребру элемента

Нижние индексы:

Производная по времени х,у Производные по координатам к Сеточные функции

О Начальные значения г, у Номер интервала сетки л'Л Положение разрыва функции с Координаты центра тяжести элемента оо условия в набегающем потоке;

К, Ь Значения на границе ячейки (где функция терпит разрыв) слева и справа boundary Значения на непротекаемой границе Верхние индексы:

А Численные значения потоков " Значения на границе ячейки (где функция терпит разрыв) слева и справа

Номер базисной функции тонн Точное решение

Осреднение по Роу п Номер итерации по времени

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

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

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

- удобство работы с неструктурированными сетками;

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

- адаптация путем локальной модификации набора базисных функций;

- автоматическая адаптация сетки к решению.

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

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

3. Метод описан, протестирован на многочисленных примерах и подготовлен к практическому применению.

ЗАКЛЮЧЕНИЕ ПО ГЛАВЕ V

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

Результаты расчетов турбулентного вязкого обтекания плоской пластины при числах f\ 7

Рейнольдса Re=10 -10 хорошо согласуются с эмпирическими данными Никурадзе, как по профилю скорости, так и по величине коэффициента местного трения. Ошибка в последнем случае составляет около 1%.

Результаты расчета распределения коэффициента давления при турбулентном обтекании кругового цилиндра при числе Рейнольдса Re=3.6-106 удовлетворительно согласуются с экспериментальными данными Achenbach. Результаты расчета турбулентного трансзвукового обтекания профиля RAE2822 на типовых режимах, используемых для верификации численных методов (CASE-6 и CASE-10) также удовлетворительно согласуются с экспериментальными данными по характеру распределения давления и положению скачка уплотнения.

Применение схемы повышенного порядка в области пограничного слоя при расчете обтекания профиля RAE2822 (CASE-6) обеспечивает получение решения с достаточной точностью на довольно грубой сетке при меньшем числе неизвестных задачи, чем в случае решения этой же задачи с использованием схемы второго порядка при сравнимой точности решения. Этот пример свидетельствует о высоких потенциальных возможностях рассматриваемой схемы с точки зрения адаптации к решению путем изменения порядка точности (p-refmement).

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

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

1. М.А.Брутян, С.В.Ляпунов. О второй вариации функционала Рябушинского. ДАН СССР, т.258, №4, 1981

2. Н.А.Владимирова, В.В.Вышинский, С.В.Ляпунов, Я.М.Серебрийский Об ускорении сходимости методов расчета плоского и пространственного трансзвукового обтекания тел в неограниченном потоке. Ученые записки ЦАГИ, т. 16, №4, 1985

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

4. А.В.Волков, В.С.Сакович, А.М.Сорокин, Ю.Б.Лифшиц, С.В.Ляпунов. Применение неструктурированных сеток в вычислительной аэродинамике. Материалы IX школы-семинара "Аэродинамика летательных аппаратов", ЦАГИ, 1998.

5. А.В.Волков, С.В.Ляпунов. Применение ориентированных схем второго порядка для решения двумерных уравнений Эйлера на неструктурированных сетках. Материалы X школы-семинара "Аэродинамика летательных аппаратов", ЦАГИ, 1999.

6. А.В.Волков, С.В.Ляпунов. Исследование схемы аппроксимации высокого порядка уравнений газовой динамики. Материалы XI школы-семинара "Аэродинамика летательных аппаратов", ЦАГИ, 2000.

7. А.В.Волков, С.В.Ляпунов. Исследование схемы аппроксимации высокого порядка уравнений газовой динамики. Современные проблемы аэрокосмической науки и техники. Международная научно-техническая конференция молодых ученых и специалистов, 2000.

8. А.В.Волков, С.В.Ляпунов. Применение схемы высокого порядка точности для решения двумерных уравнений Эйлера и Навье-Стокса на неструктурированных сетках. Материалы XII школы-семинара "Аэродинамика летательных аппаратов", ЦАГИ, 2001.

9. А.В.Волков, С.В.Ляпунов. Численный метод высокого порядка точности для решения системы уравнений Навье-Стокса и Рейнольдса. Материалы XIII школы-семинара "Аэродинамика летательных аппаратов", ЦАГИ, 2002.

10. А.В.Волков, С.В.Ляпунов. Разработка новой монотонной схемы высокого порядка точности на компактном шаблоне для решения уравнений Эйлера и

11. Навье-Стокса. Материалы XVI школы-семинара "Аэродинамика летательных аппаратов", ЦАГИ, 2003.

12. А.В.Волков, С.В.Ляпунов, А.Н.Храбров. Влияние установившегося вращения на аэродинамические характеристики профиля при наличии отрыва потока. Ученые записки ЦАГИ, t.XXXIV, №3-4,2003

13. А.В.Волков, С.В.Ляпунов. Применение метода конечных элементов к расчету вязких турбулентных течений на неструктурированных адаптивных сетках. Материалы XV школы-семинара "Аэродинамика летательных аппаратов", ЦАГИ, 2004.

14. А.В.Волков, С.В.Ляпунов. Разработка новой монотонной схемы высокого порядка точности на компактном шаблоне для решения уравнений Эйлера и Навье-Стокса. Материалы XVI школы-семинара "Аэродинамика летательных аппаратов", ЦАГИ, 2005.

15. А.В.Волков, С.В.Ляпунов. О решении акустической задачи методом Галеркина с разрывными базисными функциями. Материалы XVII школы-семинара "Аэродинамика летательных аппаратов", ЦАГИ, 2006.

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

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

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

19. С.К.Годунов. Разностный метод численного расчета разрывных решений уравнений гидродинамики. Математический сборник, т.47(89), №3, 1959.

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

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

22. С.В.Ляпунов. Неплоские крылья минимального индуктивного сопротивления. Изв.РАН, МЖГ, №2,1993

23. С.В.Ляпунов. Ускоренный метод решения уравнений Эйлера в задаче о трансзвуковом обтекании профиля. Мат.моделирование, N4, 1991

24. Э.Митчелл, Р.Уэйт. Метод конечных элементов для уравнений с частными производными, М., "Мир", 1981.

25. С.Г.Михлин, Х.Л.Смолицкий. Приближенные методы решения дифференциальных и интегральных уравнений. М., Наука, 1965

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

27. П.Роуч. Вычислительная гидродинамика. М., Мир, 1980.

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

29. Г.Шлихтинг. Теория пограничного слоя. М., Наука, 1974.

30. E.Achenbach. Distribution of Local Pressure and Skin Friction Around a Circular Cylinder in Cross-Flow up to Re=5*106. J. Fluid Mech., Vol. 34, part 4, 1968.

31. Test Cases for Inviscid Flowfield Methods, AGARD AR-211, June, 1986.

32. H.L.Atkins, C.W.Shu. Quadrature-free implementation of Discontinuous Galerkin methods for hyperbolic equations, Technical Report 96-51, ICASE, 1996,

33. T.J.Barth. 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, 1992, pp.6-1 6-61

34. F.Bassi, S.Rebay. High-Order Accurate Discontinuous Finite Element Solution of the 2D Euler Equations. J. of Сотр. Phys. 138, 251-285, 1997

35. F.Bassi, S.Rebay. High-Order Accurate Discontinuous Finite Element Method for the Numerical Solution of the Compressible Navier-Stokes Equations. J. of Сотр. Phys. 131, 267-279, 1997

36. V.D.Chin, D.W.Peters, F.W.Spaid, RJ.McGhee. Flowfield Measurements About A Multi-Element Airfoil At High Reynolds Numbers. AIAA-93-3137.

37. B.Cockburn, C.W.Shu. TVB Runge-Kutta Local Projection Discontinuous Galerkin Finite Element Method for Scalar Conservation Laws II: General Framework. Math. Сотр. 52:411-435,1989.

38. D.Cockburn, C.W.Shu. The Runge-Kutta Discontinuous Galerkin Finite Element Method for Conservation Laws V: Multidimensional Systems, J. Сотр. Phys., 141:199-224, 1998e

39. D.Cockburn, C.W.Shu. TVB Runge-Kutta Local Projection Discontinuous Galerkin Finite Element Method for Scalar Conservation Laws II: General Framework, Math. Comp., 52:411-435, 1989

40. D.Cockburn, S.Hou, C.W.Shu. TVB Runge-Kutta Local Projection Discontinuous Galerkin Finite Element Method for Conservation Laws IV: The Multidimensional Case, Math. Comp., 54:545-581, 1990

41. D.Cockburn, S.Y.Lin, C.W.Shu. TVB Runge-Kutta Local Projection Discontinuous Galerkin Finite Element Method for Conservation Laws III: One Dimensional Systems, J. Comp. Phys., 84:90-113, 1989

42. P.Collela, P.Woodward. The Piecewise-parabolic Method (PPM) for Gasdynamic Equations. J.Comput.Phys. 54, pp. 174-201, 1984.

43. Computational Aerodynamics based on the Euler Equations. Ed. by J.F.Sloof, W.Schmidt, AGARD-AG-325, 1994.

44. P.H.Cook, M.A.McDonald, M.C.P.Firmin. Aerofoil RAE 2822 Pressure distributions, and boundary layer and wake measurements. AGARD-AR-138, 1979.

45. M.Coutanceau, R.Bouard, Experimental Determination of the Viscous Flow in a Wake of a Circular Cylinder in Uniform Translation. Part I. Steady Flow. J.Fluid Mech. 79 pp.231-256,1977.

46. C.A.J.Fletcher. Computational Techniques for Fluid Dynamics, Springer-Verlag, 1988.

47. A.Harten, B.Engquist, S.Osher, S.Chakravarthy. Uniformly High Order Accurate Essentially Nonoscillatory Schemes III, J.Comput.Phys. 71, pp.231-303, 1987

48. A.Harten, S.Osher. Uniformly High Order Accurate Nonoscillatory Schemes I, SLAM J.Num.Anal. 24, pp.279-309, 1987

49. A.Harten. High Resolution Schemes for Hyperbolic Conservation Laws. J. Comp. Phys., 49, No.3, 357-393, 1983.

50. G.Hauke, T.J.R.Hughes. A Comparative Study of Different Sets of Variables for Solving Compressible and Incompressible Flows. Comp. Meth. Appl. Mech. Engrg. 153,1-44, 1998.

51. T.J.R Hughes, F.Brooks. A Multidimensional Upwind Scheme with no Crosswind Diffusion, AMD vol 34, Finite Element Methods for Convection Dominated Flows, ed. T.J.R.Hughes, ASME New York, 1979

52. A.Jameson, W.Schmidt, E.Turkel. Numerical Solution of the Euler Equations by Finite-Volume Methods Using Runge-Kutta Time-Stepping Schemes. AIAA 811259,1981.

53. C.Johnson, J.Pitkaranta. An Analysis of the Discontinuous Galerkin Method for a Scalar Hyperbolic Equation. Math. Comp., 46:1-26, 1986

54. C.Johnson. Finite Element Methods for Flow Problems. In Special Course on Unstructured Grid Methods for Advection Dominated Flows, AGARD-R-787, 1992

55. P.D.Lax, B.Wenroff. Difference Schemes for Hyperbolic Equations with High Order of Accuracy, Comm. Pure Appl. Math., 17,No.3, 381-398, 1964.

56. P.LeSaint, P.A.Raviart. On a Finite Element Method for Solving the Neutron Transport Equation, C. de Boor, editor, Mathematical Aspects of Final Elements in Partial Differential Equations, pp. 89-145, Academic Press, 1974

57. RJ.LeVeque. Numerical methods for Conservation Laws. Lectures in Mathemetics, Birkhauser Verlag, Basel-Boston-Berlin, 1992.

58. Lomtev, G.E.Karniadakis. Simulations of Viscous Supersonic Flows on Unstructured hp-meshes. AIAA-97-0754, 1997.

59. S.V.Lyapunov, A.V.Wolkov. 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.

60. 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

61. S.V.Lyapunov, A.V.Wolkov. A separated flow calculation about airfoils and high-lift systems on basis of viscous-inviscid interaction methods. EUROMECH, 3rd European Fluid Mechanics Conf., 1997

62. S.V.Lyapunov, A.V.Wolkov. Application of the Viscous-Invisid Interaction Model to Calculation of Two-Dimensional Separated Flows. TsAGI Journal, vol.2,1,1996

63. S.V.Lyapunov, A.V.Wolkov. Calculation of separated flows about airfoils and high-lift systems using viscous-inviscid interaction approach. Proc. Of the 5th russian-chinise simposium on aerodynamics and flight dynamics, 1997

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

65. S.V.Lyapunov. Accelerated Method of the Euler Equation Solution in Transonic Airfoil Flow Problem. ICAS Proceedings, ICAS-92-4.2.3,1992

66. S.V.Lyapunov. Convergence acceleration and wave drag determination in transonic airfoil calculations. ICAS Proc. \ ICAS-90-6.9.2, 1990

67. L.Martinelly. Calculations of Viscous Flows with a Multigrid Method. PhD Thesis, Dept. of Mechanical and Aerospace Engineering, Princeton University, 1987.

68. D.J.Mavriplis, A.Jameson. Multigrid Solution of the Navier-Stokes Equations on Triangular Meshes. AIAA J., 28, No.8,1415,1990.

69. N.B.Petrovskaya, A.V.Wolkov, S.V.Lyapunov. Modification of basis functions in high order discontinuous Galerkin schemes for advection equations. University of Birmingham, School of mathematics, Preprint 2006/26, 2006.

70. N.B.Petrovskaya, A.V.Wolkov, S.V.Lyapunov. Modification of basis functions in high order discontinuous Galerkin schemes for advection equations. Applied Mathemetical Modelling, Elsevier, 2006.

71. W.H.Reed, T.R.Hill. Triangular Mesh Methods for the Neutron Transport Equation. TR LA-UR-73-479, Los Alamos Scientific Laboratory, 1973

72. P.Roe. Approximate Riemann Solvers, Parametric Vectors and Difference Schemes. J. of Comp. Phys., vol.43, pp.357-372, 1981.

73. Y.Saad, M.H.Schultz. GMRES: A Generalized Minimal Residual Algorithm for Solving Nonsymmetric Linear Systems. SIAM J.Stat.Comput. vol.7, No.3, 1986.

74. Y.Saad. Iterative Methods for Sparse Linear Systems. PWS Publishing Company, Boston, 1996.

75. 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.

76. F.W.Spaid, F.T.Lynch. High Reynolds Number, Multi-Element Airfoil Flowfield Measurements. AIAA-96-0682.

77. P.R.Spalart, S.R.Allmaras. A One-Equation Turbulent Model for Aerodynamic Flows. AIAA-92-0439.

78. J.L. Steger, R.F. Wanning. Flux Vector Splitting of the Invicsid Gasdynamic Equations with Application to Finite-Difference Methods, J.Comp.Phys, vol.40, pp.263-293, April, 1981.

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

80. E.Turkel. Improving the Accuracy of Central Difference Schemes. In 11th International Conference on Numerical Methods in Fluid Dynamics, SpringerVerlag, Lecture Notes in Physics, vol.323, pp.586-591, 1988.

81. C.P.van Dam. The aerodynamic design of multi-element high-lift systems for transport airplanes. Progress in Aerospace Sciences 38 (2002) 101-144.

82. B. van Leer. Towards the Ultimate Conservative Difference Scheme I, The Quest of Monotonicity, Lect.Notes.Phys. 18, pp.163-168,1973.

83. B. van Leer. Towards the Ultimate Conservative Difference Scheme II, Monotonicity and Conservation combinrd in a Second Order Scheme, J.Comput.Phys. 14, pp.361-370, 1974

84. B. van Leer. Towards the Ultimate Conservative Difference Scheme III, Upstream-Centered Finite-Difference Schemes for Ideal Compressible Flow, J.Comput.Phys. 23, pp.263-275, 1977

85. B. van Leer. Towards the Ultimate Conservative Difference Scheme IV, A New Approach to Numerical Convection, J.Comput.Phys. 23, pp.276-299, 1977

86. B. van Leer. Towards the Ultimate Conservative Difference Scheme V, A Second Order Sequel to Godunov's Method, J.Comput.Phys. 32, pp.101-136, 1979

87. B.van Leer. Flux Vector Splitting for the Euler Equations. Lecture Notes in Physics, vol. 170, pp.501-512, 1982.

88. C.P.van Dam. The aerodynamic design of multi-element high-lifit systems for transport airplanes. Progress in Aerospace Sciences 38 (2002) 101-144.

89. N.A.Vladimirova, V.V.Vyshinsky, S.V.Lyapunov, Ya.M.Serebriisky Convergence acceleration of computational methods for two- and three-dimensional flows about bodies. Fluid Mechanics, v.2, 1986

90. T.C.Warburton, I.Lomtev, R.M.Kirby, G.E.Karniadakis. A Discontinuous Galerkin Method for the Navier-Stokes Equations in Hybrid Grids. In M.Hafez and J.C.Heirich, ed., 10th Interantional Conference on Finite Elements in Fluids, Tucson, 1998.

91. L.Wigton, N.Yu, D.Young. GMRES Acceleration of Computational Fluid Dynamic Codes, AIAA 85-1494, 1985.

92. H.C.Yee. A Class of High-Resolution Explicit and Implicit Shock-Capturing Method, von Karman Institute Lecture Series 1989-04 (NASA TM-101088), 1989.

93. R.J.Zwaan, LANN wing, Pitching Oscillation, AGARD Report 702, Compendium of Unsteady Aerodynamic Measurements Addendum No 1, Data Set 9, May 1985