CC BY
165
Том XL
УЧЕНЫЕ ЗАПИСКИ ЦАГИ 2009
№ б
УДК 532.516.2 519.63:517.96
ОСОБЕННОСТИ ПРИМЕНЕНИЯ МЕТОДА ГАЛЕРКИНА К РЕШЕНИЮ ПРОСТРАНСТВЕННЫХ УРАВНЕНИЙ НАВЬЕ — СТОКСА НА НЕСТРУКТУРИРОВАННЫХ ГЕКСАЭДРАЛЬНЫХ СЕТКАХ
А. В. ВОЛКОВ
Метод Галеркина с разрывными базисными функциями (РМГ), ранее исследованный в работах [1—3], был адаптирован к решению пространственных уравнений Эйлера и Навье — Стокса на основе использования неструктурированных гексаэдральных сеток. Реализованный алгоритм позволяет выполнять расчеты с порядком точности 4 и учитывать кривизну обтекаемой границы. В качестве метода ускорения решения используется гибридный подход, сочетающий в себе p-многосеточный метод и традиционный агломерационный h-многосеточный метод. Тестирование порядка точности алгоритма и оценка потребных компьютерных ресурсов были осуществлены в процессе решения различных тестовых задач.
В работе приводятся примеры решения задач о невязком обтекании цилиндра, ламинарном течении около пластины, пространственном вязком течении внутри изогнутой трубки, пространственном турбулентном обтекании изолированного крыла. Рассмотрено также решение задачи о распространении акустической сферической волны, описываемой в рамках линеаризованных уравнений Эйлера. Результаты расчетов и компьютерные ресурсы сопоставлены с результатами, получаемыми методом конечного объема.
Ключевые слова: метод Галеркина с разрывными базисными функциями, метод конечного объема, схема высокого порядка точности.
Повышение точности расчета обтекания тел сложных геометрических форм требует использования очень мелких расчетных сеток. Убедительная демонстрация этого факта приводится в трудах известного форума, посвященного точному расчету сопротивления крейсерской конфигурации летательного аппарата (Drag Prediction Workshops — DPW [4]). Здесь отмечается, что даже сетки с количеством ячеек более чем 20 миллионов узлов оказываются недостаточными для описания таких относительно простых геометрий, как «крыло + фюзеляж». Определяющим воздействием на результат являются не только тип и густота расчетных сеток, но и используемая модель турбулентности. В связи с этим, как правило, исследуется влияние различных параметров сеток и параметров моделей турбулентности на увеличение точности расчетов. При этом вопрос о порядке точности схем не обсуждается, так как все промышленные расчетные коды базируются на конечно-объемных методах второго порядка точности.
Вопрос о применимости расчетной схемы высокого порядка точности к решению задач внешней аэродинамики обсуждался, например, в работе [5]. Здесь изучалась возможность создания надежного, высокоточного алгоритма решения уравнений Рейнольдса при обтекании тел сложной геометрической формы. В качестве примера рассматривался расчет самолета Boeing 747-200 на базе блочно-структурированных сеток.
Одним из наиболее перспективных подходов к высокоточной аппроксимации как на структурированных, так и неструктурированных сетках является метод Галеркина с разрывными базисными функциями. В англоязычной литературе метод имеет название Discontinuous Galerkin Method (DGM) [6, 7]. В последние годы этот метод вызывает большой интерес многих исследователей вследствие его общности, гибкости и надежной теоретической обоснованности. Поэтому
неудивительно, что метод был интенсивно исследован применительно к решению пространственных течений. Например, на структурированных сетках РМГ был использован в работах [8—10], а на тетраэдальных неструктурированных сетках он был применен в работе [11]. Конечно, и другие методы высокого порядка точности могут быть использованы в качестве альтернативы методу конечного объема второго порядка точности. В работе [12] дан обзор различных подходов к аппроксимации высокого порядка точности.
Как было показано ранее [13], при расчете пространственных течений неструктурированные гексаэдральные сетки имеют значительное преимущество по сравнению с сетками, образованными тетраэдрами. Поэтому РМГ был адаптирован к решению уравнений Эйлера и Навье — Стокса именно на неструктурированных гексаэдральных сетках.
Настоящая работа представляет результаты дальнейшего развития работ [1—3]. Метод Га-леркина с разрывными базисными функциями впервые был применен к решению пространственных уравнений Навье — Стокса с использованием неструктурированных гексаэдральных сеток. Алгоритм расчетной схемы позволяет выполнять расчеты с порядком точности 4. С целью достижения высоких порядков точности в расчетную схему была включена возможность учета кривизны обтекаемой поверхности. Был внедрен оригинальный многосеточный метод ускорения решения (р-многосеточный метод в комбинации с традиционным агломерационным Ь-много-сеточным подходом). В работе на примере решения различных тестовых задач приводятся результаты тестирования порядка точности полученного алгоритма и оценка потребных компьютерных ресурсов. Рассмотрены решения задач о невязком обтекании цилиндра, ламинарном течении около пластины, пространственном вязком течении внутри изогнутой трубки, пространственном турбулентном обтекании изолированного крыла. Приведено также решение задачи о распространении акустической сферической волны, описываемой в рамках линеаризованных уравнений Эйлера.
Результаты расчетов и компьютерные ресурсы сопоставлены с результатами промышленного кода РГКЕ™/Неха (компании КИМЕСЛ Гп1;., Бельгия), основанного на широко используемом в настоящее время методе конечного объема второго порядка точности.
1. РМГ и решаемые уравнения. Рассмотрим систему уравнений Навье — Стокса, записанную в консервативной форме:
поток и ^(и, Уи) — вязкий поток; 8 — источниковый член, возникающий при использовании
модели турбулентности. Система уравнений (1) решается численно на основе использования метода Галеркина с разрывными базисными функциями (РМГ) и решение в каждой ячейке хранится в примитивных переменных О = (, и, V, м>, р). При этом величины давления р и полной энергии
Е связаны уравнением состояния.
Уравнения Навье — Стокса могут быть переписаны в примитивных переменных следующим образом:
(1)
Здесь и — вектор консервативных переменных и = (р, ри, pv, р^, рЕ), Е(и) — невязкий
(2)
где матрица
есть Якобиан преобразования от консервативных переменных к примитивным.
В каждой ячейке сетки определяются локальные полиномиальные базисные функции:
о( X) = 2иу ()фу (х),
і=1
где и у (у) — коэффициенты разложения, которые должны быть определены в процессе решения; Ку — количество базисных функций в ячейке. Величина Ку связана с максимальной степенью базисного полинома К следующим образом:
кг =
(к + 1)(к + 2 )(к + 3)
6
(5)
В настоящей работе полиномиальный порядок К принимал значения К =0, 1, 2, 3, и таким образом число базисных функций принимало значения К^ = 1, 4, 10, 20. Для К =3 использовался
следующий набор базисных функций:
І о 1 2 3 4 5 6 7 8 9
ф І 1 х у г х2 у2 2 г ху хг уг
і іо 11 12 13 14 15 16 17 18 19
ф і х3 у3 3 г х2 у х2 г у 2 х у2 г г 2х г 2 у хуг
Отметим, что базисные функции были нормализованы следующим способом:
Фі
(х - х0)
Ф8
(х - Хо )( - го )
Ф9 =
( - Л )г - го )
Ф19 =
( - хо )(у -уо )( - го)
(6)
У'г
х у г
ны
Здесь х0, у0, г0 — координаты центральной точки рассматриваемого гексаэдра, а величи-X, Ну, Н2 определяют размер ячейки вдоль соответствующих осей. Было показано, что именно такая нормализация обеспечивает лучшую точность и сходимость особенно в случае рассмотрения вязких задач с сильно вытянутыми ячейками сетки.
Система сеточных уравнений для коэффициентов и у (У) формулы (4) получается в соответствии со стандартной процедурой метода конечного элемента Галеркина, когда требуется ортогональность невязки решаемых уравнений (левая часть системы (2)) каждой базисной функции, используемой в реконструкции решения. Это требование ортогональности формулируется через условие равенства нулю интеграла от произведения системы решаемых уравнений на каждую из
базисных функций фг- ( = 1,..., К^). После интегрирования по частям ]
: имеем:
-^ф-ГО—О = - £ Фг ( -Е ))Е + ^Уф,- ( - Е)П + ^ф^—П. (7)
Здесь —X — элемент площади, ориентированной в направлении нормали п = [пх, пу, п2), и
—П — элемент объема ячейки.
Уравнение (7) состоит из объемных интегралов и поверхностных интегралов по границе ячейки. Значения всех зависимых переменных терпят разрыв на границе элементов и ключевую роль играют правила вычисления значений переменных и потоков на этих границах.
Как и в методе конечного объема (МКО), в РМГ величина невязкого потока на грани между двумя ячейками определяется как результат решения задачи Римана о распаде произвольного
И
х
разрыва. В настоящей работе используется приближенная, линеаризованная методика решения задачи Римана, предложенная Роу:
Ерань = 1 (г( + ЕЯ ) -)А|(и( - иЯ), (8)
здесь верхние символы «Ь» и «Я» означают, что соответствующие величины вычисляются соответственно с левой и с правой стороны интерфейса, и матрица А является матрицей Якоби невязкого потока:
А = —. (9)
эи
Вязкие потоки определяются через градиенты примитивных переменных Еу = г (оН, уд Н),
У^ (ддн д—Н д—Н ^ г г.
где У—Н = | ——, —--, —;— |. Градиенты примитивных переменных могут быть найдены путем
дх ду дх
непосредственного дифференцирования формулы (4) для реконструкции решения в ячейке. Однако, как было показано в работе [7], такой способ вычисления градиентов является причиной отсутствия аппроксимации. Поэтому в РМГ для вычисления вязких потоков градиенты примитивных переменных также представляются в виде линейной комбинации базисных функций:
эо Кт
-д-( х) = X§',у(У)Фу(х). (10)
дХ, у=1
Здесь ' = 1, 2, 3 соответствует х, у и х координатам. После умножения уравнения (10) на тестовые функции (для метода Галеркина тестовые функции являются базисными функциями) и интегрирования по частям получаем следующую систему линейных уравнений относительно коэффициентов разложения Ц' у:
КГ д
|2&-,уфуф*—О+£ ф^° п-х-1 ~ф^о—О=0; к=1,..., Кт. (11)
О у =1 X о '
Решение этой системы определяет градиенты примитивных переменных (10) и дает возможность вычислить вязкие потоки в ячейке.
В РМГ как градиенты примитивных переменных, так и значения потоков имеют разрывы на границах ячеек. Однако эти значения требуются при вычислении контурных интегралов и они не могут быть получены осреднением. Правила аппроксимации [7] требуют переменного выбора левой или правой ячейки для реконструкции решения на границе. Так, если при вычислении вязких потоков в (7) значения примитивных переменных и их градиентов были выбраны из левой
ячейки
Е^ = Еу(-Ь, У-Ь), (12)
то в формуле (11) при вычислении контурных интегралов значения примитивных переменных следует брать из правой ячейки
-| ь = -Я. (13)
Наконец, система нелинейных уравнений для коэффициентов иг- (у) формулы (4) получается из уравнения (7) в предположении малой вариации якобиана Г (3) внутри ячейки:
—У
Здесь М — матрица масс, ' = 1, ., Ку.
£ ф ( - Еу ))Х + |оУфг ( - Еу ))О + |офг Ъ—П
(14)
2. Вычисление интегралов и учет кривизны. Высокий порядок аппроксимации РМГ предполагает точное вычисление объемных и поверхностных интегралов в системе уравнений (14). Для вычисления интегралов используются квадратурные формулы Гаусса. Вычисление объемного интеграла внутри гексаэдра произвольной формы выполняется в параметрическом пространстве, в котором каждый гексаэдральный элемент имеет форму единичного куба. Это триквадратичное изопараметрическое преобразование построено на базе 20 узлов, из которых 8 — это вершины куба, а 12 являются центрами всех ребер куба. Наличие центральных узлов на ребрах позволяет учитывать кривизну обтекаемой границы. При таком подходе поверхностные и объемные интегралы содержат нелинейные якобианы преобразования, что накладывает более высокие требования к квадратурным формулам интегрирования и требует большего количества гауссовых точек. Таким образом, вычисление объемных интегралов потребовало применение квадратур, точных для полиномов порядка 3К.
Отметим, что объем вычислений может быть сокращен при использовании безквадратурно-го подхода, предложенного в [14]. Этот подход был использован только на грубых уровнях многосеточного решателя.
Общеизвестно, что получение высокого порядка точности требует учета кривизны обтекаемой границы. Применения 20-точечного преобразования для граничных интегралов можно избежать при расчете течений невязкой жидкости. При этом можно использовать оригинальный подход, предложенный в работе [15]. В невязком течении граничное условие на поверхности обтекаемого тела формулируется в виде требования равенства нулю компоненты скорости, нормальной к телу. В [15] показано, что для корректного учета кривизны в невязком течении достаточно использовать только «правильное» направление нормалей, ориентированных строго нормально к истинной криволинейной границе в каждой гауссовой точке. При этом само интегрирование выполняется на гексаэдре без учета его искривленности у границы.
Правильное направление нормалей к обтекаемой границе может быть получено либо в процессе построения биквадратичного преобразования искривленного поверхностного четырехугольника в плоский квадрат, либо вычислением нормалей к оригинальной поверхности, заданной в системе автоматизированного проектирования.
На рис. 1 представлена аппроксимация передней кромки сечения крыла ЬАМК [21]. Сравниваются полигональное представление контура с биквадратичной интерполяцией поверхности и реальной геометрией.
0 0.01 0.02 х
Рис. 1. Передняя кромка крыла ЬАМЫ". Сравнение биквадратичной интерполяции поверхности с полигональной границей сетки и реальной поверхностью
1СГ2
1СГ6
А-----А------А Полигональная граница
ООО Биквадратичная интерполяция
10“'
102
103
Число граней цилиндра
Рис. 2. Ошибка в определении нормали к поверхности цилиндра как функция числа граничных точек
На рис. 2 представлена величина ошибки в вычислении нормали к поверхности цилиндра в случае полигонального представления границы и биквадратичной интерполяции. Здесь величи-
линдра. Буквой N обозначена вычисленная нормаль в точке /, малой пг- — точная нормаль и М — общее количество контрольных точек на поверхности цилиндра. Видно, что биквадратич-ная интерполяция обеспечивает более точное описание поверхности.
3. Решение системы сеточных уравнений и многосеточный метод ускорения решения. Решение системы сеточных уравнений (14) осуществляется явным методом интегрирования по времени. Поиск стационарного решения позволяет использовать локальный шаг по времени, зависящий от размера ячейки, что заметно ускоряет процесс сходимости — стремление невязки (правой части уравнения (7)) к машинному нулю. Используется 5-шаговый метод Рунге — Кутта, обеспечивающий наибольшую стабильность метода, или, другими словами, быстрое убывание невязки в процессе итераций по времени.
Применение явного рунге-куттовского метода поиска решения практически не требует дополнительных ресурсов памяти, однако, несмотря на локальный шаг по времени и наличие многошаговой схемы, количество итераций, требующихся для достижения стационарного решения, слишком велико. Значительное ускорение решения достигается за счет применения многосеточного метода.
Многосеточный метод является эффективным подходом к ускорению решения задач математической физики. Основная идея метода состоит в быстрой передаче информации между различными частями расчетной сетки. Для этой цели в классическом многосеточном методе (в англоязычной литературе его называют h-multigrid) решение определяющих уравнений выполняется на наборе последовательно измельченных сеток. Решение на грубой сетке используется для определения решения на мелкой сетке, и наоборот, решение на мелкой сетке используется как начальное приближение для решения на грубой сетке. При этом управляющие уравнения на грубой сетке модифицируются с целью учета неточности аппроксимации на крупных ячейках. На каждом сеточном уровне метод Рунге — Кутта (или другой метод решения) обеспечивает эффектив-
представлена в зависимости от числа разбиений поверхности ци-
ное подавление высокочастотных ошибок. Но высокочастотные ошибки на сетке с крупными ячейками являются низкочастотными ошибками для мелкой сетки. Таким образом, малое количество переменных на грубой сетке дает возможность быстро разрешить все низкочастотные особенности течения. Перевод решения, полученного на грубой сетке, к решению на мелкой сетке позволяет ускорить передачу информации и обеспечить коррекцию низкочастотной составляющей основного решения.
В методе конечного элемента решение в каждом элементе является линейной комбинацией базисных функций. Эта комбинация образует полиномиальную реконструкцию с максимальной степенью К. Коэффициенты перед базисными функциями являются искомыми степенями свободы. В настоящей работе используются иерархические базисные функции. Это означает, что коэффициенты перед базисными функциями в разложении решения имеют ясный математический смысл, а именно: они отражают величину осредненного решения в ячейке, величину градиента решения по трем направлениям, вторую производную решения и т. д. Таким образом, при данном выборе базисных функций решение представляется в виде тейлоровского разложения решения в ячейке. Или, другими словами, коэффициенты в реконструкции решения отражают вклад различных гармоник в решение.
Идеи р-многосеточного метода (p-multigrid — в англоязычной литературе [16—18]) для конечно-элементного подхода аналогичны идеям классического многосеточного метода, описанного выше. Уменьшение или увеличение максимального порядка (К) полиномов базисных функций или, другими словами, изменение набора базисных функций в элементе аналогично укрупнению или измельчению ячеек в сетке. Использование базисных функций низкого порядка (К =0) позволяет быстро разрешить низкочастотную компоненту решения и обеспечить коррекцию общего решения при максимальном наборе базисных функций.
В настоящей работе наибольший выигрыш в скорости сходимости был получен в результате комбинирования традиционного многосеточного метода и р-многосеточного метода. В основе традиционного многосеточного метода был использован агломерационный многосеточный метод [19], реализованный в решателе промышленного кода [22]. Этот метод применялся при решении с использованием простейшей схемы аппроксимации (К =0 — одна кусочно-постоянная функция). Затем на самой мелкой сетке выполнялось постепенное наращивание количества используемых базисных функций, т. е. был реализован р-многосеточный метод.
4. Тестовые расчеты. Свойства реализованной схемы РМГ, достигающей четвертый порядок точности (К =3), изучены в процессе решения серии тестовых задач, охватывающих пространственные невязкие, вязкие и турбулентные течения. Результаты были сопоставлены с данными, получаемыми МКО [22].
В МКО второго порядка точности для ускорения решения использовался агломерационный трехуровневый многосеточный метод. В РМГ (К =1) был задействован четырехуровневый метод. При этом первые три уровня являлись традиционным многосеточным методом, где использовались кусочно-постоянные базисные функции (К =0). На четвертом уровне аппроксимация систем уравнений осуществлялась на мелкой сетке, но на базе кусочно-линейных функций (К =1). Решение для РМГ (К =2) находилось с использованием пятого уровня, а шестой уровень задействовался для РМГ (К =3). Во всех случаях в качестве сглаживающей процедуры для высокочастотной составляющей ошибки использовался 5-шаговый метод Рунге — Кутта.
Невязкое обтекание цилиндра при числе М = 0.15. Невязкое потенциальное течение около кругового цилиндра является известным тестовым случаем, позволяющим определить порядок точности численной схемы. Расчет выполняется при малом числе Маха набегающего потока М = 0.15, когда влияние сжимаемости пренебрежимо мало. Все численные эксперименты были выполнены на последовательности из четырех аналитических сеток, размер которых варьировался в диапазоне от 16 X 4 до 128 X 32, путем удвоения количества точек по каждому из направлений (рис. 3). N концентрических окружностей содержит М равномерно распределенных узлов в плоскости X — У. Радиус каждой из окружностей определяется следующим соотношением:
, ] = 0,..., N.
(15)
XaV" /X \ / \ \ \ / \ \ \ \ / \ VV 1 '^"-аГ4 V 1 Ту/ ' х х X / \ / / \ г "V ЧУ\\ /А/ / \/ \ \ \ \ / / / \у К\\\тт////Л, - Ас А ХА \ \ / X А Зг^ / Уг/УхК, \ А-"'"
\ \ \ А \ / \ \ / \ \ Г-ЛГ \ \ / 4 / \ \ / \ /А V \ \—А- \ >AfAxM А ААУ/ А'\ А- V / / ааДа \ \ А/ Ах // 1 \ \ \ А /А ^ X / \ \ X А ХА / / \ \ ХА х
\ \Х\\\ \Д4~ AXaXI ^ АС A АХ \ \ \ ЛХ-г ^ААAAXVvYuT /А A AXXXW\ 11 §11 ХХаЖ / АА / / / гтЧ-А / /V /// ГГ 7/7а//аа А\ Т/аааХ/ ХШаАаААА.АА^ Вш- ITAAXAA 1 1 \ - > ч ьч
Рис. 3. Сетки размером 16 х 4, 32 х 8, 64 х 16 и 128 х 32 ячеек около цилиндра
а) б)
Рис. 4. Распределение давления на поверхности цилиндра:
а — численное решение для РМГ К =1, 2, 3 и МКО на грубой сетке 16 X 4; б — численное решение для схем с равным количеством
степеней свободы
Порядки точности МКО и РМГ К =1, 2, 3 схем
Сетка Количе- ство элемен- тов МКО РМГ К — 1 РМГ К — 2 РМГ К —3
Сх ошибка энтропии По- ря- док Сх ошибка энтро- пии По- ря- док Сх ошибка энтро- пии По- ря- док Сх ошибка энтро- пии По- ря- док
16 X 4 64 1.18е-1 1.03е-2 7.47е-2 1.00е-2 8.94е-3 5.55е-4 2.66е-3 1.13е-4
32 X 8 256 2.73е-2 3.11е-3 1.73 7.34е-3 2.73е-3 1.88 7.01е-4 7.79е-5 3.25 5.83е-5 5.26е-6 4.42
64 X 16 1024 1.4е-3 6.85е-4 2.18 2.75е-3 6.05е-4 2.17 2.37е-4 7.50е-6 2.95 6.92е-6 3.66е-7 3.85
128 X 32 4096 9.05е-4 1.3 6е-4 2.33 1.56е-3 1.18е-4 2.36 3.74е-5 8.93е-7 3.07 1.95е-6 2.12е-8 4.11
Радиус цилиндра равен 0.5. Параметр а в соотношении (15) определялся из условия, что радиус наибольшей окружности, ограничивающий расчетную область, равен 20. На внешней границе накладывалось условие дальнего поля, основанное на инвариантах Римана.
На рис. 4 представлен расчетный коэффициент распределения давления в сравнении с аналитическим несжимаемым решением (точечная линия). На рис. 4, а изображены эпюры распределения давления вдоль поверхности цилиндра, полученные на самой грубой расчетной сетке (16 X 4) методом РМГ второго (К —1), третьего (К —2) и четвертого (К —3) порядков точности. Здесь же представлены результаты расчетов по МКО [22]. На рис. 4, б приводится сравнение эпюр, полученных на различных сетках при условии использования приблизительно равного количества неизвестных, или, другими словами, числа степеней свободы задачи ^, которое определяется произведением числа ячеек на количество переменных в ячейке. Представленные результаты демонстрируют, что РМГ обеспечивает лучшую точность даже в случае использования эквивалентного числа ^.
Порядок точности численной схемы оценивался по величине ошибки энтропии:
е =_^р!- 1
^энт ' *• •
../р!
Так как рассматриваемое течение изоэнтропическое, то при точном решении энтропия должна быть равна нулю. В табл. 1 представлены: размер сетки, общее количество элементов в ней, — норма энтропийной ошибки, вычисленной по каждому элементу сетки, величина со-
противления цилиндра и полученный порядок точности численной схемы, который оказался близким к теоретическому значению О(Ик+1). Порядок точности рассчитывался по величине
убывания ошибки при последовательном переходе от крупной сетки (/ — 1) к более мелкой (/) с использованием следующего соотношения:
Порядок — 2-
е
Вязкое течение около цилиндра при числах М = 0.15, Яе = 40. Расчет вязкого ламинарного течения был выполнен на последовательности сеток, описанных в предыдущем разделе, при Яе = 40. Порядок точности был оценен на основе величины полного сопротивления цилиндра. Анализ асимптотического поведения величины полного сопротивления был выполнен в зависимости от осредненного размера сеточной ячейки, который рассчитывался по формуле:
И = 1/л/^.
Отметим, что при таком определении осредненный шаг сетки зависит не только от ее размера, но и от метода аппроксимации. Для метода конечного объема, где используется только од-
на переменная в ячейке, на равномерной прямоугольной сетке усредненныи размер ячеики идентичен действительному шагу сетки.
Асимптотически предельная величина сопротивления с * = 1.5685 получена путем экстраполяции значений с двух наиболее мелких сеток к сетке с бесконечно малым шагом. Предпола-■ - ^(И—,)
— с * + О (
х
имеем:
Порядок — 2-
4"1 - -
Результаты расчетов, представленные в табл. 2, демонстрируют ожидаемый порядок точности всех рассмотренных численных схем.
Таблица 2
Коэффициент сопротивления и оценка порядка точности схем в задаче о вязком обтекании цилиндра
Яе = 40 для МКО и для РМГ К =1, 2, 3
Сетка Количество элементов МКО РМГ К — 1 РМГ К —2 РМГ К —3
Сх Порядок Сх Порядок Сх Порядок Сх Порядок
16 X 4 64 0.403598 0.298054 0.190749 0.163027
32 X 8 256 0.271965 1.10 0.191675 2.02 0.161390 2.90 0.156464 3.99
64 X 16 1024 0.180503 2.30 0.166170 1.90 0.157401 2.05 0.156829 4.00
128 X 32 4096 0.161464 2.46 0.159011 2.11 0.156790 3.11
Ламинарное обтекание плоской пластины при числах М = 0.35, Ке = 76 000. Расчетная область, показанная на рис. 5, а, имела следующие размеры: -0.24 < х < 0.4; 0 < у < 0.05;
0 < 2 < 0.00625. Пластина располагалась в диапазоне 0 < х < 0.4. Численные исследования были выполнены на четырех сгенерированных сетках, представленных на рис. 5. Это двумерные сетки,
Рис. 5. Расчетная область и фрагменты сеток: а — сетка 1; б — сетки 2, 3 и 4
Параметры сеток для расчета ламинарного обтекания пластины
Сетка Количество ячеек Дистанция до первого слоя у+ Коэффициент расширения слоев Кол-во сеточных линий в у-направлении
1 6739 110-5 1.02 25
2 1707 64 110-5 1.20 13
3 680 64 210-5 2.00 8
4 281 32 410-5 2.50 5
имеющие только одну ячейку в г-направлении. Некоторые параметры сеток, представленные в табл. 3, показывают количество ячеек в сетке, расстояние от поверхности стенки до первого слоя ячеек, коэффициент увеличения расстояния между последовательными рядами сеточных точек и общее количество узлов в у-направлении. Количество элементов в сетке 1 приблизительно в 4 раза больше, чем в сетке 2, в то время как сетки 3 и 4 имеют соответственно в 10 и 20 раз меньшее количество сеточных ячеек. Такой выбор сеток обеспечивает возможность проводить сравнения различных численных схем в условиях близкого количества степеней свободы.
Некоторые результаты расчетов представлены на рис. 6 и 7 в виде профилей скорости в сравнении с аналитическим решением Блазиуса, полученным для несжимаемого течения. Все
Рис. 6. Сравнение численных схем на различных сетках в условиях одинакового количества степеней
свободы:
а — ^-компонента скорости; б — Г-компонента скорости
Рис. 7. Эффект увеличения порядка точности расчетной схемы от К =2 до 3 для Г-компоненты скорости
Рис. 8. Сравнение аналитического решения с решением РМГ К = 3, полученным на грубой сетке 4 для коэффициента локального трения
зависимости представлены комбинацией линий и маркеров. Каждая линия построена через множество точек и в точном соответствии с формулой для реконструкции решения в ячейке. Маркеры же на каждой линии расположены приблизительно в серединах ячеек. Таким образом, количество маркеров на линии отражает густоту сетки.
^-компонента скорости показана на рис. 6 для всех исследуемых схем аппроксимации. Все схемы демонстрируют хорошее согласование с аналитическим решением. Г-компонента скорости, представленная на рис. 6, б, более чувствительна к качеству сетки. Результаты, полученные с РМГ К =2 на сетке 3 и К =3 на сетке 4, находятся в лучшем согласовании с аналитическим
решением, чем результаты, полученные традиционным МКО на самой мелкой сетке 1, несмотря на эквивалентное количество задействованных степеней свободы.
На рис. 7 показано влияние порядка точности схемы на Г-компоненту решения. Расчет на самой грубой сетке 4 выполнен с РМГ К =2 и К =3. Здесь наблюдается заметное улучшение согласования, несмотря на экстремально малое количество узлов расчетной сетки, расположенных поперек пограничного слоя.
Сравнение распределения локального коэффициента трения, полученного на самой грубой сетке 4, с аналитическим решением Блазиуса показано на рис. 8. Видно, что, несмотря на использование экстремально грубой сетки, удалось получить хорошее согласование расчетов с точным решением.
Пространственное ламинарное течение в изогнутой трубке. Этот пространственный тестовый случай для ламинарного течения внутри изогнутой под 90° трубки с постоянным сечением был экспериментально исследован [20] лазеро-доплеровским измерителем скорости с водой в качестве рабочей жидкости. Исследования проводились при трех различных числах Рейнольдса, включая ламинарный случай при Яе = 500, который и рассмотрен в настоящей работе. Течение, установившееся в трубке, определяемое балансом вязких и невязких сил, характеризуется наличием пары вихрей, закрученных в противоположные стороны, образуемых вниз по потоку после прохождения трубочного изгиба. Имеются экспериментально измеренные профили скорости в пяти сечениях вдоль трубки: сечение 1 отстоит на 0.58 диаметра вверх по потоку от изгиба; сечение 2 — 30° по изгибу; сечение 3 — 60° по изгибу; сечение 4 — 75° по изгибу; сечение 5 — один диаметр вниз по потоку от изгиба.
Пять неструктурированных гексаэдральных сеток было сгенерировано с использованием генератора сеток НЕХРКЕ88™ [13]. Количество ячеек в сетках варьировалось от мелкой до самой грубой в диапазоне 140 000 ^2500. Некоторые фрагменты сеток представлены на рис. 9. Для сетки 5 пунктирной линией нанесено изображение поверхности трубки. Здесь особенно наглядно видна разница между сеточной и реальной границей. Учет кривизны обтекаемой поверхности позволил учесть эту разницу.
Расчеты были выполнены для всех сеток с использованием различных схем аппроксимации. Профили скорости в различных сечениях трубки представлены на рис. 10. Здесь результаты, полученные на самой грубой сетке со схемами РМГ К =2 и 3, сопоставлены с экспериментальными данными. Видно, что расчеты по схеме РМГ (К =3) имеют хорошее согласование с экспериментальными профилями скорости. Кроме того, на этом рисунке демонстрируется эффект от увеличения порядка точности схемы.
Рис. 10. Профиль скорости в четырех сечениях трубки. Сопоставление эксперимента и численных результатов,
полученных с РМГ К =2 и 3 на грубой сетке 5
Преимущества схемы четвертого порядка точности также иллюстрируются на рис. 11. Здесь профили скорости, рассчитанные по схеме РМГ К =3 на сетке в 2500 ячеек, сравниваются с результатами схемы МКО [22], полученными на густой сетке, содержащей 140 000 элементов. Отметим, что расчеты, выполненные по традиционной схеме МКО на сетке с числом элементов 62 000 весьма далеки от экспериментальных данных. Таким образом, результаты со схемой МКО были получены с удовлетворительной точностью лишь на сетке в 140 000 узлов.
Сравнение вычислительных затрат приведено в табл. 4. Видно, что схема РМГ четвертого порядка точности расходует в 3 раза больше времени по сравнению со схемой МКО при условии использования эквивалентного количества степеней свободы (сетка 62 000). Принимая во внимание большую точность результатов РМГ, правильнее было бы сопоставлять общее время расчетных схем, требуемое для получения результатов с одинаковой точностью. Из таблицы видно примерно одинаковое время расчета для новой схемы на сетке в 2500 ячеек и для схемы МКО на сетке в 140 000 ячеек, что ставит под сомнение пользу от применения схемы высокого порядка точности в данном численном эксперименте.
Отметим, однако, что настоящая реализация РМГ схемы далека от оптимальной и здесь еще имеются потенциальные возможности сокращения расчетного времени. В частности, одним из ключевых моментов является использование оптимизированных квадратурных правил нахождения объемных и поверхностных интегралов. Это является темой дальнейших исследований, связанных с оптимизацией применения новой схемы аппроксимации.
Рис. 11. Сопоставление эксперимента и численных результатов, полученных с РМГ К =3 на грубой сетке 5 и МКО
на сетке 1
Т аблица 4
Оценка времени СРи для различных схем
Метод Количество ячеек Количество степеней свободы ЭД Время на 50 МО итераций, с Время СРи на степень свободы Отношение времен
РМГ К =3 6-уровневый 2 502 50 040 1185 4.74-10-4 3.00
РМГ К =2 5-уровневый 6 393 63 930 826 2.58-10-4 1.63
РМГ К = 1 4-уровневый 14 265 57 060 439 1.54-10-4 0.99
МКО 3-уровневый 62 689 62 689 496 1.58-10-4 1.00
МКО 3-уровневый 139 902 139 902 1115 1.59-10-4
Расчет пространственного турбулентного обтекания крыла ЬЛММ (М = 0.82,
Яе = 7.3 • 106, а = 0.6°). Расчет пространственного турбулентного трансзвукового обтекания крыла ЬАМК [21] был выполнен с использованием схемы второго порядка точности РМГ К = 1. Результаты были сопоставлены с экспериментом и расчетом, выполненным на более густой сетке по традиционной схеме МКО [22].
Все расчеты выполнялись на неструктурированных гексаэдральных сетках, сгенерированных кодом НЕХРКЕ88™ компании КЦМЕСА їй! Для схемы РМГ использовалась сетка с количеством ячеек 190 213. Результаты с использованием МКО были получены на более густой сетке, содержащей 625 076 ячеек и обеспечивающей использование эквивалентного количества переменных. В обоих случаях первый слой узлов был расположен на расстоянии у + ~ 1 от тела. Рис. 12 демонстрирует фрагменты более грубой сетки, используемой для анализа новой схемы.
Все расчеты были выполнены при М = 0.82, а = 0.6°, Яе = 7.3 -106, что соответствовало экспериментальным условиям. Использовалась модель турбулентности Спаларта — Алмараса.
Сравнение распределений давления в сечении крыла г = 0.475 для двух схем дискретизаций представлено на рис. 13. Видно, что, несмотря на использование более грубой расчетной сетки, результаты, полученные с РМГ-схемой, обеспечивают лучшее положение скачка уплотнения.
Рис. 12. Фрагменты неструктурированной гексаэдральной сетки вокруг крыла ЬАМЫ
О О О Эксперимент [21], г = 0.475
Рис. 13. Распределение давления в сечении крыла ЬАМЫ при трансзвуковом турбулентном
обтекании:
-----— МКО на сетке 625 076 ячеек;---------РМГ на сетке 190 213
Задача о распространении сферической акустической волны. Решение задач аэроакустики требует применения численных схем высокого порядка точности, обладающих минимальными дисперсионными и диффузионными свойствами. Исследуемая схема обладает такими свойствами и рассматривается для последующей промышленной реализации. В частности, свойства схемы были изучены в классической задаче о распространении акустического импульса, имеющего сравнительно малое начальное возмущение.
Численная задача рассматривается в кубической расчетной области [0 < х < 1]; [0 < у < 1];
[0 < г < 1]. Имеются следующие начальные условия при t = 0:
и (х, у, г, 0) = 0, V(х, у, г, 0) = 0, w(х, у, г,0) = 0,
Р = Р
Я2 Я2
г 2 1 + е 2 Г0 , Р Рю 1 + -у 2 ^ с2
^ '0
г0 = 0.02, е = 0.001,
где Я2 =(х - 0.5)2 +(у - 0.5)2 +(г - 0.5)2.
Решение линеаризованной задачи в этом случае описывается выражением:
г Г (г-ct )2 (г+^)2 ^ Г (г-^ )2 (г+С)2 ^ \
Рточн рю 1 + £ 1 2 2 (N1 + г02 2г 2 (N1 + г02
V V V ) )
где с — скорость звука. Решение сеточных уравнений получалось с использованием метода Рунге — Кутта пятого порядка точности с глобальным шагом по времени. Представленные результаты соответствуют времени 0.4/с, т. е. когда передний фронт возмущения еще не дошел до
границы расчетной области.
На рис. 14, а, б приводится сравнение результатов для РМГ К = 1, 2, 3 схем и традиционной схемы МКО [22]. Здесь же изображено точное решение этой задачи. Случай, когда исследуемые схемы использовались на индивидуальных сетках, подобранных таким образом, чтобы задейст-
Рис. 14. Сравнение точного и численного решений для акустического импульса: а — эквивалентное число степеней свободы; б — эквивалентная точность решений
вованное количество степеней свободы было эквивалентным, представлен на рис. 14, а. Здесь видно, что схемы с большим порядком точности обеспечивают лучшее согласование с точным решением. Однако эти схемы потребляют и большие компьютерные ресурсы. Поэтому сравнение потребных ресурсов осуществлено для случая, когда обеспечивается одинаковая точность расчетов. На рис. 14, б отражены результаты для РМГ К=2 и 3, полученные на индивидуальных сетках. Видно, что точность двух расчетов совпадает. Сравнение компьютерных ресурсов для этого случая приведено в табл. 5. Схема четвертого порядка точности обеспечила хороший результат в 5.6 раза быстрее, использовав при этом в 1.4 раза меньшую память. Отметим, что получить такую же точность результатов со схемой РМГ К =1 и тем более с использованием традиционной конечно-объемной схемы не представилось возможным из-за чрезмерно большой потребной компьютерной памяти.
Т аблица 5
Сравнение компьютерных ресурсов в случае одинаковой точности расчетов для схем РМГ К=2 и 3
Схема Сетка Количество степеней свободы Используемая память, Kb Соотношение памяти Время CPU Соотношение времени
РМГ K = 2 60 x 60 x 60 2 160 000 743 1.4 13ч 48мин 5.6
РМГ K =3 32 x 32 x 32 655 360 528 1.0 2ч 28мин 1.0
Заключение. Метод Галеркина с разрывными базисными функциями применен к решению уравнений Эйлера и Навье — Стокса на основе использования неструктурированных гексаэд-ральных сеток. Рассмотренные тестовые случаи демонстрируют хорошее согласование с теоретическими или экспериментальными результатами при использовании схемы четвертого порядка точности даже на экстремально грубых сетках.
Выполнены сравнения компьютерных затрат новой схемы и схемы, использующей метод конечного объема второго порядка точности. В случае использования эквивалентного числа степеней свободы в дискретной задаче схема высокого порядка точности требует больших ресурсов.
Исследования показали, что при решении некоторых задач в условиях решения с одинаковой точностью результатов время расчета по схемам РМГ и МКО может быть близким. Такой вывод делает пока неочевидным заключение о необходимости использования схем высокого порядка точности. Только в акустической задаче о распространении сферической волны преимущества от использования схем высокого порядка были очевидными.
Отметим, однако, что практическая реализация исследуемой схемы имеет множество возможностей существенного сокращения расчетного времени. Именно решение этих задач об оптимальной реализации разрывного метода Галеркина откроет перспективы его широкого использования. Вопросы оптимальной реализации и робастности схемы являются предметом будущих исследований.
Автор выражает признательность президенту компании NUMECA проф. Ч. Хиршу (Ch. Hirsch) за предоставление промышленной программы метода конечного объема и помощь в проведении настоящих исследований, а также С. В. Ляпунову за полезные дискуссии.
Работа частично выполнена за счет средств гранта РФФИ №09-01-00243а.
ЛИТЕРАТУРА
1. Волков А. В., Ляпунов С. В. Исследование эффективности использования численных схем высокого порядка точности для решения уравнений Навье — Стокса и Рейнольдса на неструктурированных адаптивных сетках // ЖВМиМФ. 2006. Т. 46, № 10.
2. Волков А. В., Ляпунов С. В. Применение конечно-элементного метода Галеркина с разрывными базисными функциями к решению уравнений Рейнольдса на неструктурированных адаптивных сетках // Ученые записки ЦАГИ. 2007. Т. XXXVIII, № 3—4.
3. Wolkov A., Hirsch Ch., Leonard B. Discontinuous Galerkin method on unstructured hexahedral grids for 3D Euler and Navier-Stokes equations / 18th AIAA Computational Fluid Dynamics Conference // AIAA Paper 2007-4078. 2007.
4. Morrison J., Hemsch M. Statistical analysis of CFD solutions from the 3rd AIAA Drag Prediction Workshop 3rd AIAA APA Drag Prediction Workshop, 2006. http://aaac.larc.nasa.gov/tsab/cfdlarc/aiaa-dpw/Workshop3/presentations/index.html.
5. Ladeinde F., A l ab i K., Safta C., C a i X. and Johnson F. The first high-order CFD simulation of aircraft: challenges and opportunities / 44th AIAA Aerospace Sciences Meeting and Exhibit / AIAA Paper 2006-1526, 2006.
6. Bassi F. and Rebay S. High-order accurate discontinuous finite element solution of the 2D Euler Equations // J. Comput. Phys., 1997. V. 138.
7. Cockburn B., Karniadakis G. and Shu C.-W. Discontinuous Galerkin methods: theory, computation and applications / Lecture Notes in Computational Science. — New York: Springer-Verlag, 1999.
8. Bassi F.Crivellini A. Rebay S.Savini M. Discontinuous Galerkin solutions of the reynolds-averaged Navier — Stokes and k-w turbulence model equations // Comput. Fluids. 2005, V. 34 (4—5).
9. Bassi F., Crivellini A., Di Pietro D. A., Rebay S. A high-order discontinuous Galerkin solver for 3D aerodynamic turbulent flows // Proceedings ECCOMAS CFD 2006 Conference, P. Wesseling, E. Onate and J. Periaux (Eds), Sept. 2006.
10. Nguyen N. C., Persson P-O., P e r a i r e J. RANS solutions using high-order discontinuous Galerkin methods // AIAA Paper 2007-0914. 2007.
11.Luo H.Baum J. D.Lohner R. A fast, p-multigrid discontinuous Galerkin method for compressible flows at all speeds // 44th AIAA Aerospace Sciences Meeting and Exhibit, 9—12, January 2006. Reno, Nevada; AIAA Paper 2006-110. 2006.
12. Wang Z. J. High-order methods for the Euler and Navier — Stokes equations on unstructured grids // Progress in Aerospace Sciences. 2007, V. 43.
13. Delanaye M., Patel A., Leonard B., Hirsch Ch. Automatic unstructured hexahedral grid generation and flow solution // Proc. ECCOMAS CFD-2001 Conference. Swansea, Wales, UK, Sept. 2001.
14. Atkins H. L., Shu C. W. Quadrature free implementation of discontinuous Galerkin method for hyperbolic equations // AIAA J. 1998. V. 36, N 5.
15. Krivodonova L., Berger M. High-order implementation of solid wall boundary conditions in curved geometries // J. of Comput. Phys. 2006. V. 211, N 2.
16. Helenbrook B.,Mavriplis D. J., Atkins H. L. Analysis of «p»-multigrid for continuous and discontinuous finite element discretizations / Proceedings of the 16th AIAA Computational Fluid Dynamics Conference // AIAA 2003-3989. 2003.
17. Nastase C. R., Mavriplis D. J. Discontinuous Galerkin methods using an hp-multigrid solver for inviscid compressible flows on three-dimensional unstructured meshes // 44th AIAA Aerospace Sciences Meeting and Exhibit, 9 —12, January 2006, Reno, Nevada; AIAA Paper 2006-107. 2006.
18. Luo H.Baum J. D.Lohner R. A p-multigrid discontinuous Galerkin method for the Euler equations on unstructured grids // J. Comput. Phys., 2006. V. 211.
19. Mavriplis D. J. and Venkatakrishnan V. Agglomeration multigrid for two dimensional viscous flows // Computers and Fluids. 1995. V. 24, N. 5.
20. E n a y e t M., Gibson M., Taylor A., Yianneskis M. Laser Doppler measurements of laminar and turbulent flow in a pipe bend // NASA Contract Report CR-3551, 1982.
21. Muller 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.
22. www.numeca.be.
PyKonucb nocmynuna 16/12009 г.
