Известия РАН. Механика жидкости и газа, 2021, № 5, стр. 25-33

МЕТОДИКА РАСЧЕТА

Для расчета параметров интересующих нас тестовых течений (пограничного слоя пластины и слоя смешения двух потоков) были составлены программы на языке FORTRAN. В этих программах использовались вычислительные методики [11, 12] 43-летней давности, обеспечивающие всего лишь первый порядок аппроксимации. Тем не менее в данной работе высокую точность расчетных характеристик удалось получить за счет использования конечно-разностных сеток с мелкими ячейками, что ранее было невозможно в связи с недостаточным ресурсом вычислительной техники того времени.

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

Использовалась декартова система координат (x,y) с осью “x” вдоль пластины. Применялась неявная конечно-разностная схема первого порядка точности. Аппроксимация производных первого порядка производилась при помощи односторонних разностей, ориентированных “по потоку”. Узлы конечно-разностной сетки сгущались у пластины. Решение конечно-разностных уравнений осуществлялось методом прогонки. Счет осуществлялся, как упоминалось выше, маршевым способом.

В случае расчета пограничного слоя пластины, кроме модели турбулентности, необходимо дополнительно задать числа Прандтля (ламинарное Pr_l и турбулентное Pr_t), а также зависимость ламинарной вязкости μ_l от температуры T. Расчеты проводились как для случая обтекания пластины воздухом, так и для случая обтекания пластины азотом.

В случае азота задавалось Pr_l = 0.71 (как в [9]), Pr_t = 0.85 (как в [9]), а для ламинарной вязкости, как и в [9], использовалось следующее соотношение [13]:

В случае воздуха задавалось Pr_l = 0.71 (как в [8]), Pr_t = 0.85 (как для азота), а для ламинарной вязкости использовалась следующая зависимость от температуры [8]:

Коэффициенты теплопроводности (как воздуха, так и азота) определялись по справочнику [14].

Методика расчета слоя смешения основывалась на поиске автомодельного решения для упомянутой выше системы уравнений пограничного слоя и, в основном, использовала приемы из работы [12]. Рассматривалась автомодельная задача о сжимаемом слое смешения. Применялась декартова система координат (x, y) с осью “x”, располагающейся вдоль границы, разделяющей два потока до их смешения. Использовалась автомодельная переменная ξ = y/x. Турбулентная вязкость предполагалась много больше ламинарной.

При расчете автомодельного слоя смешения, кроме модели турбулентности, необходимо дополнительно задать турбулентное число Прандтля Pr_t и отношение удельных теплоемкостей k. В [15] отмечается, что в условиях свободной турбулентности (струи, следы за телом) чаще всего полагают Pr_t = 0.7. Это значение использовалось и в нашем случае. Отношение удельных теплоемкостей полагалось равным k = 1.4.

РЕЗУЛЬТАТЫ РАСЧЕТОВ ПОГРАНИЧНОГО СЛОЯ ПЛАСТИНЫ

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

Рассмотрим гиперзвуковой режим обтекания пластины азотом (M = 7.87, T_w∕T_r = 0.48, см. [8]). Здесь T_w – температура пластины, T_r – температура восстановления [9]. На рис. 1 сплошной линией представлена зависимость коэффициента трения C_f от числа Рейнольдса Re_ϑ, посчитанного по толщине потери импульса. Здесь применялась модель турбулентности SA, а в качестве начальных профилей были использованы приближенные соотношения для пограничного слоя пластины.

Видно, что в процессе счета имел место процесс, внешне похожий на турбулизацию ламинарного пограничного слоя с переходом при Re_ϑ = 1000. Смоделируем более ранний “переход” – при Re_ϑ = 500. Для этого введем дополнительное локальное слагаемое C_turb в источниковый член модели SA (аналог так называемого “trip” – члена модели SA, см. [3]). В приближении пограничного слоя это слагаемое имеет следующий вид:

Здесь $\tilde {\nu }$ – преобразованный коэффициент турбулентной вязкости, введенный в [3], для которого и сформулирована модель SA. Видно, что это слагаемое заметно отличается от нуля только в некоторой окрестности Re_ϑ = 500, поэтому воздействие этого слагаемого носит локальный характер.

На рис. 1 пунктирной линией представлена зависимость коэффициента трения C_f от числа Рейнольдса Re_ϑ для случая более раннего “перехода” (при Re_ϑ = 500), смоделированного описанным выше способом. Видно, что внесение дополнительного возмущения уже при Re_ϑ > 2000 практически не влияет на зависимость C_f (Re_ϑ). Учитывая эту особенность, во всех проведенных в этой работе расчетах использовались в качестве начальных профилей грубые оценки пограничного слоя пластины, и контролировалась достаточная удаленность финального значения Re_ϑ (ради достижения которого и проводился расчет) от “переходного” значения Re_ϑ.

