Теплофизика высоких температур, 2022, T. 60, № 1, стр. 63-75

ВВЕДЕНИЕ

Полет самолетов происходит в условиях атмосферной турбулентности. Возможны также ситуации, когда они пересекают следы или реактивные струи от пролетевших ранее летательных аппаратов. Попадание турбулентных вихрей в воздухозаборник (ВЗ) самолета может существенно повлиять на его работу, увеличить турбулентные пульсации на входе в двигатель, в частности пульсации давления, а также вызвать отрыв потока либо привести к помпажу ВЗ. Экспериментально создать набегающий поток с заданными параметрами турбулентности очень сложно и затратно, а в случае высоких скоростей скорее всего невозможно. Выходом может быть расчетное исследование таких течений. Известно, что при численном моделировании течений с лидирующей ролью турбулентных эффектов предпочтительнее использовать вихреразрешающие подходы, которые позволяют явным образом описать нестационарный характер течения, обусловленный турбулентными вихрями, и получить пульсации всех параметров течения. Кроме того, они позволяют задать “синтетическую” турбулентность, которая будет имитировать турбулентность набегающего потока. С помощью методов RANS нельзя решить перечисленные задачи.

С учетом того, что в настоящее время активно развивается сверхзвуковая авиация, ведутся работы над проектами сверхзвуковых пассажирских самолетов исследование влияния турбулентности набегающего потока на работу сверхзвукового ВЗ в особенности актуально. В настоящей работе для решения данной задачи используется комбинированный RANS/ILES-метод высокого разрешения [1], который ранее с успехом применялся для расчета турбулентных течений в широком диапазоне скоростей, включая сверхзвуковые. Анализ доступной литературы показал, что в силу сложности задачи число публикаций, в которых представлены результаты расчетов течений в сверхзвуковых ВЗ вихреразрешающими методами, невелико. Некоторые из них приведены в обзорах литературы в статьях [2, 3]. Работ, в которых численно исследовано влияние турбулентности набегающего потока на течение в сверхзвуковых ВЗ, авторам настоящей статьи найти не удалось. Однако в работе [4] с помощью метода [1] и его усовершенствованного варианта [5] показано, что даже относительно небольшая неоднородность в температуре набегающего потока может оказывать существенное влияние на течение в сверхзвуковом ВЗ. Таким образом, исследование влияния атмосферной турбулентности на течение в сверхзвуковых ВЗ является практически значимой задачей. Учет турбулентного контента набегающего потока потребует сетки с мелкими ячейками, необходимыми для его разрешения, в области внешнего течения. Избежать чрезмерного измельчения сетки позволит использование методов высокого разрешения, к которым относится и RANS/ILES [1].

В рамках данной работы с помощью метода [1] исследовано влияние характеристик атмосферной турбулентности на течение в сверхзвуковом ВЗ при различном дросселировании. Геометрия и режимные параметры ВЗ взяты из экспериментальной работы [6]. В настоящей статье исследуются диссипативные свойства разностной схемы, использованной в [1]. Анализ влияния характеристик турбулентности набегающего потока на течение в сверхзвуковом ВЗ предполагается выполнить в следующих публикациях. Мотивация для выделения в отдельную статью следующая. Монотонные схемы с разностями против потока, которые используются в ILES, позволяют проводить расчеты сверхзвуковых течений со скачками уплотнения. Однако они имеют довольно высокий уровень схемной вязкости, что может приводить к затуханию высокочастотных колебаний и недостаточно корректному описанию этой составляющей турбулентного течения. Таким образом, оценка диссипативности метода [1] в широком диапазоне скоростей и сравнение его по этой характеристике с другими видами методов LES, реализованных в коммерческих программных комплексах, является важной для практических приложений RANS/ILES-метода задачей. В данной статье эта оценка произведена на примере моделирования распада однородной изотропной турбулентности – часто используемой задачи для тестирования численных методов. Это редко реализующееся в природе явление позволяет наглядно показать, насколько диссипативны используемые схема дискретизации конвективных членов и подсеточная модель. Например, в [7] с помощью данной задачи исследована диссипативность различных схем, реализованных в программном комплексе ЛОГОС. А в [8, 9] на примере распада однородной изотропной турбулентности производится оценка диссипативности модификации схемы MP. Следует отметить, что задача о распаде изотропной турбулентности обычно рассматривается для малых скоростей. В настоящей работе исследованы случаи с около- и сверхзвуковыми скоростями, что позволило оценить диссипативные свойства схемы [1] в режиме ILES для широкого диапазона скоростей, чего не было сделано для других численных схем.

