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

Журнал вычислительной математики и математической физики, 2023, T. 63, № 4, стр. 533-547

Простой критерий оценки сеточной детализации для rans методов

А. А. Гаврилов 1, *, А. А. Дектерев 1, 2, **, А. В. Шебелев 2, ***

1 Институт теплофизики им. С.С. Кутателадзе СО РАН
630090 Новосибирск, пр-т акад. Лаврентьева, 1, Россия

2 Сибирский федеральный университет
660041 Красноярск, пр-т Свободный, 79, Россия

* E-mail: gavand@yandex.ru
** E-mail: dekterev@mail.ru
*** E-mail: aleksandr-shebelev@mail.ru

Поступила в редакцию 20.07.2022
После доработки 10.10.2022
Принята к публикации 15.12.2022

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

Аннотация

Предложен простой критерий оценки детализации сетки для получения сеточно-независимого решения турбулентного течения в рамках статистического подхода моделирования турбулентности. Критерий выведен на основе апостериорной оценки локальной ошибки интерполяции поля кинетической энергии турбулентности. Хорошая детализация сетки должна приводить к малым значениями ошибки интерполяции сеточного поля турбулентной энергии. Использование уравнения переноса для кинетической энергии турбулентности и условия реализуемости позволили свести оценку относительной ошибки интерполяции к явной формуле для оценки максимального шага сетки, необходимого для получения сеточно-независимого решения. В настоящей работе предложенный критерий применен к стационарной задаче о течении за уступом и к нестационарной задаче об обтекании полукругового профиля, установленного под нулевым углом атаки при числе Рейнольдса Re = 45 000. В результате численного исследования показано, что введенный критерий позволяет правильно оценить необходимую детализацию сетки для получения сеточно-независимого решения вдали от стенки. Критерий может быть использован как для оценки полученного решения на сеточную независимость, так и для адаптации расчетной сетки. Библ. 10. Фиг. 12. Табл. 1.

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

1. ВВЕДЕНИЕ

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

В инженерной практике к конечно-разностной или конечно-объемной сетке предъявляются дополнительные требования, часто противоречащие основному критерию качества сетки (инженерные критерии):

1. Степень детализации не должна ухудшать итерационную сходимость решения.

2. Сетка должна описывать характерные особенности исходной геометрии.

3. Сетка должна быть экономной. Объем данных на сетке должен “помещаться” в выделенные на задачу вычислительные ресурсы. Расчет на сетке должен выполняться за приемлемое время.

4. Процесс построения сетки также должен быть “экономным”. Построение сетки должно выполняться за приемлемое время.

Наиболее популярным методом в инженерной практике для моделирования турбулентных течений является подход с использованием осреднения по Рейнольдсу (Reynolds-averaged Navier–Stokes, RANS). Нестационарный вариант подхода (Unsteady Reynolds-averaged Navier–Stokes, URANS) позволяет непосредственно разрешить только крупномасштабные когерентные вихревые структуры, в то время как мелкомасштабные турбулентные пульсации моделируются с помощью выбранной модели турбулентности. Как правило, при решении уравнений Рейнольдса применяются конечно-объемные или конечно-разностные схемы второго порядка аппроксимации по пространству и времени. Для аппроксимации конвективных слагаемых используются противопоточные схемы, которые в областях больших градиентов снижают порядок аппроксимации для получения гладкого решения.

Ошибки дискретизации уменьшаются при увеличении степени детализации сетки. И можно ожидать, что при некоторой степени детализации сетки дальнейшее уменьшение шага сетки не приведет к заметному изменению численного решения. Определить достаточную детализацию сетки и, более того, оценить ошибки дискретизации и порядок аппроксимации численного метода можно с помощью расчетов на последовательности сеток с разной детализацией. При исследовании сеточной сходимости нет необходимости определять индекс сеточной сходимости (grid convergence index) или порядок аппроксимации численного метода, достаточно показать монотонную сеточную сходимость интересующей величины на последовательности сеток (см. [1]). Количество сеток должно быть не меньше трех. Простейший способ оценки сходимости состоит в использовании метода Ричардсона, позволяющего сделать оценку ошибки дискретизации путем сравнения решений, полученных на последовательности сеток с разной детализацией.

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

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

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

