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

Вычислительные технологии

Том 5, № 4, 2000

ЛОКАЛЬНАЯ МНОГОСЕТОЧНАЯ ТЕХНОЛОГИЯ РЕШЕНИЯ ЗАДАЧ ТРАНСЗВУКОВОЙ АЭРОДИНАМИКИ *

В. Б. КАРАМЫШЕВ Институт вычислительных технологий СО РАН Новосибирск, Россия e-mail: kary@net.ict.nsc.ru О. В. КУКАРЦЕВА Новосибирский государственный университет, Россия

This paper presents an effective, fast converging method for calculating 3-D stationary transonic flows. One of its features is the combination of a local multigrid technique with the algorithm of accelerating internal iterative processes by using the method of least squares.

Одной из проблем компьютерного моделирования стационарных трансзвуковых течений идеального газа является медленная скорость сходимости итерационных методов решения систем сеточных уравнений. Смешанный тип уравнений газовой динамики (эллиптический — в дозвуковых зонах, гиперболический — в сверхзвуковых) и обширные пространственные расчетные области приводят к необходимости решения плохо обусловленных систем конечно-разностных уравнений большой размерности. Одним из эффективных методов ускорения сходимости релаксационных процессов является многосеточная технология, которая в настоящее время находит все более широкое применение в вычислительной аэродинамике. Многосеточная техника ускорения итераций была предложена Р. Федорен-ко [7], затем, в работе [14], она была распространена на методы расчета трансзвуковых потенциальных течений. В наши дни многосеточная технология широко используется в алгоритмах численного решения уравнений Эйлера и Навье — Стокса. На ее основе, в сочетании с техникой блочной генерации сеток [5], разработаны эффективные локальные многосеточные методы (например, [4]).

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

* Работа выполнена по молодежному гранту СО РАН, при финансовой поддержке Российского фонда фундаментальных исследований, грант №99-01-00619, и докладывалась на Научной сессии Президиума СО РАН 15.04.99.

© В. Б. Карамышев, О. В. Кукарцева, 2000.

1. Исходные уравнения

Уравнение для потенциала скоростей в консервативной форме имеет вид

д д д дХ (рфх) + ^ (рфу) + ^ (рф*) = 0, (1.1)

где р — плотность, и = фХ, V = фу, г = фг — декартовы компоненты скорости. Плотность р и давление р определяются из интеграла Бернулли и условия изэнтропичности течения:

р

1 + ^мс (1 - фХ - Ф2 - Ф2)

7-1

Р

7 МС

1 + I-1 мС(1 - фХ - Ф2 - Ф2)

2

7-1

Здесь Мс — число Маха набегающего потока, а 7 — показатель адиабаты. Уравнение (1.1) записано относительно безразмерных величин

р

Р

и

Рс

и и

с

V

и с

г

го и

с

Р=

Ри

рси С

X У

х = У У = У

г

7'

где I — характерный масштаб (например, наибольшая полуось обтекаемого эллипсоида), р, и, V, и т. д. — размерные величины. Введем преобразование координат:

С = С(х'У,г)'

П = п(х, У, г),

С = С (х'У,г)-

В криволинейных координатах уравнение (1.1) примет вид:

Ш)+Ш)+Ш=0.

Здесь [/, V, Ж — контрвариантные составляющие скорости:

и = Лф? + Лфп + ^5фС, V = Лф? + Лфп + А6фс, Ж = А5 ф? + Абфп + Аэфс,

где

А1 = & + С2 + £, а2 = п2 + п2 + п2,

Аз = С2 + Су2 + С2, А = схпх + Су Пу + ^ п^

А5 схсх + СУ Су + С-г С

Аб = ПхСх + Пу Су + п* С*

(1.2)

а J = Сх(пуСг - п*Су) + п Х(С^ Су Су Сг ) + Сх(Суп* - С*пу) — якобиан преобразования.

1

1

V

г

Воспользуемся следующими очевидными тождествами:

Сх = J (Уп^С - УС^) ' Пх = J (Ус*? — У?^С) ' Сх = J (У? ¿п - Уп*?)' Су = J (х ^п — Хп ^), Пу = J (ж? — 2? ), Су = J (Хп 2? — ж? ^п),

& = J (жп Ус — ЖС Уп) ' П^ = J (жс У? — ж? УС) ' С* = J (ж? Уп — жп У?) • Тогда декартовы компоненты скорости будут вычисляться по формулам:

фх = Схф? + Пх Фп + СхФс'

фу = Суф? + Пу фп + Су фС > ф = &ф? + п* фп + С* фс •

2. Конечно-разностная схема

Для решения уравнения (1.2) воспользуемся консервативной конечно-разностной схемой

¿? (ри) + </ (р1) + ¿с (рЖ) = 0, (2.1)

где

Р

Р J'

¿?(ри) = Р-+1 3— р*-1 31 1

3+2 1

7-2 7-2 ] ! 7-2