В ходе расчетов получены энергетические спектры поля скоростей при различных значениях турбулентного числа Маха и зависимости энергии турбулентности от времени. Для несжимаемого случая произведено сравнение с соответствующими экспериментальными данными [10] и с результатами DNS (Direct Numerical Simulation) из работы [11]. Проанализировано влияние различных параметров и настроек метода [1] на его диссипативность. Полученные результаты сравниваются с результатами расчетов с помощью коммерческого программного комплекса ANSYS CFX с использованием различных подсеточных моделей турбулентности.

ПОСТАНОВКА ЗАДАЧИ

Моделирование производилось в кубической расчетной области на равномерной сетке. По всем пространственным координатам использовалось периодическое граничное условие. Расчет методом [1] в режиме ILES выполнялся на трех сетках, содержащих 46 × 46 × 46 (сетка 1), 64 × 64 × 64 (сетка 2), 90 × 90 × 90 ячеек (сетка 3). Расчеты c помощью CFX производились только на сетке 1. В качестве начального условия использовалось искусственное поле скорости, сгенерированное по заданному пространственному энергетическому спектру. Моделирование проводилось для трех режимов:

1) несжимаемого, соответствующего значению турбулентного числа Маха M_t$ \ll $ 1 (режим 1);

2) слабо сжимаемого с M_t = 0.3 (режим 2);

3) сжимаемого с M_t = 1.2 (режим 3).

Турбулентное число Маха определяется по соотношению

где $\bar {c}$ – осредненная по объему скорость звука; $\overline {u_{1}^{2}} ,~\overline {u_{2}^{2}} ,~\overline {u_{3}^{2}} $ – осредненные по объему квадраты компонент скорости.

Приведенные режимы различаются энергетическим спектром, используемым для генерации начального поля скорости. Параметры начального поля для каждого случая приведены в табл. 1. Расчет производился в течение определенного времени. Для несжимаемого случая время для удобства сравнения результатов с экспериментальными данными [10] выражено в физических единицах, а в остальных случаях для сопоставления с данными DNS [11] оно отнесено к временнóму масштабу начального поля. В заданные моменты времени производилось сохранение поля скорости и построение его пространственного энергетического спектра. Для режима 1 спектр сравнивался с результатами эксперимента [10]. Для всех режимов получены зависимости кинетической энергии турбулентности от времени, для режима 2 проведено сравнение с результатами DNS [11].

Известно, что спектр поля турбулентности состоит из трех интервалов, расположенных в порядке увеличения волнового числа k: генерационного интервала, содержащего наиболее крупные вихри и большую часть кинетической энергии турбулентности, инерционного интервала и диссипативного интервала, в котором кинетическая энергия наиболее мелких вихрей переходит в тепло [12]. Для правильного описания эволюции вихревых структур численный метод должен адекватно моделировать инерционный интервал энергетического спектра, в котором происходит процесс переноса энергии от более крупных вихрей к более мелким. Угол наклона кривой зависимости энергии от волнового числа E = f(k) на этом участке спектра определяется законом Колмогорова, согласно которому он должен соответствовать степенному закону ~k^–5/3. Задача о распаде однородной изотропной турбулентности удобна для проверки этой особенности. Поэтому угол наклона энергетических спектров, полученных в результате расчетов, сравнивается с законом Колмогорова в интервале волновых чисел, соответствующих инерционному.

