1. На главную
  2. Электронные версии
  3. Журнал вычислительной математики и математической физики
  4. T. 61, № 1, 2021

Журнал вычислительной математики и математической физики, 2021, T. 61, № 1, стр. 3-19

EBR схемы с криволинейными реконструкциями переменных вблизи обтекаемых тел

А. П. Дубень 1, Т. К. Козубская 1, П. В. Родионов 1, *, В. О. Цветкова 1

1 ИПМ им. М.В. Келдыша РАН
125047 Москва, Миусская пл., 4, Россия

* E-mail: rodionov.cs@gmail.com

Поступила в редакцию 13.05.2020
После доработки 05.06.2020
Принята к публикации 18.09.2020

Полный текст (PDF)

Аннотация

Работа посвящена развитию вершинно-центрированных EBR схем повышенной точности для расчета задач газовой динамики на неструктурированных сетках. Предлагается оснастить данные схемы возможностью использования квазиодномерного криволинейного шаблона для реконструкции переменных в областях структурированной или полуструктурированной анизотропной сетки вблизи обтекаемого тела. В двумерном случае использование криволинейных реконструкций приводит к естественной трансформации EBR схемы в структурированный конечно-объемный метод на структурированной сетке. В трехмерном случае для реализации криволинейных реконструкций переменных в призматических слоях разработан оригинальный алгоритм поиска точек шаблона и соответствующих метрических коэффициентов. Эффект от использования криволинейных реконструкций в EBR схемах при решении задач внешнего обтекания демонстрируется на известной тестовой задаче о течении вокруг аэродинамического профиля NACA0012, рассмотренной в двумерной и трехмерной постановках. Валидация нового алгоритма проводится путем сравнения с известными экспериментальными данными, а также результатами расчетов других авторов. Возможность использования криволинейных реконструкций в EBR схеме приводит к улучшению устойчивости метода и повышению точности численных результатов. Библ. 14. Фиг. 14. Табл. 4.

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

ВВЕДЕНИЕ

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

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

В работе рассматриваются алгоритмы, основанные на использовании вершинно-центрированных EBR (Edge-Based Reconstruction) схем для неструктурированных сеток [1]. Повышенная точность данных схем обеспечивается за счет реконструкций потоковых переменных на расширенных реберно-ориентированных шаблонах, а их экономность – за счет квазиодномерной природы такого подхода. EBR схемы, изначально разработанные для произвольных тетраэдральных сеток, допускают обобщение на гибридные неструктурированные сетки [2]. Однако на криволинейной структурированной сетке при большой степени анизотропии сеточных элементов использование EBR схем в их оригинальной формулировке, подразумевающей прямолинейную реконструкцию (т.е. реконструкцию на прямолинейных шаблонах), может приводить как к снижению точности расчета, так и к развитию неустойчивости. В настоящей работе исследуются причины таких негативных явлений. Для преодоления возникающих трудностей предлагается использование криволинейных шаблонов для реконструкции переменных (иначе говоря, криволинейных реконструкций), учитывающих структуру сетки в пристеночной области. Ранее аналогичный подход был реализован для метода коррекции потока [3], [4].

Целью настоящей работы является реализация техники криволинейных реконструкций для их использования в EBR схемах в двумерном случае, а главное, разработке нового эффективного алгоритма построения криволинейных реконструкций для полуструктурированных сеток в трехмерном случае. Эффект от использования криволинейных реконструкций демонстрируется на численном решении одной из классических валидационных задач об обтекании профиля NACA0012 [5], рассмотренной в двумерной и трехмерной постановках.

1. МАТЕМАТИЧЕСКАЯ МОДЕЛЬ

Рассмотрим систему уравнений Навье–Стокса для сжимаемого газа

(1)
$\frac{{\partial {\mathbf{Q}}}}{{\partial t}} + \nabla \cdot F({\mathbf{Q}}) = \nabla \cdot {{F}_{{v}}}({\mathbf{Q}},\nabla {\mathbf{Q}}),$
записанную относительно консервативных переменных, где
${\mathbf{Q}} = \left( \begin{gathered} \rho \\ \rho u \\ E \\ \end{gathered} \right),\quad F = \left( \begin{gathered} \rho u \\ \rho uu + pI \\ (E + p)u \\ \end{gathered} \right),\quad {{F}_{{v}}} = \left( \begin{gathered} 0 \\ \sigma \\ \sigma u - q \\ \end{gathered} \right).$
Здесь $\rho $ – плотность, $u$ – вектор скорости, $p$ – давление, $E$ – полная энергия, $I$ – единичная матрица. Тензор напряжений и вектор теплового потока определяются как $\sigma = \mu (\nabla u + {{(\nabla u)}^{{\text{т}}}} - $ $ - \;(2{\text{/}}3)(\nabla \cdot u)I)$ и $q = - k\nabla T$ соответственно, где $\mu $ – коэффициент динамической вязкости, $k$ – коэффициент теплопроводности, $T$ – температура. Коэффициент динамической вязкости будем определять согласно закону Сазерленда.