2. МАТЕМАТИЧЕСКАЯ ФОРМУЛИРОВКА КРИТЕРИЯ

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

$\rho \frac{{dk}}{{dt}} = \nabla \cdot \left[ {\rho \left( {\nu + {{\nu }_{t}}} \right)\nabla k} \right] + {{P}_{k}} - \rho \varepsilon ,$
где Pk – скорость генерации турбулентной энергии за счет скорости деформации среднего течения, ε – скорость вязкой диссипации турбулентной энергии, ν и νt – кинематические молекулярная и турбулентная вязкости соответственно, ρ – плотность среды, d/dt – субстанциональная производная. С использованием гипотезы Буссинеска для тензора рейнольдсовых напряжений генерация турбулентной энергии описывается формулой
${{P}_{k}} = \rho {{\nu }_{t}}2{{S}_{{ij}}}{{S}_{{ij}}} = \rho {{\nu }_{t}}{{S}^{2}},$
где Sij и $S = {{\left( {2{{S}_{{ij}}}{{S}_{{ij}}}} \right)}^{{1{\text{/}}2}}}$ – компоненты и модуль тензора скоростей деформации среднего течения.

Для дискретизации исходной системы уравнений применяется метод конечного объема второго порядка точности по пространству. Аппроксимация производных в центре контрольного объема находится на основе значений функции в центрах контрольных объемов. Непрерывная аппроксимация (n + 1)-го порядка точности полевой функции f в окрестностях контрольного объема P строится с помощью разложения в ряд Тейлора относительно центра контрольного объема:

$\hat {f}\left( {\mathbf{x}} \right) = {{f}_{P}} + {{\left( {\nabla f} \right)}_{P}}\left( {{\mathbf{x}} - {{{\mathbf{x}}}_{P}}} \right) + \ldots + \frac{1}{{n!}}{{\left( {{{\nabla }^{n}}f} \right)}_{P}}{{\left( {{\mathbf{x}} - {{{\mathbf{x}}}_{P}}} \right)}^{n}}.$

Оценка ошибки дискретизации для схем второго порядка определяется интегрированием ошибки дискретизации по контрольному объему и для случая изотропной сетки с сеточным шагом Δ имеет вид

$e\left( f \right) = \frac{{{{{{\Delta }}}^{2}}}}{{24}}\left| {\frac{{{{\partial }^{2}}}}{{\partial {{x}_{i}}\partial {{x}_{i}}}}f} \right| = \frac{{{{{{\Delta }}}^{2}}}}{{24}}\left| {{{\nabla }^{2}}f} \right|.$
Потребуем, чтобы разностное решение на данной сетке удовлетворяло условию малости относительной ошибки дискретизации
$e\left( f \right) < \delta \left| f \right|,$
где δ – некоторая малая величина. В результате можно сформулировать условие на шаг сетки:

${{{{\Delta }}}^{2}} < 24\delta \frac{{\left| f \right|}}{{\left| {{{\nabla }^{2}}f} \right|}}.$

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

$\rho \left( {\nu + {{\nu }_{t}}} \right){{\nabla }^{2}}k + {{P}_{k}} \approx 0.$
В данном случае детализация сетки должна быть достаточной для описания диффузионного переноса турбулентной энергии из области ее генерации. Условия реализуемости для тензора рейнольдсовых напряжений (см. [3]) позволяют сделать оценку сверху для производства турбулентной энергии:
${{P}_{k}} \leqslant 0.3\rho kS.$
В результате для модуля лапласиана турбулентной энергии получаем оценку
$\left| {{{\nabla }^{2}}k} \right| \leqslant \frac{{0.3kS}}{{\nu + {{\nu }_{t}}}}.$
Требование малости относительной ошибки локальной интерполяции приводит к неравенству

$\frac{{e\left( k \right)}}{k} \leqslant \frac{{{{{{\Delta }}}^{2}}}}{{24}}\frac{{0.3S}}{{\nu + {{\nu }_{t}}}} < \delta .$

Окончательно для шага сетки в областях с доминированием сдвига средней скорости получаем неравенство

