CC BY
20
УДК 517.958:533.7
ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ЛАМИНАРНО-ТУРБУЛЕНТНОГО ПЕРЕХОДА В ТЕЧЕНИИ ЗА ОБРАТНЫМ УСТУПОМ
Т. Г. Елизарова*^, П. Н. Никольский
(.кафедра математики) E-mail: telizar@afrodita.phys.msu.ru
На примере задачи о течении вязжого сжимаемого газа в жанале с внезапным расширением демонстрируется возможность численного моделирования ламинар-но-турбулентного перехода в рамжах жвазигазодинамичесжих уравнений.
Введение
В связи с развитием вычислительной техники перспективным подходом к численному исследованию турбулентных течений становится метод прямого численного моделирования, основанный на уравнениях Эйлера или Навье-Стокса без использования полуэмпирических моделей турбулентности [1]. При этом роль турбулентной диссипации на масштабах, меньших шага пространственной сетки, выполняют искусственные регуляризаторы. Диееипа-тивный механизм, присущий этим регуляризаторам, выполняет роль фильтра, который сглаживает под-сеточные пульсации газодинамических величин. Однако детальное исследование диееипативных характеристик искусственных регуляризаторов, используемых в нелинейных и нестационарных численных моделях, является сложной математической задачей и, как правило, не проводится.
В настоящей работе для прямого численного моделирования турбулентного течения вязкого сжимаемого газа впервые используются квазигазодинамические (КГД) уравнения [2-4]. Упомянутые уравнения отличаются от уравнений Навье-Стокса и Эйлера дополнительными дивергентными слагаемыми. Их появление связано с введением дополнительного по сравнению с уравнениями Навье-Стокса сглаживания по времени при определении газодинамических величин — плотности, скорости и температуры среды. Эти дополнительные слагаемые выполняют роль регуляризаторов и обеспечивают устойчивость и точность численных алгоритмов, построенных на их основе. Для КГД уравнений доказана теорема о неубывании полной термодинамической энтропии, что подтверждает диееипативный характер соответствующих регуляризаторов.
Численный алгоритм расчета дозвуковых течений вязкого газа, основанный на КГД уравнениях, изложен в [5]. На его основе проводится численное моделирование течения в плоском двумерном канале
с внезапным расширением, где с увеличением числа Рейнольдса наблюдается переход от ламинарного течения к турбулентному. Результаты расчетов сопоставлены с данными эксперимента [6], который для указанной задачи считается классическим и до настоящего времени используется для тестирования численных алгоритмов (см., напр., [7, 8] и цит. лит.). В указанном эксперименте воспроизведены условия, при которых течение за уступом с достаточной точностью можно рассматривать как двумерное.
Постановка задачи и метод численного
решения
Система КГД уравнений в обычных обозначениях имеет вид
др
д(ри)
at
■diV/m =0,
■ div(jm <g> и) + Vp = div П,
(1) (2)
d_ at
p(y + £) +div i™(y
dlvq ■
= div(II •«).
(3)
Здесь вектор плотности потока массы, тензор вязких напряжений и вектор теплового потока вычисляются как
jm = p(u-w), w = -[div(pu®u) + Vp].
Р
П = II A S + TU (
г/
q= —xVT — три
p(u-V)u + Vp
(и • V)p + тр divu (и- V)e+p(u-
(4)
(5)
(6)
Щга — тензор вязких напряжений Навье-Стокса, т — параметр временного сглаживания, который
Институт математического моделирования РАН.
имеет размерность времени и определяется в виде
V р Sc
h
а-
(7)
где г/ = щ{Т/То)ш — коэффициент динамической вязкости, Бс — число Шмидта, /г — минимальный пространственный масштаб, на котором разрешается структура течения. В данной задаче это шаг пространственной сетки, с= \fyRT — скорость звука, 0 <а < 1 — численный коэффициент. Величина /г/с соответствует времени, за которое возмущения пересекают разностную ячейку.
Рассматривается плоское двумерное течение вязкого газа в канале длины и высоты 2Н с обратным уступом высоты Н. Во входном сечении задан профиль Пуазейля со средней по сечению скоростью 11о- Течение характеризуется безразмерными числами Рейнольдса и Маха
Re =
P0U0H
Ма =
ип
(8)
% с0
где ро и 7q — плотность и температура газа во входном сечении канала. Твердые стенки канала предполагаются адиабатическими с условиями прилипания и непротекания для скорости, на выходной границе ставятся условия сноса потока. В качестве начального условия используется течение Пуазейля. Численный алгоритм, который представляет собой явную по времени разностную схему с аппроксимацией второго порядка точности для всех пространственных производных, детально описан в работе [5].
В расчетах использовалась равномерная пространственная сетка по обоим направлениям. Шаг по времени выбирался из условия устойчивости вида dt = ßh/c, где коэффициент 0 < ß < 1 подбирался в процессе вычислений.
Результаты численных экспериментов
Расчеты проведены для дозвукового течения воздушного потока при нормальном давлении. Молекулярные параметры составляют 7= 1.4, Рг = 0.737, ш = 0.74, Sc = 0.746.
В таблице приведены основные характеристики расчетов при Ма = 0.1. Здесь L — длина расчетной области, Tf — безразмерное время окончания расчета. На рис. 1 рассчитанные картины течений сопоставлены с экспериментальными данными [6] по длине отрывной зоны Ls. Согласно [6] длина отрывной зоны почти линейно увеличивается с ростом числа Рейнольдса*^ до Reo = 1200, а затем уменьшается в области 1200 < Reo < 6600. Уменьшение длины отрывной зоны связывается с возникновением нестационарности в течении за уступом.
Число Рейнольдса в [6] вычисляется по формуле (8), где вместо высоты уступа Н выбрано значение = 2Н.
Параметры и результаты расчетов
Расчеты из работы [5]
Ls 18 17 16 15 14 13 12 11 10 9 8 7 6 5
Re h а ß dt L Ls [6] Ls Tf
100* 0.05 0.5 0.2 0.0010 10 5.0 5.0 81
200* 0.05 0.5 0.2 0.0010 14 8.5 8.4 204
300 0.05 0.5 0.3 0.0015 20 11.3 10.4 300
400* 0.05 0.5 0.2 0.0010 20 14.2 12.7 626
600 0.05 0.3 0.5 0.0025 30 17.0 16.0 400
0.5 0.3 0.0015 15.5
700 0.05 0.3 0.5 0.0025 30 15.5 15.7 400
0.5 0.3 0.0015 17.0
1000 0.03 0.3 0.5 0.0015 30 13.8 10.6 400
0.5 0.3 0.0009 11.4
2000 0.05 0.3 0.5 0.0025 30 9.0 9.1 400
0.5 0.3 0.0015 10.9
2800 0.05 0.3 0.5 0.0025 30 6.1 7.2 400
0.5 0.3 0.0015 8.8
3500 0.05 0.3 0.5 0.0025 30 8.0 7.8 400
0.5 0.3 0.0015 9.0
0
Рис. 1. Сопоставление рассчитанной длины отрывной зоны (пунктир, маркерами обозначены проведенные расчеты) с данными эксперимента [5] (сплошная кривая)
Для ламинарных течений (Ие <600) возвратное течение за уступом является стационарным, и длина отрывной зоны хорошо соответствует данным эксперимента и известным расчетам. При Re > 600 в поле течения возникают осцилляции и оно приобретает нестационарный, турбулентный характер. В этом случае длина отрывной зоны определяется по осредненному полю скорости. Видно уменьшение размера отрывной зоны с ростом Re и достаточно хорошее соответствие расчета данным эксперимента. Заметим, что немонотонный характер
500 1000 1500 2000 2500 3000 3500 Re
а = 0.5
•О
8 ВМУ. Физика. Астрономия. № 4
уменьшения Ь8 в эксперименте может быть связан с возникновением в течении трехмерных вихревых структур, появление которых в рассматриваемом численном расчете невозможно. Согласно многочисленным экспериментальным данным и численным расчетам ([7, 8] и библиографии к ним), при числах Ке > 4000 размер отрывной зоны практически постоянен и составляет 5 < Ь8 < 8 в зависимости от геометрии задачи. В частности, согласно [8], Ь8 = 5.5 при Ке = 3700.
На рис. 2 и 3 представлены эволюции во времени продольной компоненты скорости для Ке = 300 и Ке = 600 в точке с координатами (8,1). Лами-
их
0 50 100 150 200 250 t
Рис. 2. Временная эволюция скорости для Не = 300
их
1.6 г
1.4 h
0.2
0 ..........................................
0 50 100 150 200 250 300 350 400 t
Рис. 3. Временная эволюция скорости для Re = 600
нарное течение при Ке = 300 в процессе установления выходит на стационарный режим. Р£е = 600 соответствует переходному режиму, для которого видно зарождение осцилляций скорости. На рис. 4 приведен фрагмент временной эволюции горизонтальной компоненты скорости для турбулентного режима течения (1£е = 1000). Здесь результирующее течение является существенно нестационарным. Размер отрывной зоны определяется после осреднения поля скоростей на интервале времени в пределах 300 < At < 400. Осредненное течение мало зависит от интервала осреднения при At>5. На рис. 5 и 6 приведены мгновенная и осредненная картины течения для Ке = 1000. Эти картины качественно соответствуют результатам [8], полученным для Ке = 3700.
их
1.0 г
0.8 -
0.6 -
0.4
0 _1_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I
300 320 340 360 380 400 t
Рис. 4. Фрагмент временной эволюции скорости для Не = 1000
В отличие от эксперимента [6], где число Маха составляет ~ 0.001, в данном расчете Ма = 0.1. Однако эффекты, связанные со сжимаемостью среды, имеют порядок 0(Ма2) и их влияние на размер отрывной зоны можно считать незначительным. Для ламинарных течений это предположение подтверждается результатами [5]. В этой же работе показана слабая зависимость длины отрывной зоны от шага сетки Н. В численных экспериментах показано, что для турбулентных течений величина Ь8 также мало зависит от Н при Н < 0.1.
При вычислении времени сглаживания т (7) имеется один свободный параметр — численный коэффициент 0 < а < 1. Согласно практике предыдущих [3, 5, 9] и выполненных в настоящей работе расчетов для ламинарных течений размер отрывной зоны мало зависит от величины а. При расчетах
5 10 15
Рис. 5. Мгновенные линии тока, Re = 1000, фрагмент
20 х
^тш и 1 1 1 1 I I I I I I I I I I I | |
их -0.4 -0.2 0 0.2 0.2 0.2 0.2 1.0 1.2 1.4
5 10 15
Рис. 6. Линии тока для осредненного течения, Re = 1000, фрагмент
20 х
турбулентных режимов наблюдается зависимость от а.
Можно ожидать, что с уменьшением шага сетки влияние параметра а на осредненное течение будет уменьшаться, так как все большая часть энергетического спектра и вихревых структур будет разрешаться пространственной сеткой [1].
Близкие результаты по составлению с экспериментом размеров отрывной зоны были получены при моделировании ламинарных и турбулентных течений вязкой несжимаемой жидкости за обратным уступом в рамках квазигидродинамических уравнений [9].
Таким образом, в настоящей работе продемонстрирована возможность использования КГД алгоритма для прямого численного моделирования ла-минарно-турбулентного перехода в течении за обратным уступом. В отличие от известных методов, в которых для расчета ламинарных и турбулентных течений используются разные численные модели (см., напр., [1, 6-8]), КГД алгоритм позволяет единообразно рассчитывать оба типа течений и получить переход от ламинарного режима к турбулентному при числах Рейнольдса, соответствующих экспериментальным данным.
Несмотря на ограничения, связанные с применением двумерной пространственной модели, результаты численного расчета хорошо соответствуют данным эксперимента. Использованный численный алгоритм естественным образом обобщается на пространственный случай и может быть эффективно
реализован на многопроцессорных вычислительных системах, что позволит в дальнейшем моделировать турбулентные течения в трехмерной постановке.
Литература
1. Белоцерковский О.М., Опарин A.M. Численный эксперимент в турбулентности. М., 2001.
2. Шеретов Ю.В. Математическое моделирование течений жидкости и газа на основе квазигидродинамических и квазигазодинамических уравнений. Тверь, 2000.
3. Елизарова Т.Г., Соколова М.Е., Шеретов Ю.В. // Журн. вычисл. матем. и матем. физики. 2005. 45, № 3. С. 545.
4. Елизарова Т.Г. Квазигазодинамические уравнения и методы расчета вязких течений. М., 2007.
5. Елизарова Т.Г., Соколова М.Е. // Вестн. Моск. ун-та. Физ. Астрон. 2005. № 3. С. 6 (Moscow University Phys. Bull. 2005. N 3. P. 6).
6. Armaly B.F., Durst F., Pereira J.C.F., Schonung В. // J. Fluid Mech. 1983. 127. P. 473.
7. Le H., Moin P., Kim J. 11 J. Fluid Mech. 1997. 330. P. 349.
8. Wee D., Tongxun Yi., Annaswamy A., Ghoniem A. // Phys. Fluids. 2004. 16, N 9. P. 3361.
9. Елизарова Т.Г., Калачинская И.С., Шеретов Ю.В., Шилъников E.B. // Труды факультета ВМиК МГУ. Прикладная математика и информатика. 2003. № 14. С. 85.
Поступила в редакцию 07.06.06
9 ВМУ. Физика. Астрономия. № 4