В настоящей работе используется также система уравнений Навье–Стокса (1), осредненных по Рейнольдсу (RANS, Reynolds-Averaged Navier–Stokes). Входящий в систему RANS уравнений тензор рейнольдсовых напряжений замкнем с помощью линейной модели турбулентности Спаларта–Аллмараса (SA) [6], записанной относительно модифицированной турбулентной вязкости. При этом общий вид вязких потоков ${{F}_{{v}}}$ останется неизменным с точностью до турбулентной вязкости, которая добавляется к динамической.

2. ЧИСЛЕННЫЙ МЕТОД НА ОСНОВЕ ОРИГИНАЛЬНОЙ EBR СХЕМЫ

Для численного решения системы уравнений (1) на произвольной вычислительной сетке построим схему с определением переменных в сеточных узлах. Далее такие схемы будем называть вершинно-центрированными. Вокруг каждого узла определим медианные ячейки, для которых, согласно конечно-объемному подходу, сформулируем разностные аналоги законов сохранения. Примеры медианных ячеек для двумерной сетки изображены на фиг. 1. Определив сеточную функцию ${{{\mathbf{Q}}}_{i}}$ как интегральное среднее функции ${\mathbf{Q}}$ по построенной вокруг узла $i$ ячейке и использовав формулу Остроградского–Гаусса, перепишем систему (1) в векторно-матричном виде

${{V}_{i}}\frac{{d{{{\mathbf{Q}}}_{i}}}}{{dt}} + \sum\limits_{j \in {{N}_{1}}(i)} {{{F}_{{ij}}}{{s}_{{ij}}}} = {{F}_{{i,{v}}}},$
где ${{V}_{i}}$ – объем ячейки, соответствующей узлу $i$, ${{F}_{{ij}}}$ – интегральное среднее функции $Fn$ по разделяющей узлы $i$ и $j$ грани ячеек, площадь которой равна ${{s}_{{ij}}}$, $n$ – вектор единичной нормали, ${{N}_{1}}(i)$ – множество соседей первого порядка узла $i$, ${{F}_{{i,{v}}}}$ – интеграл функции вязкого потока ${{F}_{{v}}}$ по ячейке, отвечающей узлу $i$. Для вычисления конвективных потоков ${{F}_{{ij}}}$ будем использовать метод Роу приближенного решения задачи о распаде разрыва:
${{F}_{{ij}}} = \frac{1}{2}(F_{{ij}}^{R} + F_{{ij}}^{L}) - \frac{1}{2}\left| {{{A}_{{ij}}}} \right|({\mathbf{Q}}_{{ij}}^{R} - {\mathbf{Q}}_{{ij}}^{L}).$
Значения ${\mathbf{Q}}_{{ij}}^{{L/R}}$ слева и справа от интерфейса определим при помощи квазиодномерных реконструкций $R_{{ij}}^{{L/R}}\{ {\mathbf{Q}}\} $, определенных на шаблонах, точки которых принадлежат прямой, содержащей ребро $ij$. Значения потоков $F_{{ij}}^{{L/R}}$ будем полагать равными $F(R_{{ij}}^{{L/R}}\{ {\mathbf{Q}}\} ){{n}_{{ij}}}$ либо $R_{{ij}}^{{L/R}}\{ F{{n}_{{ij}}}\} $ в зависимости от выбранного типа реконструкций [7]. Здесь ${{n}_{{ij}}}$ – интегральное среднее вектора n по общей грани между ячейками, соответствующими узлам $i$ и $j$, $\left| {{{A}_{{ij}}}} \right| = \left| {\frac{{dF{\mathbf{n}}}}{{d{\mathbf{Q}}}}({{{\mathbf{Q}}}_{{ij}}})} \right| = {{S}_{{ij}}}\left| {{{\Lambda }_{{ij}}}} \right|S_{{ij}}^{{ - 1}}$, ${{{\mathbf{Q}}}_{{ij}}}$ – среднее по Роу, вычисленное по значениям ${\mathbf{Q}}_{{ij}}^{{L/R}}$, ${{S}_{{ij}}}$ и ${{\Lambda }_{{ij}}}$ – соответствующие ${{{\mathbf{Q}}}_{{ij}}}$ матрицы собственных векторов и собственных значений оператора $\frac{{dF{\mathbf{n}}}}{{d{\mathbf{Q}}}}$. Построенную таким образом схему, использующую реберно-ориентированные реконструкции переменных, можно в широком смысле отнести к классу EBR схем. В более узком смысле, согласно работе [1], в EBR схемах используются такие реконструкции, которые на трансляционно-инвариантных (переходящих в себя при переносе на вектор любого сеточного ребра) сетках приводят к трансформации данного метода в конечно-разностную схему высокого порядка. При этом схема данного семейства называется EBR n схемой, если в линейном случае ее порядок на трансляционно-инвариантных сетках равен n.

Фиг. 1.

Шаблон схемы EBR5 для ребра $ij$ на неструктурированной треугольной сетке.

