Монотонизация метода конечного элемента в задачах газовой динамики Текст научной статьи по специальности «Физика»

Том XL

УЧЕНЫЕ ЗАПИСКИ ЦАГИ

2009

№ 4

УДК 533.6.011.35:629.7.025.73 519.6

МОНОТОНИЗАЦИЯ МЕТОДА КОНЕЧНОГО ЭЛЕМЕНТА В ЗАДАЧАХ ГАЗОВОЙ ДИНАМИКИ

А. В. ВОЛКОВ, С. В. ЛЯПУНОВ

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

Ключевые слова: метод Галеркина с разрывными базисными функциями, метод конечного объема, схема высокого порядка точности, ограничители градиентов решения.

Растущие потребности высокоточного компьютерного моделирования приводят к необходимости использования расчетных сеток с очень большим числом ячеек. Сокращение объема расчетных сеток при неизменной точности получаемых численных решений возможно при повышении порядка точности численной схемы. В традиционном методе конечного объема, как правило, ограничиваются вторым порядком точности. Увеличение порядка точности приводит к расширению шаблона аппроксимации и усложнению расчетной схемы. Одним из известных альтернативных способов построения численных схем высокого порядка точности является метод Галеркина с разрывными базисными функциями (РМГ — разрывный метод Галеркина, в англоязычной литературе — DG — Discontinuous Galerkin) (например, [1, 2]). Метод является результатом комбинации метода конечного элемента и конечного объема, объединяя их преимущества в одном подходе. В частности, РМГ дает возможность локального увеличения или уменьшения порядка точности схемы.

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

В методе конечного объема для монотонизации решения используют ограничители градиентов решения (которые в дальнейшем в соответствии с англоязычным термином называются лимитерами, а сама процедура — лимитированием). Это дополнительный коэффициент X в формуле реконструкции решения в ячейке:

U(x) = U0 + XVU(x -x0).

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

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

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

к

и(х) = ^изФз (х)

з=1

где коэффициенты и^ — являются искомыми переменными задачи, К — общее количество

линейно-независимых базисных функций.

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

полиномов Лежандра в промежутке от —1 до 1 ф = 1, Ф2 = х, фз =—(3х2 — 1) и т. д. и равняется

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

к

и(х) = и ^изфз (х).

з=2

При этом происходит ограничение вклада высоких гармоник в решение, вплоть до использования кусочно-постоянного восполнения решения (при X = 0) с первым порядком точности. Впервые такой подход был применен в работе [1] при использовании явной схемы интегрирования по времени Рунге — Кутта.

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

С целью устранения сингулярности матрицы Якоби в настоящей работе предложена так называемая РМГ-Х (лямбда) схема. Ее основная идея состоит в использовании специфических базисных функций, которые получаются путем комбинации с коэффициентом X базисных функций высоких гармоник и кусочно-постоянных базисных функций с коэффициентом (1 — X). Схема может быть представлена следующим образом:

РМГ (К, X) = ХРМГ (К) + (1 — X) РМГ (0).

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

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

В работах [6, 7] рассматривается принципиально другой способ монотонизации решения, основанный на локальном введении искусственной вязкости в численную схему. Величина

дополнительной вязкости, добавляемой к уравнениям в «проблемных» ячейках, порядка

к (1 -X)

что обеспечивает получение гладкого решения. Здесь к — величина, отражающая

О

К

размер ячейки. Отметим, что протяженность скачка уплотнения при таком способе монотонизации решения меньше размера одной ячейки. Более того, ширина скачка уменьшается при увеличении порядка точности схемы как 1/ К. Подчеркивая особенность данного подхода, авторы используют термин «подсеточные лимитеры». В отличие от методов локального понижения порядка точности этот подход не требует сгущения сетки в областях особенностей течения. Однако особенность его реализации предполагает решение дополнительного уравнения для величины искусственной вязкости с дополнительным набором эмпирических констант [7]. Кроме того, здесь также необходимы методы определения областей сингулярностей течения.

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

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

1. Описание схемы РМГ (К, X). Рассмотрим уравнения законов сохранения в консервативной форме:

^^(К= (1.1)

ОТ

Здесь К = К (и) — невязкий поток, К = К (и, Уи) — вязкий поток.

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

Р2 (0,1) и Р3 (0,0). Любой треугольный элемент физического пространства преобразуется

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

Для реализации РМГ-Х схемы каждый сеточный элемент разбивается на множество под-элементов Е = и и • • • % , и с каждым подэлементом ассоциируется только одна узловая точка.

Здесь N — общее количество подэлементов. Конкретные разбиения треугольных элементов на подэлементы для случаев кусочно-линейной и кусочно-квадратичной реконструкций приведены на рис. 1. Здесь же приведены узловые точки, использованные при выборе базисных функций. В случае, когда максимальный порядок полинома в элементе ограничен первым порядком (кусочно-линейная реконструкция), треугольник разбивается на два треугольных и один четырехугольный элемент (N = 3, рис. 1, а) и в данном элементе вводятся три линейные базисные функции:

¥1 = X ¥2 = ^ ¥з =1 - х - У.

2

2

6

3

13 4 5

4

а)

б)

Рис. 1. Разбиение треугольника в параметрическом пространстве для РМГ-Х схемы:

а — К = 1; б — К = 2

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

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

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

Базисные функции фг- схемы РМГ(К, X) определяются как взвешенная сумма кусочно-линейных (или квадратичных) и кусочно-постоянных базисных функций:

Решение в каждом элементе Е представляется как линейная комбинация базисных функций:

порядка базисного полинома.

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

у = -х + 2 х2, у2 = -у + 2 у2, Уз = 1 - Зх - 3 у + 2 х2 + 2у2 + 4ху,

у = 4х - 4х2 - 4ху, = 4ху, у6 = 4у - 4у2 - 4ху.

Ф/ = ХУ/ +(-Х)Ф/

(1.2)

количество базисных функций в зависимости от К — максимального

/ = 1,2, к, Кг.

Здесь £ — скалярное произведение длины ребра подэлемента и нормали п (пх, пу, п2)

к ребру и О — площадь подэлемента.

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

Креб d £ = 1 (к £ Я ) £ -2 Л| (и *-и Я

(1.4)

Здесь верхние символы «£» и «Я» означают, что соответствующая величина решения вычислена в левой или правой стороне ребра, и матрица Л есть якобиан невязких потоков:

Л =

э( • & £)

эи

(1.5)

Вязкие потоки зависят от градиентов независимых переменных К = К (и, Уи). В методе РМГ градиент решения также представляется как линейная комбинация базисных функций:

К/

эи

(, X) = Е()Ф] (х),

ох> ]=1

(1.6)

где 1 = 1, 2 соответствует координатам х и у. После умножения уравнения (1.6) на тестовую функцию и интегрирования по частям по каждому из подэлементов, получается следующая система для определения коэффициентов разложения gi в (1.6):

N

Е

вк =1

Г КГ д )

| Еgl,]Ь£Щи(• &£)- | и&О :

V Ови ]=1

£ви

к = 1, к, Кг,

Ови

дфк

Эх,-

= 0,

(1.7)

где в, — ,-я компонента единичного вектора.

Величина вязкого потока на границах элементов в (1.3) определяется из левой по отношению к рассматриваемому ребру ячейки:

(1.8)

реб = К (', Уи'

а значение переменной и в выражении (1.7) определяется из правой ячейки:

и1реб ^ .

(1.9)

Такое чередование при определении этих величин на ребрах ячеек было предложено в работе [1].

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

■ = -М

-1

N

е £ % ( -к )&£- i уф] (к-к )&о- £ о

к £Е

О Е

О Е

(1.10)

где М — массовая матрица:

М/,, =

{ДЕ Ф/ Ф 0

2. Численная оценка порядка точности схемы в задачах конвекции-диффузии. Оценки порядка точности схемы были выполнены на примерах решения задач конвекции и конвекции-диффузии, имеющих гладкое решение. Уравнение конвекции имело вид:

— + 0.625соз Г 5 (х + 0.1)) — = 0. дх П ду

(2.1)

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

(2.2)

1.25п(-0.3 + у - 0.12581П Г5(х + 0.1)]) Другая задача отличалась присутствием диффузионного члена в правой части уравнения:

ди г_/ Л1\-|ди д2и ..

— + 0.625собI 5(х + 0.1)1— = |1—-. (2.3)

дх -1 ду ду

-2

Коэффициент вязкости при этом был равен Д = 1 -10 . Точное решение, накладываемое в овия на границе втек