¿п(р1) = рк+1 ^4+1 —рк-1 2

1

2 Л+4

2 Л-2

¿с ) = р;+ 1 — — р;_ 12

'+ 2

2 г- 2 1 1- 2

Заметим, что без нарушения общности можно положить ¿С = ¿П = ¿С = 1-

Компоненты скорости [/3+1, 2, Ж+1 в нецелых узлах вычисляются по формулам:

1

1

31 = Л,7+2 (ф3 + 1 — ф7 ) + 2(А4>3+^пф3 + 1 + ¿пф7 ) + 2 (А5ф7+1 + А5,7¿Сф7 )

2

14+1 = ^2^+2 (фк+1 — фк) + ^А^^? фк+1 + ¿? фк) + ^(Аб^+^с фк+1 + ¿с фк)'

1 = Аз,1+ 2 (фг+1 — фг) + 1(А5)|+1£?фг+1 + А5Д¿?фг) + ^(Аб^+^пфг+1 + А6,г£пфг),

где ¿?фз = (фз+1 — фз-1)/2, ¿пфк = (фк+1 — фк-1)/2, ¿сфг = (фг+1 — фг-1)/2 — центральные аппроксимации второго порядка производных потенциала вдоль направлений С, П, С -

Для стабилизации схемы в сверхзвуковых областях, где уравнение (1-2) имеет гиперболический тип, значения плотности р* на гранях расчетной ячейки задаются выражениями

Р,*+1 = (1 — 32 )р7-+1 + 32

1 + в

1 — в

2 "р7+г + 2— р7-1+3г

1

р1+1 =(1 - ик+1К+2 + Vк+2

Р*+1 = (1 - V¿+2 )р+2 + V1+2

1 + в

'рк+т +

1-в

~ рк—1+3т

1+ в 2

-Р1 +п +

1-в

2

■рг—1+зп

В (2.2) параметры г, т, п определяются следующим образом:

(2.2)

0 при 1 > 0,

1 при й+1 < 0,

т

0 при йк+1 > 0,

1 при Ук+1 < 0,

п

0 при \¥г+1 > 0,

1 при Шг+1 < 0.

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

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

4+2

при и+1 > 0

шах[М| - 1, 0]С шах[Л^?2+1 - 1, 0]С при и+1 < 0,

где М = С\[Л1 /и, с2 = р7/р. Параметр С подбирается эмпирическим путем (С ~ 1). Значения Vk+1, 1 определяются аналогичным образом.

2

2

г

3. Базовый релаксационный метод

Запишем конечно-разностную задачу в виде

Ьнфн = А, (3.1)

где Ь — нелинейный разностный оператор. Для ее решения воспользуемся неявным итерационным процессом

[а1 - N] Дфь. = ЬфП - А,

ФП+1 = ФП + шДф^, (3.2)

где а и ш — релаксационные параметры, а N — линейный оператор, аппроксимирующий Ь^. В настоящей работе Ь получался из N путем "замораживания" значений плотности. Отметим, что при а = 0 мы получим упрощенный вариант метода Ньютона.

Линейная система уравнений решалась методом верхней релаксации (если течение дозвуковое) либо методом типа Гаусса — Зейделя (при сверхкритических режимах обтекания).

Замечание 3.1. Использование приближенной факторизации неявного оператора (АБ1, попеременно-треугольного методов и т. д.) упрощает реализацию базового итерационного процесса. Так, расщепляя стабилизирующий оператор по пространственным направлениям, мы получим известную схему Стегера [6] для уравнения (1.2), реализуемую скалярными прогонками. Применение АБТ метода в многосеточном алгоритме изложено, например,

в [14]. Предобусловленные методы, основанные на приближенной ЬИ факторизации неявного оператора, рассмотрены в [8,15]. Однако в этих случаях, из-за ошибок приближенной факторизации, неявный оператор не будет приближаться к Ь^ по мере уменьшения а. Возникает проблема выбора оптимальных параметров а и ш , которая сама по себе является далеко не очевидной.

4. Алгоритм ускорения итераций

Для ускорения сходимости внутренних итераций в методе (3.2) использовался алгоритм А. Слепцова [13]. Суть алгоритма заключается в построении приближения для вектора погрешности в подпространстве векторов невязок. Изначально метод ускорения сходимости был предложен для линейных итераций. В дальнейшем он был обобщен на итерационные процессы численного решения широкого класса задач аэрогидродинамики [10].

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

Ах = Ь (4.1)

итерационным методом вида

хк+1 = Тхк + g, (4.2)

где Т — сжимающий оператор.

Шаг 1. Начиная с некоторого х0 вычисляем (т + 2) приближения х (я = 1, ... , те + 2) методом (4.2).

Шаг 2. Вычисляем приращения ъп = хп+1 — хп (п =1, ... , т +1).

Шаг 3. Строим ут+2 — приближение погрешности ут+2 = хт+2 — х (х — точное решение (4.1)):

т

ут+2 = ^ а^+Ъ i=1

Коэффициенты а,1 находим из системы уравнений

т

— Zi+l, ъ,- — Zj+l)аi = (ът+1, ъ,- — ъ,-+1), ] = 1, ... , т.

