Турбулентная вязкость и крупномасштабная завихренность в моделировании большими вихрями Текст научной статьи по специальности «Физика»

Вычислительные технологии

Том 8, № 1, 2003

ТУРБУЛЕНТНАЯ ВЯЗКОСТЬ И КРУПНОМАСШТАБНАЯ ЗАВИХРЕННОСТЬ В МОДЕЛИРОВАНИИ БОЛЬШИМИ ВИХРЯМИ

А.М. Бллонишников Санкт-Петербургский государственный инженерно-экономический

университет, Россия e-mail: [email protected], [email protected]

The new parametrization of turbulent viscosity coefficient taking into account the large-scale vorticity for turbulent flow incompressible liquid is offered. This parametrization has only one empirical coefficient as well-known Smagorinsky parametrization. The model qualitatively agrees with other multiparameter models, which give a decrease of turbulent viscosity coefficient in the presence of large-scale vorticity.

Введение

Подсеточное моделирование было введено как попытка преодолеть трудности прямого численного моделирования развитой турбулентности на основе исходных уравнений Навье — Стокса, поскольку современным компьютерам не хватает быстродействия для расчета динамики вихрей [1]. В подсеточных моделях выбирают некоторый масштаб h (шаг разностной сетки по пространству) и рассматривают вихри, диаметр которых превосходит h. Из-за нелинейности уравнений Навье — Стокса, описывающих движение несжимаемой жидкости, возникают существенные трудности. В уравнения, описывающие крупные вихри, входят подсеточные напряжения Рейнольдса, которые определяются скоростями мелкомасштабных вихрей, чей диаметр меньше h. Подсеточные напряжения Рейнольдса зависят от масштаба h и градиентов (и/или более высоких производных) скоростей больших вихрей и описывают переход энергии от крупных вихрей к более мелким с последующим переходом в тепло. Масштаб h не должен превышать диаметра энергосодержащих вихрей [2].

В идеальном случае выражение для подсеточных напряжений Рейнольдса желательно выводить аналитически из уравнений для скоростей мелкомасштабных вихрей с последующим осреднением произведения компонент мелкомасштабной скорости, которое определяет напряжения Рейнольдса по объему кубической ячейки с длиной ребра h. В этом случае подсеточное моделирование было бы свободно от эмпирических констант и такую модель можно было бы применять для всех типов течений, включая внутренние течения, свободную турбулентность, ламинарно-турбулентный переход. К настоящему моменту исследователи далеки от решения этой задачи. Различные подходы к подсеточному моделированию отражены, например, в работе [3].

© А. М. Балонишников, 2003.

Простейшей и наиболее используемой параметризацией является следующее выражение Смагоринского для подсеточного коэффициента турбулентной вязкости [4], V = еЬ Б, величина которого определяет напряжения Рейнольдса, где с — некоторая константа; Н — шаг разностной сетки по пространству при дискретизации уравнений гидродинамики; Б — инвариант тензора деформации крупномасштабной скорости. Величина Б вычисляется по формуле

Б = л/2Бц Бцг, Бц = + д3

где Ц (г = 1, 2, 3) — компонента крупномасштабной скорости. В действительности коэффициент турбулентной вязкости может зависеть и от антисимметричной части Пц = 1/2(д Ц — дц Ц) тензора градиента скорости. Имеется простая связь этой антисимметричной части с завихренностью крупномасштабного движения Пг : (П = гоЮ) : Пг = П^, где Е^к — тензор Леви-Чевитта. В работе [5] показано, что для несжимаемой жидкости имеются всего пять независимых инвариантов симметричной и антисимметричной частей тензора градиента деформации крупномасштабной скорости: Тг{82}, Тг{П2}, Тг{83}, Тт{П2Б}, Тг{П282}, где Тг — след матрицы. Получить зависимость коэффициента турбулентной вязкости от этих пяти инвариантов эмпирическим путем весьма затруднительно [6]. Зная коэффициент турбулентной вязкости как функцию производных крупномасштабной скорости, можно определить подсеточные напряжения Рейнольдса и тем самым получить замкнутое описание динамики крупномасштабных вихрей.

1. Вывод выражения для коэффициента подсеточной турбулентной вязкости

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