${{{{\Delta }}}^{2}} < 24\delta \left( {\frac{{\nu + {{\nu }_{t}}}}{{0.3S}}} \right) = {{C}_{1}}\left( {\frac{{\nu + {{\nu }_{t}}}}{{0.3S}}} \right),\quad {{C}_{1}} = 24\delta .$
Это неравенство можно переписать для локального сдвигового числа Рейнольдса:

${\text{R}}{{{\text{e}}}_{m}} = \frac{{{{{{\Delta }}}^{2}}S}}{{\left( {\nu + {{\nu }_{t}}} \right)}} < {{C}_{2}},\quad {{C}_{2}} = 80\delta .$

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

$\left( {\nu + {{\nu }_{t}}} \right){{\nabla }^{2}}k - \varepsilon \approx 0,\quad \left| {{{\nabla }^{2}}k} \right| \leqslant \frac{\varepsilon }{{\nu + {{\nu }_{t}}}}.$

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

$\frac{{e\left( k \right)}}{k} \leqslant \frac{{{{{{\Delta }}}^{2}}}}{{24}}\frac{1}{{\left( {\nu + {{\nu }_{t}}} \right)}}\frac{\varepsilon }{k} < \delta ,$
и ограничение для шага сетки в области с доминированием диссипации турбулентной энергии
${{{{\Delta }}}^{2}} < 24\delta \left( {\nu + {{\nu }_{t}}} \right){{T}_{t}},$
где ${{T}_{t}} = k{\text{/}}\varepsilon $ – временной масштаб турбулентности.

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

${{\Delta }} < {{C}_{3}}{{\left[ {\left( {\nu + {{\nu }_{t}}} \right){\text{min}}\left( {\frac{1}{{0.3S}},{{T}_{t}}} \right)} \right]}^{{1{\text{/}}2}}},\quad {{C}_{3}} = \sqrt {24\delta } .$
Введем сеточный критерий Qm как отношение шага сетки к линейному некоторому масштабу Lm:
${{Q}_{m}} = {{\Delta /}}{{L}_{m}},$
${{L}_{m}} = {{\left[ {\left( {\nu + {{\nu }_{t}}} \right){\text{min}}\left( {\frac{1}{{0.3S}},{{T}_{t}}} \right)} \right]}^{{1{\text{/}}2}}}.$
Для оценки значений коэффициентов зададим относительную ошибку δ ~ 5%, и в результате имеем следующие значения коэффициентов:
${{C}_{1}} \approx 1.2,\quad {{C}_{2}} \approx 4,\quad {{C}_{3}} \approx 1.1.$
Окончательно критерий для получения сеточно-независимого критерия формулируется в виде неравенства

${{Q}_{m}} = {{\Delta /}}{{L}_{m}} < 1.$

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

${{\nu }_{t}} = {{c}_{\mu }}\frac{{{{k}^{2}}}}{\varepsilon },\quad {{T}_{t}} = \frac{k}{\varepsilon },\quad {{c}_{\mu }} = 0.09,$
в результате получаем оценку для линейного масштаба Lm:
${{L}_{m}} \approx {{\left( {{{\nu }_{t}}{{T}_{t}}} \right)}^{{1{\text{/}}2}}} = {{c}_{\mu }}^{{1{\text{/}}2}}\frac{{{{k}^{{3{\text{/}}2}}}}}{\varepsilon } = 0.3{{L}_{t}},$
где ${{L}_{t}} = {{k}^{{3{\text{/}}2}}}{\text{/}}\varepsilon $ – масштаб турбулентности или характерный размер энергонесущих вихрей. Критерий для RANS моделирования дает ${{L}_{t}} > 3.3{{\Delta }}$, т.е. не менее трех шагов сетки должно приходиться на размер крупного турбулентного вихря. При использовании методов моделирования крупных вихрей (Large Eddy Simulation, LES) для моделирования турбулентного течения вдали от стенки требуется, чтобы шаг сетки удовлетворял условию (см. [3])
${{L}_{t}} > 10--4{{\Delta }}{\text{.}}$
Таким образом, сетка для LES моделирования должна быть в 3–4 раза детальнее (в каждом направлении) по сравнению с сеткой для RANS моделирования.