i=1

Шаг 4- Поправляем хт+2:

хт+2 + Ут+2 ^ хт+2.

Если критерии сходимости не выполнены, берем хт+2 в качестве начального приближения в шаге 1.

В настоящей работе все расчеты проводились при т = 2.

Замечание 4.1. Алгоритм ускорения легко применяется к уже запрограммированным итерационным методам.

Замечание 4.2. Алгоритм проще в реализации, а во многих случаях более эффективен, чем широко используемый обобщенный метод минимальных невязок [12] (см. [10]).

Замечание 4.3. Метод ускорения может быть применен к нелинейным релаксационным методам (см., например, [10,11]), поэтому возможно использование описанного выше алгоритма для ускорения сходимости базовой релаксационной схемы (3.2).

5. Многосеточный метод на двух сетках

Известно, что релаксационные методы довольно быстро подавляют высокочастотные составляющие ошибки приближенного решения. Однако решение в низкочастотной части спектра сходится медленно. Очевидно, что низкочастотная компонента на мелкой сетке становится высокочастотной на грубой. Таким образом, после малого числа итераций на мелкой сетке ошибка будет иметь только низкочастотные компоненты, которые будут эффективно подавляться на более грубой сетке. Этот принцип лежит в основе многосеточного метода Р. Федоренко [7].

Для двух сеток алгоритм выглядит следующим образом.

Пусть требуется найти решение системы конечно-разностных уравнений (3.1). Шаг 1. Некоторым итерационным методом (в нашем случае — методом (3.2)) находится фн — приближение для решения системы (5.1).

Замечание 5.1. Дальнейшая стратегия заключается в нахождении приближения для погрешности сн из проекции условия Ьн(фн + сн) = ¡н на иерархически вложенную сетку. Шаг 2. На сетке с удвоенными шагами решается (как можно точнее) система уравнений

¿2Н02Н = Я {-Ьнфн + ¡н) + ¿2н (Яфн) , (5.1)

где Я : дн^д2н, где дн — сеточная функция. В результате получаем ф2н — приближение для ф2Н.

Шаг 3. Вычисляем сН — приближение для сн:

Сн = Р (ф2н - Яфн) ,

где Р : д2н^дн .

Шаг 4- Поправляем фн по формуле

фн + Сн ^ фн.

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

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

Замечание 5.3. Формально алгоритм легко обобщается на сколь угодно большое число вложенных сеток.

6. Многосеточные технологии и их приложения 6.1. Тестовые задачи

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

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

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

Рис. 1. Пример расчетной области и эллипсоидальной сетки.

6.2. V циклы