Опишем метод построения квазиодномерных реконструкций, используемых в оригинальной формулировке EBR схем [1], на примере схемы EBR5 в двумерной постановке (фиг. 1). Рассмотрим ребро $ij$, в середине которого необходимо реконструировать значение функции ${\mathbf{Q}}$, и для каждого из его узлов построим множества топологических соседей первого и второго порядка. Обозначим точку пересечения луча $ji$ с множеством граней, все узлы которых являются соседями второго порядка узла $i$, индексом –2, а точку пересечения данного луча с множеством граней, все узлы которых являются соседями первого порядка узла $i$, индексом –1. В случае неединственности первой точки, индексом –2 обозначим точку, наиболее удаленную от узла $i$. Аналогично для луча $ij$ и узла $j$ получим точки с индексами 3 и 2 соответственно. Значения функции ${\mathbf{Q}}$ в точках {–2, –1, 2, 3} определим при помощи линейной интерполяции по соответствующим пересекаемым лучом граням. Если присвоить узлу $i$ индекс 0, а узлу $j$ – индекс 1, то операторы реконструкции функции ${\mathbf{Q}}$ в терминах разделенных разностей

$\Delta _{m}^{L}\{ {\mathbf{Q}}\} = \frac{{{{{\mathbf{Q}}}_{{m + 1}}} - {{{\mathbf{Q}}}_{m}}}}{{\left| {{{{\mathbf{r}}}_{{m + 1}}} - {{{\mathbf{r}}}_{m}}} \right|}},\quad \Delta _{m}^{R}\{ {\mathbf{Q}}\} = \Delta _{{ - m}}^{L}\{ {\mathbf{Q}}\} $
могут быть записаны как
(2)
$\begin{gathered} R_{{ij}}^{L}\{ {\mathbf{Q}}\} = {{{\mathbf{Q}}}_{i}} + \frac{{\left| {{{{\mathbf{r}}}_{i}} - {{{\mathbf{r}}}_{j}}} \right|}}{2}\sum\limits_m {{{\beta }_{m}}} \Delta _{m}^{L}\{ {\mathbf{Q}}\} , \\ R_{{ij}}^{R}\{ {\mathbf{Q}}\} = {{{\mathbf{Q}}}_{j}} - \frac{{\left| {{{{\mathbf{r}}}_{i}} - {{{\mathbf{r}}}_{j}}} \right|}}{2}\sum\limits_m {{{\beta }_{m}}} \Delta _{m}^{R}\{ {\mathbf{Q}}\} , \\ \end{gathered} $
где ${{\beta }_{{ - 2}}}$ = –1/15, ${{\beta }_{{ - 1}}}$ = 11/30, ${{\beta }_{0}}$ = 4/5, ${{\beta }_{1}}$ = –1/10 [1]. В схеме EBR3, для определения которой нужен более короткий шаблон, данные коэффициенты принимают значения ${{\beta }_{{ - 2}}}$ = 0, ${{\beta }_{{ - 1}}}$ = 1/3, ${{\beta }_{0}}$ = 2/3, ${{\beta }_{1}}$ = 0.

Как уже отмечено выше, для линеаризованных уравнений на трансляционно-инвариантных сетках теоретический порядок точности схем EBR5 и EBR3 равен пятому и третьему соответственно. В произвольном случае численный порядок точности схем EBR3 и EBR5 в зависимости от качества используемой неструктурированной сетки может варьироваться от 5/4 [8] до 3 [1].

Используемый в работе численный метод решения системы уравнений (1) имеет в основе своей схему EBR5 или EBR3 для аппроксимации конвективных потоков. Аппроксимация же вязких потоков проводится при помощи метода Галeркина с кусочно-линейными базисными функциями (с диагонализированной матрицей масс). Для интегрирования по времени применяется неявная схема первого порядка с линеаризацией системы сеточных уравнений по Ньютону. Решение системы линейных алгебраических уравнений в рамках одной итерации по методу Ньютона осуществляется с помощью метода бисопряженных градиентов с ILU0 предобуславливателем.

3. EBR СХЕМА С КРИВОЛИНЕЙНЫМИ РЕКОНСТРУКЦИЯМИ В СТРУКТУРИРОВАННЫХ ОБЛАСТЯХ ДВУМЕРНОЙ СЕТКИ ВБЛИЗИ ОБТЕКАЕМОГО ТЕЛА

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

Фиг. 2.

Шаблоны схем EBR5 и EBR5 IJK (EBR5 SS) для ребра $ij$ на анизотропной структурированной сетке вблизи обтекаемого тела.

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

