CC BY
248
Том XXXVIII
УЧЕНЫЕ ЗАПИСКИ ЦАГИ 200 7
№ 3 — 4
УДК 533.6.011.35:629.7.025.73 533.6.011.32:629.7.025.3 519.63:517.96
ПРИМЕНЕНИЕ КОНЕЧНО-ЭЛЕМЕНТНОГО МЕТОДА ГАЛЕРКИНА С РАЗРЫВНЫМИ БАЗИСНЫМИ ФУНКЦИЯМИ К РЕШЕНИЮ УРАВНЕНИЙ РЕЙНОЛЬДСА НА НЕСТРУКТУРИРОВАННЫХ
АДАПТИВНЫХ СЕТКАХ
А. В. ВОЛКОВ, С. В. ЛЯПУНОВ
Конечно-элементный метод Галеркина с разрывными базисными функциями (DG — Discontinuous Galerkin) применен к решению уравнений Рейнольдса (уравнения Навье —
Стокса и модель турбулентности Спаларта — Аллмараса) с использованием адаптивных неструктурированных сеток. Приводится описание численного метода. Рассмотрены примеры расчета вязкого трансзвукового обтекания изолированного профиля и дозвукового обтекания трехэлементной конфигурации. Решения этих практически важных задач получены с использованием технологии адаптивных сеток.
Несмотря на большой прогресс в разработке численных методов решения уравнений Рейнольдса, расчет обтекания тел сложной геометрии все еще сталкивается со значительными трудностями. Частным случаем таких течений является обтекание взлетно-посадочной механизации, которое характеризуется наличием отрывов потока, взаимодействием пограничных слоев, сходящих с различных элементов механизации, возможным наличием местных сверхзвуковых зон и пр. Сложные геометрические формы обуславливают применение неструктурированных сеток, алгоритмы построения которых позволяют сформировать сетку с необходимыми
свойствами вокруг обтекаемых тел в соответствии с особенностями течения. Создание
эффективного метода расчета требует оптимизации распределения узлов, т. е. адаптации
расчетной сетки.
Одним из наиболее распространенных методов решения уравнений Навье — Стокса является метод конечного объема, который используется в целом ряде коммерческих программ. Неизвестными при этом являются средние значения параметров потока в контрольных объемах, и применяется полиномиальная аппроксимация численного решения в этих объемах с использованием некоторого набора ячеек (шаблона). В частности, второй порядок точности на гладких решениях обеспечивается линейной аппроксимацией решения внутри контрольного объема. Пространственная аппроксимация является результатом интегрирования потоков вдоль границ контрольных объемов. Расчет потоков на границе между элементами, где решение, вообще говоря, терпит разрыв, осуществляется путем точного (метод Годунова) или приближенного (методы Роу, Ошера и пр.) решения задачи Римана о распаде произвольного разрыва. Метод конечного объема хорошо работает на сетках, близких к равномерным, и в этом случае обеспечивает надежную сходимость даже с использованием простых явных схем интегрирования по времени при небольших затратах дополнительной памяти. Для ускорения процесса получения решения
на густых равномерных сетках используют многосеточную технологию.
Однако использование равномерных сеток в задачах со сложной структурой течения ограничено из-за необходимости построения слишком больших расчетных сеток и, как следствие, слишком значительных затрат компьютерных ресурсов. Например, расчетная сетка
вокруг взлетно-посадочной механизации должна сгущаться не только вблизи твердой границы для адекватного расчета пограничного слоя, но и в зонах предполагаемого прохождения вязких следов. Таким образом, ввиду невозможности точного предсказания поведения следа за элементами механизации равномерно густые сетки вблизи сложных геометрических форм взлетно-посадочных систем получаются неоправданно громоздкими. Уменьшить количество узлов сетки при сохранении точности решения можно путем адаптации ее к решению. Становится очевидным, что для создания эффективного метода требуется единое рассмотрение процессов генерации, адаптации сетки и решения уравнений на ней. При этом возможно использовать не только механизм измельчения (укрупнения) сетки, но и механизм повышения (понижения) порядка точности аппроксимации.
Расчетные сетки, адаптированные к особенностям течения при больших числах Рейнольдса, являются, как правило, неравномерными и содержат сильно вытянутые ячейки. В таких сетках крупные ячейки правильной формы могут соседствовать с мелкими сильно вытянутыми ячейками, что может, в частности, приводить к снижению точности определения градиентов параметров потока. Это может приводить к отсутствию сходимости не только явных итерационных методов, но и неявных ньютоновских алгоритмов. Особенно затруднительным на таких сетках является использование методов высокого порядка точности, кроме того, повышение порядка конечно-объемного метода приводит к расширению шаблона аппроксимации.
Альтернативу конечно-объемным методам составляют методы конечного элемента. Метод, основанный на конечно-элементном подходе, может использовать высокий порядок аппроксимации в ячейках при ограниченном шаблоне, условие равномерности сетки имеет меньшее влияние на точность. В каждом элементе расчетной области решение представляется в виде линейной комбинации непрерывных базисных функций. Система дискретных уравнений формируется путем умножения дифференциальных уравнений на пробные функции с последующим интегрированием по частям по области расчета (слабая формулировка). В схеме Галеркина пространства пробных и базисных функций совпадают. Метод конечного элемента хорошо обоснован для решения эллиптических задач. Применение его к решению гиперболических задач, в частности
в форме законов сохранения (уравнения Эйлера, Навье — Стокса), требует его модификации [1]. Один из способов такой модификации реализован в методе SUPG (Streamwise Upwind Petrov Galerkin) путем добавления слагаемых, обеспечивающих устойчивость метода [2].
Другой конечно-элементный подход к численному решению уравнений в форме законов сохранения предложен сравнительно недавно [3, 4] и получил название в литературе DG (Discontinuous Galerkin). В отличие от SUPG базисные функции здесь терпят разрыв на границах элементов. С целью обеспечения устойчивости метода потоки на границах определяют, как и в методе конечного объема, с использованием приближенных решений задачи Римана о распаде разрыва.
В качестве базисных функций обычно выбирают полиномы. Порядок точности решения возрастает с увеличением порядка используемых базисных полиномов. В отличие от метода конечного объема повышение порядка точности не требует расширения шаблона.
Использование схемы с высоким порядком точности для решения уравнений Навье — Стокса может существенно сократить количество узлов расчетной сетки, особенно в областях гладкости решения. Вышеописанная схема позволяет локально изменять порядок точности, что дает возможность осуществлять адаптацию к решению не только путем изменения размеров ячеек сетки (й-refinement), но и путем изменения порядка полиномов, аппроксимирующих решение
(p-refinement). Например, в областях, содержащих ударные волны и пограничные слои (с толщиной порядка размера ячейки), можно использовать полиномы низкого порядка и сгущать сетку, и, напротив, в областях гладкости решение можно получать с высоким порядком полиномов.
В данной работе дается описание метода Галеркина с разрывными базисными функциями применительно к решению уравнений Навье — Стокса и Рейнольдса. Иллюстрируется возможность получения решения уравнений движения вязкого газа с высоким порядком
точности. С использованием технологии адаптивных сеток выполнены решения задач, имеющих практически важное значение.
1. Реализация метода Галеркина с разрывными базисными функциями при решении плоских уравнений Навье — Стокса с моделью турбулентности Спаларта — Аллмараса.
В консервативной форме система уравнений Навье — Стокса с моделью турбулентности Спаларта — Аллмараса [5] записывается в следующем виде:
д(м + V- Е1 (м/) = V- Еу (м, Ум) + S (м, Ум). (1)
Здесь вектор консервативных переменных w = (р, ри, ру, ре, ру)Т состоит из значений
плотности р, компонент скорости и, V, внутренней энергии е и величины дополнительной турбулентной вязкости V. Значения функций Е1 и Еу для невязких и вязких потоков, а также различные замыкающие соотношения можно найти, например, в работе [6]. Источниковый член S присутствует в последнем уравнении, моделирующем турбулентную вязкость [5].
Пусть вся расчетная область разбита на неперекрывающиеся элементы, вообще говоря, произвольной формы. Численное решение для каждой переменной рассматриваемой системы (1) в каждом элементе представляется как линейная комбинация базисных функций: ф, (х, у):
к (х, у, г) = 2 и, (X)ф (х, у), г =1, ..., КсГ = (К + 02К + 2).
г=1 2
КсГ-1
Система базисных функций ф, (х, у) = ф, (р (х, у), ц (х, у)) = 2 АпР“ ЦР, а + р = п
п=0
определена для стандартного треугольника, вершины которого расположены в точках р (0,0), (1, 0), Р, (0,1) параметрического пространствар, ц. Любой треугольный элемент с вершинами Р1 (х1, у1), Р2 (х2, у2 ), Р3 (х3, у3 ) в физическом пространстве х, у преобразуется в стандартный треугольник через линейное преобразование:
| х = (х1 - х3 ) р + (х2 - х3 ) ц + х3,
[у = (у1 -у3)Р + (у2 -у3)ц + у3.
Коэффициенты Ап подбираются из условия ортогональности функций. Ксг — количество базисных функций, которое меняется от ячейки к ячейке и для различных уравнений и определяется максимальной степенью базисных полиномов К. Коэффициенты разложения и, (х)
являются искомыми неизвестными.
В соответствии со стандартным методом конечного элемента Галеркина система сеточных уравнений относительно и, (х) получается из условия ортогональности уравнения (1) к каждой
из пробных функций, выбранных из того же функционального пространства, что и базисные функции. Это условие ортогональности записывается в так называемой «слабой форме», что соответствует умножению (1) на каждую из пробных функций и интегрированию по частям по каждому элементу О:
4-1 ф/МкйО + 0ф/ (Е - Еу ) - Пйъ - | Уф,- (Е - Еу )йО - | ф, БйО = 0. (2)
О дО О О
Здесь п — вектор нормали к границе. В отличие от стандартного метода Галеркина в
настоящем методе базисные функции на границах элементов терпят разрыв. Это, с одной
стороны, обуславливает компактность шаблона аппроксимации (в случае решения уравнений Эйлера шаблон включает в себя только текущий и соседние с ним элементы), а с другой стороны, приводит
к необходимости расчета интегралов невязких ^ (и к ) и вязких ¥у (и к ) потоков по границе 50 каждого элемента.
Значение невязких потоков ^ (и к) на границах между контрольными объемами
(интерфейсами) определяется на основе идеи, используемой в конечно-объемных ориентированных схемах, в которых с целью стабилизации решения уравнения переноса используется точное (метод Годунова) или приближенное (методы Роу [7], Ошера [8] и др.) решение задачи о распаде произвольного разрыва (задача Римана):
Здесь индексы Ь и Я соответствуют значениям потока справа и слева от рассматриваемого ребра; Б — диссипативное слагаемое, определяемое в результате решения задачи Римана.
Расчет вязких потоков (wh, ) требует определения не только зависимых переменных,
но и их градиентов. В данной работе в соответствии с [9] для этого вводятся дополнительные независимые переменные g и к, понижающие порядок дифференциального уравнения:
Здесь индексы х и у обозначают производные величин в соответствующих направлениях. При этом градиенты переменных в каждом контрольном элементе также представляются в виде линейной комбинации базисных функций фг- (х, у):
тестовые функции фу и интегрирования по частям в каждом из контрольных объемов имеем:
Необходимые в (2) значения вязких потоков и необходимые в (5) значения зависимых переменных на границах элементов вычисляются следующим образом:
Такое определение потоков и зависимых переменных обеспечивает компактность шаблона, отсутствие четно/нечетного расщепления решения и максимально возможный порядок точности
На ребрах, образующих обтекаемое тело, компоненты скорости считаются нулевыми (условие прилипания), в то время как оставшиеся переменные рассчитываются по текущим значениям коэффициентов разложения в приграничных ячейках. Используется также условие теплоизолированной стенки, т.е. поток энергии на обтекаемой границе считается нулевым.
Граничные условия на внешней границе формируются на основе одномерного анализа характеристик уравнений Эйлера. При этом используются инварианты Римана, тангенциальная составляющая скорости и энтропия, которые определяются либо из внутренних ячеек, либо считаются заданными в невозмущенном потоке в зависимости от втекания или вытекания сверхзвукового или дозвукового потока. Величина турбулентной вязкости на обтекаемом теле
(3)
8 = {рх, (р«)х, (ру)х, (ре)х , (р*)х}, и = {ру, (р«)у , (ру)у, (ре)у, (р*)у }. (4)
Ко/
Ко/
8и (х y, *)=2 (* )ф (х у ) (х y, *)=2 н (* )ф (х у )
і=і
і=і
где Ог- (X), И1 (X) — искомые коэффициенты разложения. После умножения соотношений (4) на
] П j ту п X
О дО О
О - I ф^к • пхёс + whdО = 0,
О
О
(5)
(6)
О ( ИК+1) [9]
принимается равной нулю, в то время как на бесконечности эта величина пропорциональна коэффициенту молекулярной вязкости vOT (в выполненных расчетах коэффициент
пропорциональности был ра-
вен 5).
Уравнения (2), (5), дополненные условиями определения потоков на границах контрольных объемов (3), (6), образуют систему обыкновенных нелинейных дифференциальных уравнений относительно коэффициентов разложения Uj (t).
Интегралы системы уравнений (2), (5) оцениваются на основе квадратур Гаусса. В случае одномерного интегрирования количество квадратурных точек равно K + 1. Кратные интегралы внутри контрольных объемов определяются либо с использованием суперпозиции одномерных Гауссовых узлов, либо на основе специальных кубатурных формул.
В настоящей работе было рассмотрено два подхода к приближенному решению задачи Римана, а именно схемы Роу [7] и Ошера [8]. Метод Роу основан на точном решении задачи Римана для специальным образом линеаризованной системы уравнений. Решение состоит из движущихся разрывов, отделенных друг от друга областями с постоянными значениями параметров. Метод обеспечивает точное выполнение соотношений Гюгонио на разрывах. В качестве недостатка схемы Роу отмечают отсутствие дифференцируемости диссипативного члена в некоторых точках течения вследствие наличия в соотношениях модуля от собственных значений матрицы потоков. В схеме Ошера решение задачи Римана аппроксимируется комбинацией простых волн Римана. Эта схема строится на основе полностью дифференцируемых функций, однако, в сравнении со схемой Роу требует примерно на 30% больше численных операций. Исследования, проведенные авторами настоящей работы, показали полную идентичность расчетов, выполненных по схеме Галеркина с разрывными базисными функциями с использованием обоих методов приближенного решения задачи Римана. Наличие недифференцируемого диссипативного члена схемы Роу
не оказало существенного влияния на процессы сходимости решения. Результаты расчетов, приведенные ниже, получены с использованием подхода Роу.
Полученная в результате система обыкновенных нелинейных дифференциальных уравнений решалась с помощью обратного метода Эйлера с линеаризацией на верхнем временном слое и с локальным шагом по времени в следующей форме:
“ J (W"h )]AWh = RES (W"h ). (7)
Здесь RES обозначает вектор невязки (уравнение (1)); J(w) = dRES — матрица Якоби;
dw
AWh = wn+1 - wn; n — порядковый номер текущей итерации; C — симметричная блочно-диагональная матрица; коэффициенты блока для элемента Q равны Cj = J фгфу-^О; At — шаг по
О
д CFL/in
времени, который определялся следующим образом: At =-----------, где lin — радиус вписанной
Vmax
окружности, Vmax = max (|Vn - а|, |Vn|, |Vn + a|) — максимальная характеристическая скорость;
CFL — число Куранта, которое выбиралось постоянным для всех элементов. Это означает, что шаг по времени - локальный и различный для разных элементов. Бесконечное число Куранта соответствует чистому методу Ньютона для стационарных уравнений.
Система линейных уравнений (7) решалась безматричным (matrix-free) методом GMRES с предобусловливателем в виде неполного LU разложения. Матрица Якоби формировалась на основе аналитического дифференцирования вязких и невязких потоков с численным дифференцированием диссипативного слагаемого Роу.
Программа решения системы нелинейных уравнений сконструирована из подпрограмм библиотеки PETSc версии 2.1.5, бесплатно распространяемой через интернет [10]. Эта
библиотека дает возможность оперировать с большими массивами данных, расположенных на различных вычислительных машинах, и производить решение больших систем уравнений с использованием множества процессоров. Таким образом, разработанный вычислительный код использует технику параллельных вычислений.
2. Использование адаптивных неструктурированных сеток для расчета вязких турбулентных течений. На адаптивных сетках может быть достигнуто наиболее оптимальное распределение узлов в расчетной сетке. Первый расчет выполняется на очень грубой начальной сетке, которая может быть получена автоматически. Затем, в соответствии с полученным решением, выполняется процедура адаптации, которая, с одной стороны, изменяет положение существующих узлов и их связей, с другой — добавляет (удаляет) определенное количество узлов. Процесс решение — адаптация — решение повторяется до достижения заданного малого изменения рассчитываемых характеристик течения от итерации к итерации на последовательности адаптивных сеток.
В настоящей работе применяется метод и программа адаптации сетки к решению, который был разработан А. Мартыновым и С. Медведевым (ИПМ им. Келдыша) [11]. Процесс адптации основан на эвристическом индикаторе ошибки, зависящем от линейной комбинации первой и второй производных местного числа Маха вдоль ребра. Вся сетка разделена на макроячейки, дробление (укрупнение) которых осуществляется изотропно. Узлы на ребрах макроячеек распределяются с учетом анизотропных свойств течения.
2.1. Трансзвуковое обтекание профиля NACA 64A010. Ниже представлены некоторые результаты расчета вязкого трансзвукового обтекания изолированного профиля NACA 64A010
при M = 0.8, Re = 2-106, а = 6.2° [12]. Этот режим характеризуется мощным скачком
уплотнения
с последующим развитием отрыва турбулентного пограничного слоя. На рис. 1 приведено значение коэффициента подъемной силы в зависимости от числа степеней свободы в адаптивной неструктурированной сетке. Видно, что результат стабилизируется уже после 6 итераций адаптации.
Фрагменты самой мелкой сетки (10-я итерация адаптации) приведены на рис. 2. Эта сетка содержит 18 506 узлов и 36 604 элементов. Видно, что расположение сеточных узлов в точности повторяет особенности течения (скачок уплотнения, пограничный слой). Расстояние от тела до первого
+2
слоя сеточных узлов, выраженное в долях y , варьируется в пределах от 1 до 2, достигая значения 7 лишь в ячейках, примыкающих к задней кромке профиля, что говорит о достаточно точном описании течения в пристеночной части пограничного слоя.
С целью монотонности решения в окрестности скачка уплотнения порядок аппроксимации искомых функций в ячейках, попадающих в рассматриваемую окрестность, понижается до первого, что соответствует кусочно-постоянной реконструкции решения. На рис. 3 черным цветом помечены ячейки, где порядок схемы был понижен. Таким образом, продемонстрирована линия скачка уплотнения. Видно, что скачок уплотнения имеет форму греческой буквы Я, что хорошо соответствует картине течения, наблюдаемой в эксперименте.
Су
50000 100000
NDOF
Рис. 1. Значение коэффициента подъемной силы в зависимости от числа степеней свободы в адаптивной неструктурированной сетке. Обтекание изолированного профиля NACA 64A010 при
М= 0.8, Re = 2-106, а = 6.2°
Рис. 2. Фрагменты адаптивной расчетной сетки для профиля ЫАСА 64А010 после десяти итераций адаптации
Рис. 4. Распределение коэффициента давления на поверхности профиля NACA 64A010 при M = 0.8,
Re = 2-106, а = 6.2°:
Рис. 3. Линия скачка уплотнения (черным цветом изображены ячейки с кусочнопостоянным восполнением решения)
- эксперимент [12]; Сурасч = 0.637;
cy = 0.62
уэксп
Распределение коэффициента давления на поверхности профиля в сравнении с экспериментом [12] представлено на рис. 4. Некоторое рассогласование расчетных и экспериментальных значений объясняется несовершенством используемой модели турбулентности и, вероятно, существенно пространственным характером турбулентного обтекания в эксперименте.
2.2. Расчет обтекания трехэлементной взлетно-посадочной механизации на последовательности адаптивных неструктурированных сеток. Расчет обтекания трехэлементной взлетно-посадочной механизации McDonnell Douglas Corporation (MCD) [13] также был выполнен
на последовательностях адаптивных неструктурированных сеток. Фрагменты полученной адаптивной расчетной сетки представлены на рис. 5. Расчеты выполнены при числах M = 0.2 и
Re = 9 • 106, что соответствует экспериментальным условиям [13].
На рис. 6 приведены профили давления торможения в следе над закрылком (х/с = 0.975)
на углах атаки а = 8, 16, 21° в сравнении с результатами эксперимента. Проведенное сравнение демонстрирует, что образованные в процессе расчета адаптивные сетки позволяют достаточно точно разрешить след, образуемый различными элементами механизации.
ч/
о ОС II а 43 II 3
Сечение А Сечение А
0.2 0.2
/; 0.1 0.1
(V I
. , 1 ■ 1 , 1 .
0 0 5
V
0 с.„ 0 5
Р1
Рис. 7. Зависимость коэффициента подъемной силы от угла атаки для трехэлементного профиля MCD:
-------расчет;
А А А — эксперимент [13]
Рис. 6. Профили давления торможения в следе над закрылком (х/с = 0.975) для трехэлементного профиля MCD на
углах атаки а = 8, 16, 21°:
- расчет; ° ° ° — эксперимент [13]
Рис. 5. Фрагменты адаптивной расчетной сетки, состоящей только из треугольных элементов, для трехэлементного
профиля MCD после восьми итераций адаптации:
М = 0.2, а = 24°, Яе = 9-106
Рис. 8. Поле местных чисел Маха вблизи носовой части предкрылка для трехэлементного профиля MCD:
M = 0.2, а = 24° (критический угол), Re = 9 • 106
Представляющая наибольший практический интерес зависимость величины подъемной силы от угла атаки представлена на рис. 7. Видно, что расчет достаточно точно предсказывает критический угол атаки и величину максимальной подъемной силы. Несмотря на то, что число Маха набегающего потока невелико, обтекание данной конфигурации характеризуется наличием небольшой местной сверхзвуковой зоны, замыкающейся скачком уплотнения, на верхней поверхности предкрылка. Использование адаптивных сеток позволило обнаружить эту важную особенность течения (рис. 8).
Работа частично выполнена при поддержке Российского фонда фундаментальных исследований (проект № 06-01-00283-а).
ЛИТЕРАТУРА
1. Venkatakrishnan V., Allmaras S. R., Kamenetskij D. S., Johnson F. T.
Higher order schemes for the compressible Navier — Stokes equations // AIAA 2003-3987.
2. Hughes T. J. R., Brooks A. A multidimensional upwind scheme with no crosswind diffusion, in Finite element methods for convection dominated flows. — New York: ASME. 1979.
3. Lesaint P. and R a v i a r t P. A. On a finite element method to solve the neutron transport equation / Mathematical Aspects of Finite Elements in Partial Differential Equations, edited by C. de Boor. — New York: Academic Press. 1974.
4. Cokburn B. Discontinuous Galerkin methods for convection-dominated problems /
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.
5. Spalart P. R. and Allmaras S. R. A one-equation turbulent model for aerodynamic flows // AIAA-92-0439.
6. B a s s i F. and R e b a y S. A high-order accurate discontinuous finite element method for numerical solution of the compressible Navier — Stokes equations // J. of Comp. Phys. 1977.
V. 131.
7. Roe P. L. Approximate Riemann Solvers, Parameter Vectors, and Difference Schemes //
J. of Comp. Phys. 1981. V. 43.
8. Osher S., Solomon F. Upwind Difference Schemes for Hyperbolic Systems of Conservation Laws // Mathematic of Computation. April 1982. V. 38, N 158.
9. Cockburn B., Hou S., Shu C. TVB Runge-Kutta local proj ection discontinuous Galerkin finite element method for conservation laws IV: the multidimensional case // Math. Comp.
1990, N 54.
10. http://www.mcs.anl.gov/petsc.
11. Martynov A. A., Medvedev S. Yu. A robust method of anisotropic grid generator. Workshop «GRID GENERATION: THEORY AND APPLICATIONS» Moscow,
June 24 —28. 2002. Computing Center Russian Academy of Sciences,
http://www.ccas.ru/gridgen/ggta02/.index.html.
12. Johnson D. A., Bachalo W. D. Transonic flow past a symmetrical airfoil — invis-cid and turbulent flow properties // AIAA J. 1980. V. 18.
13. Chin V. D., Peters D. W., Spaid F. W., McGhee R. J. Flowfield measurements about a multi-element airfoil at high Reynolds numbers // AIAA-93-3137.
PyKonucb nocmynuna 11/IV 2006 г.