ПАРАМЕТРЫ НАЧАЛЬНОГО ПОЛЯ СКОРОСТЕЙ

Для определения параметров начального поля скоростей по пространственному энергетическому спектру, заданному в виде зависимости энергии от волнового числа E(k), использовались следующие величины и соотношения, справедливые для случая однородной изотропной турбулентности [12].

Скорость диссипации кинетической энергии турбулентности

где ν – кинематическая вязкость.

Интегральный линейный масштаб турбулентности

где ${{k}_{t}} = \int_0^\infty {E\left( k \right)dk} $ – кинетическая энергия турбулентности.

Связь масштаба Тейлора λ со скоростью диссипации:

где ${{u}_{{{\text{rms}}}}} = \sqrt {\frac{{2{{k}_{t}}}}{3}} $ – среднеквадратическое значение пульсаций скорости; ${{\bar {\rho }}}$ и ${{\bar {\mu }}}$ – осредненные по объему значения плотности и динамической вязкости.

Интегральный временной масштаб

Число Рейнольдса, определенное по масштабу Тейлора:

Газ рассматривался как совершенный, коэффициент адиабаты γ = 1.4, газовая постоянная R = = 287 Дж/(кг К). Динамическая вязкость µ в зависимости от температуры T определялась по формуле Сазерленда

где µ₀ – контрольная вязкость при некоторой контрольной температуре T₀.

Начальные давление, температура и плотность задавались постоянными во всей расчетной области. Пульсации скорости продуцировались по заданному для каждого расчетного режима пространственному энергетическому спектру. Использованный метод генерации пульсаций скорости основан на задании в волновом пространстве случайно направленных фурье-образов скорости, модуль которых определяется по энергетическому спектру, и на дальнейшем применении обратного преобразования Фурье для получения поля скорости. Данный метод аналогичен использованному в работе [13].

Значения и порядок задания параметров начальных условий различны для каждого из трех режимов расчета.

Несжимаемый режим. Начальная температура T_н и начальное давление p_н равнялись 288 К и 10⁵ Па соответственно. Начальное поле скоростей генерировалось по энергетическому спектру E(k), полученному в эксперименте [10]. Значения контрольных параметров в формуле Сазерленда составляли T₀ = 288 К, µ₀ = 18.27 × 10^–6 Па с. Размер расчетной области равнялся 0.5 м. Значения чисел Маха, Рейнольдса и турбулентных масштабов вычислялись с помощью соответствующих соотношений на основе заданного распределения E(k). Их значения приведены в табл. 1. Сохранение поля скоростей и расчет его энергетического спектра производились, как и в эксперименте, в моменты времени 0.28 и 0.66 с.

Слабо сжимаемый режим. Турбулентное число Маха в начальный момент времени на данном режиме было принято равным M_t = 0.3. Значения остальных параметров начального поля соответствовали использованным в работах [8, 9, 11] для расчетных случаев с этим турбулентным числом Маха: плотность в начальный момент времени принималась равной ρ_н = 1 кг/м³, число Рейнольдса (4) составляло Re_λ = 72, размер расчетной области равнялся 2π м. Начальный спектр определялся в зависимости от волнового числа k следующим образом:

где A – константа, равная A = 1.3 × 10^–4; k₀ – волновое число, определяющее положение максимума, k₀ = 8.

По формуле (1) находилась средняя по объему скорость звука. Затем по ее значению вычислялась начальная температура T_н. Эта же температура принималась за контрольную температуру T₀ в зависимости Сазерленда (5). Из формул (2) и (3) путем исключения скорости диссипации выражается масштаб Тейлора:

Далее из формул (4) и (5) определялась начальная динамическая вязкость µ_н и контрольная вязкость µ₀. Значения всех параметров начального поля скорости для данного случая представлены в табл. 1. Сохранение поля скоростей и расчет пространственного спектра производились в моменты времени t/τ = 1.85 и t/τ = 7.0.

Сжимаемый режим. Начальные температура T_н и давление p_н полагались равными 288 К и 10⁵ Па. Турбулентное число Маха M_t = 1.2. Размер расчетной области составлял 2π м. Значения контрольных параметров в формуле Сазерленда: T₀ = 288 К, µ₀ = = 18.27 × 10^–6 Па с. Начальный спектр задавался зависимостью (6) c k₀ = 4. Значение константы A вычислялось по формулам (1) и (6) из условия обеспечения заданного значения числа Маха M_t:

Все прочие параметры определялись так же, как и в несжимаемом случае. Их значения представлены в табл. 1. Сохранение поля скоростей и расчет энергетического спектра производились в моменты времени t/τ = 1.94 и t/τ = 7.0.

ПАРАМЕТРЫ ВАРИАНТОВ РАСЧЕТОВ

В RANS/ILES-методе [1] конвективные потоки на гранях расчетных ячеек в уравнениях Навье–Стокса вычислялись с помощью схемы Роу:

Здесь f(q_L), f(q_R) − векторы конвективных членов уравнений с левой и правой стороны грани ячейки соответственно, |A| – “модуль” матрицы Якоби. Высокое разрешение метода достигалось за счет применения монотонной противопоточной схемы девятого порядка МР9 [14] для вычисления параметров q_R и q_L на гранях ячеек.

С помощью коэффициента α ≤ 1 регулируется уровень схемной вязкости. При α < 1 получается комбинация противопоточной и центрально-разностной схем с уменьшенным уровнем вязкости. Для несжимаемого и слабо сжимаемого случаев расчеты проводились для разных значений α с целью исследования его влияния на диссипативность метода. Для сжимаемого случая расчеты методом [1] проводились только при α = 1: при этом значении схема обладает свойством монотонности, что важно в случае сверхзвуковых течений. Также было оценено влияние на диссипативность шага по времени. Значения параметров, при которых проводились расчеты RANS/ILES-методом, приведены в табл. 2.

В программном комплексе ANSYS CFX реализован метод LES с явными подсеточными моделями. Моделирование распада однородной изотропной турбулентности проводилось с использованием моделей Смагоринского и WALE. Для аппроксимации конвективных членов в несжимаемом и слабо сжимаемом режимах использовалась центрально-разностная схема второго порядка. Однако данная схема не обладает свойством монотонности и не может быть применена для расчета течений со скачками уплотнения. В сжимаемом режиме помимо нее при расчетах использовалась схема высокого разрешения, которая, согласно документации CFX, может использоваться для течений с большими градиентами величин. Все расчеты в CFX были проведены на сетке 1, их параметры представлены в табл. 3. Далее для обозначения центрально-разностной схемы используется аббревиатура CD, а для схемы высокого разрешения ‒ HR.

РЕЗУЛЬТАТЫ РАСЧЕТОВ

Несжимаемый режим. На рис. 1 представлены пространственные энергетические спектры полей скорости, полученные в настоящих расчетах, спектр синтетического поля, используемого в качестве начального условия в расчетах, а также спектры [10]. На каждом из графиков спектры, полученные в результате расчетов, представлены в моменты времени t = 0.28 и 0.66 с. Здесь и далее на аналогичных рисунках большему моменту времени соответствует меньшая площадь под графиком спектра. На горизонтальной оси графиков отложены значения волнового числа k, которое обратно пропорционально масштабу пульсаций. Для заданной сетки максимальное из разрешимых волновых чисел ограничивается ее шагом, а минимальное – максимальным линейным размером расчетной области. На всех графиках вертикальная пунктирная линия обозначает правую границу интервала разрешимых на данной сетке волновых чисел. Наклонная пунктирная линия на всех рисунках показывает угол наклона инерционного интервала, соответствующий степенному закону Колмогорова.