и(х у) = егТ£ = ■

качестве граничного условия на границе втекания, представляется в виде функции ошибок:

у - 0.6 - 0.125б1п Г5(х + 0.1)]

^4ц(х + 0.1)

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

Рис. 2. Ошибки схемы аппроксимации на последовательности измельченных сеток в зависимости

от параметра X:

а — решение уравнения конвекции; б — решение уравнения конвекции-диффузии

ности называется величина q в асимптотической оценке нормы ошибки численного решения вида О (к11). Здесь ошибка в норме ¿2 представлена в зависимости от числа степеней свободы (неизвестных) задачи Поскольку величина ^ пропорциональна к-2, порядок точности схемы равен удвоенному тангенсу угла наклона зависимостей, приведенных на рис. 2. В отличие от традиционной схемы РМГ при определенном выборе X удалось получить не целый порядок точности. Особенно ясно это наблюдается при решении конвективной задачи (рис. 2, а). Проведенный ряд численных экспериментов позволил вывести простую эмпирическую формулу для порядка

точности в зависимости от выбора параметра X. Порядок точности схемы РМГ (К, Х = 1 -кК -1+р) есть q = К + р, при этом величина р может изменяться от нуля до единицы. Так, например, при РМГ (К = 1, Х = 1 - к12) получен порядок 1.5, а при РМГ (К = 2, Х = 1 - к3/2) — схема с порядком 2.5.

В задаче с диффузией эти коэффициенты не обеспечили определенного порядка точности. Из рис. 2, б видно, что высокий порядок на грубых сетках переходит в низкий на более мелких сетках. Однако получение промежуточных порядков не является целью настоящей работы. Важно, что сконструирована схема, в которой плавное изменение параметра X приводит не к дискретному, а к плавному изменению порядка ошибки от монотонной схемы первого порядка точности при X = 0 до схемы с порядком К +1 при X =1.

3. Сенсоры сингулярности решения и определение величины X. В основополагающей работе [1], посвященной РМГ, лимитирование решения осуществлялось при использовании явной схемы интегрирования по времени Рунге — Кутта. Алгоритм выбора лимитеров, как и в методе конечного объема, основан на сравнении решения в рассматриваемой ячейке и соседних с ней. Применение лимитеров приводит к ослаблению вклада высоких гармоник в решение вплоть до их обнуления в ячейках с локальными экстремумами решения. Однако такой алгоритм приводит к нежелательному снижению порядка точности схемы в некоторых областях с гладким решением. Сужение области применения лимитеров возможно при использовании дополнительного сенсора сингулярностей, который бы более четко указывал зоны возможных немонотонностей решения, где и происходит лимитирование. Ниже приводится сравнение двух алгоритмов выбора сенсора сингулярностей решения.

Первый подход базируется на оценке величины разрыва решения на границах между элементами [3—5]. Так в [3] демонстрируется, что в случае гладкого решения величина разрыва

решения между элементами в РМГ пропорциональна О(кК+1), где К — максимальная степень

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

Таким образом, сенсор ^ опирается на вычисление безразмерной величины разрыва решения на границах элемента:

£ (( - ипЬг| ) I

5 -—-. (3.1)

1 ||ЭО|| кК+1 ^ }

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

2 + 2

Величина энтальпии Н = уЕ — (у — 1)—^— в каждой ячейке может быть представлена как

Кг

линейная комбинация базисных функций: Н = ^ а1 фг-. Коэффициенты этого разложения а1 мо-

г=1

гут быть вычислены следующим образом:

аг = М-1 Л Н ф ^П.

Используя значения энтальпии на границе рассматриваемой ячейки, можно получить аль-

К/

тернативное разложение: Н = ^ Ь1 фг-. Здесь коэффициенты этого разложения могут быть найде-

г=1

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

К/ (а — Ь Л2

«2 = (3.2)

1=1 У а1 )

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

1, если s2 < 50 — к,

Х =

— ^0 )

1 + 81П—---

2к

^ ' ' 2 '"' если 50 — к<52 <50 + к, (3.3)

2

\ / 0, если s2 > 50 + к.

Здесь = ^ «2, ^0 и к — некоторые настраиваемые параметры, которые, как показала практика, достаточно консервативны по отношению как к выбору решаемой задачи, так и к типу расчетной сетки.

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