Сначала рассмотрим криволинейные реконструкции в EBR схемах и методы построения соответствующих им криволинейных шаблонов для двумерных задач аэродинамического обтекания, где в областях пограничных слоев, как правило, используются слои структурированных сеток, состоящих из трапециевидных элементов (фиг. 2). Учитывая структурированность данных сеток, в качестве точек шаблона криволинейной реконструкции будем выбирать не точки пересечения прямой $ij$ с соответствующими изоповерхностями порядка соседства (точки {–2, –1, 2, 3} на фиг. 2), а то или иное количество узлов справа и слева от узлов $i$ и $j$, являющихся их структурными соседями (узлы {–2 S, –1 S, 2 S, 3 S} на фиг. 2). Видно, что данный выбор обеспечивает примерно равный шаг между точками шаблона реконструкции, а также малую вариацию значений реконструируемой функции в пристеночной области. Выбор такого криволинейного шаблона вместо прямолинейного повышает устойчивость итоговой схемы и снижает ограничения на допустимую геометрию сетки при использовании тех же формул (2) для вычисления коэффициентов реконструкций. Стоит также отметить, что областях с резкими изменениями направлений сеточных линий, например вблизи острого угла геометрии, криволинейные реконструкции могут терять указанные свойства, и даже наоборот приводить к повышению ошибки и неустойчивости счета. Для устранения данного недостатка достаточно задать ограничение на максимальный угол между прямой $ij$ и направлениями криволинейного шаблона реконструкции, при нарушении которого будет происходить переключение на оригинальную схему EBR.

Алгоритм перехода от криволинейных 5-точечных шаблонов реконструкций схемы EBR5 вблизи обтекаемого тела к стандартным прямолинейным реконструкциям при удалении от него удобно организовать поэтапно на основе анализа локальной структурированности сетки. Так, при отсутствии структурных соседей второго порядка, но при наличии соответствующих соседей первого порядка можно использовать укороченный криволинейный шаблон реконструкции, соответствующий схеме EBR3, а переходить к прямолинейным шаблонам реконструкции – только в случае отсутствия структурных соседей первого порядка. Отметим, что такой подход к построению шаблонов в структурированной подобласти сетки применим вдоль сеточных линий как в тангенциальном, так и нормальном направлении относительно обтекаемого тела. Описанный метод выбора шаблонов реконструкций уже использовался в работе [9]. В дальнейшем при описании численных результатов двумерные EBR схемы, явно использующие ijk-топологию при построении криволинейных шаблонов вблизи тела, будем обозначать “EBR IJK”.

Сформулируем несколько иной вариант алгоритма построения шаблонов криволинейных реконструкций, преимущество которого заключается в простоте обобщения на трехмерный случай. Согласно этому алгоритму каждому узлу приграничной сетки сопоставляется уровень удаленности от поверхности обтекаемого тела, определяемый как минимальное число сеточных ребер, которыми этот узел может быть соединен с поверхностью. Для узлов на поверхности тела уровень удаленности полагается равным нулю. Если ребро ij лежит на изолинии уровня удаленности, то шаблон реконструкции будем составлять из узлов требуемого порядка соседства, имеющих тот же уровень удаленности. Формулы (2) для вычисления коэффициентов реконструкций, как и ранее, оставим неизменными. Нетрудно видеть, что построенные таким образом шаблоны реконструкции вдоль сеточных линий тангенциального направления будут совпадать с шаблонами реконструкции, построенными с использованием автоматического поиска структурированности сетки. Реконструкции же для ребер, расположенных в нормальном направлении по отношению к телу и вершины которых имеют разный уровень, будем проводить согласно оригинальной EBR схеме, т.е. при помощи прямолинейных шаблонов. Заметим, что в силу топологии приграничной структурированной сетки использование прямолинейных или криволинейных реконструкций в нормальном направлении не приводит к сколько-нибудь существенной разнице в результатах.

4. EBR СХЕМА С КРИВОЛИНЕЙНЫМИ РЕКОНСТРУКЦИЯМИ В ПРИЗМАТИЧЕСКИХ СЛОЯХ СЕТКИ ВБЛИЗИ ОБТЕКАЕМОГО ТЕЛА

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

Фиг. 3.

Шаблоны схемы EBR5 SS для ребра $ij$ на призматической сетке вблизи обтекаемого тела.

Для определения шаблонов криволинейных реконструкций в тангенциальном направлении построим следующий алгоритм. На стадии препроцессинга по аналогии с двумерным случаем определим уровни удаленности сеточных узлов от поверхности обтекаемого тела. Как и в двумерном случае, численный поток между узлами, имеющими разные уровни удаленности, будем определять согласно оригинальной EBR схеме, использующей прямолинейные реконструкции. Реконструкцию переменных на ребре ij, соединяющем соседние узлы i и j с одним уровнем удаленности, будем проводить по следующему алгоритму, иллюстрация к которому представлена на фиг. 4.

Фиг. 4.

Алгоритм нахождения точек криволинейного шаблона реконструкции переменных для схемы EBR5 SS.

1. Для каждого из узлов ребра $ij$ построим множества узлов, являющихся их соседями 1-го и 2-го порядка.

2. Из этих множеств исключим узлы, уровень удаленности которых отличается от уровня удаленности узлов i и j.

3. Каждому из построенных множеств узлов сопоставим множество сеточных ребер, оба узла которых принадлежат соответствующему множеству узлов.

4. Зададим проекционную плоскость, проходящую через ребро $ij$, посредством вектора $\vec {P}$, являющегося полусуммой нормалей к граням, инцидентным ребру $ij$ и принадлежащим изоповерхности уровня удаленности узлов i и j. Если инцидентная грань только одна, вектор $\vec {P}$ положим равным нормали к ней, а если такие грани отсутствуют, реконструкцию вдоль ребра $ij$ будем проводить согласно оригинальной схеме EBR с прямолинейными реконструкциями, пропуская последующие шаги данного алгоритма.

