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

А. М. Балонишников, А. В. Копыльцов

ГЛАВНЫЙ ИНВАРИАНТ ДЛЯ ПОДСЕТОЧНЫХ МОДЕЛЕЙ ТУРБУЛЕНТНОСТИ

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

Введение

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

настоящему моменту исследователи далеки от решения этой задачи. Различные подходы к подсеточному моделированию отражены, например, в работе [3]. Простейшей и наиболее используемой параметризацией является следующее выражение Смагоринского для подсеточного коэффициента турбулентной вязкости vt [4], v = ch2S, величина которого определяет напряжения Рейнольдса, где с — некоторая константа, h — шаг разностной сетки по пространству при дискретизации уравнений гидродинамики, S — инвариант тензора деформации крупномасштабной скорости:

S = /," • S,= WU, +д,и,).

где U (i = 1,2,3) — компонента крупномасштабной скорости. В действительности коэффициент турбулентной вязкости может зависеть и от антисимметричной части Q, = -2(5U, — дjUi) тензора градиента скорости. Имеется простая связь этой антисимметричной части с завихренностью крупномасштабного движения Qi (Q = rotU): Qi = 2sijkQjk, где g,,k — тензор

Леви—Чевитта. В работе [5] показано, что всего для несжимаемой жидкости имеется пять независимых инвариантов симметричной и антисимметричной части тензора градиента деформации крупномасштабной скорости:

Tr{S2}, Tr{Q2}, Tr{S3}, Tr{Q2S}, Tr{Q2S2},

где Tr — след матрицы. Получить зависимость коэффициента турбулентной вязкости от этих пяти инвариантов эмпирическим путем весьма затруднительно [6].

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

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

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