Менее популярный критерий оценки сеточной детализации для LES основан на локальной величине скорости сдвига течения (см. [4]):

${{\Delta }} < {{L}_{S}} = \sqrt {\frac{\varepsilon }{{{{{\left( {{{S}_{{ij}}}{{S}_{{ij}}}} \right)}}^{{3{\text{/}}2}}}}}} = {{2}^{{3{\text{/}}4}}}\sqrt {\frac{\varepsilon }{{{{S}^{3}}}}} .$
В случае равновесного сдвигового течения ${{\nu }_{t}}{{S}^{2}} \approx \varepsilon $ сдвиговый масштаб турбулентности LS оказывается практически равен введенному масштабу Lm:
${{L}_{m}} \approx 1.08{{L}_{S}}.$
В сделанных допущениях оказывается, что сеточные критерии для RANS и для LES моделирования совпадают.

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

3. ТЕСТОВЫЕ РАСЧЕТЫ

В качестве тестовых задач выбраны установившееся развитое турбулентное течение в плоском канале, течение за уступом и нестационарное обтекание полусферического тела. Использовались две низкорейнольдсовые RANS модели турбулентности, основанные на гипотезе вихревой вязкости: k-ω SST с поправкой на кривизну линий тока (см. [5]) и модель с эллиптической релаксацией и k–ε–ζ–f (см. [6]). Расчеты выполнялись с использованием программы σFlow, численный алгоритм которой базируется на методе конечного объема для пространственных уравнений гидродинамики на неструктурированной сетке. Программа вычислительной гидродинамики σFlow использовалась для решения широкого класса задач турбулентных течений (см. [7]).

3.1. Пристеночный пограничный слой

Оценим требование к сеточному шагу для моделирования пристеночного сдвигового турбулентного течения без использования пристеночных функций с разрешением вязкого подслоя. В вязком подслое, пренебрегая турбулентной вязкостью, получаем значение линейного масштаба ${{L}_{m}} = {{\left( {\frac{\nu }{{0.3S}}} \right)}^{{1{\text{/}}2}}}$ в безразмерных вязких масштабах $L_{m}^{ + } \approx 1.8$. Неравенство для введенного критерия Qm < 1 требует, чтобы шаг сетки по направлению к стенке был меньше двух вязких масштабов: Δy+ < 2, а положение первого пристеночного узла порядка единицы: y + (1) ~ 1. В используемом методе контрольного объема пристеночный шаг сетки в 2 раза больше расстояния до первого узла.

Для оценки влияния шага сетки на разрешение пристеночного слоя выполнены расчеты развитого установившегося турбулентного течения в плоском канале при вязком числе Reτ = 2000. Расчеты проводились на трех сетках с одинаковым коэффициентом сгущения узлов к стенке, равным 1.05, и разным шагом сетки по направлению к стенке. Первый пристеночный узел имел безразмерные координаты y + (1) = 0.36, 1.25 и 2.4. Отличие расчетного значения вязкого трения, полученного на более грубой сетке, от значения, полученного на самой детальной сетке, растет с увеличением пристеночного шага сетки и составляет значения 1 и 3% для модели k–ε–ζ–f и 5 и 7% для модели k–ω SST.

Введенный критерий совпадает с общепринятым требованием к моделированию турбулентного пристеночного слоя, которое заключается в расположении первого пристеночного узла сетки в вязком подслое y + (1) ≤ 1 (см. [8]). При нарушении критерия детализации сетки для разрешения пристеночных градиентов оказывается недостаточной для определения интегрального сопротивления.

3.2. Моделирование стационарного течения за уступом

Рассматривается плоскопараллельное турбулентное течение за обратным уступом. Расчеты выполнены в плоской постановке. Высота уступа – h = 1 м, число Рейнольдса ${\text{Re}} = \rho Uh{\text{/}}\mu $ = = 34 000. Высота расчетной области – 9h, длина – 60h.

На входе ставится условие полностью развитого турбулентного течения в плоском канале. На выходе заданы неотражающие граничные условия, на поверхности обтекаемого тела – условия прилипания. Детализация сетки вблизи стенок достаточна для полного разрешения пристеночного турбулентного слоя без введения пристеночных функций (для пристеночных узлов y+ ~ 0.5).