Как известно, многосеточная стратегия определяется, в основном, двумя факторами: последовательностью вложенных сеток и техникой циклического пересчета на ней искомых величин. Число вложенных сеток зависит от того, насколько мы в состоянии аккуратно решить вспомогательную задачу о построении приближения для вектора погрешности. Как уже упоминалось во введении, при численном моделировании пространственных стационарных трансзвуковых течений мы сталкиваемся с проблемой решения плохо обусловленной системы сеточных уравнений большой размерности. Однако, как легко видеть, уже для двух сеток вспомогательная задача сводится к решению системы в 8 раз меньшей размерности. Если при решении этой задачи возникают сложности, то вводим дополнительную вложенную сетку и рекуррентно повторяем процедуру, описанную в разделе 5. В данной работе, для случая двух и трех вложенных сеток, использовались схемы циклического пересчета, изображенные на рис. 2.

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

Замечание 6.2. Как показали численные эксперименты, при расчетах дозвуковых течений достаточно делать лишь одну внутреннюю итерацию, а алгоритм ускорения применять непосредственно к базовому релаксационному методу (см. замечание 4.3). При сверхкритическом режиме использование более 16 (ускоренных) внутренних итераций не приводит к существенному повышению скорости сходимости.

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

Рис. 2. Число внутренних итераций в V циклах.

Рис. 3. Распределение вычислительных затрат в V цикле на двух сетках.

Рис. 4. Распределение вычислительных затрат в V цикле на трех сетках.

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

Замечание 6.3. Хотя это так, использование большого числа сеток требует повышенного внимания к точности вычислений, так как многократно повторяющаяся процедура интерполирования с грубой сетки на более мелкую приводит к накоплению ошибок округления. Как показали предварительные исследования, на большом числе сеток целесообразно использовать немонотонные циклы. Однако в данной работе этот вопрос не рассматривается.

Замечание 6.4. Предложенная структура V циклов является далеко не единственной. Некоторые распространенные циклы описаны в монографии [3].

6.3. Некоторые результаты расчетов

Вначале рассмотрим результаты численного моделирования дозвуковых течений. Для их расчета использовалась последовательность вложенных сеток, содержащих 38 х 38 х 38, 19 х 19 х 19 и 10 х 10 х 10 узлов. В качестве начального приближения задавался равномерный набегающий поток. Релаксационные параметры полагались а = 0, а и = 1. Критерий сходимости — С-норма невязки меньше 10-6. Это условие гарантировало получение численного решения с точностью до 1 % (рис. 5) при расчете гипозвукового течения (Мж = 0.01). Оценка погрешности осуществлялась по известному точному решению (см.,

напр., [1]).

Об эффективности предложенной методики можно судить по данным, представленным на рис. 6. Ориентиром для сравнения служит время получения решения на самой мелкой сетке методом верхней релаксации (203 итерации). Для получения того же решения многосеточным методом понадобилось выполнение 16 (6) V циклов на двух сетках и 10 (4) — на трех. В скобках указано число V циклов при использовании алгоритма ускорения внутренних итераций. Аналогичное соотношение временных затрат было получено при расчете дозвуковых течений (вплоть до сверхкритических) при различных углах атаки.

Рис. 5. Погрешность приближенных решений.

Рис. 6. Относительные временные затраты на расчет дозвукового течения.

Рис. 7. Расчет трансзвукового обтекания эллипсоида при = 0.85. Распределение давления.

Рис. 8. Относительные временные затраты на расчет сверхкритического течения.

При моделировании трансзвуковых течений преимущества многосеточной технологии становятся более очевидными. На рис. 7 представлен расчет сверхкритического обтекания эллипсоида с использованием тех же сеток, что и в предыдущей задаче, а на рис. 8 — временные затраты на его получение. Релаксационный параметр а устремлялся к нулю по мере сходимости итераций. Для получения численного решения на одной сетке требовалось более 300 итераций по методу (3.2) (без ускорения сходимости внутренних итераций), на двух сетках — 18(10) V циклов, а на трех — 11(6) V циклов.

6.4. Локальный многосеточный метод

Численное моделирование течений сложной структуры требует привлечения адаптивных сеток. В настоящее время все более широкое распространение получают алгоритмы блочной адаптации, когда решение строится на последовательности локально вложенных, почти равномерных сеток. Технология блочной адаптации достаточно хорошо разработана, на ее основе построены эффективные локальные многосеточные методы (см., например, [4,5,9]).

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

Шаг 1. В расчетной области строим равномерную эллипсоидальную сетку 02,н (20 х 20 х 20).

