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

Том ХЫ

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

№ 3

УДК 519.63

МЕТОДЫ РЕШЕНИЯ СЕТОЧНЫХ УРАВНЕНИЙ КОНЕЧНО-ЭЛЕМЕНТНОЙ АППРОКСИМАЦИИ ПРОСТРАНСТВЕННЫХ ТЕЧЕНИЙ

А. В. ВОЛКОВ

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

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

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

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

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

Основу реализованного подхода составил так называемый полиномиальный многосеточный метод. Традиционный многосеточный метод был предложен в 1961 г. в работе [3] и в настоящее

время получил широкое распространение для ускорения численного решения систем уравнений, получающихся при аппроксимации МКО. Этот метод показал свою эффективность при решении задач с использованием мелких сеток. В методе же конечного элемента, и в частности в РМГ, высокая точность достигается не за счет измельчения сетки, а за счет увеличения числа линейно независимых полиномов (базисных функций), используемых в ячейке, и соответствующего увеличения числа степеней свободы задачи. В последнее время философия построения многосеточного метода была пересмотрена [4]. В целом ряде работ [5 — 8] метод был адаптирован к использованию в конечно-элементном методе аппроксимации течений. По аналогии с решением задачи на последовательности измельченных сеток в методе конечного элемента решается задача на фиксированной сетке, но с использованием различного набора базисных функций. Модифицировалось и название метода — полиномиальный, или p-многосеточный метод. Латинская буква <ф» отражает теперь полиномиальный характер метода, а традиционный многосеточный метод с недавнего времени предваряется буквой «Ь», отражающей необходимость выполнения расчетов на последовательности сеток с различным шагом Ь. В настоящей работе р-многосеточный метод впервые применен для решения пространственных уравнений Эйлера и Навье — Стокса.

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

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

1. Явный алгоритм поиска решений. Система сеточных уравнений, аппроксимирующих уравнения Эйлера или Навье — Стокса, может быть записана в виде [1]:

М ^ + Я(я) = 0. (1.1)

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

Явный метод решения основывается на оценке вектора невязки на текущем временном слое п:

(я”+1 - яп) / ч

М^-----—----+ Я(яп ) = 0. (1.2)

Простейший алгоритм решения, применяемый в современных СРБ-пакетах, как правило, базируется на методе Рунге — Кутта и выглядит следующим образом:

я(0) = Я”

Я(1) = Я(0)-ахЫ М-1Я ((0))

.................. (1.3)

Я(к) = ч(0) _ак А( м-1Я -1))

Яй+1 = я( к).

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

* СБЬ/

Величина шага по времени А! = ——— также ограничивается условием устойчивости схемы и

определяется допустимым числом Куранта СБЬ, характерным размером сеточной ячейки I и характерной скоростью распространения возмущения —. Например, пятишаговая схема с коэффициентами а! = 1/5, а2 = 1/4, аз = 1/3, а 4 = 1/2, а5 = 1 обеспечивает пятый порядок точности по

времени при относительно невысоком максимальном числе СБЬ= 1. Другая схема (а1 = 0.0695,

а2 = 0.1602, а3 = 0.2898, а4 = 0.5060, а5 = 1) обладает лишь первым порядком точности по

времени, но при этом допускает использование более высоких чисел Куранта (СБЬ= 1.5—2.0).

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

я(0) = Я”,

г-1

.О') =

:^(аг, я(к)+Р,к А! М-1Я ((к))), , = 1,..., т, (1.4)

к=0

я»+1 = я( т).

С целью минимизации возникновения в решении нефизичных осцилляций было разработано семейство так называемых оптимальных ТУБ-схем Рунге — Кутта [9], т. е. схем с наибольшим допустимым числом Куранта, обеспечивающим при этом сохранение по времени полной вариации решения:

dя (х)

ТУ (я (х ))=/ -яХ2* = ТУ (я0 (х )).

Например, оптимальная ТУБ-схема второго и третьего порядка точности, впервые предложенная в [10], имеет следующие коэффициенты (см. табл. 1).

Т аблица 1

Коэффициенты схемы Рунге — Кутта

Порядок агк Рг'к максимальный СЕЬ

2 1 1 1

1/2 1/2 0 1/2