Расчеты выполнялись на равномерной сетке c выделением пристеночного сеточного слоя. Вне пристеночного слоя сетка квадратная. Для оценки сеточного критерия рассматривается влияние детализации сетки непосредственно вниз по течению за уступом в области слоя смешения и рециркуляционного течения. Прямоугольная область, в которой варьируется сеточная детализация, находится на расстоянии 0.2h от уступа и имеет высоту 2 h и длину 15h. Выполнены стационарные расчеты на пяти сетках, отличающиеся детализацией в выбранной области. Примеры сеток показаны на фиг. 1. Размеры шагов сетки Δ в выделенной области представлены в табл. 1. Оси x и y направлены вдоль и поперек течения соответственно.

Фиг. 1.

Течение за уступом. Сеток М2 (а) и M5 (б).

Таблица 1.  

Течение за уступом. Параметры расчетных сеток

Сетка Шаг сетки, Δ Количество ячеек, тыс.
M1 0.0145 372
M2 0.029 305
M3 0.059 280
M4 0.118 274
M5 0.23 272

Для анализа течения было выбрано две линии. Первая линия перпендикулярна направлению течения и расположена на расстоянии 1h от уступа по течению: x = 1h. Вторая линия расположена в сдвиговой области течения, идет вдоль течения от точки отрыва потока: y = 1h. Линии проходят через область максимальной генерации турбулентности в сдвиговом течении.

Для количественной оценки сходимости решения выполнен расчет изменения решения при удвоении сетки вдоль двух выбранных линий. Относительное изменение сеточного распределения некоторого поля f на сетке с шагом Δ при измельчении сетки в два раза рассчитывается в норме L2 следующим образом:

${{E}_{{{\Delta }}}}\left( f \right) = \left\| {{{f}_{{{\Delta }}}} - {{f}_{{{{\Delta /}}2}}}} \right\|{\text{/}}\left\| {{{f}_{{{{\Delta /}}2}}}} \right\|.$

Изменение решения вдоль линии, перпендикулярной течению, при измельчении сетки показано на фиг. 2 для распределений турбулентной энергии k и продольной U и поперечной компонент V средней скорости, полученных с использованием модели k–ε–ζ–f. Аналогичные результаты получаются и для второй RANS модели турбулентности. Решение монотонно сходится при увеличении детализации сетки, и распределения на самых детальных сетках отличаются менее, чем на 5%. Сходимость решения имеет первый порядок по Δ.

Фиг. 2.

Течение за уступом. Изменение решения при измельчении сетки для модели k–ε–ζ–f. Линии 1 и 2 соответствуют первому и второму порядку по шагу сетки соответственно.

На фиг. 3 приведены графики турбулентной энергии и сеточного критерия вдоль линии y = 1h. Решение на грубых сетках демонстрирует заниженное значение турбулентной энергии и существенное смещение максимума по потоку. На этих сетках наблюдаются обширные области со значением сеточного критерия Qm больше единицы. Решения на сетках М1 и М2 отличаются незначительно, при этом на сетках М1 и М2 есть небольшие области с Qm > 1.

Фиг. 3.

Распределение энергии турбулентных пульсаций k (а) и критерия Qm (б) вдоль линии y = 1h. Модель k‒ε–ζ–f.

Построенные на фиг. 4 распределения коэффициента трения вдоль нижней стенки и распределение энергии турбулентных пульсаций вдоль линии x = 1h показывают сходимость решения на сетках М1 и М2 и хорошее согласование с экспериментальными данными (см. [9]).

Фиг. 4.

Распределение коэффициента трения Cf вдоль нижней стенки (а) и распределение энергии турбулентных пульсаций k вдоль линии x = 1 h (б). Модель k–ε–ζ–f. Символы – эксперимент (см. [9]).

3.3. Моделирование нестационарного обтекания полуцилиндрического тела

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

