CC BY
4«Труды МАИ». Выпуск № 82
www.mai.ru/science/trudy/
Вопросы моделирования турбулентности для расчета сверхзвуковых
высокотемпературных струй
Кравчук М. О.*, Кудимов Н. Ф.**, Сафронов А. В.***
Центральный научно исследовательский институт машиностроения, ЦНИИмаш, ул. Пионерская, 4, Королев, Московская область, 141070, Россия *e-mail: xxxmaxxxx3@yandex. ru **e-mail: nfkudim @tsiimash. ru ***e-mail:. avsafron@tsniimash.ru
Аннотация
В работе, на основе осредненных уравнений Навье-Стокса с модифицированной моделью турбулентной вязкости SST, представлены результаты расчетов сверхзвуковых холодных и горячих турбулентных струй. Валидация численной методики проведена на экспериментальных данных. Даны рекомендации по применению моделей турбулентности для горячих сверхзвуковых струй.
Ключевые слова: сверхзвуковые струи, турбулентность, изобарические струи, начальный участок струи.
1. Введение
Сверхзвуковые струйные течения широко встречаются в авиационной и ракетной технике. Правильное моделирование таких течений является сложной задачей, включающей физически корректное моделирование процессов турбулентности и разрешение сложной ударно-волновой структуры. Наиболее чувствительны к параметрам турбулентности, являются сверхзвуковые изобарические струи, поэтому в данной статье рассматриваются проблемы моделирования таких течений методом RANS с различными моделями турбулентности и методом LES. Приведено сравнение расчетов с экспериментом. Даны рекомендации по применению моделей турбулентности для горячих сверхзвуковых струй.
2. Постановка задачи
Расчетная область представляла собой цилиндр. Диаметр основания цилиндра равен 40Ra (где Ra - радиус среза сопла), длина цилиндра - 100 Ra. Расчет проводился от среза сопла для холодной и горячей струи. На всех границах расчетной области - давление и температура окружающей среды (в случае сверхзвукового течения данные и параметры потока «сносятся» из расчетной области).
Расчетная сетка обеспечивала порядка 20 ячеек на радиус среза сопла. Для метода RANS была использована расчетная сетка в осесимметричной постановке (двумерная), общее число расчетных ячеек порядка 200 тыс. Время счёта на 8 логических ядрах - 3 часа. Для метода расчета LES использовалась трёхмерная
сетка, число расчетных ячеек порядка 5млн. Время счёта на 8 логических ядрах - 14 суток. Использовался расчетный пакет Ansys Fluent V.15.0, для модификаций стандартных моделей турбулентности использовались user-defined function -подпрограммы, написанные на языке С++, которые компилировались в библиотеки и подключались к основной программе.
При описании турбулентных течений в рамках осредненных по Рейнольдсу/Фавру уравнений Навье-Стокса, для замыкания системы уравнений использовалась двухпараметрическая модель Ментера SST[1] c учетом сжимаемости Вилкокса[2].
Модель турбулентности SST[2] является двухпараметрической и предполагает решение уравнений для энергии турбулентности k и удельной скорости диссипации a:
+ p0Mj -[p + Oafy)
dx . j
YP -вра2 + 2( 1 - F1)
pOa2 dk da ;
a dx . dx jj
где
P = T..
V dx .
F1 = tanh(arg'4 ),arg1 =
= mm max
slk 500v'] 4P°a2k
* ' Л ' Л
в ad d2a ) CDk d2
'2
CD = max
1р°а2 dk da 1 20
a dx ■ dx ■'
dx ■ dx ■ j J
2
F2 = tanh(arg2 ),arg2 = max 2
Константы в уравнениях определяются следующим образом:
ф = Fфф + (1 - Fl )ф
= в1 -ас1к
2
ул=-±- - °п ,а1л = 0.85,а л = 0.5,в = 0.075
'1 в* Г* к1 ®1 1 2
Г2 =в* 2 = 10 а т2 = 0 856в = 00828
в* = 0.09( 1 +1.5^(М*)),к = 0.42^ = 0.31
F(Mt) = (Щ2 -М1)ЩЩ -мш)
М2 = 22 ,М{0 = 0.25,а = л[уКГ
турбулентная вязкость связана с энергией турбулентности и скоростью
ралк
диссипации следующим соотношением: п =_1_ где s = /2ГХ7.
* max(ala,SF2) , V Ч Ч
Для численного решения уравнений применяется TVD схема 2-ого порядка аппроксимации по времени и пространству.
На срезе сопла задаются кинетическая энергия пульсаций скорости (которая является физической величиной) и масштаб турбулентности. Согласно [3, 4, 5], для сверхзвуковых струй эти величины можно принять следующими:
ка = 3 (аиа )2 ,а = 0.01 - 0.05, (1)
где а - интенсивность турбулентного потока.
£а = ЬЯа ,Ь = 0.01 -0.05. (2)
На внешней границе свободной струи, в случае отсутствия искусственной турбулизации потока, кинетическая энергия пульсаций ке и турбулентная вязкость уе принимаются малыми величинами: ке/иа2~10-6, уе/(иа^а)~10-6, а скорости диссипации вычисляются по этим уровням с помощью соотношений (1). Практически
полагается, что турбулентная вязкость на внешней границе струи соответствует ламинарной [3].
Принято: на срезе сопла а=0.02, Ь=0,01, на внешней границе свободной струи:
£е/иа2=10"6, КДВДа)=Ю-6.
3. Экспериментальные данные
Наиболее известные систематические измерения влияния температуры на распределение осевой скорости сверхзвуковых изобарических струй даны в работе [6]. В этих испытаниях рассматривались изобарические горячие струи воздуха при па=Ра/Ре=1, числе Маха на срезе сопла Ма=2 и температуре в камере, которая изменялась в диапазоне 300-1370К.
Для получения общей картины сравнения расчетов с экспериментом, экспериментальные данные работы [6] были обобщены в работе[3] следующими зависимостями:
и/иа=1-ехр(-0.216//-5.84//5), при/<4 , (3)
/=Х/Хст, Хст=ХС300(0.64+0.36Ки), (4)
здесьХст - длина начального участка «горячих» струй; Кн=СреГе/(СраГ0) - энтальпийный фактор;
Сре,Сра -теплоемкости газов внешней среды и на срезе сопла, соответственно;
Те,Т0 - температуры внешней среды и в камере, соответственно.
^е300 - длина начального участка "холодных" струй (Те = Т0), которая
вычисляется с учетом влияния показателя адиабаты на срезе сопла у для горячих
2
струй по критерию ?Ма :
Хс300/Яа=0.44+8.97(?Ма2)а45 , (5)
здесь Яа - радиус среза сопла, иа - скорость на срезе сопла. Характер зависимости длины начального участка струи от температуры приведен в таблице 1.
Таблица 1.
Характер зависимости длины начального участка сверхзвуковой струи от
температуры Т0 при Сре=Сра, Те=300К.
Т0, К 300 600 1000 1500 2500 3000
КН=СреТе/(СраТ0) 1 0.5 0.3 0.2 0.12 0.1
0.64+0.36КН 1 0.82 0.748 0.712 0.6832 0.676
4. Результаты расчетов 4.1. Холодные струи (То=300К)
На рисунках представлены результаты расчётов изобарических сверхзвуковых осесимметричных струй с применением модели турбулентности к-юББТ, в сравнении с зависимостями(1-3) . Принятый критерий, по которому можно судить о длине начального участка струи (Х075 ) - расстояние от среза сопла на котором скорость вдоль оси струи падает на 25% по сравнению со скоростьюна срезе (и/Иа=0.75).
На рисунке 1 приведена зависимость длины начального участка холодной струи, отнесенная к радиусу сопла (Х/Яа), от числа Маха.
Рисунок 1. Длина начального участка изобарической холодной струи
Как видно из графика, модель турбулентности без поправок Вилкокса на сжимаемость занижает дальнобойность струй (длину начального участка), а турбулентная интенсивностьа потока оказывает незначительное влияние на длину начального участка изобарической струи, показатель 0.02 для этого параметра является оптимальным. В дальнейшем, для всех расчетов примем это значение параметра а.
На рисунке 2 по оси абсцисс отложены расстояния от среза сопла, отнесенные к длине начального участка ХСТ, на оси ординат - скорость на оси струи, отнесенная к скорости на срезе сопла. Данные приведены для холодных струй. Экспериментальные данные построены по обобщающей зависимости из [3]:
их/иха=1-ехр (-0.357//-3.5//2); (6)
и/иа
О 0.5 1 1.5 2 2.5 3 3.5 4
■ аппроксимация эксп. данных М= 2 М=3 М=4
Рисунок 2. Распределение скорости по оси холодной струи, модель ББТс поправкой на сжимаемость. Мю=0.25, а=0.02.
При низких числах Маха, расчетные данные хуже согласовываются с экспериментальными - они немного занижают скорость на дальнем расстоянии от среза сопла.
4.2. Горячие струи
Рассмотрим график зависимости длины начального участка горячей струи (Т0>500К), отнесенной к радиусу сопла (Х/Яа), от температуры.
Рисунок 3. Длина начального участка горячей изобарической струи.
Далее, рассмотрим график профиля скорости по оси струи, отнесенный к
начальному участку.
Рисунок 4. Распределение скорости на оси горячей струи. Модель ББТ поправкой на
сжимаемость. Мг0=0.25.
Из рис. 4 и рис. 5 видно, что для расчета горячих струй, стандартная модель турбулентности ББТк-ю не подходит. В этом случае возможны следующие подходы: вводить дополнительные поправки на температуру (например, [7], учитывающие градиент полной температуры) или модифицировать модель, поменяв пороговое значение числа Маха турбулентных пульсаций Мю (для холодных струй Мю=0.25).
Путем проведения методических расчетов были установлены оптимальные значения параметра Мю, которые необходимо задавать для моделирования горячих струй, они сведены в таблице ниже.
Таблица 2.
Рекомендуемое пороговое значение числа Маха турбулентных пульсаций в зависимости от температуры при Ма=2, в сравнении с аппроксимацией (3-5).
То, К 300 600 1000 1500 2500 3000
Мю 0.25 0.26 0.31 0.33 0.35 0.36
Сравним экспериментальную и расчетную (при значениях Мю указанных в таблице) зависимость дальнобойности струи от температуры.
Рисунок 5. Длина начального участка изобарической струи.
Рисунок 6. Распределение скорости на оси струи. Модель ББТ с поправкой на
сжимаемость.
5. Поправки на температуру Абдоль-Хамида.
Еще один способ модифицировать модель турбулентности для расчета горячих струй - ввести поправки на температуру. В стандартной двухпараметрической модели Ментера 8БТ[1] в формуле удельной скорости диссипации:
а
е
в* к
коэффициент в принимается константой в = С = 0.09.
Абдоль-Хамид предположил, что коэффициент Сзависит от градиента полной температуры с помощью следующей зависимости:
С^ = 0.09
1 +
Тг
0.041 (Мг)
где
ч=
УТ
0
( и л к / /Б V У
Т
0
и значение числа Маха турбулентных пульсаций:
М, ="' *%
Ниже представлен график профиля струи для модели 88Тк-ю,с измененным Мю и график профиля струи для модели турбулентности с поправками Абдоль-
Хамида.
*
Рисунок 7.
Из рис. 7 видно, что модель турбулентности с поправками Абдоль-Хамида лучше согласовывается с экспериментальными данными, чем модель ББТк-ю.
6. LES
Идея LES (Large Eddy Simulation) состоит в разрешении турбулентных вихрей, чей размер больше размера ячейки расчетной сетки и моделировании мелких. При такой постановке, необходимо создавать расчетную сетку с такими ячейками, которых достаточно что бы разрешить большую часть энергетического спектра турбулентности.
При расчёте использовалась подсеточная модель вихревой вязкости WALE. В модели WALE, вихревая вязкость описывается соотношением:
,2 ^f u = р2- -у ;
j Р +7 )
где ,s и Sijd в модели WALE, определены, как:
(
Ls = min
kd,CwV/3 I;
w
V
где k - константа Кармена, d - расстояние до стенки.
ed 1 (-2,-2 ^ 1 о 2 dui
SU = 2(giJ + gU J-sV2k' giJ =8X7
В ANSYS Fluente, по умолчанию, задана константа Cw=0.325, как наиболее подходящая для большинства задач. В модели WALE, для ламинарной подвижной среды, турбулентная вязкость приравнивается к нулю. Это позволяет правильно моделировать ламинарные зоны. Модель Смагоринского-Лилли, например, не дает значение турбулентной вязкости, равное нулю. Поэтому модель WALE предпочтительнее.
Расчет методом LES нужно проводить в нестационарной постановке. Все величины при этом пульсируют с течением времени (Рисунок 8). Метод LES имеет потенциал для получения более достоверных результатов, чем методы RANS. Однако, метод LES более требователен к вычислительным ресурсам, чем RANS метод.
Рисунок 8. Расчет методом RANS (верхняя картинка) и методом LES (нижняя картинка, мгновенное поле) горячей струи [6], Ma=2, T0=1370K. Поле числа Маха.
В рассматриваемом случае использовалась расчетная сетка с общим числом ячеек порядка 5 млн. Для экономии вычислительных ресурсов в расчетах не учитывалось течение в сопле с детальным разрешением пограничного слоя, расчеты проводились от среза сопла. В качестве расчетных схем использовались схемы второго порядка аппроксимации по времени и пространству с решением задачи распада-разрыва методом Roe, а также схемой третьего порядка MUSCL и низкодиссипативным методом Roe. Различие в применяемых подходах можно видеть на Рисунке 8 - применение низкодиссипативной схемы приводит к более детальной картине развития процессов смешения.
Roe, 2-order upwind Low-diffusion Roe, 3-order MUSCL
Рисунок 9 - Распределение завихренности в слое смешения сверхзвуковой струи [6] (M=2, T0=1370K) для различных вариантов применяемых численных схем
Сравнение полученных распределений скорости по оси струи в рамках различных методов приведено на Рисунке 10. Видно, что классический метод LESWALE хорошо согласуется с данными эксперимента для холодной струи, но для горячих струй получено существенное занижение дальнобойности струи. На Рисунке 11 представлено сравнение энергии турбулентности по оси струи методами LES и RANS. Видно, что методом LES ядро струи имеет меньшую протяженность, и начало турбулентного смешения начинается раньше. Профиль турбулентной энергии для метода LES более крутой и достигает больших значений, что оказывает влияние на профиль скорости. Отметим также, что для получения осредненной скорости потребовалось порядка 3 величин пролетного времени, в то же время для получения среднеквадратичных значений (как, например, турбулентной энергии) потребовалось 8 величин пролетного времени.
О 20 40 60 XlRa о 10 20 30 40 XJRa
• Safronov • LES • SST standart * SST без поправок «Seiner • SST + Abdol-Hamid SSTMtcHLÎ ♦ SST standart ■ LES
Рисунок 10. Сравнение скорости по оси струи различными моделями RANS и LES для холодной T0=300K, Ma=3 (слева) и горячей T0=1370K, Ma=2 (справа) струй
0.00 20.00 40.00 60.00 00.00 x/R
Рисунок 10. Сравнение турбулентной энергии по оси струи методами RANS и LES
для горячей струи Ma=3, T0=300K
7. Выводы
Путем сравнения результатов расчетов с экспериментальными данными
установлено:
1) Модель k-ю с поправками Вилкокса на сжимаемость существенно завышает дальнобойность сверхзвуковых горячих струй, а без поправок занижает.
2) Устранение этого недостатка для изобарических струй возможно с помощью изменения порогового значения числа Маха турбулентных пульсаций (Mt0), либо, вводом поправок Абдоль-Хамида, учитывающих градиент полной температуры.
3) Даны рекомендации по изменению числа Mt0 в зависимости от температуры струи при числе Маха на срезе сопла Ма=2.
4) Метод LES (WALE) для холодных струй согласуется с экспериментальными данными, а для высокотемпературных струй существенно завышает уровень турбулентности и тем самым занижает дальнобойность.
Работа выполнена при частичной поддержке гранта РФФИ № 14-08-01149.
Библиографический список
1. Menter F. «Two-Equation Eddy-Viscosity Turbulence Models for Engineering Applications.» AIAA Journal, Vol. 32, 1994, pp. 1598-1605.
2. Wilcox D. C. Turbulence Modeling for CFD, DCW Industries, Inc., La Canada, CA 91011, 3rd ed., 2006.
3. Сафронов А. В. О применимости моделей турбулентной вязкости для расчета сверхзвуковых струйных течений // Физико-химическая кинетика в газовой динамике 2012, Т.13, выпуск №1 :http://chemphys.edu.ru/issues/2012-13-1/articles/305/
4. Kenzakowski D.C. Turbulence modeling improvements for jet noise prediction using PIV datasets//AIAA-2004-2978, 10th AIAA/CEAS Aeroacoustics Conference and Exhibit, 10-13 May.- 2004.
5. Кудимов Н. Ф., Сафронов А. В., Третьякова О. Н. Численное моделирование взаимодействия многоблочных сверхзвуковых турбулентных струй с преградой // Электронный журнал «Труды МАИ», 2013, № 70: http://www.mai.ru/science/trudy/published.php?ID=44440 (дата публикации 25.11. 2013)
6. Seiner, J.M., Ponton, M.K., Jansen, B.J., and Lagen, N.T. The Effect of Temperature on Supersonic Jet Noise Emission// DGLR/AIAA Paper 92-02-046.1992.
7. Abdol-Hamid, K. S., Massey, S. J., and Elmiligui, A., "Temperature Corrected Turbulence Model for High Temperature Flow," Journal of Fluids Engineering, Vol. 126, 2004, pp. 844.