3 1 1 1

3/4 3/4 0 1/4

1/3 0 2/3 0 0 2/3

Реализация оптимальных ТУБ-схем Рунге — Кутта более высоких порядков точности требует использования сопряженных операторов невязки. Такой оператор соответствует оператору невязки с измененным знаком перед диссипативным членом. Ниже приведен пример шестишаговой ТУБ-схемы Рунге — Кутта пятого порядка точности [9] с максимально допустимым числом 7

Куранта СБЬ = -

30

я(1)= я(0) + 2 L (я(0) ) •

я(2)= 3 я (0) +1 я(1) +1 ь (я (1)),

64

.(3)-1 ч(0) - 8 ь (я(0))+8 я(1) - ^ ь (я(і))+2 я(2)+2ь (я(2)) •

ь (т)+1ч(1) - 64ь (ч(1))+1<т+1 ь с*)+2ч(3)+16 ь и

(5)= 89 537 я(0) 2 276 219

2 88°°°°

4° 32°°°°

1511

28°°

64

ь (>)

87

4°7 °23 2 88°°°°'

,(1)

4°7 °23 672 °°°

1511

12°°°'

,(2)

261

14°

ь ((2)) + -ІІ-ч(3) - 261 ь ((3)) + ±ч(4) + 8ь ((4))

2°° 14° 15 7

4_

15

,(6)--

(°)

■15ь С)-45ь (ч<11)+45 о'3’+3ь И

14 —«

15

(5)

9°

■к51 )•

(1.5)

Здесь Ь (я ) = АМ-1Я (я ) — вектор-функция невязки, а Ь (я ) = АМ-1Я (я ) — вектор-функция сопряженной невязки.

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