Моделируется обтекание полуцилиндрического тела при угле атаки, равном нулю. Число Рейнольдса задачи Re = 45 000. Для исследования сеточной сходимости нестационарного решения в ближнем следе (на расстоянии >2L от тела) построено четыре сетки с разной детализацией в области следа. Высота расчетной области Ly = 40L, ширина в третьем направлении H = 1L (сетка трехмерная). Размеры вдоль течения – 4.5L от входной границы до передней кромки тела, 14.5L от задней кромки до выходной границы.

Для всех рассмотренных сеток детализация в окрестностях тела на расстоянии 1L от тела одинаковая. Около обтекаемого тела выделены пять слоев с последовательным уменьшением размеров ячейки сетки при приближении к телу (фиг. 5а). Размер шага сетки в каждом слое постоянен и составляет 0.00245, 0.0049, 0.0098, 0.0195 и 0.039 L соответственно. Расстояние от поверхности профиля до внешней границы слоя составляет 0.05, 0.2, 0.6, 1.0 и 2.0 L соответственно. В непосредственной близости тела построен пристеночный сеточный слой толщиной 0.03L, состоящий из 28 узлов с пристеночным шагом h = 9 × 10–5L. Сеточные линии в пристеночном слое сгущаются к стенке с коэффициентом сгущения 1.10. Слой вблизи передней кромки показан на фиг. 5б. Такая детализация пристеночного слоя дает значение y+ первого пристеночного узла ~0.2. Размер ячейки внешней сетки 0.31L.

Фиг. 5.

Обтекание профиля. Сетка М2: (а) – пять слоев сетки около тела, (б) – пристеночный слой, (в) – сетка с областью следа.

Для оценки сеточного критерия рассматривается влияние детализации сетки вниз по течению за обтекаемым телом. Для этого в области следа за телом выделяется подобласть с постоянной детализацией с квадратными ячейками. Высота подобласти 4L. Сетки отличаются пространственной детализацией в этой подобласти. Четыре сетки имеют следующие шаги сетки в области следа: M1 – 0.02L, M2 – 0.04L, M3 – 0.08L, M4 – 0.16L. Самая детальная сетка – М1, самая грубая сетка – М4. На фиг. 6 показаны распределения узлов вблизи обтекаемого тела для сеток М1 и М4.

Фиг. 6.

Обтекание профиля: (а) – сетка М1, (б) – сетка М4.

Расчеты выполнены с использованием модели k–ω SST с коррекцией на кривизну линий тока и модели k–ε–ζ–f. Основные детали численной схемы следующие: конвективная схема на скорость – QUICK; метод нестационарного шага для уравнения движения – неявная направленная схема второго порядка; конвективная схема на турбулентные величины – UMIST TVD; метод нестационарного шага для турбулентных величин – неявная первого порядка; шаг по времени Δt = 0.01 L/U.

Распределения осредненного давления и вязких напряжений на стенке обтекаемого тела показаны на фиг. 7. Картины осредненного течения в виде траекторий частиц жидкости представлены на фиг. 8. Более детальную информацию о численном моделировании и экспериментальном исследовании обтекания полуцилиндрического тела можно найти в [10].

Фиг. 7.

Обтекание профиля. Распределения осредненного давления Cp (а) и вязкого трения τ (б) на поверхности для моделей k–ω SST (линия 1) и k–ε–ζ–f (линия 2). Символы – эксперимент (см. [10]).

Фиг. 8.

Обтекание профиля. Линии тока осредненного течения для моделей (а) k–ω SST с коррекцией на кривизну линий тока и (б) k–ε–ζ–f.

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

Для оценки влияния детализации сетки в следе проводится сравнение осредненной продольной компоненты скорости U, моделируемой кинетической энергии турбулентных пульсаций $k\left( {{\text{mod}}} \right) = \left\langle k \right\rangle $, разрешаемой энергии пульсаций скорости $k\left( {{\text{res}}} \right) = \frac{1}{2}\left\langle {u_{i}^{'}u_{i}^{'}} \right\rangle $ и пульсаций тензора скоростей деформации $e = \left\langle {2S_{{ij}}^{'}S_{{ij}}^{'}} \right\rangle $. Треугольные скобки означают осреднение по времени, штрих относится к разрешаемой пульсации величины. Сравниваются распределения, полученные на разных сетках, вдоль трех выбранных линий. Поперечные к потоку линии y1 и y2 расположены на расстоянии 5L и 10L за телом, линия x1 идет от задней кромки обтекаемого тела вдоль направления течения.