5. Спроектируем полученные в п. 3 множества ребер на проекционную плоскость, задаваемую вектором $\vec {P}$ и содержащую ребро $ij$, и построим на ней двумерный шаблон прямолинейной реконструкции в полном соответствии с оригинальной EBR схемой в двумерной постановке.

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

При использовании в качестве базовой схемы EBR5 в полуструктурированных областях сетки переключение между криволинейными и прямолинейными реконструкциями в неструктурированных зонах также предлагается выполнять поэтапно по аналогии с двумерным случаем. В дальнейшем при описании численных результатов трехмерные EBR схемы с криволинейными реконструкциями в полуструктурированных областях, построенными согласно вышеописанному алгоритму, будем обозначать “EBR SS” (Semi-Structured).

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

5. ЧИСЛЕННЫЕ РЕЗУЛЬТАТЫ

5.1. Физическая постановка задачи

Для оценки модификаций схем EBR, описанных в предыдущем разделе, были проведены серии расчетов задачи об обтекании аэродинамического профиля NACA0012 [5]. Постановка задачи формулируется следующим образом. Профиль NACA0012 с хордой единичной длины помещается в однородный поток воздуха с числом Маха ${{M}_{\infty }} = 0.15$ и температурой ${{T}_{\infty }} = 300$ K. Число Рейнольдса относительно длины хорды $c$ составляет ${{\operatorname{Re} }_{c}} = 6 \times {{10}^{6}}$. В настоящей работе рассмотрены углы атаки 0°, 10° и 15°.

5.2. Вычислительная постановка задачи

При проведении расчетов используется фоновый поток, характеризующийся высоким уровнем турбулентности. Это обеспечивается за счет задания входного условия ${{\nu }_{t}}{\text{/}}\nu = 1$, где ${{\nu }_{t}}$ – коэффициент турбулентной вязкости, а $\nu $ – коэффициент молекулярной вязкости.

В двумерной постановке расчетная область представляет собой квадрат $ - 500 \leqslant x{\text{/}}c$, $y{\text{/}}c \leqslant 500$, в центре которого, в точке $(0,0)$, находится передняя кромка профиля. В трехмерной постановке расчетной областью является прямоугольный параллелепипед $ - 500 \leqslant x{\text{/}}c$, $y{\text{/}}c \leqslant 500$, $0 \leqslant z{\text{/}}c \leqslant 0.25$, передняя кромка профиля в данном случае принадлежит прямой $(0,0,z)$. На поверхности профиля задаются условия прилипания $u = 0$, ${{\nu }_{t}} = 0$ и адиабатичности. На свободных границах $x{\text{/}}c = \pm 500$ и $y{\text{/}}c = \pm 500$ в случае входа держатся все параметры потока, кроме давления, которое экстраполируется, в случае выхода, наоборот, держится только давление, а основные параметры экстраполируются из внутренней области. В трехмерной постановке на границах $z{\text{/}}c = 0$ и $z{\text{/}}c = 0.25$ задаются условия периодичности.

Для проведения двумерных расчетов использовалась последовательность гибридных сеток, являющихся структурированными трапециевидными вблизи профиля и неструктурированными треугольными в остальной области (фиг. 5). Характерные параметры указанных сеток приведены в табл. 1, где $N$ соответствует общему числу узлов сетки, а ${{N}_{{{\text{surf}}}}}$ – числу узлов, принадлежащих границе обтекаемого профиля. Для проведения трехмерных расчетов использовалась аналогичная последовательность гибридных сеток, являющихся полуструктурированными призматическими вблизи профиля и неструктурированными тетраэдральными в остальной области. Сетки данной последовательности совпадали с соответствующими указанными ранее двумерными сетками на границах $z{\text{/}}c = 0$ и $z{\text{/}}c = 0.25$. Характерные параметры последовательности трехмерных сеток сведены в табл. 2, где при помощи ${{N}_{{{\text{surf}},\;z = 0}}}$ обозначено количество точек на поверхности профиля в плоскости $z = 0$. Отметим, что названия “x1”, …, “x8” введены исключительно для удобства дальнейших обозначений и не являются следствием последовательного разбиения.

Фиг. 5.

Конфигурация расчетных сеток.

Таблица 1.  

Параметры двумерных сеток

2D сетка x1 x2 x3 x4 x8 897 × 257
$N$ 55 тыс. 51 тыс. 69 тыс. 83 тыс. 84 тыс. 231 тыс.
${{N}_{{{\text{surf}}}}}$    102    162    246    442    930    513
Таблица 2.  

Параметры трехмерных сеток

3D сетка x1 x2 x3 x4 x8
$N$ 515 тыс. 801 тыс. 1.4 млн. 2.9 млн. 9.7 млн.
${{N}_{{{\text{surf}}}}}$ 3 тыс. 7 тыс. 15 тыс. 40 тыс. 176 тыс.
${{N}_{{{\text{surf}},\;z = 0}}}$    102    162      246      442      930