(

Р°= Ра

?2 Л

1 + 82 г°

г°=°.°2, в-°.°°1, В2 -(х - °.5)2 +(у - °.5)2 .

Задача имеет аналитическое решение [11]:

Р (В І) = Ра

8

2а

| J° (кВ)008 (екі)кє 4аёк

где с — скорость звука; J° (х) — функция Бес-

1п ( 2)

селя первого рода; а -

(г° )

2 '

Аналитическое решение дает возможность оценить точность расчетной схемы в любой момент времени. Разница точного и численного значений давления была получена в норме ¿2 на вложенной последовательности структурированных сеток в момент времени ^ = 0.4/с. Аппроксимация по пространству выполнена с использованием РМГ с максимальной степенью полинома базисной функции К =4, что теоретически должно обеспечивать пятый порядок точности. Расчеты выполнялись как с использованием обычной схемы (1.3), так и ТУБ-схемы Рунге — Кутта (1.5) с вычислением сопряженной невязки. Сравнительные результаты оценки порядка точности численных

Рис. 1. ¿2 норма точности решения задачи о распростра-

нении цилиндрическои волны в зависимости от густоты схем представлены в табл. 2 и графически расчетноИ сетки. Сравнение обычной (пунктирная линия) на рис. 1. Видно, что только использование и ТУБ (сплошная) Р-К схем интегрирования по времени

Сетка Число степеней свободы Обычная Р-К схема TVD Р-К схема

L2 норма ошибки Порядок точности норма ошибки Порядок точности

10x10 1 500 1.27e-1 1.43е-1

15x15 3 375 6.59e-2 1.62 7.94е-2 1.45

20x20 6 000 2.81e-2 2.97 3.5°е-2 2.84

25x25 9 375 1. 19e-2 3.85 1.48е-2 3.86

30x30 13 500 4.91e-3 4.85 5.69е-3 5.24

40x40 24 000 1.23e-3 4.83 9.4°е-4 6.26

50x50 37 500 8.51e-4 1.64 2.8°е-4 5.43

ТУБ-схемы обеспечило более высокий порядок точности, превышающий теоретический пятый на самой мелкой сетке.

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

л СБЬI

ный шаг по времени, т. е. величина А^ = ——— в разных ячейках различна и определяется ее размерами и характерной локальной скоростью распространения возмущений. Отметим, что при решении нестационарных задач шаг по времени в численной схеме один для всех ячеек и ограничивается размерами минимальной ячейки и максимальной скоростью распространения возмущений. Число Куранта в явной схеме ограничивает скорость сходимости к стационарному решению. Поэтому использование таких схем требует большого количества итераций и редко применяется для расчета вязких течений при больших числах Рейнольдса.

2. Неявный алгоритм поиска решений. Неявный алгоритм решения теоретически допускает использование бесконечно больших чисел Куранта. Один из простейших вариантов неявной схемы требует оценки невязки в уравнении (1.1) на следующем п +1 временном слое:

( - ЯП) , ч

М^------А--- + К (яп+1 ) = 0. (2.1)

После разложения невязки в ряд Тейлора

^ + К (я" ) + £ ( - я " )-°

/ пП+1 пП\

_ _чі+r (ч.)+®

поправка к вектору решения находится путем решения следующей системы линейных уравнений:

ІСБЬ 5я

или ввода соответствующих обозначений:

(q”1 -4") = -R(q") (2.2а)

Ax = b. (2.2б)

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

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

яп+1 = яп + юх; 0 <ю<1. (2.3)

Новое решение должно быть физически корректно. Например, локальные величины давления и температуры должны быть положительны. Если в какой-либо ячейке это условие не выполняется, то величина коэффициента ю уменьшается. Другое условие на коэффициент ю требует уменьшения ¿2 нормы новой невязки. Иногда для поиска оптимального коэффициента ю находят промежуточное значение ¿2 нормы невязки при половинном шаге ю = 0.5. Предполагая квадратичную зависимость ¿2 нормы невязки от параметра ю, легко найти его наилучшее значение.

Процедуру расчета начинают с малых чисел Куранта. Затем на следующем шаге по времени число Куранта определяется в соответствии с соотношением

СРЬпете = СРЬоМ2(2ю-1). (2.4)

Таким образом, предполагается, что при большом значении параметра ю >0.5 задача хорошо разрешается и число Куранта может быть увеличено. И напротив, проблемы поиска нового решения, приводящего к уменьшению ¿2 нормы новой невязки, отражены при значении параметра ю <0.5, что автоматически приводит к уменьшению числа Куранта на следующей итерации.

Решение системы линейных уравнений, как правило, осуществляется обобщенным методом минимальных невязок (ОМКЕ8 [12]). Для ускорения расчета используют предобуславлива-тель Р — легко обратимую матрицу, умножение на которую снижает жесткость решаемой линейной системы уравнений:

Р * А; АР-1 = В * I. (2.5)

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

АР-1 Рх = Ву = Ь, (2.6)

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

х = Р-1у.

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

Широкое распространение получили предобуславливатели, построенные на основе Ьи факторизации матрицы А [14]. Предобуславливатель Р = Ьи = А идеален, но его построение обходится дороже, чем непосредственное прямое решение системы (2.2). Известно, что места расположения ненулевых элементов (портрет) матриц Ь и и, повторяя внешний контур портрета матрицы А, внутри этого контура полностью заполняют матрицы Ь и и вплоть до главной диагонали. Таким образом, количество ненулевых элементов этих матриц значительно превышает эту величину в основной матрице. Как результат, полная факторизация требует огромных затрат как памяти, так и количества операций, необходимых для вычисления каждого элемента.

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

эффективный предобуславливатель. Рассматривают также заполненные варианты портретов фак-торизированных матриц 1Ьи[1], 1ЬИ[2] и т. д., портреты которых больше, чем 1Ьи[0], но существенно меньше, чем в случае полного ЬИ разложения. Существуют и более усложненные алгоритмы, когда добавление нового элемента в предобуславливатель осуществляется только при значимом его воздействии на результат [14].

Отметим, что в настоящий момент вычислитель освобожден от необходимости программировать эти алгоритмы. Можно свободно воспользоваться одним из современных пакетов программ, содержащих основной набор предобуславливателей и методов решения систем линейных уравнений. Одним из таких известных пакетов является РБТ8е [15]. Его отличительная особенность заключается в способности выполнения алгебраических операций с использованием многопроцессорной техники.

Неявный метод решения систем сеточных уравнений весьма эффективен, но он требует грандиозных затрат оперативной памяти компьютера, необходимой для хранения матрицы Якоби невязки. К значительной экономии приводит использование безматричной формы метода ОМКЕ8 [16]. Алгоритм решения требует использования не самой матрицы А, а только результата матрично-векторного произведения Ар, где р — некий вектор. Если к невязке Я добавить слагаемые, соответствующие производной решению по времени, и определить функцию расширенной невязки

Г (я )=М—¡7 я+Я (я).

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

Ар *Ь,тГ(я+ 8р)-Г(я). (2.7)

8——0 £

При этом основные затраты оперативной памяти приходятся лишь на хранение матрицы В в (2.5).

Эффективность неявного метода решения была апробирована в задаче расчета многоэлементного профиля на последовательности адаптивных сеток [1]. Фрагменты одной из сеток представлены на рис. 2. Отметим, что в данном примере гладкость изменения форм и размеров соседних ячеек неудовлетворительна с точки зрения МКО. Однако используемый в настоящей работе РМГ не теряет точности даже на таких высокоанизотропных сетках. Эта особенность РМГ в сравнении с МКО была проанализирована в работе [17]. Отметим, что решение нелинейных аппроксимирующих уравнений на таких сетках требует использования очень эффективных методов, в частности неявной схемы, которая и была успешно использована в [1].

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

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

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

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

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

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

В настоящей работе наибольший выигрыш в скорости сходимости получен при комбинировании традиционного многосеточного метода и p-многосеточного метода. В основе традиционного многосеточного метода использован агломерационный подход [18], реализованный в промышленном коде FINE™/Hexa, компании NUMECA Int [20]. Этот h-метод с агломерациями ячеек применяется только на нижних многосеточных уровнях с использованием максимально простой схемы аппроксимации (K =0 — одна кусочно-постоянная функция). Затем на самой мелкой сетке выполняется постепенное наращивание количества используемых базисных функций, т. е. используется p-многосеточный метод. Такое последовательное чередование р-уровней образует одну итерацию многосеточного метода.

3.1. Полиномиальный многосеточный метод для РМГ. В основе реализованного р-много-сеточного метода лежит известная схема полной аппроксимации (FAS) [19], применяемая к решению нелинейных задач традиционным h-многосеточным методом.

Определяется последовательность уровней многосеточного метода l = 1,...,L, где уровень l +1 в сравнении с уровнем l имеет либо большее количество ячеек, либо большее количество базисных функций (степеней свободы). При применении традиционного сеточного измельчения (h-измельчение) используется последовательность вложенных гексаэдральных сеток, генерируемая кодом HEXSTREAM (компании NUMECA Int. [20]) на основе алгоритмов агломерации. При измельчении по степеням свободы (р-измельчение) используются наборы из следующих базисных функций.

К =0 К = 1 К = 2 К =3

Базис набора К =0 дополняется функциями: Базис наборов К =0, 1 дополняется функциями: Базис наборов К =0, 1, 2 дополняется функциями:

1 х, у, г 2 2 2 х , у г , ХУ, хг, уг х3, у3, г3, 222 х у, х г, ух, 222 у г, ?х, г у, хуг

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

(х - x0 )а(у - у0 - z0 )

ф, (х) = ----- } К в ’ К,= 1,...,Kf; 0<а + р + у<К.