На фиг. 9 и 10 представлены распределения осредненной моделируемой кинетической энергии турбулентности и разрешаемой кинетической энергии турбулентности, полученные на разных сетках с использованием модели k–ω SST вдоль выделенных линий.

Фиг. 9.

Обтекание профиля. Распределения вдоль линии x1: (а) – энергия разрешаемых пульсаций, (б) – осредненная моделируемая кинетическая энергия турбулентности, (в) – осредненная x-компонента скорости, (г) – пульсации скорости деформации.

Фиг. 10.

Обтекание профиля. Распределения вдоль линии y1: (а) – энергия разрешаемых пульсаций, (б) – осредненная моделируемая кинетическая энергия турбулентности, (в) –осредненная x-компонента скорости, (г) – пульсации скорости деформации.

Приведенные графики демонстрируют монотонную сходимость решения при увеличении детализации сетки, т.е. при уменьшении сеточного шага. Качественно можно сказать, что сетка М4 имеет очень грубую детализацию, что приводит к сильному размазыванию распределений в следе и занижению разрешаемых пульсаций. Детализация сетки М3 оказывается также недостаточной для получения сеточно-независимого решения. На сетках М1 и М2 получены близкие решения, но на более грубой сетке М2 разрешаемые и моделируемые пульсации немного ниже.

Для количественной оценки сходимости решения выполнен расчет изменения решения при удвоении сетки вдоль трех выбранных линий. Изменения распределений k(mod), k(res) и e, полученные с использованием k–ε–ζ–f модели, вдоль трех пространственных линий при измельчении сетки показаны на фиг. 11.

Фиг. 11.

Обтекание профиля. Изменение решения при измельчении сетки для модели k–ε–ζ–f. Линии 1 и 2 соответствуют первому и второму порядку по шагу сетки соответственно.

Из графиков видно, что наблюдается монотонная сходимость решения при увеличении детализации сетки. Порядок сходимости по шагу сетки лежит в диапазоне от 1 до 2. Наибольший порядок сходимости имеет величина разрешаемых пульсаций скорости, которая сходится со вторым порядком по пространственному шагу на линиях y1 и y2. Порядок сходимости осредненной моделируемой энергии турбулентности составляет 1.3–2, а порядок сходимости пульсации скорости деформации 1.3–1.8. Следует отметить наличие сеточной сходимости как среднего течения, так и пульсационной составляющей течения.

На всех линиях изменение k(mod) при переходе от сетки М2 к самой детальной сетке М1 меньше 5%. Для разрешенной составляющей турбулентной энергии k(res) это изменение составляет меньше 2%. Изменение осредненной продольной компоненты скорости не обладает монотонной зависимостью от сеточного шага. Отклонение этой величины при переходе от сетки М2 к М1 для всех выбранных линий меньше 0.5%. Приведенные числа и графики позволяют нам утверждать, что детализации сеток М1 и М2 являются достаточными для получения сеточно-независимого нестационарного решения.

Рассмотрим распределение критерия Qm для разных сеток. На фиг. 12 построены распределения критерия Qm вдоль линии, параллельной набегающему потоку и расположенной в области следа. Для сеток, разрешающих нестационарные течения М1 и М2, значение критерия не превосходит значения 1, а значения сеточного сдвигового числа Rem всюду меньше 3. Для самой грубой сетки М4 критерий имеет значение в диапазоне от 1.5 до 4. Сетка М3, имеющая недостаточное сеточное разрешение, дает локальные максимумы критерия Qm на уровне 1.5, а числа Rem – на уровне 9. Отсюда можно сделать следующий вывод. Значение критерия Qm для получения сеточно-независимого RANS решения должно быть меньше единицы, а значение сеточного сдвигового числа Rem – меньше 4. Возможно наличие локальных небольших областей со значением критерия Qm, превышающего 1, не приведет к ухудшению численного решения. Однако, если размеры этих областей будут сопоставимы с характерным размером задачи, то численное решение на такой сетке однозначно будет недоразрешенным, и потребуется детализация сетки в этих областях для получения сеточно-независимого решения.