[—1, если у — 0.6 — 0.1258т [х + 0.1] < 0,

"точи = ^ (3.4)

[+1, иначе.

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

После решения вышеописанных задач конвекции в каждой ячейке расчетной сетки выполнялся расчет величины сенсора. Отметим, что для расчета сенсора «2 вычислялись полная энтальпия в уравнениях Эйлера или Навье — Стокса. При численных экспериментах с уравнением конвекции энтальпия была заменена на некоторую сохраняемую в данном решении искусственную величину, а именно, в задаче с гладким решением Н = и2 +(1 — и^^), а в задаче

со скачком Н = и 2.

Рис. 3. Пример сетки, адаптивной к решению модельной задачи

Рис. 4. Величина сенсоров сингулярности решения (а — ^ б — ^2) в задачах с гладким и разрывным решениями

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

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

Рис. 4, а соответствует сенсору 51, рис. 4, б — сенсору 52. Здесь черные линии являются результатом для РМГ (К =1), а серые для РМГ (К=2). Сплошные линии соответствуют гладкому решению, пунктирные — разрывному решению. Прямоугольные маркеры обозначают решение на равномерной сетке, а треугольные маркеры — решения, полученные на адаптивных сетках с ячейками произвольной формы.

Из представленных зависимостей видно, что сенсор 51 (см. рис. 4, а) не позволяет отличить гладкое решение на адаптивной сетке от разрывного, полученного на равномерной сетке.

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

Таким образом, сенсор 52, основанный на постоянстве определенной величины, продемонстрировал преимущество над сенсором 51 и все последующие расчеты были выполнены именно с ним.

4. Пример решения уравнения конвекции с разрывом. Уравнение конвекции с разрывом решения (2.1), (3.4) решено с использованием описанной выше процедуры лимитирования (сенсор 52) на последовательности адаптивных сеток. При вычислении величины X были выбраны следующие параметры: ¿о = -15, к =1. При рассмотрении многих задач (уравнения конвекции-диффузии, уравнения Эйлера и Навье — Стокса) эти параметры были неизменными, так как именно такой выбор обеспечил высокий порядок точности в гладких областях решения и необходимую монотонность решения в областях скачков уплотнения.

Рис. 5. Поведение решения в сечении сеткиX = 0.99 с/без использования лимитирования при К= 1, 2: а — расчет на начальной сетке — нулевая итерация адаптации; б — расчет на адаптивной сетке после пяти итераций адаптации

Результаты расчетов, полученные на начальном и конечном этапах адаптации, демонстрируются на рис. 5. Здесь сравниваются результаты расчета, полученные с использованием процедуры лимитирования (сплошная линия) и без (пунктирная линия) при РМГ K = 1 (черная линия) и РМГ K =2 (серая). Величина решения представлена в сечении расчетной сетки х = 0.99. Отметим, что лимитирование решения в процессе итераций адаптации привело к устранению осцилляций решения, однако наблюдается локальное понижение порядка точности схемы в достаточно узком слое вблизи разрыва решения, что является неотъемлемой платой за монотонность.

5. Результаты расчета обтекания изолированного профиля невязким транс- и сверхзвуковым потоком сжимаемого газа. Результаты решения уравнений Эйлера представлены для случаев обтекания транс- и сверхзвуковым потоком изолированного профиля RAE 2822 под углом атаки а = 3°. Распределение давления для различных чисел Маха набегающего потока в диапазоне от трансзвуковых значений M = 0.7 до сверхзвуковых M = 2 представлено на рис. 6. Видно, что для всех рассмотренных режимов решение не имеет каких-либо забросов.

Разработанная РМГ-Х схема изучалась также на примере расчета известного тестового случая профиля NACA 12 при числе M = 0.8 и угле атаки а = 1.25° на последовательности неструктурированных адаптивных сеток. Метод и программа адаптации сетки к решению были разработаны в [8]. Процесс адаптации основан на эвристическом индикаторе ошибки, зависящем

Рис. 6. Примеры решений уравнения Эйлера с использованием схемы РМГ (К = 2, X): распределение давления

при разных числах Маха набегающего потока

Рис. 7. Расчет на последовательности адаптивных сеток. Зависимость коэффициента подъемной силы от количества степеней свободы расчетных схем РМГ K = 1,2:

a — расчет невязкого обтекания профиля NACA 12 при M = 0.8 и угле атаки а = 1.25°; б — расчет турбулентного обтекания профиля RAE 2822 при M = 0.725, Re = 6.5 • 106

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

Величины коэффициента подъемной силы, полученные в процессе адаптации на разных сетках при К =1 и 2, показаны на рис. 7. Видно, что при превышении определенного числа степеней свободы (число ячеек X число степеней свободы в ячейке > 40 000) оба подхода дают приблизительно одинаковые результаты. Как и предполагалось, наличие протяженного скачка уплотнения в расчетной области в целом заметно снижает порядок точности решения с использованием обеих схем, и как следствие, приводит к одинаковым результатам. Однако в задачах, где сингулярность течения незначительна, а основная картина обтекания формируется плавным решением, следует ожидать преимущества от использования схем высокого порядка точности. В приведенных примерах демонстрируется лишь возможность надежного получения результата со схемой высокого порядка даже в задачах с существенными скачками уплотнения.

ООО Решение AGARD AR-138

ср ....... Решение РМГ

/ . <Т ■ О- . / Т-О-о. © "С / р 0 i.

ó 5 lili ¡ i i

] 0.5 . 1

х/с

Рис. 8. Результаты расчета с использованием РМГ (K = 1, X) схемы на последней адаптивной сетке (16-я итерация адаптации — 39 001 элемент):

профиль NACA 0012, M = 0.8, а = 1.25°

Распределения давления, полученное с РМГ K = 1 после 16 итераций адаптации, представлено на рис. 8.

6. Результаты расчета профиля при обтекании трансзвуковым турбулентным потоком. Серия расчетов была выполнена также для известного тестового случая Case-6 [9] обтекания профиля RAE 2822 при M = 0.725, Re = 6.5 106. Экспериментальный угол атаки составляет 2.92°, расчетный угол был выбран из условия равенства экспериментального и расчетного коэффициентов подъемной силы cy = 0.738. Расчеты производились на последовательности адаптивных сеток с использованием как кусочно-линейной РМГ K =1, так и кусочно-квадратичной РМГ K =2 реконструкции решения в ячейке.

Рис. 9. Финальная расчетная сетка, полученная в процессе расчета: а — РМГ К = 1; б — К =2

Рис. 10. Распределение давления: — экспериментальные результаты [9] а = 2.92°, cy = 0.738, cx = 0.0127; — расчет РМГ K =1 а = 2.485°, cy = 0.738, cx = 0.0151

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

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

Авторы выражают признательность CFD лаборатории компании Боинг и лично Форстеру Джонсону (Forrester Johnson) за полезные дискуссии при выполнении настоящих исследований.

Работа частично выполнена при финансовой поддержке Российского фонда фундаментальных исследований (проект № 06-01-00283a).

ЛИТЕРАТУРА

1. Cockburn B., Shu C. W. TVB Runge — Kutte local projection Discontinuous Galerkin methods for scalar conservation laws II/ General framework// Math. Comp. 1989. V. 52.

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

3. KrivodonovaL., Xin J.,Remacle J.-F., Chevaugeon N.,Flaherty J. E. Shock detection and limiting with Discontinuous Galerkin methods for hyperbolic conservation laws // Applied Numerical Mathematics. 2004. V. 48.

4. Luo H. and Baum J. D., Lohner R. A fast, p-multigrid Discontinuous Galerkin method for compressible flows at all speeds // AIAA 2006-110.

5. Luo H. and Baum J. D., Lohner R. A hermite WENO-based limiter for Discontinuous Galerkin method on unstructured grids // AIAA 2007-510.

6. Persson P.-O. and P e r a i r e J. Sub-cell shock capturing for Discontinuous Galerkin method // AIAA-2006-112.

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

8. Мартынов А. А., Медведев С. Ю. Надежный способ построения сеток с вытянутыми ячейками / Рабочий семинар «Построение расчетных сеток: теория и приложения». — Москва, 24—28 июня 2002 г., Вычислительный центр РАН им. А. А. Дородницына, http://www.ccas.ru/gridgen/ggta02/papers/Medvedev.pdf.

9. Cook P. H., McDonald M. A., &Firmin M. C. P. Aerofoil RAE 2822 — pressure distribution, and boundary layer and wake measurements // AGARD-AR-138. 1979.

Рукопись поступила 9/VI2008 г.