На рис. 1а представлены результаты расчетов RANS/ILES-методом при различных значениях шага по времени, а на рис. 1б ‒ при различных значениях параметра α. Здесь и далее номера вариантов расчетов разными методами – это номера вариантов из соответствующей таблицы. Уменьшение шага по времени не приводит к заметному изменению уровня численной диссипации. Уменьшение же α вызывает значительное снижение диссипативности схемы в области мелких масштабов. Тем не менее в несжимаемом режиме на высоких частотах вблизи границы интервала разрешимых на сетке масштабов используемая в методе [1] схема обладает повышенным уровнем численной диссипации. В остальной же части интервала разрешимых масштабов спектр рассчитанного поля достаточно точно согласуется с экспериментальными данными. Полученный в расчетах спектр содержит как генерационный, так и инерционный интервалы, угол наклона последнего неплохо соответствует степенному закону.

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

Результаты вариантов расчетов с помощью CFX № 3 и методом [1] № 3 представлены на рис. 1д. Расчеты по CFX проведены с использованием модели WALE при значении константы 0.5, которое установлено по умолчанию. В этом случае достигается хорошее совпадение с экспериментальными данными во всем диапазоне разрешимых на сетке волновых чисел.

Анализ представленных результатов расчетов, а также результатов тестирования различных схем дискретизации конвективных членов, проведенного в [7] для программного комплекса ЛОГОС, показал, что повышенная диссипативность в области мелких масштабов является особенностью многих противопоточных схем или представляющих собой комбинацию противопоточной и центрально-разностной схем. Сочетание центрально-разностной схемы и явной подсеточной модели позволяет добиться хорошего совпадения с экспериментом путем калибровки констант модели. Однако центрально-разностные схемы не могут быть использованы для моделирования течений со скачками уплотнения, кроме того, требуется точная калибровка подсеточной модели, которая для сложных течений может быть затруднительна.

Рис. 2 демонстрирует влияние шага расчетной сетки на спектр. Здесь представлены спектры, полученные методом [1], для вариантов №№ 3–5. При увеличении числа ячеек, что при постоянном размере расчетной области соответствует уменьшению шага, верхняя граница интервала разрешимых волновых чисел смещается вправо. Участок с повышенной диссипативностью, примыкающий к этой границе, при уменьшении шага также смещается вправо, сохраняя свою ширину.

На рис. 3 показана зависимость от времени кинетической энергии турбулентности K_t, отнесенной к ее начальному значению K_t0. График построен для расчетов с помощью CFX № 3 и методом [1] № 3, показавших наилучшее совпадение с данными эксперимента [10]. В этом режиме течения для t/τ < 1.0 при расчете методом [1] наблюдается более высокая скорость уменьшения кинетической энергии турбулентности, чем при расчете с помощью CFX. Однако далее скорость уменьшения кинетической энергии при расчете методом [1] становится меньше, чем при расчете с помощью CFX. В результате за оставшийся интервал времени разница в абсолютных величинах кинетической энергии турбулентности в этих двух расчетах уменьшается.

Дополнительной иллюстрацией полученных результатов может служить визуализация турбулентных структур с использованием изоповерхностей Q-критерия в различные моменты времени, приведенная на рис. 4, 5. Q-критерий определяется как второй инвариант тензора градиента скорости [15] и представляет собой разность норм матриц тензоров завихренности и деформации. Его положительные значения соответствуют преобладанию завихренности над деформацией. Для большей наглядности изоповерхности окрашены в цвета поля скорости. Рис. 4а, 4б демонстрируют влияние коэффициента α при расчете RANS/ILES-методом: уменьшение α заметно снижает диссипативность метода, что проявляется в большем количестве вихревых структур. На рис. 4в показаны изоповерхности Q-критерия, полученные по результатам расчета № 3 с помощью CFX. Из всех расчетов на сетке 1 в данном случае наблюдается наилучшее совпадение спектра с экспериментом, что проявляется в большей, чем на рис. 4а, 4б, наполненности области течения турбулентными структурами. На рис. 5 представлено влияние шага сетки на картину течения. Видно, что уменьшение шага приводит к существенному наполнению расчетной области мелкомасштабными структурами.