Фиг. 12.

Распределения критерия Qm вдоль выбранной линии x1.

4. ЗАКЛЮЧЕНИЕ

Предложен простой критерий оценки сеточной детализации сетки для получения сеточно-независимого численного решения турбулентного течения в рамках RANS подхода с использованием полуэмпирической модели турбулентности. Критерий применен к стационарному течению за уступом и к нестационарному отрывному турбулентному течению периодического типа. Оценка приемлемости критерия сеточной детализации основывается на исследовании сеточной сходимости течений в выделенных областях. Показано, что предложенный критерий позволяет правильно оценить необходимую детализацию сетки для получения сеточно-независимого решения вдали от стенки. Если значение критерия ~1, то измельчение сетки в 2 раза приводит к изменению численного решения меньше, чем на 5%. Тестовые расчеты были ограничены приближением плоскопараллельного турбулентного течения и сетками с квадратными ячейками в исследуемой области.

Критерий может быть использован как для оценки полученного решения на сеточную независимость, так и для локальной адаптации сетки. Поскольку предложенный критерий позволяет оценить сеточную детализацию для получения URANS решения, то критерий может быть использован и как ограничитель в моделях с частичным разрешением турбулентного спектра, например, в моделях частично осредненных уравнений Навье–Стокса, (Partially Averaged Navier–Stokes, PANS). Критерий может быть индикатором для переключения модели турбулентности с переменным разрешением в вихреразрешающую ветку при достаточной для URANS модели сеточной детализацией.

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

  1. Oberkampf W.L., Roy C.J. Verification and Validation in Scientific Computing, Cambridge Univ. Press, 2010.

  2. Habashi W.G., Dompierre J., Bourgault Y., Ait-Ali-Yahia D., Fortin M., Vallet M.G. Anisotropic mesh adaptation: towards user-independent, mesh-independent and solver-independent CFD. Part 1: general principles // Inter. J. Numer. Meth. Fluid. 2000. V. 32. P. 725–744.

  3. Pope S.B. Turbulent Flows. Cambridge Univ. Press, 2000.

  4. Picano F., Hanjalić K. Leray-α Regularization of the Smagorinsky-Closed Filtered Equations for Turbulent Jets at High Reynolds Numbers // Flow Turbulence Combust. 2012. V. 89. P. 627–650.

  5. Isaev S.A., Baranov P.A., Zhukova Yu.V., Usachov A.E., Kharchenko V.B. Correction of the shear-stress-transfer model with account of the curvature of streamlines in calculating separated flows of an incompressible viscous fluid // J. Engineer. Phys. Thermophys. 2014. V. 87. № 4. P. 1002.

  6. Hanjalić K., Popovac M., Hadžiabdić M. A robust near-wall elliptic relaxation eddy viscosity turbulence model for CFD // Inter. J. Heat Fluid Flow. 2004. V. 25. № 6. P. 1047–1051.

  7. Дектерев А.А., Гаврилов А.А., Минаков А.В. Современные возможности СFD кода SigmaFlow для решения теплофизических задач // Современная наука: исследования, идеи, результаты, технологии. 2010. Т. 2. Вып. 4. С. 117–122.

  8. Wilcox D.C. Turbulence Modeling for CFD. DCW Industries, Inc., La Canada, CA, 1993.

  9. Driver D.M., Seegmiller H.L. Features of Reattaching Turbulent Shear Layer in Divergent Channel Flow // AIAA J. 1985. V. 23. № 2. P. 163–171.

  10. Isaev S., Baranov P., Popov I., Sudakov A., Usachov A., Guvernyuk S., Sinyavin A., Chulyunin A., Mazo A., Demidov D., Dekterev A., Gavrilov A., Shebelev A. Numerical simulation and experiments on turbulent air flow around the semi-circular profile at zero angle of attack and moderate Reynolds number // Comput. Fluid. 2019. V. 188. P. 1–17. https://doi.org/10.1016/j.compfluid.2019.03.013

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

Инструменты

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

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

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