Предположим, что крупномасштабную скорость и ее тензор градиента можно считать постоянными при рассмотрении динамики мелкомасштабной компоненты скорости [7], т.е. откажемся от широко используемого предположения Праудмена — Бэтчелора [8] о линейности крупномасштабной или средней скорости (что дает вклад в спектр скорости мелкомасштабных вихрей). Тогда уравнение для Фурье-компонент мелкомасштабной скорости будет иметь следующий вид [9]:

(dt + vk2)ua(k,t) = - 2 PaßY (k)J^ Uß (p,t)u7 (k - p,t) + Paß ömUß um(k,t), (1)

p

где Pj = öij — kikj/k2, PaßY = kßPaY + kYPaß; индексы в этих формулах пробегают значения 1, 2, 3 (здесь и всюду по тексту по повторяющимся индексам подразумевается суммирование). Если жидкость или газ для простоты считать несжимаемыми (kiui = 0), то уравнение (1) можно трансформировать к двум независимым компонентам мелкомасштабной скорости, направленным вдоль единичных векторов ej и е2 (j = 1, 2, 3), ортогональных вектору k, как это было сделано для изотропной турбулентности в отсутствие градиентов крупномасштабной скорости dU [10]:

(dt + vk2)uY (k, t) = a'V + ®Yaßua(p,t)uß (q,t), (2)

p,q,p+q=k

где а™ = ejemdmUj; Ф7ав(k, p, q) = —ikmej(k)ea (p)em(q), верхние индексы пробегают значения 1, 2, а нижние — 1, 2, 3. Входящие в эти формулы векторы имеют следующие декартовы компоненты, выраженные через углы Эйлера [10]:

k = (k cos 9 cos n, k sin 9 cos n, k sin n),

где cos 9 = ki/k", sin 9 = кг/к", cos n = k"/k, sin n = k3/k, k" = yk^+kf. Если выбрать единичный вектор e = (cos 9, sin 9, 0), то единичный вектор ei(k) определится соотношением

e1 = e х k/|e х k|

или e1 = (sin 9, — cos 9, 0), а единичный вектор e2(k) — соотношением e2 = k х e1(k)/|k х e1| или

e2 = (cos 9 sin n, sin 9 sin n, — cos n).

При инверсии вектора k справедливы соотношения e1(k) = —e1(—k),e2(k) = e2(—k). Верно соотношение e1 х e2 = k/k. Связь обычных Фурье-компонент мелкомасштабной скорости с ее поляризационными компонентами выражается соотношением

u (k,t) = e?(kK(k,t), (3)

причем

k ■ eM = 0, eM(k) ■ eA(k) = ef(k)e^k) = Pij(k). (4)

С помощью замены переменных u = Bv, где матрица B составлена из собственных векторов матрицы A = {аав}, можно привести матрицу A к диагональной форме (в общем случае к жордановой матрице) [11], а само уравнение для поляризационных компонент мелкомасштабной скорости примет вид

(d + vk2)vA = JAV + (B-1)Ay Y, ф7"в(Bv)a(Bv)e, (5)

где JАм — элемент диагональной матрицы, на диагонали которой стоят собственные числа матрицы A — А1 и Л2, определяющиеся как

P P2

А1,2 = "2 Т — Q, (6)

где P = ek^dmUj + e2^дтUj;

д = (фтадН^ВД) - (ф^т^Х^дЩ).

Если разложить вИ на симметричную и антисимметричную части = +

и ввести единичный вектор п с компонентами щ = кг/к, то, используя свойства векторов еА и условие несжимаемости, можно выражения для Р и д упростить:

Р = -ПгЩБгу , (7)

Q = (e1emSmj)(e?epSip) — (e^Smj)2 + > ■ П)2. (8)

1

41

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

eA:

q=(eme2—e1em)e1 epsmj Sip+i(n ■ ^)2-

Введем антисимметричный тензор тт1 = — £1£т • Выписав члены в Q, отличные от 0 с т = I, после очевидных преобразований получим