' иа % и] ■' f р '

Здесь х0, у0, z0 — некая внутренняя точка гексаэдра, а Их, Иу и Иг — линейные размеры

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

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

Kf

Я (, х y, z) = £ и] () Фу (x, У, z),

,=1

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

КГ (К) =

(К +1)( К + 2 )( К + 3)

6

(3.1)

Решаемая система уравнений Навье — Стокса на самом верхнем (мелком) многосеточном уровне Ь может быть записана в виде:

М^ + Яь (Яь ) = 0.

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

Чь = Яь + ЛЯь .

Запишем теперь сеточные уравнения на более грубом (расположенном ниже) многосеточном уровне I = Ь -1:

М + Я, (я, ) = ¥1.

(3.2)

Появившаяся в уравнении правая часть призвана учесть ошибку, связанную с потерей точности аппроксимации при переходе на нижний уровень, где уравнения записываются либо на более грубой сетке, либо с использованием полиномов меньшей степени. Для определения вида правой части рассмотрим эту систему уравнений на двух последовательных уровнях — уровне I (уравнение (3.2)) и верхнем уровне I +1:

(3.3)

Итерационная схема должна быть построена таким образом, чтобы решение, перенесенное

Вычитая из уравнения (3.4) уравнение (3.5), получим рекуррентное соотношение для правой части:

На самом верхнем уровне уравнения не должны содержать правой части Е,=ь = 0. Здесь

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

где начальное приближение для решения на более грубом уровне является результатом его переноса с верхнего уровня:

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

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

f=1, Nf

Здесь г — индекс элемента, а п — индекс соседнего элемента, отделенного от элемента г стороной f (суммирование производится по всем сторонам элемента расчетной сетки); —

площадь соответствующей стороны.

с более верхнего уровня на нижний Чі = 111+1 [чі+! ], также являлось решением уравнения на нижнем слое. Перепишем уравнение (3.2) в виде:

(3.4)

После перенесения уравнения (3.3) на нижний уровень имеем:

(3.5)

(3.6)

!/+1 [■], 1^+1 [■] — операторы переноса решения и невязки решения с высокого многосеточного

Ч і = Я/0)+АЯ/,

ч/0) = II+1 [ч+1 ].

(3.7)

Чкорр ={Іі-1 Чі-і Іі 1 [чі ] },

где Іі-1 [■] является оператором переноса решения с нижнего на верхний многосеточный уровень. Выражения в фигурных {■} скобках подвергаются осреднению по следующему правилу:

q корр = ^оррD / , ^ г\ | Г. (38)

max (q г |, 1г [q г J, |q t _i|, sj

В выполненных расчетах коэффициент демпфирования D варьировался в диапазоне значений 0.25-0.75, s = 1 -10_20.

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

3.2. H-измельчение. В случае перевода решения с мелкой сетки на крупную используются следующие операторы:

!/+1 [q/+1J =

EH+1)Ql+1.

Е«*1 ’

1гг+1 [^+1 ]=^{+1}.

Здесь суммирование производится по всем ячейкам уровня I +1, содержащимся внутри крупной ячейки уровня I; О — объем ячейки.

Обратный перевод решения с крупной сетки на мелкую осуществляется оператором первого порядка точности:

+1 к ] = Я.

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

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

М т

q/ = Е^фу и q/-1 = ^ujlфjl в рассматриваемой ячейке. Количество различных функций

У =1 у=1

в базисах определяется в соответствии с соотношением (3.1) М = К^ (К1), т = К^ (к1 -1). Равенство решений здесь понимается в смысле выполнения интегральных условий: К - q/-1) фу^О = 0, получаемых умножением разности решений на систему базисных функций

О

фу уровня, на который переводят решение. Таким образом, операторы интерполяции \\+1 и сборки I+1 могут быть представлены в виде:

I+1 к ]=[> рр ]-1 м ,

\\+1 [+1 ] = [Mдд ] 1 Mgpql+1. (3.9)

Оператор (3.10) для перевода невязки в соотношении (3.6) может быть отличен от оператора (3.9):

3+1 [ I+1 ] = Mgp [Mpp ]_1 Rl+l. (3.10)

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

о

о

о

о

3.4. Стратегии полиномиального многосеточного подхода. Изложенный выше алгоритм многосеточного метода использует множество параметров, настраиваемых при решении конкретных задач. В первую очередь необходимо определить оптимальное количество уровней многосеточного метода. Количество уровней р-измельчения равно К — максимальной степени полиномов базисных функций. Количество уровней вложенных сеток Ь-измельчения зависит от геометрической сложности и общего количества ячеек самой мелкой сетки. Предполагается, что расчет начинается на самом мелком (верхнем) сеточном уровне, затем решение переводится на уровень ниже, и после нескольких итераций сглаживания (путем итерационного решения уравнения (3.3)) решение опять переводится на уровень ниже, такая процедура выполняется вплоть до достижения самого грубого уровня. На следующем этапе осуществляется коррекция на верхних многосеточных уровнях с использованием формулы (3.7). Весь этот полный цикл является одной многосеточной итерацией (У-цикл). Данный алгоритм допускает варьирование как типа сглаживателя (явный или неявный метод), так и количества итераций применения сглаживателя и коэффициента демпфирования в выражении (3.8). В проведенных численных экспериментах на верхнем уровне использовалась одна итерация сглаживания. По мере спуска к нижнему уровню число итераций увеличивалось на единицу.

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

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

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

