Численное моделирование турбулентного течения в канале Текст научной статьи по специальности «Математика»

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

Том 12, № 5, 2007

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ТУРБУЛЕНТНОГО ТЕЧЕНИЯ В КАНАЛЕ

У. С. АБДИБЕКОВ, Н. Б. УСЕНБАЕВ, О. Л. КАРУНА Казахский национальный университет им. аль-Фараби, Алматы

e-mail: [email protected]

Turbulent flow in a channel is studied using Large Eddy Simulation method applied to assemble averaged Navier—Stoks equations. Fourier method and Chebyshev transformation were used. Smogorinsky and dynamic models were employed as sub-grid models. Averaged and pulsatio characteristics of the flow were obtained.

Введение

В настоящее время для численного моделирования проблем турбулентности используются три подхода: прямое численное моделирование, метод осреднения по Рейнольдсу и метод крупных вихрей [1-9]. Прямое численное моделирование предполагает разрешение вихрей всех масштабов, что не реализуемо на существующих ЭВМ средней мощности. Метод Рейнольдса позволяет моделировать только очень крупные вихри, соизмеримые с физической областью, все среднемасштабные и мелкомасштабные вихри исчезают при осреднении. Наиболее продуктивным является метод крупных вихрей, так как он позволяет моделировать и крупные вихри, и вихри среднего размера, вплоть до размера ячейки расчетной сетки, а вихри мельче расчетной ячейки моделируются с помощью различных гипотез.

1. Основные уравнения

Процессы турбулентного переноса представляют собой сложное физическое явление, теоретическое изучение которого опирается на основные законы физики, и описываются уравнениями гидродинамики. Основными уравнениями гидродинамики для описания турбулентных течений являются уравнения Навье—Стокса [10]:

— + — (ии) = - — + — 9 2щ (1)

дЬ дх2 % 3 дхл И,е дх3 дх3-'

дил . ,

дХ =(2)

Система уравнений записана в декартовой системе координат хл (г = 1, 2, 3), в физическом пространстве; три компоненты скорости и и давление, р; Ь — безразмерное

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2007.

время, И,е — число Рейнольдса. В системе уравнений (1)-(2) и далее по повторяющимся индексам следует производить суммирование.

Для выделения основных энергонесущих вихрей используется следующая форма осреднения по пространству:

и (х, Ь) = / С (г, х) и (х — г, Ь) ¿г,

(3)

V

где и = («1, и2, и3), черточка сверху означает осреднение по объему; г = (г1, г2, г3) — вектор координатных точек, по которому производится интегрирование; V — объем вычисляемой области и интегрирования; С (г, х) — функция-фильтр. В качестве функции-фильтра чаще используется известный гауссовский фильтр вида

3 / 6 V/2 С (г) = С (г! = и ехр

6т2

д2

(4)

а также "фильтр-ящик":

С (г") = Д3

0,

д

|г| ^ 2,

|П| > Д, Уг = 1, 2, 3,

(5)

где Д — характерная длина фильтра, имеющая порядок размера ячейки сетки. Обычно у нее такой вид:

д = (й1№)1/3

2\1/2

д = (й2 + ^2 + ь3)

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

Применяя операцию осреднения (3) с фильтрами типа (4) и (5) к уравнениям (1) и (2), можно вывести следующую систему уравнений [1-5]:

ди, + д (щи.;)

дЬ

дх,

др 1

дх,

д2 иг

дт,

г]

И,е дх; дх;

дх;

(6)

ди, дхг

т. ит, иг..

Подсеточный член тг; отвечает за мелкомасштабную турбулентность. Для его определения используется модель Смагоринского:

5] о

тг]--3 ткк = —2^8СЯ8

¿31

где турбулентная подсеточная вязкость представляется Уяая = (СяД) (28г]8гз)

2/0* РГ \ 1/2