Помимо указанных сеток, для валидационных расчетов в настоящей работе используется двумерная структурированная сетка 897 × 257, применяющаяся в [5] для получения эталонных численных результатов. Стоит отметить, что геометрия профиля в построенных сетках и сетке 897 × 257 отличались в пределах 1%.

Распределения значений безразмерной высоты первой пристеночной ячейки ${{y}^{ + }}$, соответствующих обтеканию профиля под углами атаки 0°, 10° и 15°, для двумерных сеток изображены на фиг. 6. Аналогичные значения ${{y}^{ + }}$ для трехмерных сеток практически не отличаются от приведенных.

Фиг. 6.

Безразмерная высота первой пристеночной ячейки y+ для двумерных сеток при обтекании профиля под углами атаки 0° (а), 10° (б) и 15° (в).

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

5.3. Валидация EBR схем с криволинейными реконструкциями

Для тестирования схем EBR IJK и EBR SS использовались самые подробные двумерные и трехмерные сетки x8, а также эталонная сетка 897 × 257. Валидация данных схем и их программной реализации проводилась путем сопоставления значений безразмерных коэффициентов подъемной силы (Cl) и сопротивления формы (Cd) для различных углов атаки. Сравнение полученных результатов с данными из [5] представлено в табл. 3.

Таблица 3.  

Результаты валидации схем EBR5 IJK и EBR5 SS по коэффициентам подъемной силы (Cl) и сопротивления формы (Cd)

Схема (сетка) 0°: Cl 10°: Cl 15°: Cl 0°: Cd 10°: Cd 15°: Cd
EBR5 (3D, x8) ∼0 1.0862   0.00810 0.01234  
EBR5 SS (3D, x8) ∼0 1.0865   0.00810 0.01233  
EBR5 (2D, x8) ∼0 1.0875 1.5339 0.00811 0.01239 0.02179
EBR5 IJK (2D, x8) ∼0 1.0871 1.5345 0.00812 0.01237 0.02166
EBR5 (897 × 257) ∼0 1.0946 1.5437 0.00810 0.01264 0.02219
EBR5 IJK (897 × 257) ∼0 1.0940 1.5436 0.00810 0.01259 0.02203
CFL3D (897 × 257) ∼0 1.0909 1.5461 0.00819 0.01231 0.02124
FUN3D (897 × 257) ∼0 1.0983 1.5547 0.00812 0.01242 0.02159
NTS (897 × 257) ∼0 1.0891 1.5461 0.00813 0.01243 0.02105
JOE (897 × 257) ∼0 1.0918 1.5490 0.00812 0.01245 0.02148
SUMB (897 × 257) ∼0 1.0904 1.5446 0.00813 0.01233 0.02141
TURNS (897 × 257) ∼0 1.1000 1.5642 0.00830 0.01230 0.02140
GGNS (897 × 257) ∼0 1.0941 1.5576 0.00817 0.01225 0.02073

Видно, что в численных данных, полученных известными программными комплексами, значения Cl отличаются друг от друга на 1%, а значения Cd – на 4%. При добавлении к ним данных тестирования схем EBR5 и EBR5 IJK на сетке 897 × 257 разница по Cl остается в пределах 1%, а отклонение по Cd увеличивается до 6%. При добавлении к исходным данным результатов тестирования схем EBR5, EBR5 IJK и EBR5 SS на двумерных и трехмерных сетках ×8 максимальное отличие по Cl увеличивается до 2%, а отклонение Cd составляет 5%. Отметим также, что в рамках одной и той же сетки при использовании различных типов реконструкций разница по Cl сохраняется в пределах 0.1%, по Cd – в пределах 1%.

Сравнение численно полученных распределений коэффициентов давления (Cp) и трения (Cf) для различных углов атаки с данными [5] приведено на фиг. 7. Видно, что результаты расчетов соответствуют экспериментальным данным, а сами результаты расчетов по коэффициентам давления и трения практически совпадают друг с другом.

Фиг. 7.

Результаты валидации схем EBR5 IJK и EBR5 SS по коэффициентам давления (Cp) и трения (Cf) для углов атаки 0° (а, б), 10° (в, г) и 15° (д, е).

В итоге, несмотря на некоторое увеличение отклонений по интегральным характеристикам от соответствующих референсных значений и с учетом небольших различий в геометрии профиля между используемыми в расчетах сетками и эталонной сеткой 897 × 257, можно заключить, что полученные с помощью схем EBR IJK и EBR SS численные результаты подтверждают применимость данных схем для решения стационарных задач аэродинамического обтекания, а также корректность их программной реализации.

5.4. Сравнительный анализ результатов, полученных EBR схемами с прямолинейными и криволинейными реконструкциями

