close

Вход

Забыли?

вход по аккаунту

код для вставкиСкачать
Заметки о теории турбулентности Колмогорова
А.С. Киселев
Московская область, Химки; [email protected]
Интерес к изучению турбулентных потоков жидкостей и газов
возник у меня в конце 30-х годов. Мне сразу стало ясно, что
основным математическим аппаратом исследований призвана стать
теория случайных функций многих переменных (случайных полей),
которая в то время только зарождалась. Кроме того, вскоре мне
стало ясно, что трудно надеяться на создание замкнутой в себе
чистой теории. За отсутствием такой теории придется опираться на
гипотезы, получаемые из обработки экспериментальных данных.
А.Н. Колмогоров [1, стр. 421]
Аннотация. Рассмотрены следствия из исходных и альтернативных гипотез теории
развитой турбулентности Колмогорова; определены универсальные функции и их
взаимозависимость; предложены модели структурных функций и функции плотности
вероятности завихренности; представлен способ подсеточного моделирования для LES
(Large Eddy Simulation).
Ключевые слова. Турбулентность, теория Колмогорова, случайные поля,
изотропность, универсальные функции, скорость диссипации, завихренность,
перемежаемость, LES.
1. Теория турбулентности Колмогорова и её следствия
Исходная теория Колмогорова о структуре развитой турбулентности [2] включает
гипотезы:
К1) при достаточно больших числах Рейнольдса статистика случайного поля
w (r ) = u(r + x0 , t ) - u( x0 , t ) локально (в универсальном диапазоне r < rIE ) однородна (не
зависит от x0 и t ) и изотропна (инвариантна относительно вращений и зеркальных
отражений). Здесь и далее u – скорость, x и t – координаты и время, rIE – характерный
масштаб, разделяющий универсальный и энергетический диапазоны турбулентности.
К2) Статистика в диссипативном диапазоне r < rDI однозначно определяется n и
e , в инерционном диапазоне rDI < r < rIE – только e , где rDI – характерный масштаб,
разделяющий
диссипативный
и
инерционный
диапазоны
турбулентности;
n
–
коэффициент кинематической вязкости; e – скорость диссипации; скобки используются
для обозначения статистического осреднения.
2
1.1. Частично изотропное поле скорости
Далее, для простоты, термин «изотропное поле» используется для случайного поля,
чья статистика однородна и инвариантна относительно вращений и зеркальных
отражений; термин «частично изотропное поле» – если статистика однородна и
инвариантна только относительно вращений. Теория Колмогорова рассмотрена в этом
параграфе в расширенной постановке, а именно, в предположении о частичной
инвариантности поля скорости w (r ) = u(r + x0 , t ) - u( x0 , t ) ; следствия будут определены
для статистических моментов второго порядка.
В общем случае, моменты второго порядка поля скорости w (r ) могут быть
записаны в виде трехточечных структурных функций
Dij (r ¢, r ¢¢) º wi (r ¢) w j (r ¢¢) = ( ui (r ¢ + x0 , t ) - ui ( x0 , t ) ) ( u j (r ¢¢ + x0 , t ) - u j ( x0 , t ) ) ,
(1.1)
для которых Dij (r ¢, r ¢¢) = D ji (r ¢¢, r ¢) . Первоначально рассмотрим двухточечные моменты
Dij (r ) º wi (r ) w j (r ) = ( ui ( x ¢¢, t ) - ui ( x ¢, t ) ) ( u j ( x ¢¢, t ) - u j ( x ¢, t ) ) ,
(1.2)
которые определены для произвольных точек x ¢ и x ¢¢ с расстоянием между ними
r º x ¢¢ - x ¢ . Очевидно, что Dij (r ) = Dij (r , r ) и Dij (r ) = D ji (r ) . Выберем систему координат,
в которой r = re1 = (r , 0, 0) , где ei – орты координат. Из инвариантности Dij (r ) = Dij (r% )
относительно вращения вокруг оси e1 на 180° ( x%1 = x1 , x%2 = - x2 , x%3 = - x3 ; w% 1 = w1 ,
w% 2 = - w2 , w% 3 = - w3 ) находим
D12 (r ) = D13 (r ) = D21 (r ) = D31 (r ) = 0 ,
(1.3)
где волнистой чертой обозначены векторы в новых координатах после преобразования.
Например, равенство D12 (r ) = 0 получаем из D12 (r ) = D12 (r% ) = w% 1w% 2 = - w1w2 = - D12 (r ) .
Из инвариантности относительно вращения вокруг оси e1 на 90° ( x%1 = x1 , x%2 = x3 , x%3 = - x2 )
следует
D23 (r ) = D32 (r ) = 0 , D22 (r ) = D33 (r ) .
Таким
образом,
есть
только
две
независимые
компоненты
(1.4)
DNN (r ) = D22 (r )
и
DLL (r ) = D11 (r ) . Для произвольной системы координат, используя вращение вокруг оси
r ´ e1 на угол j = arccos(re1 / r ) , перейдем к координатам, в которых r = (r , 0, 0) и,
следовательно, все компоненты структурных функций известны. Общий вид моментов
3
Dij (r ) = DNN (r )dij +
DLL (r ) - DNN (r )
ri rj .
r2
(1.5)
выводится из уравнений преобразований моментов при этом вращении. Свертки (1.5)
Dii (r ) = DLL (r ) + 2 DNN (r ) ;
¶ 2 Dij (r )
¶ri ¶rj
¶Dij (r )
¶ri
=
ü
1 ì dr 2 DLL (r )
- 2rDNN (r ) ý rj ;
3 í
r î
dr
þ
ü
1 d ì dr 2 DLL (r )
- 2rDNN (r ) ý .
í
2
r dr î
dr
þ
=
(1.6)
В инерционном диапазоне rDI < r < rIE структурные функции не зависят от
коэффициента вязкости, и из теории размерностей следует, что
Dii (r ) = CD e
2/3
r 2/3 ,
(1.7)
где CD – универсальная константа. Из теории размерностей также следует, что
характерное время вихрей с характерным размером r – tr = e
характерное
время,
турбулентности, t IE = e
разделяющее
-1/3
инерционный
и
-1/3
r 2/3 , в частности,
энергетический
диапазоны
rIE2/3 . В безразмерной форме связь масштабов длины и времени
имеет вид: tr / th = r 2/3 , где
h º n3/ 4 e
-1/4
, th º n1/ 2 e
-1/ 2
(1.8)
колмогоровские микромасштабы длины и времени и
r º r /h.
(1.9)
Далее, из гипотезы К2 следует, что функция
f D (r , CD ) º ( e n )
-1/2
Dii (r )
(1.10)
универсальна в диапазоне r < rIE , в инерционном диапазоне
f D (r , CD ) = CD r 2/3 .
(1.11)
Дополнительно определим
R 2 (r ) º
1 ¶ 2 Dii (r )
1 d 2 dDii (r )
= 2
r
dr
2 ¶rj ¶rj
2r dr
(1.12)
и универсальную функцию
f R (r , CD ) º e
в инерционном диапазоне
-1
nR 2 ( r ) =
1 d 2 df D (r , CD )
r
,
2r 2 dr
dr
(1.13)
4
f R (r , CD ) = (5 / 9)CD r -4/3 .
(1.14)
Заменив разности
un (r ¢¢ + x0 , t ) - un (r ¢ + x0 , t ) = ( un (r ¢¢ + x0 , t ) - un ( x0 , t ) ) - ( un (r ¢ + x0 , t ) - un ( x0 , t ) )
(1.15)
в выражении
Dij (r ¢¢ - r ¢) = ( ui (r ¢¢ + x0 , t ) - ui (r ¢ + x0 , t ) ) ( u j (r ¢¢ + x0 , t ) - u j (r ¢ + x0 , t ) ) ,
(1.16)
находим
Dij (r ¢, r ¢¢) + Dij (r ¢¢, r ¢) = Dij (r ¢) + Dij (r ¢¢) - Dij (r ¢¢ - r ¢) .
(1.17)
Рассматривается только несжимаемая жидкость, в этом случае поле скорости –
соленоидальное поле: ¶w j (r ) / ¶rj = 0 , и из (1.1) следует, что
¶Dij (r ¢, r ¢¢) ¶Dij (r ¢, r ¢¢)
=
= 0.
¶ri¢
¶rj¢¢
(1.18)
Принимая во внимание последнее уравнение (1.6), определим функцию
F(r ) º
¶ 2 Dij (r )
¶ri ¶rj
=
ü
1 d ì dr 2 DLL (r )
- 2rDNN (r ) ý ,
í
2
r dr î
dr
þ
(1.19)
отсюда
r
r dDLL (r )
1
DNN (r ) - DLL (r ) = - ò F(r )r 2 dr .
2 dr
2r 0
(1.20)
Дифференцируя (1.17) с использованием (1.18), получим
¶ 2 Dij (r ¢, r ¢¢)
¶ 2 Dij (r ¢¢ - r ¢) ¶ 2 Dij (r )
==
= F(r ) ,
¶rj¢¶ri¢¢
¶rj¢¶ri¢¢
¶rj ¶ri
(1.21)
следовательно,
¶ ( ui (r ¢ + x0 , t ) - ui ( x0 , t ) ) ¶ ( u j (r ¢¢ + x0 , t ) - u j ( x0 , t ) )
¶ui ( x ¢, t ) ¶u j ( x ¢¢, t )
=
¶x¢j
¶xi¢¢
¶rj¢
¶ri¢¢
=
¶ 2 wi (r ¢) w j (r ¢¢)
¶rj¢¶ri¢¢
¶ 2 Dij (r ¢, r ¢¢)
=
= F (r ) ,
¶rj¢¶ri¢¢
(1.22)
где r ¢ = x¢ - x0 , r ¢¢ = x ¢¢ - x0 , r = x ¢¢ - x ¢ = r ¢¢ - r ¢ , а x0 – произвольная точка. Запишем
свертку (1.2) в виде
Dii (r ) = ui ( x¢, t )ui ( x¢, t ) + ui ( x¢¢, t )ui ( x ¢¢, t ) - 2 ui ( x¢, t )ui ( x¢¢, t ) ,
отсюда,
(1.23)
5
¶ui ( x¢, t ) ¶ui ( x¢¢, t )
1 ¶ 2 Dii (r ) 1 ¶ 2 Dii (r )
==
= R 2 (r ) .
¶x¢j
¶x¢¢j
2 ¶x¢j ¶x¢¢j
2 ¶rj ¶rj
(1.24)
Завихренность определяется как
æ ¶u ¶u ¶u ¶u ¶u ¶u ö
ω º Ñ´u = ç 3 - 2 ; 1 - 3 ; 2 - 1 ÷ ,
è ¶x2 ¶x3 ¶x3 ¶x1 ¶x1 ¶x2 ø
(1.25)
или
wi = eikm
¶uk
,
¶xm
(1.26)
где ненулевые компоненты
e132 = 1 , e123 = -1 , e 213 = 1 , e 231 = -1 , e321 = 1 , e312 = -1 .
(1.27)
Используя равенство
wi ( x¢, t )wi ( x¢¢, t ) =
¶ui ( x ¢, t ) ¶ui ( x¢¢, t ) ¶ui ( x¢, t ) ¶u j ( x¢¢, t )
¶x¢j
¶x¢¢j
¶x¢j
¶xi¢¢
(1.28)
и уравнения (1.22, 1.24), определим функцию
W 2 (r ) = wi ( x ¢, t )wi ( x ¢¢, t ) = R 2 (r ) - F (r ) , r = x ¢¢ - x ¢ ,
которая
характеризует
распределение
интенсивности
завихренности
(1.29)
в
каскаде
Ричардсона. Тензор скорости деформации
1 æ ¶u ( x, t ) ¶u j ( x, t ) ö
Sij ( x, t ) º ç i
+
÷.
2 çè ¶x j
¶xi ÷ø
(1.30)
Используя равенство
2 Sij ( x¢, t ) S ji ( x¢¢, t ) =
¶ui ( x¢, t ) ¶ui ( x¢¢, t ) ¶ui ( x¢, t ) ¶u j ( x¢¢, t )
+
¶x¢j
¶x¢¢j
¶x¢j
¶xi¢¢
(1.31)
и уравнения (1.22, 1.24), определим функцию
S 2 (r ) = 2 Sij ( x ¢, t ) S ji ( x ¢¢, t ) = R 2 (r ) + F(r ) , r = x ¢¢ - x ¢ ,
(1.32)
характеризующую распределение интенсивности скорости деформации. Уравнения (1.29,
1.32) можно представить в виде
S 2 (r ) - W 2 (r ) = 2F(r ) , S 2 (r ) + W 2 (r ) = 2 R 2 (r ) .
(1.33)
Из второго уравнения (1.33) следует, что в универсальном диапазоне R 2 (r ) > 0 , поскольку
S 2 (r ) > 0 и W 2 ( r ) > 0 .
6
1.2. Изотропное поле скорости
Теория Колмогорова рассматривается здесь в исходной интерпретации, а именно, в
предположении, что поле скорости w (r ) = u(r + x0 , t ) - u( x0 , t ) изотропно. Инвариантность
относительно зеркальных отражений не использовалось до сих пор, рассмотрим, что дает
это условие. Произведем преобразование координат
x% = r ¢ + r ¢¢ - x
(1.34)
которое реализуется сдвигом x (1) = x - (r ¢ + r ¢¢) / 2 , тремя зеркальными отражениями (по
одному отражению для каждой оси) x (2) = - x (1) и сдвигом x% = x (2) + (r ¢ + r ¢¢) / 2 . В новой
системе
координат
w% (r% ) = - w (r ¢ + r ¢¢ - r ) ,
инвариантность
относительно
этого
преобразованию дает
Dij (r ¢, r ¢¢) = w% i (r%¢) w% j (r%¢¢) = wi (r ¢¢) w j (r ¢) = Dij (r ¢¢, r ¢) ,
(1.35)
Dij (r ¢, r ¢¢) = (1/ 2) { Dij (r ¢) + Dij (r ¢¢) - Dij (r ¢¢ - r ¢)} .
(1.36)
и, из (1.17),
Следовательно, во-первых, достаточно использовать двухточечные статистические
моменты (1.2), моменты (1.1) находятся из (1.36). Во-вторых, дифференцируя (1.36) и
принимая во внимание (1.18) получим
2
2
¶ 2 Dij (r ¢, r ¢¢)
1 ¶ Dij (r ) 1 ¶ Dij (r )
==
= 0 , r = r ¢¢ - r ¢ ,
¶ri¢¶rj¢¢
2 ¶ri¢¶rj¢¢
2 ¶ri ¶rj
(1.37)
и, из (1.19),
F(r ) = 0 ,
(1.38)
из (1.20),
DNN (r ) = DLL (r ) +
r dDLL (r )
.
2 dr
(1.39)
То есть, для изотропного поля скорости уравнения из предыдущего параграфа (1.1-1.33)
должны использоваться при F(r ) = 0 . В частности, приведем важные следствия
изотропности:
¶Dij (r )
¶ri
=0;
¶ui ( x ¢, t ) ¶u j ( x ¢¢, t )
= 0 ; S 2 (r ) = W 2 (r ) = R 2 (r ) .
¢
¢¢
¶x j
¶xi
Полезные зависимости находим из первого уравнения (1.6) и (1.39):
(1.40)
7
Dii (r ) =
D (r ) - DLL (r )
1 dr 3 DLL (r )
; DNN (r ) = ii
.
2
r
dr
2
(1.41)
Дополнительно определим универсальные функции, которые характеризуют
распределения интенсивности завихренности
f W (r , CD ) º e
-1
nW 2 (r )
(1.42)
-1
nS 2 ( r ) .
(1.43)
и интенсивности скорости деформации
f S (r , CD ) º e
в каскаде Ричардсона.
Подчеркнем, что для изотропного поля скорости реализуется вырождение
R 2 (r ) = W 2 (r ) = S 2 (r ) , f R (r , CD ) = f W (r , CD ) = f S (r , CD ) .
(1.44)
1.2.1. Одноточечные корреляции
Используя (1.36), получим выражения для моментов производных скорости
¶un ( x ¢, t ) ¶un ( x ¢¢, t )
¶wn (r ¢) ¶wn (r ¢¢)
¶ 2 Dnn (r ¢, r ¢¢)
=
=
¶xm¢
¶xm¢¢
¶rm¢
¶rm¢¢
¶rm¢ ¶rm¢¢
1 ¶ 2 Dnn (r ) 1 ¶ 2 Dnn (r )
==
,
2 ¶rm¢ ¶rm¢¢
2 ¶rm2
(1.45)
где w (r ) = u( x + r , t ) - u( x , t ) , r ¢ = x ¢ - x0 , r ¢¢ = x¢¢ - x0 , r = r ¢¢ - r ¢ ; x0 – произвольная точка.
Здесь отсутствует суммирование по повторяющимся индексам n и m. Переходя в систему
координат, в которой вектор r параллелен оси em , получим
¶un ( x ¢, t ) ¶un ( x ¢¢, t )
1 ¶ 2 Dnn (r ) d nm d 2 DLL (r ) (1 - d nm ) d 2 DNN (r )
=
=
+
.
¶xm¢
¶xm¢¢
2 ¶rm2
2
dr 2
2
dr 2
(1.46)
DLL (r ) = (1/ 2) ( d 2 DLL (0) / dr 2 ) r 2 ; DNN (r ) = (1/ 2) ( d 2 DNN (0) / dr 2 ) r 2 ,
(1.47)
При r << 1
и, из (1.39), находим
d 2 DNN (0)
d 2 DLL (0)
=
2
.
dr 2
dr 2
Отсюда получаем выражения для одноточечных корреляций
(1.48)
8
æ ¶un ( x, t ) ö
ç
÷
è ¶xm ø
2
=
d nm d 2 DLL (0) (1 - dnm ) d 2 DNN (0) (2 - d nm ) d 2 DLL (0)
+
=
.
2
dr 2
2
dr 2
2
dr 2
(1.49)
По определению, скорость диссипации
e( x, t ) º nSij ( x, t ) S ji ( x, t ) ,
(1.50)
следовательно,
æ n d 2 dDii (r ) ö
r
e º nSij ( x , t ) S ji ( x , t ) = nS 2 (0) = lim ç 2
÷
r ® 0 2 r dr
dr ø
è
=
3n d 2 Dii (0) 15n d 2 DLL (0)
=
.
2 dr 2
2
dr 2
(1.51)
Окончательно, для изотропного поля скорости
2
d 2 Dii (0)
d 2 DNN (0)
d 2 DLL (0) 4 e
5
10
=
=
=
dr 2
dr 2
dr 2
3 n
(1.52)
и
æ ¶un ( x , t ) ö
ç
÷
è ¶xm ø
2
=
(2 - d nm ) e
.
15
n
(1.53)
1.2.2. Моменты завихренности второго порядка
Используя (1.1, 1.36), получим
wi ( x ¢, t )w j ( x ¢¢, t ) = eikl e jmn
= eikl e jmn
=-
¶ ( uk (r ¢ + x0 , t ) - uk ( x0 , t ) ) ¶ ( um (r ¢¢ + x0 , t ) - um ( x0 , t ) )
¶rl¢
¶rn¢¢
¶ 2 wk (r ¢) wm (r ¢¢)
¶ 2 Dkm (r ¢, r ¢¢)
= eikl e jmn
¶rl¢¶rn¢¢
¶rl¢¶rn¢¢
eikl e jmn ¶ 2 Dkm (r ) eikl e jmn ¶ 2 Dkm (r )
,
=
2
¶rl¢¶rn¢¢
2
¶rl ¶rn
(1.54)
где r ¢ = x¢ - x0 , r ¢¢ = x ¢¢ - x0 , r = r ¢¢ - r ¢ . Отсюда,
Wij2 (r ) º wi ( x , t )w j ( x + r , t ) =
eikl e jmn ¶ 2 Dkm (r )
.
2
¶rl ¶rn
(1.55)
Подставив Dkm (r ) из (1.5), находим
Wij2 (r ) = W 2NN (r )dij + éëW 2LL (r ) - W 2NN (r ) ùû
W 2NN (r ) = W 2LL (r ) +
r d W 2LL (r )
,
2 dr
ri rj
r2
,
(1.56)
(1.57)
9
и
W 2 (r ) º Wii2 (r ) = R 2 (r ) =
1 d 2 dDii (r )
r
.
2r 2 dr
dr
(1.58)
Окончательно, моменты завихренности второго порядка определяются (1.56, 1.58) и
r
1
1 dDii (r )
W 2 (r ) - W 2LL (r )
W (r ) = 3 ò W 2 (r )r 2 dr =
; W 2NN (r ) =
.
r 0
2r dr
2
2
LL
(1.59)
Подчеркнем, что формулы для моментов завихренности второго порядка (1.56-1.57)
идентичны
уравнениям,
получаемым
из
предположения,
что
случайное
поле
завихренности изотропно.
1.2.3. Осреднение по объёму сферы
Для сферы радиусом r , объемом Vr = (4p / 3)r 3 , поверхностью sr = 4pr 2 и
нормалью к поверхности n определим осреднение завихренности по объёму сферы
1
wr ,i º
Vr
Vr
s
eikm r
ò0 wi ( x, t )dx = Vr ò0 uk ( x, t )nm dsr .
(1.60)
Функция
W 2r º wr ,i wr ,i
(1.61)
характеризует интенсивность вихрей, локализованных в сфере, на основании этого
определим универсальную функцию
f r (r , CD ) º e
=
n
e V2
-1
nW 2r =
V V
n
e Vr2
2
ò ò W (| x¢¢ - x¢ |)dx¢dx¢¢ =
0 0
Vr Vr
òò
wi ( x¢, t )wi ( x¢¢, t ) dx ¢dx¢¢
0 0
1
Vr2
Vr Vr
òò f
W
(| x ¢¢ - x ¢ |, CD )dx¢dx¢¢ .
(1.62)
0 0
Из (1.12, 1.44) следует, что ¶ 2 Dii (r ) / (¶x¢j ¶x¢¢j ) = -¶ 2 Dii (r ) / (¶rj ¶rj ) = -2 R 2 (r ) = -2W 2 (r ) , где
r = x ¢¢ - x ¢ . Используя это уравнение и теорему о градиенте
V
s
¶f( x )
ò0 ¶xi dx =ò0 f( x )ni ds ,
(1.63)
для произвольной области объёмом V , поверхностью s и нормалью к поверхности n ,
находим
n
e V2
V V
n
ò0 ò0 W (| x¢¢ - x¢ |)dx¢dx¢¢ = - 2 e V 2
2
¶ 2 Dii (| x¢¢ - x¢ |)
ò0 ò0 ¶x¢j ¶x¢¢j dx¢dx¢¢
V V
10
=-
n
2 e V2
s s
ò ò (n¢n¢¢) Dii (| x¢¢ - x¢ |)ds¢ds¢¢ = 0 0
h2
2V 2
s s
ò ò (n¢n¢¢) f
D
(| x¢¢ - x ¢ |, CD )ds¢ds¢¢ .
(1.64)
0 0
Отсюда, для сферы,
f r (r , CD ) = -
h2
2Vr2
sr sr
ò ò (n¢n¢¢) f
D
(| x ¢¢ - x ¢ |, CD )dsr¢ dsr¢¢ .
(1.65)
0 0
Интегрируя по углам, находим
r
r
ü
8p ì
3
2
I ( x¢) º ò (n¢n¢¢) f D (| x¢¢ - x¢ |, CD )dsr¢¢ = - 2 í2ò f D (2r , CD )r dr - r ò f D (2r , CD )rdr ý . (1.66)
r î 0
0
0
þ
sr
Учитывая, что интеграл I ( x¢) не зависит от x ¢ , определяем
s
r
h2 r
9h2 I 18
9
f r (r , CD ) = - 2 ò I ( x ¢)dsr¢ = = 6 ò f D (2r , CD )r 3dr - 4
4
2Vr 0
8pr
r 0
r
r
òf
D
(2r , CD )rdr ,
(1.67)
0
или, эквивалентно,
r
r
72
36
f r (r / 2, CD ) = 6 ò f D (r , CD )r 3dr - 4 ò f D (r , CD )rdr .
r 0
r 0
(1.68)
Если универсальная функция f r (r , CD ) известна, средняя скорость диссипации
определяется из определения (1.62):
e = nW 2r f r-1 (r , CD ) .
(1.69)
В инерционном диапазоне f D (r , CD ) = CD r 2/3 , отсюда, f r (r , CD ) = (27 / 14)CD (2r ) -4/3 и
e = 4(27CD /14) -3/ 2 r 2 W3r = Cr2 D 2 W3r ,
где
D º Vr1/3 = (4p / 3)1/3 r
(1.70)
и Cr = 2(4p / 3) -1/3 (27CD / 14)-3/ 4 , например, Cr = 0.172
при
CD = 7.2 . Формула (1.69) определяет среднюю скорость диссипации неявно, поскольку
r º r / h = n -3/ 4 e
1/ 4
r.
Используя
D 2 º D 2 / h2 = D 2 e
1/2
n -3/ 2 = f r-1/ 2 (r , CD )(n -1D 2 W r ) ,
определим универсальную функцию g r (W r , CD ) в виде параметрической зависимости:
{
}
W r º n -1D 2 W r = D 2 f r1/ 2 (r , CD ) ; g r (W r , CD ) º f r (r , CD ) (1 + Cr2W r )
-1
.
(1.71)
Эта функция позволяет выразить среднюю скорость диссипации
e = g r (n -1D 2 W r , CD ) {nW 2r + Cr2 D 2 W3r }
(1.72)
явно. В диссипативном ( W r << 1 ) и инерционном ( W r >> 1 ) диапазонах g r (W r , CD ) = 1 .
Отметим, что для изотропного поля скорости уравнения (1.62-1.72) справедливы,
если заменить W r на Sr , где
11
Sr ,ij º
1
Vr
Vr
ò Sij ( x, t )dx =
0
1
2Vr
sr
ò ( u ( x, t )n
i
j
+ u j ( x , t )ni ) dsr
(1.73)
0
и
Sr2 º 2 S r ,ij S r , ji .
(1.74)
1.3. Преобразование Фурье
1.3.1. Трехмерный спектр
Функциям в декартовом пространстве r соответствуют распределения в волновом
пространстве k . Определим прямое
1
fˆ ( k ) º F { f (r )} º 2
2p
¥ ¥ ¥
ò ò ò exp(-ikr ) f (r )dr
(1.75)
-¥ -¥ -¥
и обратное
f (r ) º F
преобразования
Фурье.
-1
Вводя
{
¥ ¥ ¥
}
1
fˆ ( k ) º
ò -¥ò -¥ò exp(ikr ) fˆ ( k)d k
4p -¥
сферические
координаты
(1.76)
( r , q, j) ,
так
что
r = ( r sin q cos j, r sin q sin j, r cos q ) и k × r = kr cos q , и интегрируя по q и j , получим для
сферической функции f (r ) = f (r ) :
¥ ¥ ¥
òò
¥
sin( kr )
f (r )r 2 dr ,
k
r
0
ò exp(-ikr ) f (r )dr = 4pò
-¥ -¥ -¥
(1.77)
следовательно,
µ
2 sin kr
fˆ ( k) º F { f (r )} = ò
f (r )r 2 dr .
p 0 kr
(1.78)
Аналогично, обратное преобразование
{
}
µ
sin kr ˆ
f (r ) º F -1 fˆ ( k) = ò
f ( k) k 2 d k .
kr
0
(1.79)
Определим распределение турбулентной энергии по масштабам вихрей
kr ( x , t , r ) º (1/ 2) ( ui ( x, t ) - ui ( x, t )
)( u ( x + r , t ) i
ui ( x + r , t )
)
= (1/ 2) ( k ( x , t ) + k ( x + r , t ) ) - (1/ 4) ( ui ( x + r , t ) - ui ( x , t ) )( ui ( x + r , t ) - ui ( x, t ) )
+ (1/ 4) ( ui ( x + r , t ) - ui ( x, t )
)
2
,
(1.80)
12
где
k ( x , t ) º (1/ 2) ( ui ( x , t ) - ui ( x , t )
)( u ( x, t ) i
ui ( x, t )
)
(1.81)
– турбулентная энергия. Для локально однородного поля скорости k ( x, t ) = k ( x + r , t ) = k ,
ui ( x + r , t ) = ui ( x, t ) и
kr (r ) = kr ( x, t , r ) = k - (1/ 4) Dii (r ) .
(1.82)
Скорости ui ( x , t ) и ui ( x + r , t ) не коррелируют при r ® ¥ , следовательно, kr ( x, t , ¥) = 0 и
Dii (¥) = 4k .
(1.83)
Распределение турбулентной энергии в волновом пространстве E ( k) определим из
µ
sin kr
E ( k)d k ,
k
r
0
kr (r ) º ò
(1.84)
отсюда, при r = 0 , имеем
µ
k = ò E ( k)d k .
(1.85)
0
Из (1.84) находим:
fˆ ( k) = k -2 E ( k) и
f (r ) = kr (r ) = ( Dii (¥) - Dii (r )) / 4 , подставив эти
выражения в (1.78-1.79), получим прямое
¥
E ( k) =
k 2 sin kr
( Dii (¥) - Dii (r ) ) r 2 dr
ò
2p 0 kr
(1.86)
и обратное
µ
æ sin kr ö
Dii (r ) = 4 ò ç1 ÷ E ( k) d k
kr ø
0è
(1.87)
преобразования Фурье. Формально, преобразования (1.86-1.87) корректны, но не
оптимальны из-за явной зависимости от Dii (¥) . Дифференцируя (1.87), находим
µ
1 d 2 dDii (r )
sin kr
r
=
2
E ( k) k 2 d k .
2
ò
dr
k
r
2r dr
0
(1.88)
По сути, уравнение (1.88) связывает распределения диссипации в волновом ( 2nk 2 E ( k) ) и
физическом ( nS 2 (r ) ) пространствах и может быть использовано для определения E ( k)
вместо (1.84). В частности, из него следует
µ
n d 2 dDii (r )
r
= e .
2
dr
dr
ò 2nE ( k)k d k = lim 2r
2
0
r -> 0
Соответствующее (1.88) прямое преобразование Фурье
(1.89)
13
µ
1 ì 1 d 2 dDii (r ) ü 1 sin kr æ d 2 dDii (r ) ö
E ( k) = F í 2
r
r
dr .
ý=
4 î r dr
dr þ 2p ò0 kr çè dr
dr ÷ø
(1.90)
Асимптотически, при k << rDI-1 и больших числах Рейнольдса,
1 sin kr æ 10CD e
E ( k) =
ç
2p ò0 kr çè
9
¥
ö
r 2/3 ÷ dr = e
÷
ø
2/3
5CD
k -5/3 ,
9G(1/ 3)
2/3
где G – гамма функция. Из определения универсальной константы CE º e
(1.91)
-2/3
k5/3 E ( k) в
инерционном диапазоне находим связь констант
CD =
9G(1/ 3)
CE .
5
(1.92)
1.3.2. Одномерный спектр
Определим прямое
¥
1
fˆ ( k1 ) º F1 { f (r1 )} º
ò exp(-ik1r1 ) f (r1 )dr1
2p -¥
(1.93)
и обратное
-1
1
f (r1 ) º F
{
¥
} ò exp(ik r ) fˆ (k )d k .
fˆ ( k1 ) º
1 1
1
(1.94)
1
-¥
преобразования Фурье. Для чётных функций f (r1 ) = f (- r1 ) :
¥
1
fˆ ( k1 ) = fˆ (-k1 ) = F1 { f (r1 )} = ò cos( k1r1 ) f (r1 )dr1 ,
p0
-1
1
f (r1 ) = F
{
}
¥
fˆ ( k1 ) = 2 ò cos( k1r1 ) fˆ ( k1 )d k1 .
(1.95)
0
Введем корреляции
knr ( x, t , r1e1 ) º (1/ 2) ( un ( x, t ) - un ( x, t )
)( u ( x + r e , t ) n
un ( x + r1e1 , t )
1 1
)
= (1/ 2) ( kn ( x, t ) + kn ( x + r1e1 , t ) ) - (1/ 4) ( un ( x + r1e1 , t ) - un ( x, t ) )( un ( x + r1e1 , t ) - un ( x , t ) )
+ (1/ 4) ( un ( x + r1e1 , t ) - un ( x , t )
)
kn ( x , t ) º (1/ 2) ( un ( x, t ) - un ( x, t )
)
2
,
(1.96)
где
2
,
(1.97)
здесь нет суммирования по повторяющемуся индексу n. Для однородной турбулентности
kn ( x, t ) = kn ( x + r1e1 , t ) = kn , un ( x + r1e1 , t ) = un ( x , t ) и
14
knr (r1e1 ) = knr ( x , t , r1e1 ) = kn - (1/ 4) Dnn (r1 ) .
Скорости
un ( x , t )
un ( x + r1e1 , t )
и
не коррелируют
при
(1.98)
r1 ® ¥ ,
следовательно,
knr ( x , t , ¥) = 0 и
Dnn (¥) = 4kn .
(1.99)
Определим одномерное спектральное распределение Enn ( k1 ) из
¥
2knr (r1e1 ) = 2kn - (1/ 2) Dnn (r1 ) º ò cos( k1r1 ) Enn ( k1 )d k1 ,
(1.100)
0
отсюда, прямое преобразование Фурье
¥
1
cos( k1r1 ) ( Dnn (¥) - Dnn (r1 ) ) dr1 .
Enn ( k1 ) =
4p ò0
(1.101)
Отметим, что
¥
1
Enn (0) =
( Dnn (¥) - Dnn (r ) ) dr ,
4p ò0
(1.102)
и, из (1.100) и условия Dnn (0) = 0 ,
¥
òE
nn
( k1 )d k1 = 2kn .
(1.103)
0
Дифференцируя (1.100) получим обратное преобразование Фурье
¥
d 2 Dnn (r1 )
= 2 ò cos( k1r1 ) Enn ( k1 ) k12 d k1 ,
2
dr1
0
(1.104)
соответствующее прямое преобразование
Enn ( k1 ) =
¥
¥
d 2 Dnn (r1 )
dD (r )
1
1
cos(
k
r
)
dr
=
sin( k1r1 ) nn 1 dr1 .
11
1
2 ò
2
ò
pk1 0
dr1
pk1 0
dr1
(1.105)
На практике удобнее вместо (1.105) использовать выражения
¥
1 d kI nn ( k)
1 sin( kr ) d 2 Dnn (r )
Enn ( k) = 2
, I nn ( k) º ò
dr .
k
dk
p 0 kr
dr 2
(1.106)
Дополнительно, можно получить из (1.105) другой способ вывода уравнения (1.99)
¥ ¥
¥
ü dD11 (r1 )
1 ì sin( k1r1 )
1 dD11 (r1 )
D (¥)
2k1 = ò E11 ( k1 )d k1 = ò í ò
d k1 ý
dr1 = ò
dr1 = 11
,
p 0 î 0 k1
2 0 dr1
2
0
þ dr1
¥
(1.107)
и, из (1.48 и 1.104), следует, что для изотропной турбулентности
d 2 D11 (0) 2 e
2ò nE11 ( k1 ) k d k1 = ò nE22 ( k1 ) k d k1 = ò nE33 ( k1 ) k d k1 = n
=
.
dr12
15
0
0
0
¥
¥
2
1
¥
2
1
2
1
(1.108)
15
1.3.3. Связь одно- и трехмерных спектров
Вначале, рассмотрим связи между E ( k) , E11 ( k1 ) и E22 ( k1 ) = E33 ( k1 ) при условии,
что поле скорости частично изотропно. Из (1.105) следует, что
µ
µ
d k 2 Eii ( k)
d 2 Dii (r )
dD (r )
1
1
dr
E
= - ò r sin kr
,
(
k
)
=
sin kr ii dr .
ii
2
ò
dk
p0
dr
pk 0
dr
(1.109)
Уравнение (1.90) запишем в виде
E ( k) =
µ
µ
d 2 Dii (r )
dD (r )
1
1
r
sin
k
r
dr
+
sin kr ii dr .
2
ò
ò
2pk 0
dr
pk 0
dr
(1.110)
Из (1.109-1.110) получим
E ( k) = -
k d ( E11 ( k) + 2 E22 ( k) )
.
2
dk
(1.111)
Дифференцирование (1.20) дает
r
d 2 D22 (r1 )
d 2 D11 (r1 ) r1 d 3 D11 (r1 )
1 1
r d F (r1 )
2
=
F (r )r 2 dr - 1
.
2
2
3
3 ò
dr1
dr1
2 dr1
r1 0
2 dr1
(1.112)
Используя равенство
ì d 3 D11 (r1 ) ü
æ
d ö ì d 2 D11 (r1 ) ü
F1 ír1
=
1
+
k
ý
ý,
ç
÷ F1 í
1
dr13 þ
d k1 ø î dr12 þ
î
è
(1.113)
после применения преобразования Фурье к (1.112), находим
E22 ( k1 ) +
ìï 1 r1
k12 d E11 ( k1 )
r d F (r1 ) üï
1
= - 2 F1 í 3 ò F(r )r 2 dr + 1
ý.
2 d k1 k1
2 k1 ïî r1 0
2 dr1 ïþ
(1.114)
Для изотропного поля скорости F(r ) = 0 , следовательно, аналог (1.39) в волновом
пространстве
k 2 d æ E11 ( k) ö
E22 ( k) = ç
÷,
2 dk è k ø
(1.115)
и, из (1.111, 1.115), получаем связи
E ( k) =
k3 d æ 1 dE11 ( k) ö d æ E ( k) ö
d æ 1 dE22 ( k) ö
ç
÷,
ç
÷ = -k
ç
÷.
2 dk è k dk ø dk è k ø
dk è k dk ø
(1.116)
1.4. Связи колмогоровских констант
В
инерционном
диапазоне
связь
функций
с
физическими
определяется через набор (в общем случае не универсальных) констант
переменными
16
CE = e
-2/3
k5/3 E ( k) , CE ,11 = e
CD , NN = e
-2/3
-2/3
r -2/3 DNN (r ) , CD , LL = e
CW = (9 / 5) e
-2/3
k5/3 E11 ( k) , CЕ ,22 = e
-2/3
-2/3
k5/3 E22 ( k) , CD = e
r -2/3 DLL (r ) , CS = (9 / 5) e
r 4/3W 2 (r ) , CF = (9 / 5) e
-2/3
-2/3
-2/3
r -2/3 Dii (r ) ,
r 4/3 S 2 (r ) ,
r 4/3F (r ) .
(1.117)
Для частично изотропного поля скорости можно использовать уравнения (1.1-1.33)
при F(r ) ¹ 0 . Существуют две независимые константы, выбрав функции
Dii (r ) = CD e
2/3
r 2/3 , W 2 (r ) = (5 / 9)CW e
2/3
r -4/3 ,
(1.118)
получим, из (1.12, 1.33),
R 2 (r ) = (5 / 9)CD e
2/3
r -4/3 , S 2 (r ) = (5 / 9)(2CD - CW ) e
F(r ) = (5 / 9)(CD - CW ) e
2/3
2/3
r -4/3 ,
r -4/3
(1.119)
в инерционном диапазоне. Используя (1.119), находим
r
r d F(r1 ) ïü
(CD - CW )
ïì 1 1
F1 í 3 ò F (r )r 2 dr + 1
e
ý=2 dr1 ïþ
27
ïî r1 0
2/3
F1 {r1-4/3 } =
2(CD - CW )
e
9G(1/ 3)
2/3
k1/3
1
(1.120)
и, из зависимости (1.114),
k12 d E11 ( k1 )
(C - CW )
E22 ( k1 ) +
=- D
e
2 d k1 k1
9G(1/ 3)
2/3
k1-5/3 .
(1.121)
Из (1.6, 1.20, 1.92, 1.111, 1.121) определяем зависимости констант от CD , LL и CD , NN
CD = CD , LL + 2CD , NN , CW = -7CD , LL + 8CD , NN , CS = 9CD , LL - 4CD , NN , CF = 8CD , LL - 6CD , NN
CE ,11 =
8CD , LL
11G(1/ 3)
, CE ,22 =
8CD , NN
11G(1/ 3)
, CE =
5CD , LL + 10CD , NN
9G(1/ 3)
,
(1.122)
или от CD и CW
CD , LL =
4CD - CW
7C + CW
, CD , NN = D
, CS = 2CD - CW , CF = CD - CW ,
11
22
8C - 2CW
7C + CW
5CD
CE ,11 = D
, CE ,22 = D
, CE =
.
33G(1/ 3)
33G(1/ 3)
9G(1/ 3)
(1.123)
Используя универсальные функции (1.10, 1.13), можно записать:
Dii (r ) = ( e n )
-1/2
f D (r , CD ) , R 2 (r ) = e n -1 f R (r , CD ) , W 2 (r ) = e n -1 f R (r , CW ) ,
S 2 (r ) = e n -1 ( 2 f R (r , CD ) - f R (r , CW ) ) , F(r ) = e n -1 ( f R (r , CD ) - f R (r , CW ) ) ,
(1.124)
здесь выражения для W 2 (r ) , S 2 (r ) и F(r ) точно выполняются в диссипативном и
инерционном диапазонах и приблизительно на переходе между ними.
17
Для изотропного поля скорости связь универсальных констант:
CD , LL = (3 / 11)CD , CD , NN = (4 /11)CD , CS = CD , CF = 0 ,
CE ,11 =
6CD
8CD
5CD
, CE ,22 =
, CE =
.
33G(1/ 3)
33G(1/ 3)
9G(1/ 3)
(1.125)
Есть только одна независимая константа, например, если выбрана CE ,11 = 0.49 , то CE = 1.5 ,
CD = 7.22 , CD , LL = 1.97 , CD , NN = 2.62 и CE ,22 = 0.65 . Для изотропного поля скорости
зависимости (1.124) вырождаются:
Dii (r ) = ( e n )
-1/ 2
f D ( r , C D ) , R 2 ( r ) = W 2 ( r ) = S 2 ( r ) = e n -1 f R ( r , C D ) , F ( r ) = 0 .
(1.126)
2. Модель для структурных функций
Предлагается следующая аппроксимация продольной компоненты структурной
функции:
D11 (r ) º ( e n )
-1/ 2
D11 (r ) = fU ,11 (r ) f E ,11 (r ) ,
(2.1)
где определены универсальная (для диссипативного и инерционного диапазонов)
fU ,11 (r ) =
(1/ 15)r 2
{1 + a11 (r / rDI )2 + b11 (r / rDI )3 + (r / rDI )4 }
1/3
(2.2)
и уникальная для каждого течения (в энергетическом диапазоне)
f E ,11 (r ) =
1
{(1 + rDI2 r 2 + r 4 )3d11 / 2 + j11 (r / r11 )2 + (r / r11 )4 }
1/6
(2.3)
функции; rDI , a11 , b11 – универсальные и r11 , j11 , d11 – уникальные для каждого течения
параметры аппроксимации.
В диапазоне реализации локальной изотропности поля скорости D11 (r ) = fU ,11 (r ) .
Асимптотика при r >> rDI : fU ,11 (r ) = CD , LL r 2/3 , отсюда и из (2.2) находим характерный
масштаб, разделяющий диссипативный и инерционный диапазоны турбулентности:
rDI = (15CD , LL )3/ 4 . В диапазоне r << rDI
D11 (r ) = (1/15)r 2 - (a11rDI-2 / 45)r 4 - (b11rDI-3 / 45)r 5 + O(r 6 ) .
(2.4)
Следовательно, во-первых, предполагается, что существует локальная изотропность для
самых мелких вихрей ( r << rDI ) и, следовательно, d 2 D11 (0) / dr 2 = (2 /15) (1.52); во-
18
вторых, параметры a11 и b11 универсальны и определяют поведение структурной
функции в диссипативном диапазоне.
В диапазоне r >> rDI :
D11 (r ) =
{r
CD , LL r 2/3
6 d11
+ j11 (r / r11 ) 2 + (r / r11 )4 }
1/6
.
(2.5)
Первый член в скобках (2.5) моделирует влияние на инерционный диапазон
перемежаемости, D11 (r ) = CD , LL r 2/3-d11 при rDI << r << r11 , это выражение разъясняет смысл
параметра d11 . Из (2.5) следует
r11 = éë D11 (¥) / CD , LL ùû
3/2
(2.6)
и зависимость D11 (r ) в диапазоне r >> r11
-2
D11 (¥) - D11 (r ) j11 æ r ö
=
ç ÷ ,
D11 (¥)
6 è r11 ø
(2.7)
которая проясняет роль параметра j11 . Характерный размер, разделяющий инерционный
и энергетический диапазоны, лежит на пересечении (2.7) и D11 (r ) = CD , LL r 2/3 :
-3/4 3/2
r11, IE = j11
CD , LL D119/ 4 (¥) .
(2.8)
Аппроксимации D22 (r ) и D33 (r ) аналогичны, приведем их для D22 (r ) .
D22 (r ) º ( e n )
-1/ 2
D22 (r ) = fU ,22 (r ) f E ,22 (r ) ,
(2.9)
где определены универсальная
2
3
4
2r 2 {1 + (5a 22 / 6)(r / rDI ) + (3b22 / 4)(r / rDI ) + CD , NN / (2CD , LL )(r / rDI ) }
fU ,22 (r ) =
,
4/3
15
{1 + a (r / r )2 + b (r / r )3 + (r / r )4 }
22
DI
22
DI
(2.10)
DI
и уникальная для каждого течения
f E ,22 (r ) =
1
{(1 + r
2
DI
4 3d22 / 2
r +r )
2
+ j22 (r / r22 ) 2 + (r / r22 ) 4 }
1/6
(2.11)
функции. Функция (2.10) построена таким образом, что при a 22 = a11 , b22 = b11 и
CD , NN = (4 / 3)CD , LL следствие изотропности (1.39) выполняется строго:
fU ,22 (r ) = fU ,11 (r ) +
В диапазоне r >> rDI
r dfU ,11 (r )
.
2 dr
(2.12)
19
D22 (r ) =
CD , NN r 2/3
{r 6d22 + j22 (r / r22 )2 + (r / r22 )4 }
1/6
.
(2.13)
Отсюда, при r ® ¥ , имеем
3/2
r22 = éë D22 (¥) / CD , NN ùû ,
(2.14)
-3/ 4 3/2
9/ 4
r22, IE = j22
CD , NN D22
(¥ ) .
(2.15)
и, аналогично (2.8),
Для определения распределений энергии Enn ( k) использовались уравнения (1.106),
но, из-за конечности диапазона r0 £ r £ rN , с коррекцией
r
1 ì dDnn (r0 ) dDnn (rN ) ü 1 N sin( kr ) d 2 Dnn (r )
I nn ( k) = í
dr ,
ý+
dr þ p rò0 kr
dr 2
p î dr
(2.16)
чтобы равенство I nn (0) = 0 выполнялось строго. Также использовались логарифмическая
сетка ( ri +1 / ri = const ) и численная схема для интегрирования
ri+1
ò
ri
( sin(kr ) )i +1/ 2
sin( kr ) dF (r )
dr =
( F (ri +1 ) - F (ri ) )
dr
kr
kri +1/2
ìï sin ( k(ri +1 - ri ) / 2 ) üï ìï sin ( kri +1/ 2 ) üï
=í
ýí
ý ( F (ri +1 ) - F (ri ) ) ,
îï ( k(ri +1 - ri ) / 2 ) þï îï ( kri +1/ 2 ) þï
(2.17)
i = 0,1, 2...N , которая получается, если положить
ri +1/ 2 =
ri +1 + ri
1
, ( sin( kr ) )i +1/ 2 =
2
(ri +1 - ri )
ri +1
ò sin(kr )dr .
(2.18)
ri
2.1. Обработка эксперимента с использованием модельных функций
Продемонстрируем,
как
данная
простая
модель
описывает
выдающийся
эксперимент [3]. В [3] представлен набор из четырёх экспериментальных результатов,
каждый набор можно охарактеризовать числом Рейнольдса по тейлоровскому масштабу
( Rl ), скоростью в ядре потока (Ue) и расстоянием от стенки (y). Основные параметры
потока представлены в таблице 1. Отметим, что в [3]
Rl º
(u
1 - u1
)
2 1/ 2
l / n , где l 2 º
(u
1
(1.53), Rl = (15 / 4)1/2 D11 (¥) = 601/2 ( e n )
- u1
-1/2
k1 .
)
2
/
( ¶u1 / ¶x1 )
2
Rl
определяется, как
, следовательно, используя
20
Rl
Ue, м/с
y, мм
k1 , м2/с2
500
600
1400
Таблица 1
1450
10
100
0.435
10
515
0.1715
50
100
13
50
400
4.75
k2 , м2/с2
0.155
0.087
3.75
2.25
k3 , м2/с2
0.240
0.106
6.5
2.975
e , м2/с3
3.1
0.33
342
49
Далее, на рисунки из [3] с экспериментальными данными (черный цвет)
наносились зависимости модели (зеленый цвет), смотри рис. 1-11. Использовались
следующие параметры: CD , LL = 1.97 , CD , NN = (4 / 3)CD , LL , a11 = a 22 = a33 = 3.2 , остальные
параметры приведены в таблице 2.
Таблица 2
Rl = 500
n
b nn
1
0.2
2
0.4
3
0.1
dnn
0.038
0.046
0.032
jnn
11
12
12
Rl = 600
n
b nn
1
0.5
2
0.2
3
0.5
dnn
0.01
0.03
0.024
jnn
5.0
3.5
2.7
Rl = 1400
n
b nn
1
0.5
2
0
3
0
dnn
0.008
0.022
0.005
jnn
10
5
5
Rl = 1450
n
b nn
1
0.5
2
0
3
0.3
dnn
0.006
0.004
0.001
jnn
5.4
3.2
3.2
21
Рис. 1. Сравнение продольного и поперечного спектров, Rl = 500 и Rl = 1400 :
[3, рис. 7] – черные линии; модель – зеленые точки
22
Рис. 2. Сравнение продольного и поперечного спектров, Rl = 600 и Rl = 1450 :
[3, рис. 8] – черные линии; модель – зеленые точки
23
Рис. 3. Сравнение колмогоровского универсального масштабирования для продольного спектра:
[3, рис. 9] – черный цвет; модель при Rl = 600 и Rl = 1450 – цветные точки
Рис. 3 наглядно показывает, что поведение продольного спектра очень близко к
универсальному в диапазоне r < rIE , то есть универсальные закономерности у хаоса
(турбулентности) существуют! Отметим, что переход к чисто диссипативному диапазону
(где D11 (r ) = (1/15)r 2 ) имеет место для модельных уравнений при k1h » 1 . Небольшое
различие двух модельных кривых при k1h > 1 определяется влиянием численных ошибок
( : 10-5 ) при интегрировании.
24
Рис. 4. Сравнение диссипативного спектра, Rl = 600 :
[3, рис. 10] – черные точки; модель – зеленые точки
25
Рис. 5. Сравнение спектра в диссипативном диапазоне, Rl = 600 :
[3, рис. 11] – черные кружки; модель – зеленые точки
Наклоны кривых в диапазоне 0.2 < k1h < 1 определяются параметрами a nn , при
выбранном значении
a11 = a 22 = a33 = 3.2 ,
модель
достаточно хорошо
описывает
экспериментальные данные. Равенство этих параметров указывает на изотропность в
диссипативном диапазоне k1 > 0.2 / h .
26
Рис. 6. Сравнение продольного и поперечного спектров, Rl = 600 :
[3, рис. 12] – черные линии и точки; модель – зеленые точки
27
Рис. 7. Сравнение продольного и поперечного спектров, Rl = 1450 :
[3, рис. 13] – черные линии и точки; модель – зеленые точки
28
Рис. 8. Сравнение продольного и поперечного спектров, Rl = 500 :
[3, рис. 14] – черные линии и точки; модель – зеленые точки
29
Рис. 9. Сравнение продольного и поперечного спектров, Rl = 1400 :
[3, рис. 15] – черные линии и точки; модель – зеленые точки
Совпадение с экспериментальными данными достаточно хорошее, за одним
исключением:
в
диссипативном
диапазоне
при
Rl = 1400 .
Вероятно,
различие
определяется невозможностью, как отмечено в [3], разрешить диссипативный диапазон
при этом числе Рейнольдса.
30
Зависимости построены при CD , NN = (4 / 3)CD , LL , достаточно хорошее совпадение с
экспериментом в инерционном диапазоне (если он существует) указывает, что в этом
диапазоне турбулентность почти изотропна. Параметры b nn определяют значения
максимумов на переходе от диссипативного к инерционному диапазону; с ростом b
величина максимума уменьшается. Для изотропной турбулентности b11 = b22 = b33 , в
нашем случае параметры менялись от 0 до 0,5. Логично предположить, что если в
диссипативном и инерционном диапазонах турбулентность почти изотропна, то она почти
изотропна на переходе и различие параметров определяется сложностью точного
экспериментального определения максимумов.
Влияние перемежаемости на универсальный диапазон происходит через параметры
dnn , но установить зависимость этих параметров от Rl не представляется возможным изза малой статистики.
Важный вывод из эксперимента [3] состоит в том, что переход от инерционного к
энергетическому диапазону (2.8, 2.15) зависит не только от Rl но, существенно, от
структуры крупных вихрей, которые определяют значения jnn .
Далее
рассмотрим
сравнение
структурных
использованы средние скорости диссипации
e
функций.
На
рис. 1-9
были
из табл. 1, которые были вычислены
интегрированием спектральных распределений энергии. Но на графиках в [3] структурные
функции представлены при e% =0,26 м2/с3 для Rl = 600 и e% =40 м2/с3 для Rl = 1450 ,
полученных из оценки поведения структурных функций третьего порядка. Структурные
функции приведены на рис. 10-11, где зеленые линии соответствуют зависимостям
e
-2/3
r -2/3 Dnn (r ) от r / h = e
r / h% = e%
1/ 4
1/ 4
n3/4 r , зеленые точки – зависимостям e%
-2/3
r -2/3 Dnn (r ) от
n3/4 r . Подчеркнем, что модельные распределения, которые представлены
точками, близки распределениям из [3], но модельным распределениям в волновом
пространстве (рис. 1-9) соответствуют линии рис. 10-11.
31
Рис. 10. Сравнение структурных функций второго порядка, Rl = 600 :
[3, рис. 17] – черные кружки; модель, e =0,33 м2/с3 – сплошные зеленые линии;
модель, e =0,26 м2/с3 – зеленые точки
32
Рис. 11. Сравнение структурных функций второго порядка, Rl = 1450 :
[3, рис. 18] – черные кружки; модель, e =49 м2/с3 – сплошные зеленые линии;
модель, e =40 м2/с3 – зеленые точки
33
3. Изотропное поле завихренности
Рассмотрим теорию Колмогорова, основанную на альтернативных гипотезах:
1) случайное поле завихренности ω( x , t ) изотропно.
2) В универсальном диапазоне r < rIE статистика определяется n и
e ; в
инерционном диапазоне rDI < r < rIE – только e .
3) Поле скорости w (r ) = u(r + x0 , t ) - u( x0 , t ) изотропно для самых мелких вихрей в
каскаде Ричардсона, следовательно, e = nW 2 (0) .
Из этих гипотез следуют уравнения (1.56-1.57) для двухточечных моментов и
W 2 (r ) = n -1 e f W (r , CW ) ,
(3.1)
где f W (r , CW ) – универсальная функция, CW – универсальная константа. В инерционном
диапазоне
f W (r , CW ) = (5 / 9)CW r -4/3 .
(3.2)
3.1. Одноточечная функция плотности вероятности завихренности
Рассмотрим одноточечную функцию плотности вероятности (ФПВ) завихренности
Pω (ω) .
Из
инвариантности
относительно
вращения
вокруг
оси
e3
на
угол
% 1 = w12 + w22 , w
% 2 = 0, w
% 3 = w3 ) следует
j = arctan(w2 / w1 ) ( w
Pω ( w1 , w2 , w3 ) = Pω
(
)
w12 + w22 , 0, w3 .
(3.3)
Здесь использовались уравнения преобразования завихренности при вращении; они
просты, например, при вращении вокруг оси e3 на угол j преобразования координат:
x%1 = cos jx1 + sin jx2 , x%2 = - sin jx1 + cos jx2 , x%3 = x3 ;
(3.4)
преобразование завихренности:
% 1 = cos jw1 + sin jw2 , w
% 2 = - sin jw1 + cos jw2 , w
% 3 = w3 .
w
(3.5)
)
(
Совершив дальнейшее вращение на угол j = arctan w3 / w12 + w22 вокруг оси e2 , находим,
что трехмерная ФПВ является сферической функцией
Pω ( w1 , w2 , w3 ) = Pω
(
)
w12 + w22 + w32 , 0, 0 = Pω ( w) .
(3.6)
34
Введем одномерную ФПВ
µ µ
Pwi (wi ) º
ò ò P (ω)Õ d w
ω
n
n ¹i
-µ -µ
.
(3.7)
Из сферичности (3.6) следует
Pw1 (w) = Pw1 (-w) , Pw1 (w) = Pw2 (w) = Pw3 (w) .
(3.8)
Отметим, что из чётности следует, что моменты нечетного порядка равны нулю, это
принципиально не так для случайных переменных ¶ui / ¶x j и Sij . Связь одно- и
трехмерных ФПВ:
µ µ
Pw1 (w1 ) º
òò
-µ -µ
¥
Pw (w)d w2 d w3 = 2p ò Pw (w)wd w , Pw (w) = w1
1 dPw1 (w)
,
2pw d w
(3.9)
отсюда,
¥
w
n
¥
= ò w Pw (w)4pw d w = -2 ò w1n+1
2
n
0
dPw1 (w1 )
0
d w1
d w1
¥
= (n + 1) ò | w1 |n Pw1 (w1 )d w1 =(n + 1) | w1 |n .
(3.10)
-¥
Следовательно,
s2 º var(ω) = w2 = 3s12 ,
дисперсия
где
s12 º var(w1 ) = w12 ,
и
нормализованные моменты
M n º (w / s) n =
n +1
| w1 / s1 |n .
n/ 2
3
(3.11)
В системе координат, привязанной к центру сферы, средняя по сфере
завихренности ω r инвариантна относительно вращения, следовательно, её ФПВ тоже
сферическая функция
Pωr (ω) = Pωr (w) .
(3.12)
Далее рассмотрим модельные ФПВ завихренности, где ω означает случайную
величину завихренности либо в точке, либо осредненную по объёму сферы. Трех- и
одномерные нормальные ФПВ:
æ 3w2 ö
æ w12 ö
33/ 2
1
PN ,ω (ω) =
exp ç exp ç - ÷ ,
÷ , PN ,w1 (w1 ) =
(2p)3/2
(2p)1/ 2
è 2 ø
è 2 ø
(3.13)
для которых
M n = wn =
2 n / 2+1 G ( (n + 3) / 2 )
,
p1/2 3n /2
(3.14)
35
где ω º ω / s , w1 º w1 / s1 .
При турбулентном течении в некоторые моменты времени происходят события с
w2 >> 1 из-за перемежаемости и эти события вносят значительный вклад в статистику.
Вероятность этих событий определяется ФПВ, которая существенно отличается от
нормальной.
Следуя
Колмогорову,
предположим,
что
распределение
случайной
переменной V = ln w имеет нормальную ФПВ
æ (V - V
1
PV ( V ) =
exp ç (2p)1/2 sln
ç
2s2ln
è
)
2
ö
÷,
÷
ø
(3.15)
wn = exp ( n V + n 2sln2 / 2 ) . Из очевидного равенства
где s2ln º var(V) . Для этой ФПВ
w2 = 1 следует, что V = -s2ln и
M n = wn = exp ( n(n - 2)s2ln / 2 ) .
(3.16)
Из равенства вероятностей в пространствах w и V получаем Pω (ω)4pw2 d w = PV (V)d V ,
следовательно, Pω (ω) = PV (V) / (4pw3 ) . Используя
æ (V - V
exp ç 2sln2
ç
è
)
2
ö
æ V 2
æV V ö
÷ = exp ç 2 ÷ exp ç - 2
ç 2sln
÷
è sln ø
è
ø
ö
æ V2 ö 1
æ sln2
exp
exp
÷
ç- 2 ÷ =
ç÷
è 2
è 2sln ø w
ø
æ ln 2 w ö
ö
exp
ç - 2 ÷ (3.17)
÷
ø
è 2sln ø
находим лог-нормальную ФПВ
PK ,ω (ω) =
æ ln 2 w ö
exp(-sln2 / 2)
exp
ç- 2 ÷ .
2(2p)3/ 2 sln w4
è 2sln ø
(3.18)
Одномерная ФПВ определяется из (3.9)
exp ( 3sln2 / 2 ) ìï
æ ln(w12 / 3) + 4s2ln
PK ,w1 (w1 ) =
1
erf
í
ç
481/ 2
23/2 sln
è
îï
где
x
erf ( x) = 2p-1/ 2 ò exp(-t 2 )dt
0
–
функция
ошибок;
ö üï
÷ý ,
ø þï
erf (- x) = -erf ( x) ,
(3.19)
erf (¥) = 1 .
Особенность этой ФПВ: она практически постоянна (рис. 12) в окрестности нуля и равна
12-1/ 2 exp ( 3sln2 / 2 ) .
36
PK ,ω (ω1 )
1
0.46
0.45
0.44
0.43
-0.4
-0.2
0
0.2
0.4
ω1
Рис. 12. Лог-нормальная ФПВ PK ,w1 (w1 ) при sln = 0.557 в окрестности нуля
Согласно
экспериментам
при
w1 > s1
существует
практически
линейная
зависимость логарифма ФПВ от w1 , поэтому в качестве модельной ФПВ предлагается
PT ,w1 (w1 ) . Введем случайную переменную x с ФПВ
Px (x) =
1 derf (x)
= p-1/2 exp(-x2 ) ,
2 dx
(3.20)
Определим связь случайных переменных x и w1 как x = aw1 (b2 + a 2 w12 ) -1/ 4 , тогда
x2 é x2 + 4b2 + x 4 ù
dx
a(2b2 + a 2 w12 )
û.
2
=
, w1 = ë
2
2 2 5/4
2
2a
d w1 2(b + a w1 )
(3.21)
Из равенства PT ,w1 (w1 )d w1 = Px (x)d x находим
PT ,w1 (w1 ) =
æ
ö
a(2b2 + a 2 w12 )
a 2 w12
exp
,
ç
1/ 2
2
2 2 5/ 4
2
2 2 1/ 2 ÷
2p (b + a w1 )
è (b + a w1 ) ø
(3.22)
и PT ,w1 (0) = a / (pb)1/2 . Из (3.22) следует, что параметр
a = lim
d ln PT ,w1 (w1 )
d w1
w1 ®¥
.
(3.23)
Удобно использовать эту ФПВ 1) из-за ясности интерпретации параметра a и 2) диапазон
его изменения a » 0.9...1.2 довольно узок. Отметим, что параметр a не является
универсальным, он медленно уменьшается с ростом Rl . Моменты
¥
| w1 |
n
= 2ò w P
n
1 T ,w1
0
¥
(w1 )d w1 = 2 ò w P (x)d x = 2p
n
1 x
0
-1/ 2
¥
òw
n
1
0
exp(-x2 )d x
37
=p
-1/2
-n
a G(1/ 2 + n) + 2p
-1/ 2
éæ 2
x + 4b2 + x4
ê
a ò ç
êç
2
0 è
ë
-n
¥
ö
÷
÷
ø
n /2
ù
- x ú xn exp(-x2 )d x ,
ú
û
n
(3.24)
причём существует рекурсивная последовательность G(n + 1/ 2) = (n - 1/ 2)G(n - 1 / 2) и
G(1 / 2) = p1/ 2 . Из равенства w12 = 1 находим связь a и b :
¥
a 2 = 3 / 4 + p-1/ 2 ò é 4b2 + x4 - x2 ù x 2 exp(-x2 )d x .
ë
û
0
(3.25)
Простая аппроксимация (3.25):
b = 1.07(a 2 - 3 / 4)1/ 2 + 1.44(a 2 - 3 / 4) + 0.08(a 2 - 3 / 4)3/2 .
Зависимости b и PT ,w1 (0) от a показаны на рис. 13-14.
β
3
2
1
0
0.8
1
1.2
1.4
α
Рис. 13. Зависимость b от a ; точки – (3.25), линии – (3.26)
PT ,ω1 (0)
1.4
1
0.6
0.8
1
1.2
1.4
α
Рис. 14. Зависимость PT ,w1 (0) от a ; точки – (3.25), линии – (3.26)
Поведение PT ,w1 (w1 ) для разных a представлено на рис. 15.
(3.26)
38
PT ,ω
1
10
0
10-1
10-2
10-3
10-4
-10
-5
5
0
10
ω1
Рис. 15. PT ,w1 (w1 ) , — – a =0.9; — – a =1.0; — – a =1.1; — – a =1.2
Нормализованные моменты M3 и M4 (3.14, 3.24) показаны на рис. 16.
Mn
7
6
5
4
3
2
1
0.8
1
1.2
1.4
α
Рис. 16. Нормализованные моменты: — – M3; — – M4
Определим распределения как подобные, когда совпадают w2 и w4 , отсюда получаем
связь между sln для PK ,w1 (w1 ) и a для PT ,w1 (w1 ) (рис. 17). Её аппроксимация
sln » 0.295 + 0.106 / a + 0.2 / a 2 .
(3.27)
39
σ ln
0.7
0.6
0.5
0.8
1
1.2
α
1.4
Рис. 17. Зависимость sln и a для подобных распределений;
точки – точные уравнения, линии – (3.27)
Сравнение трех модельных ФПВ показано на рис. 18.
P
100
10-2
10-4
10-6
-10
-5
0
5
10
ω1
Рис. 18. ФПВ: — – PN ,w1 (w1 ) ; — – PK ,w1 (w1 ) при sln = 0.557
и — – PT ,w1 (w1 ) при a = 1.1
3.2. Моделирование диссипации для сферы
Покажем метод моделирования диссипации для инерционного диапазона.
Рассмотрим сферу радиуса r , объёмом Vr = (4p / 3)r 3 , поверхностью sr = 4pr 2 и нормалью
40
к поверхности n. Определим осреднение по объёму завихренности wr ,i (1.60), w2r º wr ,i wr ,i
и s2r = W 2r º w2r . Предполагая зависимость e = e(wr ) , ищем её в виде
¥
e = å an wnr ,
(3.28)
n =1
отсюда,
¥
e = å an M n srn .
(3.29)
n =1
Средняя скорость диссипации e = 4(27CW /14) -3/2 r 2 s3r (1.70), следовательно, ненулевой
коэффициент только a3 и
e = 4(27CW / 14)-3/ 2 M 3-1r 2 w3r .
(3.30)
При реализации события с завихренностью wr масштаб минимальных вихрей в каскаде
Ричардсона:
hwr = ( e / e ) h = (wr / sr )-3/ 4 M 31/ 4h .
1/ 4
(3.31)
Таким образом, при реализации случайной переменной wr величина скорости
диссипации может быть получена из (3.30), масштаб минимальных вихрей из (3.31).
Подчеркнем, что при реализации события перемежаемости с wr >> sr диссипация
реализуется в вихрях очень малого масштаба. Следовательно, в общем случае,
моделирование диссипации необходимо даже при использовании прямых численных
расчетов уравнений Навье-Стокса, в которых линейный размер разностных ячеек D
сравним с колмогоровским микромасштабом h.
4. Теория Колмогорова и LES
Рассмотрим LES модель. Разностная сетка для LES строится так, чтобы размер
любой разностной ячейки D лежал в диапазоне D < rIE . Здесь и далее характерный размер
ячейки определяется через её объём VD : D = VD1/3 ; соответствующее характерное время
tD = ( D / h)
2/3
th . В LES используются уравнения, которые формально совпадают с
уравнениями Навье-Стокса, следовательно, для них выполняются условия сохранения
массы, импульса и энергии и они инвариантны относительно основной группы
преобразований: переносов, вращений и зеркальных отражений.
41
Отличие от точных уравнений Навье-Стокса состоит в следующем. Во-первых, все
используемые гидродинамические параметры – это осредненные по времени tD
переменные, поскольку высокочастотные движения с частотой, большей tD-1
не
разрешимы в принципе. Во-вторых, масштаб минимальных вихрей в этой модели
соответствует размеру ячейки D , несмотря на то, что в действительности диссипация
происходит в колмогоровских вихрях. Поэтому прямое применение уравнений НавьеСтокса приводит к занижению диссипации в » ( h / D )
4/3
раз, и, по этой причине, для
обеспечения реальной скорости диссипации используются уравнения Навье-Стокса с
эффективным коэффициентом вязкости. Для его определения используется подсеточное
моделирование в предположении, что устанавливается равновесное распределение вихрей
каскада Ричардсона в диапазоне r < D , соответствующее мгновенным значениям скорости
диссипации.
Ранее было установлено, что равенство
S 2 (r ) = W 2 (r ) следует из условия
изотропности поля скорости, но в экспериментах и расчетах при умеренных числах
Рейнольдса наблюдается их существенное различие. Это можно показать из обработки
экспериментальных данных [3] (рис. 19-20): использовались аппроксимации D11 (r ) ,
D22 (r ) и D33 (r ) (рис. 6-7), затем вычислялись R 2 (r ) (1.12), F(r ) (1.19) , S 2 (r ) и W 2 (r )
(1.33).
6
5
4
3
2
1
0
10-1
100
Рис. 19. Rl = 600 , функции: — – e
101
-2/3
102
r 4/3W 2 ( r ) ; — – e
103
-2/3
104
r
r 4/3 R 2 ( r ) ; — – e
-2/3
r 4/3 S 2 ( r ) ;
42
6
5
4
3
2
1
0
10-1
100
101
Рис. 20. Rl = 1450 , функции: — – e
-2/3
102
r 4/3W 2 ( r ) ; — – e
103
-2/3
104
r
r 4/3 R 2 ( r ) ; — – e
-2/3
r 4/3 S 2 ( r ) ;
Очевидно, что с ростом числа Рейнольдса различие функций в инерционном диапазоне
уменьшается, но при Rl = 600 оно значительное. Поэтому желательно не использовать
гипотез, из которых следует, что функции S 2 (r ) и W 2 (r ) равны и универсальны в
инерционном диапазоне. Близко к универсальному поведение продольной структурной
функции (см. рис. 3). Можно оценить отклонения от универсальности для CD ,nn из рис. 1011 (зеленые линии). Определив
{
CD ,nn = max e
r
-2/3
}
r -2/3 Dnn (r ) , CW = -7CD ,11 + 4CD ,22 + 4CD ,33 ,
(4.1)
находим:
CD ,11 = 1.83 , CD ,22 = 2.26 , CD ,33 = 2.30 , CD = 6.39 , CW = 5.43 при Rl = 600 ;
CD ,11 = 1.88 , CD ,22 = 2.56 , CD ,33 = 2.59 , CD = 7.03 , CW = 7.44 при Rl = 1450 .
Видно, что наименьшие изменения характерны именно для константы CD ,11 , то есть её
поведение наиболее близко к универсальному.
В результате, используются следующие гипотезы теории Колмогорова:
1) случайное поле w (r ) = u(r + x0 , t ) - u( x0 , t ) частично изотропно.
2) Статистика в диссипативном диапазоне определяется n и e , в инерционном
диапазоне – e . При умеренных числах Рейнольдса не вся статистика универсальна из-за
отклонения от универсальности констант CD и CW ; близка к универсальной статистика,
которая зависит исключительно от константы CD , LL = (4CD - CW ) /11 .
43
3) Случайное поле w (r ) = u(r + x0 , t ) - u( x0 , t ) изотропно для самых мелких вихрей
в каскаде Ричардсона, следовательно, e = nS 2 (0) = nW 2 (0) .
Определяя линейную форму
P 2 (r ) º (1 - b) S 2 (r ) + bW 2 (r )
(4.2)
и используя (1.122-1.124) и уравнение
af R (r , C1 ) + (1 - a) f R (r , C2 ) = f R (r , aC1 + (1 - a)C2 ) ,
которое
строго
выполняется
в
диссипативном
и
инерционном
(4.3)
диапазонах
и
приблизительно на переходе между ними, находим
P 2 ( r ) = e n -1 f R ( r , C P ) ,
(4.4)
где
CP = (9 - 16b)CD , LL + (-4 + 12b)CD , NN .
(4.5)
Функция P 2 (r ) почти универсальна (зависит только от CD , LL ) при
b = 1/ 3 , CP = (11/ 3)CD , LL , P 2 (r ) = (2 S 2 (r ) + W 2 (r )) / 3 .
(4.6)
Рассмотрим моделирование диссипации для сферы в инерционном диапазоне.
Используя осредненные по объёму wr ,i (1.60), W 2r (1.61), S r ,ij (1.73), S r2 (1.74) определим
pr2 º (1 - b)2 S r ,ij S r , ji + bwr ,i wr ,i
P 2r º p2r = (1 - b) S r2 + bW 2r .
и
Уравнения
(1.62-1.72)
справедливы, если заменить W 2r на P 2r . Тогда, средняя скорость диссипации
e = 4(27CP / 14)-3/ 2 r 2 P 3r ,
(4.7)
и, аналогично выводу уравнения (3.30), скорость диссипации
e = 4(27CP / 14) -3/ 2 M 3-1r 2 p3r .
(4.8)
Процедура моделирования скорости диссипации для разностной ячейки аналогична
процедуре, используемой для сферы, различие состоит в том, что отсутствует
аналитическое решение (1.68). Рассмотрим ячейку объёмом VD , поверхностью sD ,
нормалью к поверхности n и определим средние по объёму завихренность и скорость
деформации
wD ,i
SD ,ij
1
º
VD
1
º
VD
VD
VD
sD
e
ò0 wi ( x, t )dx = VikmD
1
ò0 Sij ( x, t )dx = 2VD
ò u ( x , t )n
k
m
sD
ò ( u ( x, t )n
i
0
dsD ,
0
j
+ u j ( x , t )ni ) dsD
(4.9)
44
и их интенсивности
wD2 º wD ,i wD ,i , S D2 º 2S D ,ij S D , ji .
(4.10)
Введем
p2D º (1 - b)S D2 + bwD2 , P 2D º (1 - b) S D2 + b w2D
(4.11)
и универсальную функцию
fD º e
=-
-1
h2
2VD2
nP 2D =
1
VD2
VD VD
òòf
R
(| x ¢¢ - x¢ |, CP )dVD¢dVD¢¢
0 0
sD sD
ò ò (n¢n¢¢) f
D
(| x ¢¢ - x ¢ |, CP )dsD¢ dsD¢¢ ,
(4.12)
0 0
которая может быть вычислена заранее для ячейки с фиксированной геометрией. Тогда
средняя скорость диссипации
e = nP 2D f D-1 .
(4.13)
Используем простейшую (однопараметрическую) монотонную аппроксимацию для
подсеточного моделирования:
f R (r , CD ) = F0 (r / r0 ) ,
(4.14)
где
F0 ( x) = 1 - (4 / 5) x 2 + (1/ 4) x 4 , x < 1 ;
F0 ( x) = (9 / 20) x -4/3 , x > 1 .
(4.15)
r0 º r0 / h = (100CD / 81)3/ 4
(4.16)
Параметр
находится из асимптотики при r >> 1 :
f R (r , CD ) = (5 / 9)CD r -4/3 . Эта аппроксимация
достаточно точна, если размер ячейки D лежит в инерционном диапазоне. Универсальная
функция f D (r , CD ) находится из (1.13):
f D (r , CD ) = r02 {(1/ 3)(r / r0 ) 2 - (2 / 25)(r / r0 ) 4 + (1/ 84)(r / r0 )6 } , r < r0 ;
f D (r , CD ) = r02 {(81/100)(r / r0 ) 2/3 - (2 / 3) + (64 / 525)(r / r0 ) -1} , r > r0 .
(4.17)
Асимптотика при r >> 1 : f D (r , CD ) = CD r 2/3 . Универсальная функция f r (r , CD ) = F1 (r / r1 )
определяется из (1.68), где:
F1 ( x) = 1 - (6 / 25) x 2 + (9 / 280) x 4 , x < 1 ;
45
F1 ( x) = (2187 /1400) x -4/3 - (256 /175) x -3 + (3 / 4) x -4 - (2 / 35) x -6 , x > 1 ;
r1 º r1 / h = (100CD / 81)3/ 4 / 2 = r0 / 2 .
(4.18)
f r (r , CD ) = (27 / 14)CD (2r ) -4/3 .
(4.19)
Асимптотика при r >> 1 :
Были вычислены зависимости и построены аппроксимации для разностной ячейки
– правильного шестигранника: -Di / 2 £ xi £ Di / 2 , i = 1, 2,3 ; D = (D1D 2 D 3 )1/3 ; индексация
такая, что D1 £ D 2 £ D 3 . Зависимости f D = f D (D1 , D 2 , D 3 , CD ) показаны на рис. 21, где
D i º D i / h . Аппроксимация
[ F (D / D ) F (D / D ) F (D / D )]
f D (D1 , D 2 , D3 , CD ) = 1 1 02 1 2 0 1 3 0
1 + 0.07 ln [ F1 (D1 / D 0 ) / F1 (D 3 / D 0 ) ]
1/3
(4.20)
является достаточно точной в широком диапазоне параметров, где D 0 = 0.922CD3/4 h , F1 –
функция (4.18).
fD
D2=D1
D2=D3
fD
100
100
D1=D3
D1=D3
D1=0.2D3
D1=0.2D3
D1=0.04D3
D1=0.04D3
10-1
10-1
10-2
10-2
10-3
100
101
102
103 D3/h
10-3
100
101
102
103 D3/η
Рис. 21. Зависимости f D = f D ( D1 , D 2 , D 3 , CD ) при CD =7.2, точки – интегралы (4.12),
линии – аппроксимации (4.20)
При D1 >> rDI (в инерционном диапазоне)
(D 0 / D) 4/3
2187
,
1400 1 + 0.124 ln 2 (D1 / D3 )
(4.21)
1400nP D2 D 4/3
=
1 + 0.124 ln 2 (D1 / D 3 )} .
{
4/3
2187 D 0
(4.22)
f D (D1 , D 2 , D3 , CD ) =
и, из (4.13),
e = nP f
2
D
-1
D
4/3
Подставив в это уравнение D 4/3
CP n e
0 = 0.922
-1/3
, находим
46
e =
3/2
0.1762
1 + 0.124 ln 2 (D1 / D 3 )} D 2 P 3D ,
3/ 2 {
(CP / 7.2)
(4.23)
и, аналогично выводу уравнения (3.30), скорость диссипации
e=
3/ 2
0.1762
1 + 0.124ln 2 (D1 / D3 )} D 2 (p2D )3/2 .
3/ 2 {
M 3 (CP / 7.2)
(4.24)
Итак, определив эффективный коэффициент вязкости n LES º e / SD2 , находим
e = n LES SD2 ,
(4.25)
где
n LES = CD2 D 2 ( S D2 )1/ 2 ,
(4.26)
3/4
2
2
3/ 4 ì (1 - b) S + bw ü
0.176
2
D
D
CD = 1/2
1
+
0.124
ln
(
D
/
D
)
{
}
í
ý
1
3
M 3 (CP / 7.2)3/ 4
S D2
î
þ
(4.27)
в инерционном диапазоне. Напомним, что D1 / D 3 – отношение минимального к
максимальному размерам ячейки.
Для коррекции этого выражения на весь универсальный диапазон используем
(1.72), тогда
n LES = g r (n -1D 2 pD , CP ) {n + CD2 D 2 ( S D2 )1/ 2 } .
(4.28)
Для определения функции g r вначале найдем из (1.68, 1.41 и 2.2) функцию
f r (r / 2, CD ) =
36 ì
1
f (r ) - 2
2 í U ,11
r î
r
r
ò
0
fU ,11 (r )rdr -
2
r4
r
òf
U ,11
0
ü
(r )r 3dr ý ,
þ
её график приведен на рис. 22 при CD , LL = 1.97 , CD = (11/ 3)CD , LL , a11 = 3.2 и b11 = 0 .
(4.29)
47
r
4/3
7
fr ( r , CD )
6
5
4
3
2
1
0
10-4
10-2
100
102
104
106 r
Рис. 22. График f r (r , CD )
Используя f r (r , CD ) определим g r (W D , CD ) из (1.71), эта функция и её аппроксимация
{
g r (W r , CD ) = 1 - 0.355exp -0.06 éë0.1ln(W r / 138) + 0.9 ln 2 (W r / 138) ùû
2
}
представлены на рис. 23.
g r ( W, C D )
1
0.9
0.8
0.7
0.6
10-5
100
105
1010 W
Рис. 23. График g r (W, CD ) ; точки – (1.71), линии – (4.30)
Окончательно, g r , CD , n LES и e могут быть получены из (4.30, 4.27, 4.28 and 4.25).
(4.30)
48
Приложение 1. Структурные функции третьего порядка
В этом параграфе рассмотрено только изотропное поле скорости. Определим
двухточечные структурные функции
Dijk ( x¢¢ - x ¢) º (ui ( x ¢¢, t ) - ui ( x ¢, t ))(u j ( x ¢¢, t ) - u j ( x ¢, t ))(uk ( x ¢¢, t ) - uk ( x ¢, t )) .
(5.1)
Из симметрии Dijk (r ) = D jik (r ) = Dikj (r ) = Dkji (r ) и изотропности следует, что
Dijk (r ) = ( DLLL (r ) - 3DLNN (r ) )
ri rj rk
r
3
+ DLNN (r )
(d r
ij k
+ dik rj + d jk ri )
r
.
(5.2)
В системе координат, в которой r = (r , 0, 0) , ненулевые компоненты:
D111 = DLLL ; D122 = D133 = D221 = D331 = D212 = D313 = DLNN .
(5.3)
Свертки (5.2):
ri ¶Dijj (r ) 1 dr 2 ( DLLL + 2 DLNN )
Dijj (r ) = ( DLLL + 2 DLNN ) ;
= 2
;
r
r
dr
¶ri
¶Dijk (r )
¶ri
¶ 2 Dijk (r )
¶ri ¶rj
ì d ( DLLL - DLNN )
ü r r ì dD
üd
= ír
+ 2 DLLL - 8 DLNN ý j 3k + ír LNN + 4 DLNN ý jk ;
dr
î
þ r
î dr
þ r
3
ì d æ drDLLL
ö ü r ¶ Dijk (r ) 1 d ì 1 d 2 æ drDLLL
öü
= 2
- 6 DLNN ÷ ý . (5.4)
= í r2 ç
- 6 DLNN ÷ ý k4 ;
r ç
í
¶ri ¶rj ¶rk r dr î r dr è dr
øþ r
øþ
î dr è dr
В общем случае, статистические моменты третьего порядка поля скорости w (r )
должны быть представлены структурными функциями
Dijk (r ¢, r ¢¢, r ¢¢¢)
º ( ui (r ¢ + x0 , t ) - ui ( x0 , t ) ) ( u j (r ¢¢ + x0 , t ) - u j ( x0 , t ) ) ( uk (r ¢¢¢ + x0 , t ) - uk ( x0 , t ) ) .
Очевидно, что
Dijk (r ¢, r ¢¢, r ¢¢¢) = D jik (r ¢¢, r ¢, r ¢¢¢) = Dkji (r ¢¢¢, r ¢¢, r ¢)
и
(5.5)
Dijk (r ) = Dijk (r , r , r ) . Из
равенства
Dijk (r ¢¢ - r ¢)
= (ui (r ¢¢ + x0 , t ) - ui (r ¢ + x0 , t ))(u j (r ¢¢ + x0 , t ) - u j (r ¢ + x0 , t ))(uk (r ¢¢ + x0 , t ) - uk (r ¢ + x0 , t )) , (5.6)
используя процедуру вывода уравнения (1.17), находим
Dijk (r ¢¢ - r ¢) = - Dijk (r ¢) + Dijk (r ¢¢) - Dijk (r ¢, r ¢¢, r ¢¢) + Dijk (r ¢¢, r ¢, r ¢) + Dijk (r ¢, r ¢, r ¢¢) - Dijk (r ¢¢, r ¢¢, r ¢)
+ Dijk (r ¢, r ¢¢, r ¢) - Dijk (r ¢¢, r ¢, r ¢¢) .
(5.7)
49
Из инвариантности относительно преобразования x% = r ¢ + r ¢¢ - x получим
Dijk (r ¢¢, r ¢, r ¢) = - Dijk (r ¢, r ¢¢, r ¢¢) , Dijk (r ¢¢, r ¢¢, r ¢) = - Dijk (r ¢, r ¢, r ¢¢) ,
Dijk (r ¢¢, r ¢, r ¢¢) = - Dijk (r ¢, r ¢¢, r ¢) ,
(5.8)
отсюда,
- Dijk (r ¢, r ¢¢, r ¢¢) + Dijk (r ¢, r ¢, r ¢¢) + Dijk (r ¢, r ¢¢, r ¢) = (1/ 2) éë Dijk (r ¢¢ - r ¢) + Dijk (r ¢) - Dijk (r ¢¢) ùû .
(5.9)
Для несжимаемой жидкости ¶Dijk (r ¢, r ¢¢, r ¢¢) / ¶ri¢ = ¶Dijk (r ¢, r ¢¢, r ¢) / ¶rj¢¢ = ¶Dijk (r ¢, r ¢, r ¢¢) / ¶rk¢¢ = 0 .
Дифференцируя (5.9) по ri¢ , rk¢¢ и rj¢¢ , имеем
¶ 3 Dijk (r )
¶ 3 Dijk (r )
== 0 , r = r ¢¢ - r ¢ ,
¶ri¢¶rk¢¢¶rj¢¢
¶ri ¶rk ¶rj
(5.10)
отсюда и из последнего уравнения (5.4) находим 6 DLNN = d (rDLLL ) / dr + CLNN . Константа
CLNN равна нулю из-за нечётности функций DLLL (r ) и DLNN (r ) , окончательно,
DLNN =
1 drDLLL
.
6 dr
(5.11)
Приложение 1.1. Подход Кармана-Ховарта
Уравнения Навье-Стокса для несжимаемой жидкости имеют вид
¶u( x, t )
1
+ (u( x , t )Ñ)u( x , t ) + Ñp ( x , t ) = nDu( x , t ) .
¶t
r
(5.12)
Рассмотрим произвольную точку x0 и определим u0 º u( x0 , t ) и w (r , t ) º u(r + x0 , t ) - u0 .
Введем r ¢ º x¢ - x0 , u¢ º u(r ¢ + x0 , t ) и w ¢ º u¢ - u0 для точки x ¢ ; для точки x ¢¢ используем
аналогичные обозначения. Записав (5.12) для точки x ¢ , дифференцируя e jnm ¶ / ¶xm¢ и
умножая на w¢¢k , получим
¶ 2 w¢j
¶ 2ui¢wn¢
¶un¢ ¶ui¢wn¢ 1 ¶p¢
¶ 2un¢ ¶w¢j
+
+
=n
;
+ e jnm
=n
;
¶t
¶xi¢
r ¶xn¢
¶xq¢ ¶xq¢ ¶t
¶xi¢¶xm¢
¶xq¢ ¶xq¢
w¢¢k
¶w¢j
¶t
+ e jnm e klp
¶ 2 w¢j w¢¢k
¶ 3ui¢wn¢ wl¢¢
=n
.
¶xi¢¶xm¢ ¶x¢¢p
¶xq¢ ¶xq¢
(5.13)
Записав аналогичное последнему уравнению (5.13) уравнение для точки x ¢¢ , суммируя,
осредняя и используя (1.45), найдем
50
e jnm e klp
ü
¶ 2 Dnl (r , t )
¶ 2 ì 1 ¶Dnl (r , t ) ¶Dinl (r ¢, r ¢, r ¢¢, t ) ¶Dinl (r ¢¢, r ¢, r ¢¢, t )
n
+ A nl ý = 0 , (5.14)
í
¶xm¢ ¶x¢¢p î 2
¶t
¶ri¢
¶ri¢¢
¶rm ¶rm
þ
где r º r ¢¢ - r ¢ , Anl – анизотропная часть уравнений. Очевидно, что критерий изотропности
æ ¶
¶ ö
Anl º - ui0 ç
+
÷ wn¢ wl¢¢ = 0 .
¢
¢¢
¶
x
¶
x
i
i
è
ø
(5.15)
Уравнения (5.14) описывают процессы изменения статистики и можно использовать
формализм теории Колмогорова, если характерные времена этих процессов значительно
превосходит t IE . Из инвариантности относительно преобразования x% = r ¢ + r ¢¢ - x следует
¶Dijk (r ¢, r ¢¢, r ¢, t ) ¶Dijk (r ¢¢, r ¢, r ¢¢, t )
=
.
¶ri¢
¶ri¢¢
(5.16)
С учётом этого равенства, дифференцируя ¶ / ¶ri¢ (5.9), находим
¶Dijk (r ¢, r ¢, r ¢¢, t ) ¶Dijk (r ¢¢, r ¢, r ¢¢, t )
1 ¶Dijk (r , t ) 1 ¶Dijk (r ¢, t )
+
=+
,
¶ri¢
¶ri¢¢
2
¶ri
2
¶ri¢
(5.17)
и, из (5.14),
e jnm e klp
¶ 2 Dnl (r , t ) ü
¶ 2 ì ¶Dnl (r , t ) ¶Dinl (r , t )
+
- 2n
í
ý=0.
¶rm¢ ¶rp¢¢ î ¶t
¶ri
¶rm ¶rm þ
(5.18)
Поскольку выражение в скобках (5.18) – чётная функция, следовательно,
¶D jk (r , t )
¶t
+
¶Dijk (r , t )
¶ri
= 2n
¶ 2 D jk (r , t )
¶rm ¶rm
-
4CLLL e
d jk ,
3
(5.19)
где CLLL – универсальная константа.
Свертка (5.19), с использованием (5.4), дает
¶Dii 1 ¶r 2 ( DLLL + 2 DLNN ) 2n ¶ 2 ¶Dii
+ 2
= 2
r
- 4CLLL e .
¶t
¶r
¶r
r
r ¶r
(5.20)
Из (5.11) получаем
DLLL + 2 DLNN =
1 ¶r 4 DLLL
,
3r 3 ¶r
(5.21)
и, используя первое уравнение (1.41),
4
¶r 4 DLL 1 ¶r 4 DLLL
¶ 1 ¶r 3 DLL 4CLLL e r
+
= 2n r 3
,
3 ¶r
3
¶t
¶r r 2 ¶r
отсюда,
(5.22)
51
r
r
3
4CLLL e r 5
¶ 4
r 4 DLLL
3 ¶ 1 ¶r DLL
2
r
D
dr
+
=
n
r
dr
.
LL
ò0 ¶r r 2 ¶r
3
15
¶t ò0
(5.23)
Интегрируя по частям, находим
r
3
òr
0
r
¶ 1 ¶r 3 DLL
¶r 3 DLL
¶r 3 DLL
¶D
=
3
dr
r
dr = r 4 LL .
2
ò
¶r r
¶r
¶r
¶r
¶r
0
(5.24)
r
ü
¶ ì3 4
¶DLL 4
- CLLL e r .
í 4 ò r DLL dr ý + DLLL = 6n
5
¶t î r 0
¶r
þ
(5.25)
Окончательно,
Для стационарной статистики подход Кармана-Ховарта позволяет определить
моменты третьего порядка:
DLLL (r ) = 6n
dDLL (r ) 4
- CLLL e r ,
5
dr
(5.26)
DLNN (r ) определяется из (5.11). При r << h структурная функция DLL = (1/ 15) e n -1r 2 ,
следовательно, DLLL = (4 / 5)(1 - CLLL ) e r . Колмогоров предположил, что при r << h
можно пренебречь DLLL в уравнении (5.26), в результате, CLLL = 1 и, таким образом,
DLLL (r ) = -(4 / 5) e r при r >> h.
Литература
[1] Колмогоров А.Н. Избранные труды. Математика и механика. М.: Наука. 1985.
[2] Колмогоров А.Н. Локальная структура турбулентности в несжимаемой вязкой жидкости
при очень больших числах Рейнольдса // Доклады АН СССР. 1941, т.30. № 4. С. 299-303.
[3] Saddoughi S.G., Veeravalli S.V. Local isotropy in turbulent boundary layers at high Reynolds
number // J. Fluid Mech. (1994), vol. 268, pp. 333-372
1/--страниц
Пожаловаться на содержимое документа