1 (диг + ди3

8,3 2 1 дх,

дхг

— тензор скоростей деформации; 5

3

г3

0, г = з

— знак Кроне-

кера; Ся — коэффициент Смагоринского, который лежит в отрезке 0.06-0.25 [1].

1

0

а

1

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

Для применения динамической модели проводится двойное осреднение с длиной фильтра А = 2А, тогда

дщ + д (ищ) = _ др_ + д2щ _ дт, _ д (щщ - щщ) дЬ дх, дхг Ив дх, дх, дх, дх,

Уравнение (1), подвергнутое осреднению с двумя фильтрами длиной А и А соответственно, имеет следующий вид:

дщ + д (щщ) = - + ^ д2щ - дТ, дЬ дх, дхг Не дх, дх, дх

где

I г, — щг *

Из (7) и (8) следует

Тг, + и^и, и^и,, и напряжения Леонарда выражаются

г.. = т.. — Т-, , ,

Тогда и Т, имеет следующий вид:

Т,--3Ткк = — 2 (С3А) (2з,3,) ^ 5г, *

А напряжения Леонарда имеют вид

с __

Г, - 3Ькк = -2 (Сз)2 [(Д 2 (25Д,-)1/2 - (А)2 (Щз,)1/2 * (9)

Из (9) при использовании метода наименьших квадратов находится значение Сз в виде

С2 = -1 И,

2 Ыгк Ыгк

Г___2 , _ _ ) 1 /2 _ <-,

где И,

(А) (2зг,3ч) 3г, (А) (2зг,3,) / 3,

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

Для решения задачи с учетом вышепредложенных моделей турбулентности рассматриваются уравнения турбулентного движения в канале в следующем виде:

дщ + д (щ щ) = - д 2щ - дт, - , (ю)

дЬ дх, дх г Ив дх, дх, дх, дхг 1 г

где последнее слагаемое в (10) — средний градиент давления, т.е. полагается, что дви

дР

жение в канале осуществляется за счет перепада давления. Здесь задается

дх1

2

«2, где

ит — скорость трения на стенке. Число Рейнольдса определяется за счет максимальной скорости на середине канала итах и половины длины канала Н. Область канала имеет следующие размеры: Ь1 х Ь2 х Ь3 = 5Н х 2Н х 3Н. Граничные условия задаются периодическими в направлении осей х и г, а по у выполняется условие стенки. Схематически область канала представлена на рис. 1.

Для определения скорости трения ит на стенке нужно воспользоваться известными эмпирическими отношениями между обычным числом Рейнольдса Ие и числом Рей-

нольдса Иет. Согласно формуле

ит

Иет

Иет

180 приблизительно соответствует

итах Ие

Ие = 3381 [2, 6]. Для численного решения задачи применяется явный метод расщепления по физическим процессам [11]. Уравнение (10) интегрируется по времени без слагаемого, содержащего давление, т. е. первый этап состоит в вычислении, где определяется промежуточное значение скорости и*

и* = и" + Аг

1

д 2ип

Ие дх, дх,

д_ дх,

дт™

ииг - дХ

дР дх1

$м

(12)

где верхний индекс п означает определенный номер временного уровня. После вычисления (12) из требования выполнения уравнения неразрывности находится давление Значения скоростей следующего уровня п +1 определяются поправкой

ипг+1 = и* - Аг.

г г дх,

(13)

Выражение (13) дифференцируется по переменной хг:

дигг+1 ди* д 2ргг+1

дх,

дх,-

дх2

(14)

Сложив (14) для всех г = 1, 2, 3 с учетом, что должно выполняться уравнение неразрывности (11) на каждом следующем временном уровне, получаем следующее трехмерное уравнение Пуассона для давления:

+ 32Г+1 +

дх2 дх2 дх3

ди* ди* ди3 1 + тт^ + 3

х1 х

2

х

3

(15)

Рис. 1. Трехмерное течение между двумя неподвижными бесконечными пластинами

Из решения уравнения давления (15) вычисляются градиенты давления, которые затем подставляются в (13) для расчета конечного значения скорости

Чтобы найти производные в направлении осей х и г, для которых заданы периодические граничные условия, применяется спектральный метод Фурье, а по у — используется спектральный метод, основанный на полиномах Чебышева. При использовании спектрального метода Фурье применяется быстрое преобразование Фурье, позволяющее сократить количество операций для вычисления производных. Рассмотрим решение уравнения Пуассона (15)

д2р д2р д 2р

дх"! + дх\ + дх\ (16)

где / — известная функция (15). Спектральный метод Фурье позволяет представить давление и правую часть (16) в виде

Р = р(к 1, у, кг) ег(к1Х+к^, (17)

к\ кз

f = /(кьУ,кэ) вг(к1Х+к^, (18)

к-1 кз

где р, / — моды соответствующих функций, а к — частота, выражаемая к = —п

Lг

у с ( N N Л ы

чп € I ——, —--II, причем п является целочисленным, N — количество точек в

направлении соответствующей оси координат. Подставляя (17) и (18) в (16), умножая обе стороны на е-г(к1х+кз*) и учитывая ортогональность данной функции, получим

- (к? + к32) р + || = - (к2 + к|) /+ /, Р> = р(к1,Ш Ук1,кз,3. (19)

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

д'р

д— (г.,к) = ^^ О (.,т) р(г,т,к), (20)

3 т=1

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

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

др(г, I, к) р (г, 0,к) - 2р(г, 1,к)+ р(г, 2, к) -2р(г, 1,к) + 2р(г, 2, к)

др2 (Нг (I))2 (Нг (I))

д'р (г, N3, к) р (г, N2 + 1,к) - 2р(г,1з,к) + р^,! - 1,к)

(21)

дх2 (Нг !))2

= -2р (г, N2, к) + 2р (г, N2 - 1,к)

(Нг т)2 , ( )

Для достижения большей точности вычисления давления по третьей координате используется неравномерная сетка следующего вида: y (j) = cos [п (j — 1) /N2].

Система уравнений (19) с учетом (20)-(22) представляется в виде матричной формы по переменной y

A\,k3 =bU Vk1 (23)

и решается методом декомпозиции LU.

После нахождения мод давления p из (23) обратным преобразованием Фурье вычисляется давление в физическом пространстве с применением алгоритма быстрого преобразования Фурье.

Отметим, что при применении модели Смагоринского применялась демпфирующая функция Ван Дриста [7] около стенок

У+

Cs = Cs\ 1 — exp где y+ — пристеночная функция.

25

3. Численные результаты

Для решения задачи заданы вычислительная сетка размерами 25 х 25 х 25, а коэффициент CS для модели Смагоринского равным 0.13. Меньшие значения коэффициента Смагоринского приводили расчет к неустойчивости. Начальное поле задано случайным образом и доводилось до состояния, когда профиль осредненной скорости не менялся по времени, это заняло свыше 2500 безразмерных единиц времени. Далее проводилась запись мгновенных значений получаемых скоростей для обеих моделей метода крупных вихрей через каждые 0.5 единиц безразмерного времени. Статистическая обработка значений вычисляемых скоростей осуществлялась на промежутке 300 безразмерных единиц времени.

Для сравнения полученных численных расчетов использовался численный эксперимент, проведенный Дж. Кимом, П. Моином и Р. Мозером [7]. В их работе рассматривалось прямое численное моделирование в канале с разрешением сетки 192 х 129 х 160, т. е. применялся метод прямого численного моделирования. Течение происходило между двумя неподвижными пластинами (ReT = 180). Граничные условия заданы так же, как

б

а

Рис. 2. Профиль средней скорости, отнесенный к ит: а — модель Смагоринского; б —динамическая модель; о — численное решение; пунктир — и+ = у+; сплошная линия — и+ = 2.5 1пу+ +5.5

и в данной работе. Полученный результат сравнивался с экспериментальными данными Креплина и Экельмана. Получено удовлетворительное согласование.

На рис. 2 представлен профиль средней скорости верхней части канала для расчетов с применением двух разных моделей. В нижней части канала прослеживается симметричность профилей скорости для обоих расчетов. Видно, что вблизи стенки выполняется известный линейный закон, а далее выполняется логарифмический закон.

На рис. 3 представлена важная турбулентная характеристика — турбулентная интенсивность, или среднеквадратичное отклонение от соответствующих средних скоростей. Видно, что согласование с работой [6] удовлетворительное, но динамическая модель

Рис. 3. Турбулентные интенсивности, отнесенные к ит: а — модель Смагоринского; б — динамическая модель; о — игт8; □ — шгт8; Д — игт8; сплошная линия — игт8 [6]; штриховая линия — и>гт8 [6]; пунктир — -игт8 [6]

Рис. 4. Сдвиговое напряжение Рейнольдса, отнесенное к u2: а — модель Смагоринского; б — динамическая модель; —(u'1u2)t: □ — расчет; сплошная линия — расчет DNS [6]; —(u'1u2)t +

L: о — расчет; штриховая линия — расчет DNS [6]

1 d (ui)t.

Re dx2

а

а

дает более точный результат. Величины wrms и vrms модели Смагоринского принимают значения намного ниже, чем те же wrms и vrms работы [6], а динамическая модель дает более точный результат.

На рис. 4 представлены моменты второго порядка сдвигового напряжения Рейнольд-са. Согласование достаточно удовлетворительное. Однако результат, полученный по динамической модели, точнее, чем данные, полученные по модели Смагоринского. В последнем случае наблюдается сильное отклонение от сдвигового напряжения Рейнольдса работы [7] вблизи стенок.

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

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

[1] Deardorff J.W. A numerical study of three-dimensional turbulent channel flow at large Reynolds numbers // J. Fluid Mech. 1970. Vol. 41. P. 453-480.

[2] Ферцигер Дж. Х. Численное моделирование крупных вихрей для расчета турбулентных течений // Ракетная техника и космонавтика. 1977. Т. 15, № 9. С. 56-65.

[3] Lesieur M., Metais O. New trends in large-eddy simulations of turbulence // Annu. Rev. Fluid. Mech. 1996. Vol. 28. P. 45-82.

[4] Davidson L. Lecture notes on large eddy simulation. Goeteborg: Chalmers Univ. of Tech., 2000. SE 412 96. P. 1-18.

[5] Germano M., Piomelli U., Moin P., Cabot W. A dynamic subgrid-scale model eddy viscosity model // Phys. Fluids A. 1991. Vol. 3, N 7. P. 1760-1765.

[6] ДжЕнкинс Г., Влттс Д. Спектральный анализ и его приложения. М.: Мир, 1972. Т. 2. 288 с.

[7] Kim J., Moin P., Moser R. Turbulence statistics in fully developed channel flow at low Reynolds number //J. Fluid Mech. 1987. Vol. 177. P. 133-166.

[8] Moin P., Kim J. Numerical investigation of turbulent channel flow //J. Fluid Mech. 1982. Vol. 118. P. 341-377.

[9] Meinke M., Schroeder W., Krause E., Rister T. A comparison of second- and sixth-order methods for large-eddy simulations // Computer&Fluids. 2002. Vol. 31. P. 695-718.

[10] Лойцянский Л.Г. Механика жидкости и газа. М.: Наука, 1987. C. 373-382.

[11] Флетчер К. Вычислительные методы в динамике жидкости. М.: Мир, 1991. Т. 1. 502 c.

Поступила в редакцию 24 апреля 2007 г., в переработанном виде — 25 мая 2007 г.