Рассмотрим теперь результаты расчетов, проведенных с использованием схем EBR c прямолинейными и криволинейными реконструкциями на последовательностях двумерных и трехмерных гибридных неструктурированных сеток x, x2, x3, x4 и x8. Как и ранее, на двумерных сетках будем использовать схему с криволинейными реконструкциями EBR IJK, а на трехмерных – схему с криволинейными реконструкциями EBR SS. Отметим, что применение прямолинейных реконструкций при проведении расчетов на достаточно грубых сетках x1–x3, как правило, приводит к неустойчивости, возможные причины которой уже рассматривались. Для обеспечения устойчивости алгоритма в таких случаях использовалось ограничение максимального допустимого отношения длин шагов на шаблоне реконструкции и, при его несоблюдении, производилось переключение на схему EBR3 с изменением, при необходимости, коэффициента реконструкции ${{\beta }_{{ - 1}}}$: если в формулах (2) значение $\left| {{{{\mathbf{r}}}_{0}} - {{{\mathbf{r}}}_{{ - 1}}}} \right| \times {{C}_{{{\text{ratioLim}}}}}$ было меньше, чем $\left| {{{{\mathbf{r}}}_{i}} - {{{\mathbf{r}}}_{j}}} \right|$, где ${{C}_{{{\text{ratioLim}}}}}$ – глобально задаваемый ограничитель, значение ${{\beta }_{{ - 1}}}$ умножалось на коэффициент ${{C}_{{{\text{ratioLim}}}}} \times \left| {{{{\mathbf{r}}}_{0}} - {{{\mathbf{r}}}_{{ - 1}}}} \right|{\text{/}}\left| {{{{\mathbf{r}}}_{i}} - {{{\mathbf{r}}}_{j}}} \right|$. Для рассматриваемой задачи и сеток x1–x3 достаточным оказалось значение ${{C}_{{{\text{ratioLim}}}}}$, равное 20.

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

Фиг. 8.

Невязкая (Cd,nonvisc) и вязкая (Cd,visc) компоненты коэффициента сопротивления формы, полученные с помощью различных схем на последовательностях двумерных и трехмерных сеток, для угла атаки 0°.

На грубых сетках x1–x3 схема EBR5 и ее версии с криволинейными реконструкциями ожидаемо демонстрируют более высокую точность в сравнении с соответствующими схемами, основанными на схеме EBR3. При этом использование криволинейных реконструкций приводит к существенному уменьшению ошибки таким образом, что значения исследуемых контрольных величин при использовании криволинейных реконструкций даже на достаточно грубых сетках x1 (как в двумерном, так и в трехмерном случаях) оказываются уже вполне близкими к истинным значениям, полученным на эталонных сетках. Также следует отметить, что при увеличении сеточного разрешения разница между результатами, полученными схемами с прямолинейными и криволинейным реконструкциями, уменьшается, что объясняется постепенным выпрямлением шаблонов криволинейных реконструкций в пограничном слое при сгущении сеток.

Приведенные выше выводы подтверждаются аналогичными результатами для углов атаки 10° и 15° (фиг. 9 и 10). Причем стоит заметить, что при ненулевом угле атаки на грубых сетках схемы EBR3, использующие укороченный криволинейный шаблон, показывает более высокую точность, чем схема EBR5, работающая на более протяженном, но прямолинейном шаблоне.

Фиг. 9.

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

Фиг. 10.

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

На фиг. 11 для иллюстрации работы рассматриваемых EBR схем на различных двумерных сетках приведены некоторые наиболее характерные распределения коэффициентов давления и трения по поверхности аэродинамического профиля. Из приведенных графиков видно, что для всех углов атаки наибольшее отклонение от эталонных численных результатов демонстрирует схема EBR3 на сетке x1, далее следуют схемы EBR5 на той же сетке, EBR3 на сетке x2, EBR5 на сетке x2 и EBR5 IJK на сетке x1. Описанная последовательность подтверждается данными, представленными на фиг. 8, 9 и 10. Еще раз отметим, что коэффициенты давления и трения, полученные схемой EBR5 IJK с криволинейными реконструкциями, даже на самой грубой сетке x1 уже вполне хорошо согласуются с эталонными результатами, полученными на самой подробной сетке x8.

Фиг. 11.

Наиболее характерные распределения коэффициентов давления (Cp) и трения (Cf), полученные в расчетах на последовательности двумерных сеток при углах атаки 0° (а, б), 10° (в, г) и 15° (д, е).

По приведенным на фиг. 11 распределениям коэффициента трения можно оценить вариации размеров зоны рециркуляции. Соответствующие координаты точки отрыва сведены в табл. 4. Из нее следует, что без использования криволинейных реконструкций в EBR схемах на грубых сетках полученный размер рециркуляционной зоны может существенно отличаться от соответствующего физически корректного результата. Особенно это заметно для угла атаки 15°. Также использование прямолинейных реконструкций на грубых сетках может приводить к ложному отрыву, что имеет место при угле атаки 10°.

Таблица 4.  

Координата $x{\text{/}}c$ точки отрыва течения для некоторых расчетов, выполненных на последовательностях двумерных и трехмерных сеток