Слабо сжимаемый режим. На рис. 6 представлены пространственные энергетические спектры поля скорости в моменты времени t/τ = 1.85 и 7.0, а также энергетический спектр поля скорости, использованного в качестве начального условия. Как и при несжимаемом режиме, увеличение α приводит к росту диссипативности схемы в области меньших масштабов. При меньших значениях α наклон спектра в инерционном интервале ближе к значению, соответствующему степенному закону ~k^–5/3. На рис. 6в представлены спектры полей скорости, полученные на сетке 2 этим же методом с различными значениями шага по времени. Видно, что изменение шага по времени в данном режиме также не приводит к заметному изменению диссипативности метода.

Сравнение результатов, полученных с помощью CFX с моделью Смагоринского при значении константы модели 0.1 (расчет № 4) и RANS/ILES-методом (расчет № 6), показывает (рис. 6г), что спектры довольно близки как по уровню энергии E при одинаковых k, так и по углу наклона графиков энергии Е на инерционном интервале. Однако они демонстрируют несколько различную динамику изменения во времени. Так, при t/τ = 1.85 спектр поля, полученного методом [1], в области малых масштабов расположен несколько выше, чем спектр поля с использованием CFX. При t/τ = 7.0 спектр, рассчитанный RANS/ILES-методом, расположен ниже, чем спектр, полученный с помощью CFX, во всем интервале масштабов.

На рис. 7а, 7б представлена зависимость от времени кинетической энергии турбулентности, отнесенной к ее начальному значению, при расчетах на сетках 1 и 2 соответственно. Расчеты RANS/ILES-методом с α = 0.1 как на сетке 1, так и на сетке 2 демонстрируют хорошее совпадение с результатами DNS [11]. Увеличение α приводит к существенному росту диссипативности на начальном интервале времени при t/τ < 2, что проявляется в несколько заниженных значениях кинетической энергии турбулентности в сравнении с результатами DNS [11]. Кривая, полученная с использованием CFX № 4, близка к данным DNS [11] в интервале t/τ < 2. При увеличении t/τ она идет выше, чем кривые, полученные с помощью DNS [11] и методом [1] для расчета № 6.

Рис. 8 иллюстрирует влияние шага сетки на энергетический спектр поля скоростей. Расчеты выполнены методом [1] при α = 0.1. Размер расчетной области постоянен. Как и следовало ожидать, уменьшение шага ведет к увеличению максимального значения волнового числа, разрешимого на сетке, в результате чего правая граница спектра смещается в сторону больших значений волнового числа.

Сжимаемый режим. На рис. 9 показаны энергетические спектры поля скорости при t/τ = 1.94 и 7.0, полученные при расчетах разными методами. Расчет RANS/ILES-методом проводился при α = 1, в случае CFX использовалась подсеточная модель WALE со значением константы 0.5 и схема CD для конвективных членов. В обоих расчетах угол наклона инерционного интервала близок к углу, определяемому степенным законом ~k^–5/3. Однако в области больших волновых чисел при k > 13 1/м разностная схема, реализованная в методе [1], демонстрирует немного большую диссипативность в сравнении со схемой CD в расчете № 5 с помощью CFX. В расчете с помощью CFX (рис. 9б) для аппроксимации конвективных членов использовалась схема HR, пригодная для случая течений со скачками уплотнения. По сравнению с разностной схемой, реализованной в методе [1], у схемы HR диссипативность больше при k > 6 1/м.