Шаг 2. На сетке 02,н получаем приближенное решение.

Шаг 3. Выделяем (априори) подобласть 01сО0, в которой локализована особенность задачи (сверхзвуковая область и скачок уплотнения).

Шаг 4- В 01 строим сетку 0\ (30 х 19 х 38), шаги которой вдвое меньше чем у

Шаг 5. Интерполируем данные с О^ на С\.

Шаг 6. Выполняем V цикл на двух сетках: ^О(0н^ОН.

Замечание 6.5. Для сетки 02,н уравнение (5.1) записывается в виде

^2Нф2Н = ЯНБ2Ъ,

где ЯН Б2н = Я (—^ фн + /н) + ^н (Яфн) — если узел принадлежит Лои^, и ЯН52н = /2Н — в противном случае.

Шаг 7. Выделяем подобласть 02СВ1, в которой строим сетку ОН/2 (52 х 19 х 74).

Рис. 9. Топология вложенных сеток.

Шаг 8. Интерполируем данные с G]l на G2h/2.

Шаг 9. Выполняем до сходимости V циклы на трех сетках: G"h/2^Gjl^G0h^G\^G\/2.

Полученное решение сопоставимо по точности с расчетом на сетке, содержащей 80 х 80 х 80 узлов, при сокращении затрат памяти ЭВМ более чем в 5 раз. Временные затраты на решение задачи соответствуют расчету на последовательности вложенных сеток, содержащих 57 х 57 х 57, 29 х 29 х 29, 15 х 15 х 15 узлов.

Список литературы

[1] Лойцянский Д. Г. Механика жидкости и газа. Наука, М., 1987.

[2] Норри Д., де Фриз Ж. Введение в метод конечных элементов. Мир, М., 1981.

[3] Флетчер К. Вычислительные методы в динамике жидкости. Мир, М., 1987.

[4] Becker K., Rill S. MELINA — A multi-block, multi-grid 3D Euler code with local subblock technique for local mesh refinement. In "ICAS Conf.", Beijing, Sept., 1992, Paper No. 92-4.3.R.

[5] Berger M. J., Colella P. Local adaptive mesh refinement for shock hydrodynamics. J. Comput. Phys., 82, 1989, 67-84.

[6] Bridgeman J.O., Steger J.L., Caradonna F.X. A conservative finite difference algorithm for the unsteady transonic potential eqution in generalized coordinates. AIAA paper No. 82-1388, 1982.

[7] FEDORENKO R. P. The speed of convergence of one iterative process. Comp. Math. and Math. Phys., 4, 1964, 227-235.

[8] JAMESON A., YOON S. LU implicit schemes with multiple grids for the Euler equations. In "24th Aerospace Sciences Meeting", 1986, AIAA paper No. 86-0105.

[9] JOUHAUD J.-C., BORREL M. A Hierarchical adaptive mesh refinement method: application to 2D flows. In "Computational Fluid Dynamics'96", John Wiley & Sons Ltd., 1996, 268-274.

[10] KARAMYSHEV V. B., KOVENYA V. M., Sleptsov A. G., Cherny S. G. Variational method of accelerating linear iterations and its applications. Computers & Fluids, 25, No. 5, 1996, 467-484.

[11] KARAMYSHEV V., KUKARTSEVA O. On the optimization of the relaxation method for solving the three-dimensional problems of transonic aerodynamics. Russ. J. Numer. Anal. Math. Modelling, 12, No. 2, 1997, 151-162.

[12] SAAD Y. Krylov subspace method for solving large unsymmetric linear system. Math. Comput., 37, No. 155, 1981, 105-217.

[13] SLEPTSOV A. G. On the convergence acceleration of linear iterations. Russ. J. Theoretical Appl. Mech. Elsevier Science Publishers, 1, No. 1, 1991, 74-83.

[14] SOUTH J.C., Brandt A. Applications of multi-level grid method to Transonic Flow Calculations. Proc. of Workshop on Transonic Flow Problems in Turbomachinery, Monterey, 1976. Hemisphere, 1977, 180-206.

[15] Street C. L., Zang T. A., HUSSIAINI M. Y. Spectral multigrid method with applications to transonic potentional flow. J. Comput. Phys., 57, No. 1, 1985, 43-76.

Поступила в редакцию 18 июня 1999 г., в переработанном виде 5 июля 1999 г.