Схема (сетка) 10° 15°
EBR5 IJK (897 × 257) – 0.91
EBR5 IJK (2D, x8) – 0.90
EBR5 SS (3D, x8) –  
EBR5 IJK (2D, x1) – 0.89
EBR5 2D (2D, x2) 0.994 0.74
EBR3 2D (2D, x2) 0.987 0.66
EBR5 2D (2D, x1) 0.980 0.40
EBR3 2D (2D, x1) 0.935 0.32
EBR5 3D (3D, x1) –  
EBR3 3D (3D, x1) 0.973  

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

Фиг. 12.

Распределения коэффициента трения в наиболее характерных областях, получаемые на подробных сетках при углах атаки 0° и 10°.

Фиг. 13.

Распределения коэффициента трения в наиболее характерной области, получаемые на подробных сетках при угле атаки 15° (слева), и распределения коэффициента давления в наиболее характерных областях, получаемые для различных углов атаки на сетках x8 (справа).

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

Поверхностные распределения коэффициента трения для трехмерных расчетов обтекания с нулевым углом атаки по схемам EBR5 на подробных сетках изображены на фиг. 14. Из них видно, что использование прямолинейных реконструкций приводит к возникновению сеточных осцилляций, в то время как использование криволинейных реконструкций позволяет избежать данного эффекта. Аналогичная ситуация имеет место как для схем EBR3, так и для ненулевых углов атаки.

Фиг. 14.

Поверхностные распределения коэффициента трения, полученные с использованием схем EBR5 и EBR5 SS на трехмерных сетках x8, x4, x3 для угла атаки 0°.

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

Все описанные в настоящей работе расчеты выполнены с использованием программного комплекса NOISEtte [14].

ЗАКЛЮЧЕНИЕ

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

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

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

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

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

Авторы благодарят П.А. Бахвалова за ценные замечания и полезные дискуссии.

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

  1. Abalakin I.V., Bakhvalov P.A., Kozubskaya T.K. Edge-based reconstruction schemes for unstructured tetrahedral meshes // Int. J. Numer. Methods Fluids. 2016. V. 81. № 6. P. 331–356.

  2. Bakhvalov P.A., Kozubskaya T.K. On Efficient Vertex-Centered Schemes on Hybrid Unstructured Meshes // AIAA Paper 2016–2966. 22nd AIAA/CEAS Aeroacoustics Conference. 2016.

  3. Katz A.J., Work D. High-order flux correction/finite difference schemes for strand grids // J. Comput. Phys. 2015. V. 282. P. 360–380.

  4. Tong O., Katz A.J., Wissink A.M., Sitaraman J. High-order methods for three-dimensional strand-cartesian grids // AIAA Paper 2015–0835. 53rd AIAA Aerospace Sciences Meeting. 2015.

  5. NASA Langley Research Center. Turbulence Modeling Resource. 2DN00: 2D NACA 0012 Airfoil Validation Case URL: https://turbmodels.larc.nasa.gov/naca0012_val.html. SA Model Results URL: https://turbmodels.larc.nasa.gov/naca0012_val_sa.html.

  6. Spalart P.R., Allmaras S.R. A One-Equation Turbulence Model for Aerodynamic Flows // AIAA Paper 92–0439. 30th Aerospace Science Meeting. 1992.

  7. Bakhvalov P.A., Kozubskaya T.K. EBR-WENO scheme for solving gas dynamics problems with discontinuities on unstructured meshes // Comput. Fluids. 2017. V. 157. P. 312–324.

  8. Бахвалов П.А. О порядке точности реберно-ориентированных схем на сетках специального вида // Препринты ИПМ им. М.В. Келдыша. 2017. № 79. 32 с.

  9. Duben A.P., Kozubskaya T.K. On Scale-Resolving Simulation of Turbulent Flows Using Higher-Accuracy Quasi-1D Schemes on Unstructured Meshes // Progress in Hybrid RANS-LES Modelling. HRLM 2016. Notes on Numerical Fluid Mechanics and Multidisciplinary Design. V. 137. P. 169–178.

  10. Kallinderis Y., Ward S. Prismatic grid generation for three-dimensional complex geometries // AIAA J. 1993. V. 31. № 10. P. 1850–1856.

  11. Connell S., Braaten M. Semi-structured mesh generation for 3D Navier-Stokes calculations // AIAA Paper 92–0439. 12th Computational Fluid Dynamics Conference. 1995.

  12. Khawaja A., Kallinderis Y. Hybrid grid generation for turbomachinery and aerospace applications // Int. J. Numer. Methods Eng. 2000. V. 49. № 1–2. P. 145–166.

  13. Athanasiadis A.N., Deconinck H. A folding/unfolding algorithm for the construction of semi-structured layers in hybrid grid generation // Comput. Methods Appl. Mech. Eng. 2005. V. 194 № 48–49. P. 5051–5067.

  14. Gorobets A. Parallel Algorithm of the NOISEtte Code for CFD and CAA Simulations // Lobachevskii J. Math. 2018. V. 39. № 4. P. 524–532.

Дополнительные материалы отсутствуют.

Инструменты

следующая статья выпуска содержание выпуска

Журнал вычислительной математики и математической физики

Архивы выпусков Информация о журнале Отправить рукопись в журнал