4. Обтекание сферы невязким потоком. Серия расчетов обтекания сферы дозвуковым потоком сжимаемого газа при М = 0.15 была выполнена на последовательности вложенных сеток с использованием кусочно-линейного, квадратичного и кубичного полиномиального базисов. Фрагменты сеток, сгенерированных около четверти сферы, представлены на рис. 3. Элементами неструктурированных сеток являются гексаэдры. В случае расчетов с высоким порядком точности (квадратичный и кубичный полиномиальные базисы) элементы, примыкающие к поверхности обтекаемого тела, имели одну криволинейную грань, отражающие кривизну сферы.

Т-к\ \-ш ЗЗг ■ щ С Т V— л \ \~1ГТ да# \ . * ряр \ V \

- —

Рис. 3. Фрагменты вложенной последовательности сеток около сферы

5-шаговая схема Рунге — Кутта - - - 5-уровневый р-многосеточный метод А ' Локально-неявный сглаживатель

_1__________I_________I

~ 0 10000 20000 30000 40000 50000

Время СРи, с

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

....— сглаживатель на основе метода Рунге — Кутта; — А — — локальнонеявный сглаживатель

При расчете на каждой из сеток путем агломераций образовывалась система из трех вложенных сеток, решения на которых получались с использованием наименьшего кусочно-постоянного базиса. Далее уровни многосеточного алгоритма измельчались путем обогащения базиса, как было описано ранее. Таким образом, расчеты РМГ с К =3 (кусочно-кубический базис) выполнялись с использованием шестиуровневого алгоритма.

Типичная история сходимости (зависимость ¿2 нормы невязки решения от времени СРи) представлена на рис. 4 для случая кусочно-квадратичной аппроксимации решения. Рисунок демонстрирует эффект ускорения решения при использовании многосеточного алгоритма. Здесь штриховая линия представляет результат использования сглаживателя на основе метода Рунге — Кутта, а сплошная линия с маркерами — результат использования более эффективного локальнонеявного сглаживателя. Отметим, что без применения многосеточного алгоритма при кусочнокубичном базисе получение решения в разумное время становится невозможным.

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

качество расчетной схемы. Анализ результатов показал, что скорость стремления к нулю коэффициента сопротивления по мере сгущения сеток выше, чем 2K +1, что говорит о достижении теоретически обоснованной [21] суперсходимости конечно-элементного метода.

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

сферы составил величину 3.2-10-6. Получить такой низкий коэффициент сопротивления методом конечного объема чрезвычайно затруднительно. Например, использование кода FLUENT, даже на сетке в 240 000 узлов не позволило получить сопротивление меньше, чем

2.1 -10-4.

5. Применение полиномиального многосеточного метода к расчету обтекания крыла. P-многосеточный алгоритм был применен также к расчету дозвукового турбулентного обтекания крыла LANN [22]. Одна из целей расчета состояла в сравнении конечноэлементного РМГ на основе кусочно-линейных базисных функций с методом конечного объема второго порядка точности при условии расчета на сетках, обеспечивающих использование эквивалентного количества переменных.

Расчеты были выполнены на неструктурированных гексаэдральных сетках. Сетка с числом элементов 190 213 ячеек была использована для РМГ, в то время как расчеты МКО выполнялись на сетке с числом элементов 625 076. В обеих сетках первый слой точек располагался над обтекаемой поверхностью на расстоянии, обеспечивающем величину у+ «1. Некоторые фрагменты крупной сетки представлены на рис. 6.

Использовались следующие условия для набегающего потока: M = 0.4, а = 0.6°,

Re = 7.17 • 106. Все расчеты были выполнены с моделью турбулентности Спаларта — Аллмараса. Начальный коэффициент турбулентной вязкости на бесконечности был равен 5.

Решение в МКО [20] было получено с использованием трехуровневого многосеточного метода. В РМГ использовался четырехуровневый p-многосеточный метод. При этом решение на первых трех уровнях (от грубой до мелкой сеток) получалось на базе РМГ с кусочно-постоян-

Рис. 6. Фрагменты неструктурированной гексаэдральной сетки около крыла LANN

