Расчет затухания однородной и изотропной турбулентности с помощью метода моделирования больших вихрей (les) второго порядка аппроксимации Текст научной статьи по специальности «Математика»

УДК 532.517.4

И. С. Босняков, Г. Г. Судаков

Центральный аэрогидродинамический институт им. проф. H. Е. Жуковского

Расчет затухания однородной и изотропной турбулентности с помощью метода моделирования больших вихрей (LES) второго порядка аппроксимации

Исследована задача о затухании однородной и изотропной турбулентности. Для неё существуют теоретические оценки закона затухания турбулентных пульсаций со временем и аналитические зависимости для корреляционного тензора скоростей. Эти оценки используются для верификации промышленного численного метода второго порядка точности по пространству и времени. Задача о диссипации турбулентности решается с помощью метода моделирования больших вихрей. Полученные результаты попадают в общий диапазон результатов, полученных другими авторами.

Ключевые слова: однородная изотропная турбулентность, моделирование больших вихрей, вырождение турбулентности, численное моделирование.

1. Введение

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

В статье на примере решения задачи о диссипации однородной и изотропной турбулентности рассмотрено описание динамических и энергетических свойств турбулентных возмущений с помощью стандартного метода моделирования крупных вихрей LES (Large Eddy Simulation).

В работе особое внимание уделяется двум вопросам. Первый вопрос посвящён методу генерирования начального турбулентного поля. Идея метода описана в работах [5, 12, 13, 14] и заключается в том, что турбулентное поле скоростей представляется как суперпозиция мод Фурье, удовлетворяющих уравнению неразрывности, со случайными амплитудами. В настоящей работе подробно описан алгоритм, реализующий основные идеи этих подходов.

Второй вопрос касается моделирования изотропной турбулентности с помощью LES. Обсуждается вопрос о применимости к данному классу задач промышленных солверов со схемами второго порядка аппроксимации по пространству. В более ранних работах данная проблема исследовалась, как правило, с помощью специализированных исследовательских программ. Имеющаяся статистика по данному направлению представлена в работе [8]. Известны также работы, например [6], где численные методы используются как способ изучения турбулентности и численной проверки теоретических выводов. Некоторые другие работы используют теоретические соображения для настройки численного метода [6, 7, 9].

Важной проблемой для рассматриваемой задачи является влияние схемной вязкости на характер развития течения. Показано, что схемы низкого порядка аппроксимации имеют большую схемную вязкость. Но и схемы высокого порядка [8] не решают этой проблемы, так как диссипация остается достаточно высокой для больших волновых чисел.

За последнее время произошёл резкий рост производительности вычислений в связи с широким распространением компактных суперЭВМ, поэтому разумным выходом из сложившейся ситуации может оказаться использование более подробных расчётных сеток и стандартных методов второго порядка аппроксимации.

2. Постановка задачи

Рассмотрим течение несжимаемой вязкой жидкости. Пусть скорость жидкости представляется в виде

V = {V) + V',