Q = (т12#1?#2р + #3р + ^зр)т?р + ^(п ■ П)2. Нетрудно убедиться в том, что т12 = п3, т23 = п1, т13 = — п2. Поэтому Q = п3(#11#22 — #12#21) + П2(5И^33 — #13 #31) + ^2(^22^33 — #23#32) + 2^1^2(^23 #13 —

— #12#33) + 2П1П3(#12#23 — #22#13) + 2П2П3(#12#13 — #11#23) + ^(п ■ П)2. (9)

В инвариантном виде Q = — 1/2Тг[(пх8)2] + 1/4(п-П)2, где Тг[(пхЯ)2] = п^#кр£раЬпа#Ьг. Окончательно собственные числа матрицы линеаризованных уравнений для мелкомасштабной скорости (без вязких членов) получаются следующим образом:

Л1,2 = — 1 Тг(СЯ) ± ^4[Тг(С8)]2 + 1 Тг[(п х Я)2] — 4(п ■ П)2, (10)

учтем, что С г2 = 'Пг'П^ •

Согласно принципу подчинения Хакена [12], неустойчивые моды подчиняют себе устойчивые. Если имеется набор неустойчивых мод, подчиняющихся одному и тому же уравнению, то следует ожидать, что наиболее существенный вклад в подсеточные напряжения Рейнольдса дадут наиболее неустойчивые моды. Если дискриминант квадратного уравнения, определяющий собственные числа Л1 и Л2, положителен, то эти собственные числа действительны. Если оба собственных числа положительны или одно положительно, а другое отрицательно, то, по-видимому, основной вклад дадут моды с наибольшим числом Л1.

Среди направлений вектора к следует выбрать то, которое дает максимальное значение Л1 .

Для нахождения этого максимума по угловым переменным считаем, что предварительно был осуществлен переход с помощью ортогонального преобразования в систему координат, в которой тензор Б диагонален. Тогда выражения (??) и (10), определяющие величины Р и Q, существенно упрощаются:

Р = —п2#1 — п2 #2 — п3#3, (11)

Q = п3#1#2 + п2#1#3 + п?#2#3 + 1(п ■ П)2, (12)

где #1 = #11, #2 = #22, #3 = #33 — собственные числа тензора Б, > #2 > #3. Условие несжимаемости будет иметь вид

#1 + #2 + #3 = 0. (13)

Компоненты единичного вектора п удовлетворяют соотношению

п? + п2 + п3 = 1. (14)

Используя соотношения (??)-(??), выражение (10) для величины Л1 запишем как

Л1 = — 1[(2п1 + п2 — 1)#1 + (2п2 + п1 — 1)#2] + Z, (15)

где _

£ = ^[М + П - ад + (2п2 + П?)#2]2 + + 2#1) + + 2#2) - ЗД - Я,

Я = (П1П1 + П2^2 + ПзПэ)2, Пз = 1 - П2 - п2. (16)

Необходимыми условиями максимума функции А1(п) являются: дЛ1/дп1 = 0, дЛ1/дп2 = 0. Используя соотношение (15), получим систему двух алгебраических уравнений:

(2#1 + #2- (2#1 + ^ММ + ^2 - ВД + (2п2 + П2 - 1)#2] - П1#2(2#1 + #2) +

+ (П + дПз Пз )(П1^1 + П2П2 + П3П3) = 0, (17)

дп1

(2#2 + #1 )П2^ - (2#2 + #1)П2[(2П1 + п2 - 1)#1 + (2п2 + п2 - 1)#2] - П2#1 (2#2 + #1) +

+ (П + дпэ Пз )(щП1 + П2П2 + П3П3) = 0, (18)

дП2

где дп3/дп1 = ^п1/у/Г-, дп3/дп2 = 1 - П - П.

Система уравнений (??) и (??) не становится проще при переходе к угловым переменным Однако при П = 0 можно получить решение Л1тах = 2#1. В других случаях необходимо использовать приближенные численные методы, например метод Ньютона. Таким образом, основным инвариантом будет I = А1тах, если подкоренное выражение

в (15) положительно и в нем максимум берется по различным направлениям вектора к. Если подкоренное выражение отрицательно, то I = 8ир(Ке(Л1(п))} при вир(Ке(Л1)) > 0

или I = 0 при 8ир{Ке(Л1)} < 0, где Б,е — действительная часть числа. В этом случае возникновение неустойчивостей возможно лишь за счет нелинейных взаимодействий, а диссипацию будут обеспечивать вязкие члены. Физически первый случай (подкоренное выражение в (15) положительно) соответствует превалированию деформации крупномасштабным полем скорости мелкомасштабных вихрей над эффектом вращения. Второй случай (подкоренное выражение в (15) равно 0 или отрицательно) соответствует обратной ситуации. Ячейка со стороной к находится в ядре когерентного вихря [13], где использовался более грубый метод идентификации когерентных структур: инвариант д тензора dU положителен (д = 1/2(ПуП^ - Б^Sji)).

В итоге можно считать, что коэффициент подсеточной турбулентной вязкости равен

ск2\//, (19)

где с — константа, которую следует определить путем сравнения численных экспериментов с натурными экспериментами; к — шаг разностной сетки.

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

Отметим, что уменьшение коэффициента турбулентной вязкости при увеличении крупномасштабной завихренности — явление достаточно известное. Так, в [14] предложена

крупномасшабная феноменологическая модель с учетом наличия крупномасштабной завихренности с коэффициентом турбулентной вязкости:

V = ^ (к2/е)[1/(1 + 0.36кП/е)],

где V) — константа; к — кинетическая энергия турбулентности; е — удельная скорость диссипации турбулентной энергии; П = \fWijWji — модуль крупномасштабной завихренности; Шц — антисимметричная часть тензора градиента скорости dU. Похожее выражение использовалось в [15] для коэффициента турбулентной вязкости

V = см/мкть

где / — коэффициент, зависящий от числа Рейнольдса; т1 = к/е + ^/е — характерное время оборота энергосодержащих вихрей; V — кинематический коэффициент молекулярной вязкости;

1

A1 + AsTi max(S, П)

А1 = 8, А5 = -\Z3cos Ф, Ф = 1агссо8(Лх), х = 21'5 ^^^.

3 о

В работах [14, 15] V ^ 0, как 1/П, при П ^ то. В нашем подходе, как следует из (10), V пропорционален аО + VЬ2О2 — П2, если П < ЬО, где а, Ь — константы; при П > ЬО V пропорционален О, т. е. коэффициент турбулентной вязкости в первом приближении перестает зависеть от завихренности. Какая ситуация более адекватна действительности, в настоящий момент сказать трудно. Величины к, е могут, кроме того, непосредственно зависеть от крупномасштабной завихренности П через свои уравнения переноса.

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

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

Автор благодарен А. В. Копыльцову за критические замечания.

Список литературы

[1] PlOMELLI U. Large-eddy simulation: achievements and challenges // Prog. Aerosp. Sci. 1999. No. 35. P. 335-350.

[2] BALONISHNIKOV A. M. Extended local balance model of turbulence and Couette-Taylor

flow // Phys. Review E. 2000. Vol. 61, No. 2. P. 1390-1394.

[3] Voelkl T., Pullin D. I., Chang D. C. A physical-space version of the stretched-vortex subgrid-stress model for large-eddy simulation // Phys. of Fluids. 2000. Vol. 12, No. 7. P. 1810-1825.

[4] Smagorinsky J. General circulation experiments with the primitive equations // Mon. Weather Rev. 1963. No. 91. P. 99-105.

[5] Pope S. B. A more general effective viscosity hypothesis //J. Fluid. Mech. 1975. No. 72. P. 331-340.

[6] JORGEN T., Gatski T. B. General explicit algebraic stress relations and best approximation for three-dimensional flows // Intern. J. Eng. Sci. 1998. No. 36. P. 739-763.

[7] Скворцов Г. е., Тимохов Л. А. К теории турбулентности // Вестник ЛГУ. 1980. Вып. 3, №13. С. 106-110.

[8] Batchelor G. K., Proudman J. The effect of rapid distorsion on a fluid in turbulent motion // Q. J. Mech. Appl. Math. 1954. No. 7. P. 83-90.

[9] БАлонишников А.М. Однопараметрическая модель неоднородной гидродинамической турбулентности // Системы автоматизации в науке и производстве. Л.: Наука, 1984. С. 50-54.

[10] Lee J. The triad-interaction representation of homogeneous turbulence //J. Math. Phys. 1975 No. 7. P. 1359-1366.

[11] Брюно А. Д. Локальный метод нелинейного анализа дифференциальных уравнений. М.: Наука, 1979.

[12] Хакен Г. Синергетика. М.: Наука, 1980.

[13] Dubief Y., Delcayre F. On coherent-vortex identification in turbulence // J. of Turbulence. 2000. No. 1. P. 2-22 (http://jot.iop.org).

[14] Zhou Y. A phenomenological treatment of rotating turbulence // Phys. of Fluids. 1995. No. 7. P. 2094-2106.

[15] Merci B., De Langhe C., Vierendeels J., Dick E. A quasi-Realizable Cubic Low-Reynolds Eddy-Viscosity Turbulence Model with a New Dissipation Rate Equation // Flow, Turbulence and Combustion. 2001. No. 64. P. 133-160.

Поступила в редакцию 16 июля 2002 г.