^---- РМГ(К= 1)

Число степеней свободы

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

Рис. 7. Фрагменты агломерированных сеток около крыла LANN

200 400 600 800

Итерации многосеточного метода

Рис. 8. История сходимости для МКО и РМГ

Рис. 9. Сравнение распределений давления в разных сечениях крыла LANN при расчете МКО на сетке 625 076 ячеек и

РМГ на сетке 190 213 ячеек

ными функциями. Затем на самой мелкой сетке решение на базе РМГ получалось с использованием кусочно-линейных базисных функций — четвертый многосеточный уровень. Фрагменты агломерированных сеток, использованных на начальных уровнях многосеточного подхода, представлены на рис. 7. Типичная история сходимости L2 нормы невязки по итерациям для МКО и РМГ представлена на рис. 8. Использование более грубой сетки в РМГ обеспечило меньшее количество многосеточных итераций и сравнимое расчетное время (~26 ч на персональном компьютере 2 ГГц). Это, в частности, свидетельствует об эффективности использованного многосеточного метода в РМГ.

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

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

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

Работа частично выполнена за счет средств гранта РФФИ № 09-01-00243а.

ЛИТЕРАТУРА

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

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

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

4. Rönquist E. V., Patera A. T. Spectral element multigrid. I. Formulation and numerical results // J. of Scientific Computing. 1987. V. 2. N 4.

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

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

7. Luo H., Baum J. D., Löhner R. A fast, p-multigrid discontinuous Galerkin method for compressible flows at all speeds // AIAA Paper 2006-110, 2006.

8. F i d k o w s k i K. J. O l i v e r T. A., L u J., D ar m o f a l D. L. p-Multigrid solution of high-order discontinuous Galerkin discretizations of the compressible Navier — Stokes equations // J. of Comp. Phys. 2005. V. 207, p. 92—113.

9. Shu C. W., Osher S. Efficient implementation of essentially non-oscillatory shock-capturing schemes// J. of Comp. Phys. 1988. V. 77, p. 439—471.

10. Co с kb urn B. Discontinuous Galerkin methods for convection-dominated problems, in High-Order Methods for Computational Physics, edited by T.Barth and H.Deconik, Lecture Notes in Computational Science and Engineering, V. 9 — Berlin: Springer Verlag, 1999.

11. Tam C. K. W., Webb J. C. Dispersion-relation-preserving finite difference schemes for computational acoustics // J. of Comp. Phys. 1993. V. 107, p. 262—281.

12. Saad Y. and Schultz M. GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems // SIAM J. Sci. Statist. Comput. 1986. V. 7, p. 856—869.

13. Saad Y. Iterative methods for sparse linear systems. Books — 1-st and 2-nd SIAM edition. Имеется в открытом доступе: http://www-users.cs.umn.edu/~saad/books.html

14. Chapman A., Saad Y., Wigton L. High-order ILU preconditioners for CFD problems // International Journal for numerical methods in fluids. 2000. V. 33, p. 767—788.

15. http://www.mcs.anl.gov/petsc/petsc-as/

16. Wigton L., Young D. GMRES acceleration of computational fluid dynamics codes // AIAA Paper 85-1494, 1985.

17. Petrovskaya N. B., Wo lko v A. V. The issues of the solution approximation in higher-order schemes on distorted grids // International Journal of Computational Methods. 2007. V. 4, N 2.

18. Mavriplis D. J. and Venkatakrishnan V. Agglomeration multigrid for twodimensional viscous flows // Computers and Fluids. 1995. V. 24, N 5, p. 553—570.

19. Brandt A. Multilevel adaptive computations in fluid dynamics // AIAA J. 1980. V. 18,

N 10.

20. www.numeca.be

21. B a r t h T. J. A posteriori estimation and mesh adaptivity for finite volume and finite element method // Published in Springer series Lecture Notes in Computational Science and Engineering (LNCSE). 2004. V. 41.

22. Müller U. R., Schulze B., Henke H. Computation of transonic steady and unsteady flow about LANN wing. Validation of CFD codes and assessment of turbulence models // ECARP Report 58. Vieweg, 1996, p. 479—500.

Рукопись поступила 7/VII2009 г.