С целью проверить применимость используемого в данной работе приближения пограничного слоя к описанию гиперзвукового пограничного слоя пластины, были проведены расчеты для режима с азотом (M = 7.87, T_w/T_r = 0.48), ранее рассчитанного как с помощью DNS-подхода [8], так и с помощью RANS-подхода [9]. В последнем случае использовалась модель турбулентности SA.

Расчеты в приближении пограничного слоя (BL-расчеты) проводились на разных сетках. Количество узлов N в сечении пограничного слоя варьировалось от 250 до 1000. Сравнение некоторых характеристик, соответствующих числу Рейнольдса, посчитанному по толщине потери импульса, Re_ϑ = 8440 (BL-расчеты и RANS-расчет [9]), приведены ниже в табл. 1.

Здесь (см. определения в [9]): Re_ϑ = ρ_∞ × u_∞ × ϑ/μ_∞, Re_τ = ρ_w × u_τ × δ/μ_w, Re_δ2 = ρ_∞ × u_∞ × ϑ/μ_w, ϑ – толщина потери импульса, H = δ*/ϑ – формпараметр, δ* – толщина вытеснения пограничного слоя, δ – толщина пограничного слоя (99% от скорости набегающего потока), Z_τ = ν_w/u_τ, u_τ – скорость трения, RA = 2 × C_h/C_f, C_h = q_w/[ρ_∞C_pU_∞(T_r – T_w)], B_q = q_w/(ρ_wC_pu_τT_w).

Видно слабое влияние числа узлов N на точность BL-расчета. Так, изменение числа узлов вдвое (от 500 до 1000) привело к изменению коэффициента трения всего на 0.125%. Полученная высокая точность BL-расчета (0.125%) позволяет судить о том, велика ли погрешность, обусловленная использованием приближения пограничного слоя. То есть, как сильно различаются полученные результаты BL-расчета от результатов RANS-расчета [9]. Видно, что использование приближения пограничного слоя привело к довольно слабым изменениям характеристик от их значений при RANS-расчете в [9]. Так, коэффициенты трения в BL-расчете и RANS-расчете отличаются тоже на 0.125%. Полученный результат демонстрирует применимость приближения пограничного слоя к описанию гиперзвукового пограничного слоя пластины.

Результаты расчета того же режима, но для несколько большего значения Re_ϑ = 9714 (см. [8]), представлены в табл. 2.

В табл. 2 также приведены результаты расчетов при использовании модели турбулентности Nut-92. Сравним с DNS-данными результаты для коэффициента сопротивления C_f при использовании различных моделей турбулентности. Видно, что модель SA завышает C_f на 1.4%, а модель Nut-92 занижает C_f на 14.2%.

Результаты расчетов остальных режимов (с воздухом в качестве рабочего газа) представлены в табл. 3–6.

Относительная погрешность (от DNS-расчетов [8]) в предсказании коэффициента трения моделями SA и Nut-92 для рассмотренных режимов (см. табл. 2–6) приведена в табл. 7.

Из анализа этой таблицы видно, что модель Nut-92 систематически занижает значения коэффициента трения, причем максимальная погрешность (20.8%) реализуется для режима с максимальным значением числа Маха (М = 13.64). Модель SA точнее, ее максимальная погрешность составляет всего 6.9%.

Повысить точность модели Nut-92 можно путем введения корректирующего множителя в источниковый член I_compr , отвечающий за учет скоростной сжимаемости (см. [4]):

Здесь (см. [4]) C₅ = 3.5, ν_t – коэффициент турбулентной вязкости, a – скорость звука,

Заметим, что в рассматриваемом нами частном случае

Предлагаемый корректирующий множитель I_cor имеет следующий вид:

Здесь ν – коэффициент ламинарной вязкости, M_t – аналог турбулентного числа Маха (см. [9]):

Относительная погрешность в предсказании коэффициента трения модифицированной таким образом модели Nut-92m для рассмотренных режимов (см. табл. 2–6) приведена в табл. 7 в нижней строке. Видно, что в результате модификации погрешность удалось понизить с 20.8 до 3.6%.

Заметим, что корректирующий множитель обращается в единицу в случае несжимаемого газа (M_t = 0). Корректирующий множитель обращается в единицу и в случае сжимаемого слоя смешения, когда ν_t$ \gg $ ν. Корректирующий множитель менее единицы в случае пограничного слоя.