где величина в скобках есть среднее значение скорости, а штрихованная величина — ее пульсационная составляющая. По определению (V') = 0.

Пусть

«=

Величина д называется уровнем турбулентности. Турбулентная энергия определяется следующим образом:

' Л = ^

к = 2 •

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

йк

£ = - т •

При больших числах Рейнольдса в турбулентном потоке присутствуют пульсации с масштабами от очень больших до очень малых. На больших масштабах порядок изменения средней скорости равен величине пульсаций. Этот размер обозначается через Ь* и называется масштабом турбулентности. В соответствии с моделью Ричардсона (см. [14]) турбулентная энергия переходит от пульсаций с большими масштабами к пульсациям с мелкими масштабами практически без диссипации. Другими словами, имеется поток энергии от крупно- к мелкомасштабным пульсациям. Далее этот поток диссипирует в тепло на самых мелких масштабах г ~ Л (масштаб диссипации), где существенную роль играет вязкость жидкости.

Колмогоров и Обухов предложили модель однородной и изотропной турбулентности на основе теории подобия (см. [14]). В соответствии с предложением Ричардсона единственной размерной константой, определяющей свойства турбулентного движения на масштабах Л << г << Ь^ является скорость диссипации турбулентной энергии е.

Вводится понятие волнового числа пульсаций О ~ 1/г, и пусть Е (О) (Ю. есть кинетическая энергия единицы массы жидкости, заключенная в пульсациях со значением волновых чисел [О, О + ^О]. Величина Е (О) называется спектральной плотностью энергии. Тогда на основании теории подобия и размерностей может быть получен закон Колмогорова-Обухова:

Е(О) = се2/3О-5/3, 1/и < О < те, с и 1,5. (1)

Здесь с — эмпирическая постоянная. Закон Колмогорова-Обухова не описывает поведение турбулентности на масштабах, больших Ь*, однако на основании анализа экспериментальных данных был получен ряд других законов распределения, описывающих турбулентность на больших масштабах. Наиболее часто используются распределения Драйдена и Кармана.

Распределение Драйдена имеет вид

8 2г (О^)4

Е(О) = ■

1 + (О!*)

2 3"

Распределение Кармана имеет вид

*(О) =

55 2 (1,339О£*)4

1 + (1, 339О£*)2

17/6 '

Легко видеть, что распределение Кармана удовлетворяет закону Колмогорова-Обухова при & ^ ж, в то время как распределение Драйдена не удовлетворяет этому закону. На рис. 1 приведены графики этих распределений.

Следуя [14, 15], можно ввести корреляционные тензоры:

Ьк = (Уг(гг)ук(г2)), Вгк = ((Уг(г\) - Ук(Г2)) (ь^Г 1) - Ьк(Г2))) , 2 2

Вгк = 3Ч Ьгк - 2Ъгк, 1,к = 1, 2, 3.

Не останавливаясь на деталях вычисления, которые находятся в [15], можно выписать следующий результат:

В,(г) = 4 / ЕК) (l - *

_ -

0

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

Вц(г) = С ■ г2/3, (2)

где С & 7,2е2/3.

Полученное поведение корреляционного тензора Вц (г) справедливо только в инерционном интервале, когда инерционные силы преобладают над вязкими силами. В задаче о затухании однородной и изотропной турбулентности такое поведение может наблюдаться, когда время t не слишком велико.

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

q2 = const ■ (t - to)-p , (3)

где to — начало отсчета времени, р — положительное число.

Теоретическая оценка степени затухания пульсаций зависит от предположений, накладываемых на свойства турбулентности. Предположим, что интеграл Лойцянского (см. [15])

в любой момент времени конечен и отличен от нуля (М. Д. Миллионщиков, 1939 г.). При больших временах t , когда мелкомасштабные турбулентные пульсации уже затухли и остались только большие вихри, подверженные ламинарной диссипации, есть теоретическая оценка р = 5/2 (вязкий интервал).

Имеются оценки для показателя степени в инерционном интервале времени. Оценка Колмогорова основывается на существовании и постоянстве интеграла Лойцянского. В этом случае р = 10/7. Такой же показатель степени будет, если принять гипотезу Кармана об автомодельности турбулентности на определенных масштабах и принять числа Рейнольдса большими.

Существует оценка, основанная на инварианте Биргхофа [4]. Полагая, что при k ^ 0 Е(к) ~ k2, получается значение р = 6/5 (большие длины волн). Исследования с помощью теории ренормализации групп Яхотом и Орзагом [10] дают степень р = 1.47. Оценка Яхота [11], приведенная для случая, когда вязкость v ^ 0 и рассматриваемое время t << L2/v ^ ж, даёт оценку р = 10/7.

Проводились также экспериментальные исследования затухания турбулентности. Для случая с малыми числами Рейнольдса эксперименты Бэтчелора и Таунсэнда [1] подтверждают предположение Лойцянского и дают р ~ 5/2 для вязкого интервала времени. В то же время представленные в работе [2] экспериментальные результаты, которые соответствуют инерционному интервалу, имеют разброс для величины р ~ 1.0^1.37.

зультаты в диапазоне р ~ 1.17^2.0 (см. обзор в [8]).

Таким образом, возникает задача оценки скорости затухания однородной и изотропной турбулентности для больших времен (3) как в инерционном, так и на вязком интервалах времени, а также оценки для корреляционного тензора (2) в инерционном интервале.

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

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

3. Метод генерации турбулентного поля

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

Энергетическая функция спектральной плотности для этой модели имеет вид [3]

27^ [1 + (а^О)2]17/6'

Где £ — масштаб, а2 = (и'2 + и'2 + и>'2') — среднеквадратическое отклонение, О = (О х, Оу, Ог) — вектор пространственной частоты, а^ ~ 1.339 — постоянная Кармана. Волновое число Оо, соответствующее масштабу турбулентности Ь, определяется из вышеприведенной формулы следующим образом:

ОоЬ = 1.

Функция Е(О) задаёт величину компонент случайного спектрального вектора скорости порыва ветра. Корреляция и дисперсия этих компонент выражается корреляционным

тензором Sik (Q). Он определяется по формуле [14]:

Е(П) {Q2ôik - QQk)

Sik (Q) =

4^

Q4

Следует подчеркнуть, что в О-прострапстве корреляция между частотами отсутствует (следствие однородности и изотропности турбулентности), но имеются корреляции между фурье-компонентами скорости в отдельно взятой точке О-прострапства (следствие уравнения неразрывности).

Наиболее общим способом задания реализаций является представление скоростей ветра в виде рядов Фурье по дискретному набору частот:

М К N

Шу (х,у,г) = ^ ^ ^Кгкп соё(2^-(Охтх+ОукУ+Огпг))+В3ткп ёт(2к{Пхтх+Пуку+Охпг)),

т=1к=1п=1

где А3ткп) — независимые по индексам т, к, п гауссовские случайные величины,