(д, + v 2 )иа (k^ t) = — iPa№ (к и в (p^ t )ur (к — t) + Ра/Зд mU рит (к t )• (1)

Р

где Р, = 8, — кк/ к2, PaY = крРаг + кгРар, индексы в этих формулах пробегают значения 1, 2, 3 (здесь и всюду по тексту по повторяющимся индексам подра-

зумевается суммирование). Если жидкость или газ для простоты считать несжимаемыми ktUi = 0 , то можно трансформировать уравнение (1) к двум независимым компонентам мелкомасштабной скорости, направленным вдоль единичных векторов s1 и s. (j = 1, 2, 3), ортогональным вектору k, как это было

сделано для изотропной турбулентности в отсутствие градиентов крупномасштабной скорости dU [10]:

(dt + V 2K(k, t) = amu ™ + £Фаиа(р, t)u ß(q, t), (2)

p,q, p+q=k

где a™ = sr.sßmdmU., Фга (k,p,q) = -ikmsr. (k)s“(p)sm (q), верхние индексы пробегают значения 1, 2, а нижние — значения 1, 2, 3.

Входящие в эти формулы вектора имеют следующие декартовы компоненты, выраженные через углы Эйлера [10]:

k = (kcos#cosn, ksin#cosn, ksinn),

где cos# = kj/k”, sin# = k2/k”, cosn = k”/k, sinn = k3/k,k” =y¡k¡ + k22.

Если выбрать единичный вектор e = (cos#,sin#,0), тогда единичный вектор s1 (k) определяется соотношением

s1 = eхk/|exk |,

или

s1 = (sin#,-cos#,0), а единичный вектор s2 (k) определяется соотношением

s2 = k x s'(k)/1 k x s1 |,

или

s 2 = (cos#sinn, sin#sinn,-cosn).

При инверсии вектора k справедливы соотношения

sl(k) = -sl(-k), s2(k) = s2(-k).

Верно соотношение s1 xs2 = k /k . Связь обычных Фурье-компонент мелкомасштабной скорости с ее поляризационными компонентами выражается соотношением

u.(k,t) = s™(k)u™(k,t), (3)

причем

k s™ = 0, s™ (k) -s\k) = 5^, s™ (k )s™ (k) = p(k). (4)

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

(dt +V >д = J + (В -уг^ФГСФ(Ву)а(ВуУ, (5)

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

(6)

Если разложить dU на симметричную и антисимметричную части дmUJ = + 0т и ввести единичный вектор п с компонентами п1 = к1 / к, то,

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

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

Окончательно имеем следующее выражение для собственных чисел матрицы линеаризированных уравнений для мелкомасштабной скорости (без вязких членов):

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

используя свойства векторов е и условие несжимаемости, можно выражения для Р и Q упростить:

(8)

(7)

могательные вектора ел :

Нетрудно убедиться, что т12 = п3, т23 = п1, т13 = —п2. Поэтому

■■■¡11 .и. С иИ о.. .

1]к J кр раЬ а т

Л,2 =— -2 Тг (С8) ±4±[Тг (С8)]2 + ± Тг[(п х 8)2] — |(п • О)2 , (10)

где Сц = n1nJ.

ющий собственные числа Л и Я2, положителен, то эти собственные числа действительны. Если оба собственных числа положительны, или одно положительно, а другое отрицательно, то, по-видимому, основной вклад дадут моды с наибольшим числом Л .

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

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

Р = -«251 - «252 - «3253, (11)

Q = «32 5152 + «22 5153 + «2 52 53 + -4(п - О)2, (12)

где 51 = 511, 52 = 522, 53 = 533 — собственные числа тензора 8, 51 > 52 > 53. Условие несжимаемости будет иметь вид:

51 + 52 + 53 = 0. (13)

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

«2 + «2 + «32 = 1. (14)

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

Л=- -2 [(2«2 + «2 -1)5! + (2«2 + «1 -1)52 ] + Ъ, (15)

где

Ъ = д/-4[(2«2 + «2 -1)51 + (2«2 + «2)52]2 + [«5 (52 + 251) + «251 (51 + 252) - 5152 ] - Я ,

я = («1о1 + «2о2 + «3о3)2, «3 = ±д/1 -«2-«2. (16)

Необходимыми условиями максимума функции Л(п) являются

Л = 0; ^ = 0.

д« д«2

Используя соотношение (15), получим систему двух алгебраических уравнений:

(251 + 52 )п1г - (251 + 52)«[(2«2 + «22 -1)51 + (2«2 + «2 -1)52 ] - «152 (251 + 52) +

д«

+ (О, + О 3)(«1О1 + «2 О 2 + «3О 3) = 0, (17)

д«1

(252 + 51 )«22 - (252 + 51 )«2 [(2«1 + «22 -1)51 + (2«22 + «I -1)52 ] - «251 (252 + 51) +

д«

+ (О2 + —^ О3)(«1О1 + «2О2 + «3О3) = 0, (18)

д«2

д«3 / /1 2 ~ «3 / /1 2 3”

где —- = + «1/л/1 - «1 - «2, —- = + «2 / л/1 - «1 - «2 .

д«1 д«2

Система уравнений (17) и (18) не становится проще при переходе к угловым переменным Г), @ . Однако в случае О = 0 можно получить решение Л тах = 251. В других случаях необходимо использовать приближенные численные методы, например метод Ньютона.

Таким образом, соответствующий основной инвариант будет I = Лтах, если подкоренное выражение в (15) положительно, где максимум берется по различным направлениям вектора к. Если подкоренное выражение отрицательно, то I = 8ир{Яе(Л (п))} при 8ир(Яе(Л)) > 0 или I = 0 при 8ир{Яе(Л)} < 0 , где Яе — действительная часть числа. В этом случае возникновение неустойчивостей возможно лишь за счет нелинейных взаимодействий, а диссипацию будут обеспечивать вязкие члены. Физически первый случай (подкоренное выражение в (15) положительно) соответствует превалированию деформации крупномасштабным полем скорости мелкомасштабных вихрей над эффектом вращения. Второй случай (подкоренное выражение в (15) равно 0 или отрицательно) соответствует обратной ситуации. Во втором случае ячейка со стороной Ь находится в ядре когерентного вихря [13], где использовался более грубый метод индентификации когерентных структур: инвариант q тензора ^ положителен, где q = у (О^.О^ - 5^5 ).

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

к, = ек2 л/7, (19)

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

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

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

К = к0 (к2 /е)[1/(1 + 0,36Ю/е)],

где к0 — константа, к — кинетическая энергия турбулентности, е — удельная скорость диссипации турбулентной энергии, О = — модуль круп-

номасштабной завихренности, Ж1у — антисимметричная часть тензора градиента скорости dU. Похожее выражение использовалось в [15] для коэффициента турбулентной вязкости у,:

к = ем/М,

где fß — коэффициент, зависящий от числа Рейнольдса, т1 — характерное

время оборота энергосодержащих вихрей: т1 = k/s + Vv/s, v — кинематический коэффициент молекулярной вязкости,

1

с =---------------------,

м A1 + Ат1 max(S, Q)

где A1 = 8, As =-J3 cos Ф, Ф = -jarccos(V6x), x = 215 ,J j к .

S

В этих двух работах vt ^ 0 как 1/ Q при Q^œ. В нашем подходе, как следует из (10), vt пропорционален aS + ->/ b2 S2 -Q2 , если Q < bS , где a, b — константы; при Q > bS vt пропорционален S, т. е. коэффициент турбулентной вязкости в первом приближении перестает зависеть от завихренности.

Какая ситуация более адекватна действительности, в настоящий момент сказать трудно. Величины к, s могут, кроме того, непосредственно зависеть от крупномасштабной завихренности Q через свои уравнения переноса.

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

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

БИБЛИОГРАФИЧЕСКИЕ ССЫЛКИ

1. Piomelli U. Large-eddy simulation: Achievements and challenges // Prog. Aerosp. Sci. 1999. Vol. 35. Р. 335.

2. Balonishnikov A. M. Extended local balance model of turbulence and Couette-Taylor flow // Physical Review E. 2000. Vol. 61. № 2. Р. 1390-1394.

3. Voelkl T., Pullin D. I., Chang D. C. A physical-space version of the stretched-vortex sub-grid-stress model for large-eddy simulation // Physics of Fluids. 2000. Vol. 12. № 7. Р. 1825.

4. Smagorinsky J. General circulation experiments with the primitive equations // Mon. Weather Rev. 1963. Vol. 91. Р. 99-105.

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

6. Jorgen T., Gatski T. B. General explicit algebraic stress relations and best approximation for three-dimensional flows // Int. J. Eng. Science. 1998. Vol. 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. Vol. 7. P. 83-90.

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

10. Lee J. The triad-interaction representation of homogeneous turbulence // J. Math. Phys. 1975. № 7. Р. 1359-1366.

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

12. Хакен Г. Синергетика. — М., 1980.

13. Dubief Y., Delcayre F. On coherent-vortex identification in turbulence // Journal of Turbulence. 2000. Vol. 1. P. 2-22.

14. Zhou Y. A phenomenological treatment of rotating turbulence // Physics of Fluids. 1995. Vol. 7. Р. 2094.

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. Vol. 64. P. 133-160.

А. Balonishnikov, A. Kopyltsov THE MAIN INVARIANT FOR SUB-GRID MODELS OF TURBULENCE

A new expression for the sub-grid coefficient of turbulent viscosity based on a simplified analysis of equations for small-scale components of the incompressible turbulent flow and the dimensional analysis has been obtained. This expression in contrast with the Smagorinsky model is determined by an invariant of the tensorgradient of the large-scale velocity, which contains explicitly the large-scale vortices, which diminishes the coefficient of turbulent viscosity.