Также исследовано влияние на уровень диссипации схемы HR подсеточных моделей с различными константами. На рис. 9в показаны результаты расчетов с использованием модели WALE при значениях константы 0.3 и 0.7, а на рис. 9г – модели Смагоринского со значениями константы 0.03 и 0.35. Видно, что схема HR во всех случаях имеет более высокий уровень численной диссипации по сравнению со схемой [1]. При этом изменение подсеточной модели и уменьшение подсеточной вязкости путем калибровки констант моделей не ведет к существенному уменьшению диссипативности.

На рис. 10 показана зависимость от времени кинетической энергии турбулентности, отнесенной к ее начальному значению. При расчете с помощью CFX с использованием схемы HR наблюдается более быстрый спад K_t по сравнению с RANS/ILES-методом. Особенно это заметно при t/τ < 1. Однако к t/τ = 7.0 значения кинетической энергии в обоих расчетах становятся одинаковыми. При расчете с помощью CFX со схемой CD значения K_t близки к результатам расчета с помощью метода [1]. При t/τ < 1 значения K_t при расчете в CFX уменьшаются немного быстрее, чем при расчете RANS/ILES-методом. На оставшемся промежутке наблюдается обратная картина: уже при расчете методом [1] кинетическая энергия убывает быстрее.

На рис. 11 представлена визуализация турбулентных структур в различные моменты времени с помощью изоповерхностей Q-критерия, окрашенных в цвета числа Маха. Поле, полученное в расчете с помощью CFX со схемой HR, содержит меньшее количество мелкомасштабных вихревых структур, чем поле, полученное в расчете методом [1]. Это согласуется с формой энергетических спектров, представленных на рис. 9б: спектр, полученный с помощью CFX, в области больших волновых чисел лежит ниже, чем спектр, полученный RANS/ILES-методом.

ЗАКЛЮЧЕНИЕ

Проведено исследование диффузионных свойств RANS/ILES-метода в режиме ILES на примере моделирования распада однородной изотропной турбулентности в широком диапазоне чисел Маха, включая сверхзвуковые. Рассмотрено влияние настроек схемы, шага по времени и размеров ячеек расчетной сетки на энергетические спектры турбулентности.

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

Также выполнены расчеты с помощью программного комплекса ANSYS CFX с использованием различных подсеточных моделей. В несжимаемом режиме при использовании схемы CD с подсеточной моделью WALE после подбора константы модели удалось получить лучшее по сравнению с RANS/ILES-методом совпадение с экспериментом. В слабо сжимаемом режиме результаты, полученные с помощью CFX, близки к результатам расчетов RANS/ILES-методом. При больших числах Маха (M_t = 1.2) расчеты с помощью CFX проводились с использованием двух различных схем для аппроксимации конвективных членов: CD и HR. Результаты расчетов продемонстрировали повышенный уровень диссипативности второй из них, ориентированной на расчет течений с большими градиентами, в области больших волновых чисел. Кроме того, смена подсеточной модели и уменьшение подсеточной вязкости путем калибровки константы при использовании схемы HR не приводили к заметному уменьшению диссипативности в области больших волновых чисел. В свою очередь схема CD при соответствующем подборе константы в подсеточной модели WALE показала меньшую диссипативность в области больших волновых чисел, чем схема, реализованная в RANS/ILES-методе. Однако схема CD не обладает свойством монотонности, что делает невозможным ее применение для течений со скачками уплотнения. В этом случае при расчетах с помощью CFX необходимо использовать схему HR, обладающую повышенным уровнем численной диссипации. Это осложняет практическое использование программного комплекса. Разностная схема, реализованная в RANS/ILES-методе, имеет бо́льшую универсальность и относительно низкий уровень численной диффузии в широком диапазоне скоростей. Важным для практики ее достоинством является возможность расчета с высокой точностью сверхзвуковых течений со скачками уплотнения без какой-либо модификации или дополнительной настройки метода.

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