описываемые корреляционной матрицей (],к = 1, 2, 3):

M

дз \k

mkn mkn

= M

R3 r k

Sjk (Qxm^ Qyk ,

В качестве примера на рис. 2 представлено поле ш-компоненты скорости для одной из реализаций.

Рис. 2. w-компонента скорости. Параметры турбулентности: q = 1 м/с, L = 20 м

4. Численный метод

В настоящей работе для расчета затухания однородной и изотропной турбулентности использован метод конечных объемов второго порядка аппроксимации как по пространству, так и по времени (промышленный пакет ANSYS CFX). Расчет производился для случая несжимаемой вязкой жидкости. В качестве модели турбулентности использована стандартная модель LES Смагоринского с константой Смагоринского Cs = 0.1.

5. Результаты расчета

Рассмотрим физическую область 0 < xi < 10 м, i = 1, 2, 3, с разбиением 200 х 200 х 200. Ей сопоставим область волнового пространства с размерами 1/10 < ki < 1/0.05, г = 1, 2, 3, с

разбиением на 200 ячеек в каждом направлении. В качестве граничных условий использованы граничные условия: циклические и тина Opening, с заданным давлением на участках втекания жидкости. В качестве начальных условий задается ноле скоростей, сгенерированное по алгоритму п. 3 с параметрами турбулентности q = 1 м/с, L = 20 м.

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

Для регистрации турбулентных пульсаций и определения статистических свойств течения были выбраны 10 контрольных точек в середине области, расположенных вдоль оси х с шагом Ах = 0.1 м. В качестве среды был выбран воздух постоянной плотности при 25 °С.

Далее производится расчет затухания турбулентности и оценка уровня турбулентности q и корреляционного тензора Вц(г) в каждый момент времени t. Средние величины при этом понимаются как средние (в окрестности момента времени t) в промежутке времени At, малого по сравнению с характерным временем затухания Т ~ 300 с, но большим по сравнению с характерным временем пульсаций скорости. Для данной задачи был выбран интервал At = 20 с. Следует напомнить, что процесс затухания турбулентных пульсаций скорости для данной задачи происходит в два этапа: первый этап инерционный интервал t < 35 с. Второй этап — вязкое затухание t > 35 с. Закон затухания (t — t0)-p справедлив как в инерционном, так и на вязком интервале, когда инерционные силы малы по сравнению с вязкими силами, но на каждом интервале имеют место разные показатели степени. Это утверждение можно проиллюстрировать расчетными данными, приведенными на рис. 3 6. На рис. 5 данные рис. 4 приведены в логарифмических переменных, что позволяет определить степень в законе затухания турбулентных пульсаций для вязкого интервала времени.

1 ------------------:------------------;------------------1------------------1-------------------:------------------:•'

0.5 - (■: ...;

1 _i_i_i_i_i_i

О 50 100 150 200 250 300

ВД

Рис. 3. Мгновенные показания датчиков скорости

Поскольку под сеточная вязкость моделирует именно инерционные силы, степень затухания в численном методе на вязком интервале не должна зависеть от константы Смагоринского. Этот факт подтверждается расчетами при разных значениях константы Смагоринского: результаты расчетов для значений констант Смагоринского Cs = 0.1 (стандартная) и Cs = 0.25 дают очень близкие результаты для закона затухания на больших временах. Прямая линия на рис. 5 построена методом наименьших квадратов. Старший коэффициент определяет степень р в законе затухания турбулентности. По результатам расчетов получено значение для р = 2.6, близкое к теоретическому значению р = 5/2, предполагающему сохранение интеграла Лойцянского. Близкий результат дает также эксперимент [1|.

Инерционный интервал

Вязкая диффузия

1 —

О 50 100 150 200 250 3Û0

t(C]

Рис. 4. Уровень турбулентных пульсаций скорости по времени: 1 инерционный интервал 0 < t < 35 с, 2 — вязкая диффузия t > 35 с

-data

I

V?

i i

■°g(t)

Рис. 5. Затухание турбулентных пульсаций скорости по времени в логарифмических переменных: р/2 = 1.3 для стадии вязкого затухания

о

Ж

m о

- -2 -3

у = - 0j.51*x - (j.49 —i_i_ -data 1 —linear-

е—--- _ _

________1_______ i _______L. ^ ■ ^ ы i i

1 1.2 1.4 1.6 1.В 2 2.2 2.4 2.6 2.8 3

iog(t)

Рис. 6. Затухание турбулентных пульсаций скорости по времени в логарифмических переменных: р/2 = 0.51 для инерционной фазы

На рис. 6 аналогичные данные приведены для инерционного интервала времени. Здесь расчетная оценка для показателя степени р = 1 в законе затухания соответствует экспериментальным значениям р ~ 1 ^ 1.37, как указано в п. 2, но несколько ниже теоретических оценок р ~ 1.2^1.47.

На рис. 7 дано сравнение зависимости корреляционного тензора Вц (г) от расстояния при £ = 9.4 с внутри инерционного интервала времени с теоретическими значениями (2).

6. Обсуждение результатов

Таким образом, численные результаты моделирования затухания однородной и изотропной турбулентности, полученные в настоящей работе, достаточно хорошо коррелируют с теоретическими оценками характеристик данного течения как в инерционном, так и на вязком интервале времени. Показано также, что численное моделирование указывает на справедливость стсисннох'о закона затухания турбулентной энергии с показателем степени р = 2.6, близким к теоретическому значению р = 5/2, для вязкого интервала времени, и р = 1 — дЛЯ инерционного интервала времени, что также удовлетворительно согласуется с теоретическими (р ~ 1.2^1.47) и экспериментальными (р ~ 1^1.37) оценками. Численные расчеты данной задачи разных авторов с помощью LES для показателя р также дают результаты в диапазоне р ~ 1.17^2.0 (см. обзор в [8]). Обращает на себя внимание большой разброс для оценок показателей степени, полученных с помощью расчетных методов. Такой разброс можно получить, если не делить весь временной интервал на две зоны и предположить, что степенной закон с единым показателем степени справедлив на всем временном интервале. Тогда для показателя степени получится некоторое среднее значение, зависящее от величины временного интервала и заключенное между предельными значениями 1 < р < 2.5. Для временного интервала Т = 300 с, принятого в настоящей задаче такая оценка дает величину р = 1.44.

Следует отметить также, что если на вязком интервале показатель степени в законе затухания не зависит от константы Смагоринского, то на инерционном интервале такая зависимость есть. Пересчет данной задачи с увеличенной константой Смагоринского Cs = 0.25 вместо стандартной Cs = 0.1 дает р = 1.4 для величины показателя степени в инерционном интервале вместо р = 1.0.

Работа была выполнена при поддержке гранта РФФИ 13-08-00346 и при финансовой поддержке Министерства образования и науки в рамках договора № 700013728 от 21.11.2012 «Разработка моделирующего комплекса реалистичного восприятия оператором (летчиком) сложных режимов полета и оценки его психофизиологического состояния» по 218 постановлению правительства РФ.

Литература

1. Batchelor G.K., Townsend A.A. Decay of turbulence in the final period // Proc. Roy. Soc.

- 1948. - A194. - N 1039. - P. 527-543.

2. Comte-Bellot G., Corrsin S. The use of a contraction to improve the isotropv of grid-generated turbulence //J. Fluid Mech. - 1966. - V. 25. - N 4. - P. 657-682.

3. Etkin B. Dynamics of Atmospheric Flight. — New York: John Wiley k, Sons Inc., 1972.

4. Birkhof G. Fourier synthesis of homogeneous turbulence // Commun. Pure Appl. Math. —

1954. _ у. 7. _ p. и) п.

5. Kraichnan R.H. Diffusion by a random velocity field // Physics of Fluids. — 1970. — V. 13, N 1.

6. Wang L., Chen S., Brasseur J.G., Wyngaard J.C. Examination of hypotheses in the Kolmogorov refined turbulence theory through high-resolution simulations. Part 1. Velocity field // J. Fluid Mech. - 1996. - V. 309. - P. 113-156.

7. Muschinski A. A similarity theory of locally homogeneous and isotropic turbulence generated by a Smagorinskv-tvpe LES //J. Fluid Mech. — 1996. — V. 325. — P. 239 260.

8. Thornber В., Drikakis D. Large-eddv simulation of isotropic homogeneous decaying turbulence 11 ECCOMAS CFD, TU Delft. - 2006.

9. Thornber В., Mosedale A., Drikakis D. On the implicit large eddy simulation of homogeneous decaying turbulence //J. Comput. Phvs. — 2007. — V. 226. — P. 1902-1929 (doi:10.1016/j.jcp.2007.06.030).

10. Yakhot V., Orszag S.A. Renormalization group analysis of turbulence // Phvs. Rev. Lett.

- 1986. - V. 57. - P. 1722-1725.

11. Yakhot V. Decay of three-dimensional turbulence at high Reynolds numbers // J. Fluid Mech. - 2004. - V. 505. - P. 87-91.

12. Wang S.-T., Frost W. Atmospheric turbulence simulation techniques with application to flight analysis 11 NASA CR 3309. - 1980.

13. Адамьяп Д.Ю., Стрелец M.X., Травин А.К. Эффективный метод генерации синтетической турбулентности на входных границах LES области в рамках комбинированных RANS-LES подходов к расчету турбулентных течений // Математическое моделирование. - 2011. - Т. 23, № 7. - С. 3-19.

14. Ландау Л.Д., Лифшиц Е.М. Теоретическая физика. Т. 6. Гидродинамика. — М.: Наука, 1988.

15. Монин А.С., Яглом A.M. Статистическая гидромеханика. Ч. 2. Механика турбулентности. — М. : Наука, 1967.

Поступила в редакцию 12.11.2013.