РЕЗУЛЬТАТЫ РАСЧЕТОВ СЛОЯ СМЕШЕНИЯ

Рассмотрим смешение двух сверхзвуковых турбулентных потоков с одинаковыми параметрами статического давления и статической температуры [10]. Числа Маха для смешивающихся потоков были равны M_min = 1.3 и M_max = 3.3. Параметр спутности был равен m = 1.3/3.3 = 0.394, а конвективное число Маха M_c = (U_max – U_min)/(2a) = 1 [10].

Результаты расчетов автомодельной задачи для этого режима с использованием рассматриваемых моделей турбулентности приведены в табл. 8. Как упоминалось выше, модифицированная модель Nut-92m в случае автомодельного слоя смешения идентична модели Nut-92.

Здесь δ_ω – толщина слоя смешения, определенная с использованием автомодельной переменной ξ = y/x и максимального градиента скорости [10], Δξ – шаг расчетной сетки. Видно, что толщина слоя смешения изменилась не более чем на 0.2% при уменьшении вдвое шага сетки.

Эффект скоростной сжимаемости принято оценивать с помощью относительной толщины слоя смешения ${{\bar {\delta }}_{\omega }}({{M}_{c}}) = {{\delta }_{\omega }}\left( {{{M}_{c}}} \right){\text{/}}{{\delta }_{\omega }}\left( 0 \right)$. Если при прежнем параметре спутности m = 1.3/3.3 изменять температуру смешиваемых потоков, то таким образом можно изменять число Маха конвективного, в том числе и до малого значения M_c = 3 × 10^–3, практически соответствующему случаю несжимаемой среды M_c = 0. Результаты таких расчетов приведены в табл. 9.

Видно, что толщина слоя смешения и в этом случае изменилась не более чем на прежние 0.2% при уменьшении вдвое шага сетки.

Заметим, что в случае смешения потоков несжимаемой жидкости с одинаковой плотностью имеет место следующая оценка [16]:

В рассматриваемом случае (m = 0.394) получаем следующее значение δ_ω(0) = 0.074. Это значение превышает как оценку по модели SA (на 15%, см. табл. 9), так и оценку по модели Nut-92 (на 12%, см. табл. 9). То есть, модели SA и Nut-92 показали в несжимаемом случае примерно одинаковую точность.

В сжимаемом случае (M_c = 1) значение DNS-расчета [10] (${{\bar {\delta }}_{\omega }}(1) = 0.332$) оказалось ниже как расчетного значения по модели Nut-92 (0.02872/0.06527 = 0.440, на 33% по сравнению с DNS-расчетом, см. табл. 8), так и расчетного значения по модели SA (0.06122/0.06265 = 0.977, на 194% по сравнению с DNS-расчетом, см. табл. 8). То есть, модель Nut-92 показала в сжимаемом случае удовлетворительную точность, а модель SA – неудовлетворительную.

Расчеты, аналогичные описанным выше, были проведены и для других значений конвективного числа Маха M_c. Результаты таких расчетов, относящиеся к относительной толщине слоя смешения ${{\bar {\delta }}_{\omega }}\left( {{{M}_{c}}} \right)$, приведены на рис. 2 сплошной и штриховой линиями для моделей Nut-92 и SA соответственно. Здесь же приведены и известные результаты различных трехмерных DNS-расчетов. Эти данные заимствованы, в основном, из работы [10] (см. рис. 3 в указанной работе).

Если учитывать данные всех приведенных на рис. 2 DNS-расчетов, то результаты модели Nut-92 окажутся внутри разброса DNS-данных. Что касается результатов модели SA, то они располагаются гораздо выше DNS-данных.

ЗАКЛЮЧЕНИЕ

1. Модель Nut-92 хорошо учитывает эффект скоростной сжимаемости в слое смешения в диапазоне M_c < 1.5. Модель SA эффект скоростной сжимаемости в слое смешения в диапазоне чисел 0.5 < M_c < 1.5 учитывает неудовлетворительно.

2. В сверхзвуковом пограничном слое (M = 2.5) все модели показали хорошую точность (расхождение с DNS-расчетом по коэффициенту трения не более 2.8%).

3. В гиперзвуковых пограничных слоях (13.64 > M > 5.84) эффект скоростной сжимаемости лучше учитывает модель SA (расхождение с DNS-расчетами по коэффициенту трения не более 6.9%), чем модель Nut-92 (расхождение не более 20.8%). В целом модель SA обеспечила хорошую точность, а модель Nut-92 – удовлетворительную. Модифицированная модель Nut-92m по точности описания коэффициента трения оказалась несколько точнее модели SA (расхождение не более 3.6%).