Об одном эффективном методе вычисления звуковых полей, генерируемых турбулентными потоками жидкости Текст научной статьи по специальности «Математика»

УДК 532.5; 535.4

ОБ ОДНОМ ЭФФЕКТИВНОМ МЕТОДЕ ВЫЧИСЛЕНИЯ ЗВУКОВЫХ ПОЛЕЙ, ГЕНЕРИРУЕМЫХ ТУРБУЛЕНТНЫМИ ПОТОКАМИ ЖИДКОСТИ

© 2009 г. М.А. Сумбатян1, Э.А. Троян

2

1Южный федеральный университет, ул. Мильчакова, 8а, Ростов н/Д, 344090, dnjme@math. sfedu.ru

2Научно-исследовательский институт механики и прикладной математики Южного федерального университета, пр. Стачки, 200/1, г. Ростов н/Д, 344090, niimpm@ns.math.rsu.ru

1Southern Federal University, Milchakov St., 8a, Rostov-on-Don, 344090, dnjme@math. sfedu.ru

2Research Institute of Mechanics and Applied Mathematics of Southern Federal University, Stachki Ave, 200/1, Rostov-on-Don, 344090, niimpm@ns.math.rsu.ru

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

Ключевые слова: турбулентный поток, неявная схема, функция Грина, звуковые поля.

We develop a numerical calculation of the oscillating components of turbulent velocity field in a channel of constant width, on the basis of the Navier-Stokes equations. On the basis of this velocity field there are calculated generated sound fields.

Keywords: turbulent flow, implicit scheme, Green's function, sound waves.

Введение

Математическая постановка задачи

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

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

Рассмотрим классическую задачу об однородном турбулентном потоке несжимаемой вязкой жидкости в канале постоянной ширины (двумерный случай). Уравнения движения Навье-Стокса в терминах пары завихренность - функция тока имеют следующий вид [1, 2]:

— = -и—1- - v—1- + vA^, <^ = A\ . dt 5x dy

(1)

При этом компоненты вектора скорости частиц жидкости u(x,y,t), v(x,y,t) связаны с функцией тока у(х,у,$ и завихренностью ^(х,у,/) следующим образом:

_ d\ _ d\ _ du dv dy dx dy dx

(2)

Для определенности будем считать, что ось x направлена вдоль канала (—»<X<<»); поперечная вертикальная координата изменяется на интервале -Л<у<Л. При этом ширина канала равна 2Л.

Если средний расход жидкости по сечению н

Q = |и(х, у)йу = у |у=н - у |у=_н и не зависит от х и ^ -н

то средняя скорость потока Um=Q/(2h) является некоторой известной постоянной величиной. Кроме того, выделение простейшего решения, обеспечивающего заданный расход и условия прилипания, позволяет пе-

реписать уравнения (1), (2) для новых функций, удовлетворяющих однородным краевым условиям на стенках канала:

3

У = У0 + т б

(

У

3h3

тл С = С°-ТЯУг -

~ ' ^ дх дх

3Q

д£_дС°___

ду ду 2h

о 3Q u = u0 + ^

4h

2 h3 1 - hr

(

V = V

(3)

дДу0 i2 0

дГ

У

у=± h

ду ду

= 0,

у=±h

f = -

о 3Q u0 + — 4h

у2 Т д^0 - V0 3Q"

h2 J дх _ ду 2h3 _

Для больших чисел Рейнольдса Re = U„fi/\ течение является турбулентным, и точный расчет потока, обеспечивающий корректное вычисление осциллирующих турбулентных составляющих скорости, становится затруднительным. Между тем информация об осциллирующих компонентах скорости является ключевой именно в аэроакустике [9]. В связи с этим к решению нелинейной системы дифференциальных уравнений (1), (2) применялись различные прямые численные методы, эффективные как в ламинарном, так и в турбулентном течениях.

Одним из известных методов численного решения является метод, основанный на классических итерациях по времени. При этом на каждом шаге итерации задача сводится к некоторой линейной эллиптической проблеме по пространственным переменным. Цель данной работы - построение замкнутого решения данной эллиптической задачи для канала постоянной ширины. Метод построения основан на использовании функции Грина, которая выписывается явно в квадратурах.

Задача (1), (2) по переменной t является классической задачей Коши, если считать известными начальные условия для скоростей. Обычно в развитом потоке такие условия редко бывают заданы a priori, однако численные эксперименты показывают, что выбор тех или иных начальных условий не меняет качественного характера течения. Известно (см., например, [1, 2]), что устойчивое численное решение можно получить различными итерационными процессами по времени, если использовать конечно-разностную аппроксимацию производной по времени в варианте «назад». Простейшей является схема Эйлера (э^0/dt)n )Д, где т - шаг по времени.

Подстановка этого соотношения в (3) с выбором функции f содержащей нелинейные члены на предыдущем временном слое, а члена со старшей 4-й производной от функции у - на новом временном слое, приводит к эллиптической краевой задаче 4-го порядка относительно функции тока

V = gn-1, gn-1 = C-i + fi. (4)

Такая схема обладает линейной скоростью сходимости с уменьшением временного шага т. Квадратичную скорость сходимости по т дает применение схемы Кранка-Николсона вместо схемы Эйлера [1]. Это приводит к эллиптической задаче, лишь незначитель-

но отличающейся от (4): -У^Д2у° = g*-1.

£-1=^0—1.

Наконец, еще более высокий порядок точности итераций достигается применением схемы Кранка-Николсона для вязкого члена, Адамса - нелинейного [6]:

. 0 У^ а 2 0 **

ДУ0 - — Д У0 = Я„-ь

(5)

80-1 = ^0—1 +уЛ£0_1 ] + ^(3/о-1 - /о-2 ).

Заметим, что все три варианта итерационного метода можно записать единым образом

ДУо0 -£ДУ0 = q,

У

у=± h

ду ду

= 0, (6)

у=± h

где е > 0 - некий малый параметр при старших производных; q - некая функция, известная из предыдущих итераций.

Построение функции Грина для завихренности

В случае канала постоянной толщины уравнение (5) допускает явное точное решение. Для этого сначала построим функцию Грина, т.е. решение краевой задачи

ЛО — е Л2О = д(£ — х)д(р — у)

Gl = —

= 0,

(7)

q=±h

где G=G(£¡,ц,xy); оператор Лапласа берется по переменным (£п); 8(х) - дельта-функция Дирака.

Применение преобразования Фурье по переменной ¿(^ ^ а) сводит (7) к линейному обыкновенному дифференциальному уравнению 4-го порядка с постоянными коэффициентами (образы Фурье обозначены знаком тильда)

й 4° (1 + 2еа 2 )+ (а2 + еа 4 )<5 = -е'ахд{т1 — у),

d-

dq

G\

_cG q=±h дq

= 0.

(8)

q=±h

Частное решение неоднородного уравнения (8) легко строится методом разделения переменных в рядах:

(9)

Gp =

(

h

i

i

Л

ßl +«12 ßl +«22

7ппу . jnnq sin —— sin -

2h

2h

- 1 — 1 П Т7 п кт

у = у + И, 1=1 + И, а1 =а, а2 =у!а + /е, рт = — .

После суммирования некоторых табличных рядов (9) может быть переписано в виде

сй[(2й — 1 — у| )а1 ]— сЩц + у)а1 ]

a1sh (lha1)

(10)

3

о

сН[(2Н - \\ - У)я2 ]- сН[(\ + у)«2 ] а2 ^ (2На2 )

Добавляя к (10) общее решение однородного уравнения, получаем

О = О р + С^Н [а (\ - Н)]+ С^Н [а (\ + н)]+

+ С3 sН [а2 (\ - Н)]+С4 sН [а2 (\ + Н)].

Неизвестные постоянные С1 - С4 должны быть найдены из 4 граничных условий (8). Заметим, что частное решение (9), (10) специально строилось так, чтобы автоматически удовлетворялись условия

Ор (\ = +Н) = 0 . Кроме того, специальный выбор

структуры решения однородного уравнения приводит к тому, что в возникающей линейной алгебраической системе размерности 4*4 для С1 - С4 из 16 элементов матрицы 4 оказываются нулевыми. Кроме того, из 4 элементов вектора правой части 2 оказываются нулевыми. В результате все 4 неизвестных коэффициента легко выписываются в явном виде, что после некоторых преобразований приводит к следующему окончательному представлению для образа Фурье функции Грина: О(\, х, у)= е,ахИ,

сН[(2Н - \ - У }х1 ] - сН[(\ + у)а ]

H V У )=■

2a1sh (2h а1)

+

+

сН[(2Н - \\ - у|)«2]- сН[(\ + у] 2а 2 яН (2На2 ) + |[а2ь'Н(2На1) - а^Н(2На2 )] х

х [Оо (у)Оо (-\) + Оо (-у)Оо (\)] +

[а2^Н (2На1 )сН(2На2 )- а^Н (2На2 )сН(2На1)] х х[Оо (у)Оо (\) + Оо (-у)Оо (-\)]}-^ ,

М)

О0{у)= Н(Н + у)а1 ] [(Н + у)а2 ]

яН (2На1) яН (2На2 ) = (а12 + а|)£Н (2На1 )уН (2На2 ) + + 2а1а2 [1 - сН(2На1 )сН(2На2 )]. Очевидно, сама функция Грина получается обращением Фурье выражения 1

О; д, х,у) = — | ео8[а(х - 4)]и(\, у)ёа. ^ о

В результате точное решение уравнения (6) выписывается в квадратурах:

ад Н

уП(ху)= I IО(д,\x,у)чп-\(д\)л&\.

Построение функции Грина для пары функция тока - завихренность

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

системы эллиптических уравнений 2-го порядка для функций С и у:

С(п) -еАС(и) = g(«-1), Ау(и) =С(И), V(n\x,0) = 0,

Я,,/«) Л|//(и)

y(n)(x,b) = ß , — (x,b) = д^_(x,b) = 0 , (11)

су

ду

(£,у)(о,\) = (£,у)(М), ^(о,*) = д£?)(Ь\),

дд дд

где по горизонтальной координате можно поставить условие периодичности с периодом Ь, если этот период достаточно велик: Ь /Ь >>1. Заметим, что в соотношениях (11) положено Ь=2Л.

Оба эллиптических уравнения 2-го порядка в (11) для функций £ и у решаются с граничными условиями Дирихле. Для у они прямо следуют из условия прилипания. Для £ применяется классическое условие Тома [9], которое переформулирует условия для нормальной производной функции у в некоторые условия для функции £.

Сначала сконструируем функции Грина для уравнений, которым удовлетворяют функции £ и у в рассматриваемом канале постоянной ширины, с нулевыми граничными условиями на боковых поверхностях, с условием периодичности вдоль переменной х. В целях сокращения записей продемонстрируем этот вывод на более простом случае функции тока: АОу = 8(4 - х) 8\ - у), Оу (д,о) = Оу (д, ь) = о,

дОу дОу

Оу(о,\) = Оу(Ь,\) , —у (о,\) = —у (Ь,\). (12)

од од

Эта функция тока может быть построена разложением в ряд Фурье по синусам вдоль переменной \, что автоматически удовлетворяет нулевым граничным условиям на продольных сторонах канала. Каждый коэффициент Фурье удовлетворяет обыкновенному дифференциальному уравнению 2-го порядка по переменной д. Его периодическое решение имеет

вид Gw = £

sin(bmv)sin(bmy) ebmS-X + ebm(L-S-x|)

жт

1 - ebmL

bm =

жт b

Как только функция Грина для у построена, соответствующее решение уравнения Пуассона (12) может быть записано в явном виде с использованием

Lb

теоремы Грина: у (х,у) = ||£(и)(4,\) .

оо

Ключевым моментом рассуждений является тот факт, что функция Грина допускает разделение переменных. Это позволяет записать соотношения, связывающие функции £ и у , для каждой компоненты Фурье отдельно:

у(т](х)=- п 1 -Ь Ь, )\еЬт (4-х|-ь) + е-Ьи|д-х|£т (№,

2Ьт (1 - е т )о

ад

у(п)(х,у) = 2ути)(х)^т(Ьту),

(13)

т=1

где

да

т

СП0)(х) = 2 Jsin(bmq)C(o)dq -

b

(14)

это опять-таки соответствующая Фурье-компонента для завихренности.

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

V(n)(x, y). Предположим, что вдоль координаты x

взято N узлов сетки, вдоль y - M = 2k . Если величина M - целая степень двойки, то для построения быстрых вычислений можно применить БПФ. Это позволяет вычислить решение в M точках y , равномерно распределенных на отрезке (0, b). Для каждой из M компонент Фурье (m = 1,2,...,M) интегрирование в (14) в силу разделения переменных требует линейного числа операций для каждого Xj, j = 1,2,...,N . Это - O(MN) арифметических операций.

Любопытно сравнить эту оценку, например, с аналогичной оценкой по методу конечных элементов (МКЭ). Последний сводит задачу к линейной алгебраической системе с матрицей размерности MN х MN, содержащей по 5 ненулевых элементов в каждой строке. Итерационные методы решения таких систем, например метод сопряженных градиентов, требуют O(h_1) итераций. Поскольку каждый шаг итерационного процесса - O(MN) арифметических операций, всего при использовании МКЭ требуется порядка O(M 2 N ) операций. В довольно естественном случае, когда N = O(M), число необходимых

операций в МКЭ растет кубически с ростом размерности M increasing, а метод, предложенный здесь

для вычисления v(n) (x, y), является квадратичным по

М.

Похожие результаты легко получаются в применении к функции завихренности. Если мы конструируем функцию Грина Gf , такую, что

Gc - sAGc = 8(£ - x) 8(n - y), Gc (|,0) = Gf (t, b) = 0,

Gc(0,q) = Gc(M), —f (0,7) = —f (L,V) ,

o^ ot,

то легко получить разложением в ряд Фурье по синусам следующее представление

Gc=-Z

sin(bmq)sin(bmy) eimS-x| + eк(L-S-x|)

m=1

к =, К +1.

bsX„

1 - e

1nnL

Если функция Грина построена, то решение первого уравнения в (11) строится с применением интегральной формулы Грина

х) =

2Km (1 - e ~ÄmL )

(15)

х i[eÄm+ eg(i°-1)(4)d4 +

b„

L\p (S,0) - (-1)mC(S,b)]x

bim (1 - e )0 x[eK х|-L) + eх|] dS-[y(o), g(о) ]

¿У l0), gl0) b)sin(bn

y).

m=1

Легко видеть, что количество необходимых арифметических операций в (15) совпадает с аналогичным значением для формулы (13), а именно О(МИ). Единственное отличие состоит в том, что для того, чтобы здесь найти нелинейный член 8<Мп'1)(х) для каждого т и всех х^, требуется вычислить соответствующие коэффициенты Фурье, при использовании подходящей формы БПФ понадобится О(1о§М) операций. Таким образом, можно заключить, что полное число арифметических операций в этом случае имеет порядок О(МИ 1о§М), т.е. чуть более высокий, чем полученная выше квадратичная оценка.

Расчет генерируемых звуковых полей по полю скоростей

Для расчета структуры звукового поля воспользуемся аэроакустической теорией Лайтхилла [9], согласно которой акустическое давление p в среде можно с высокой степенью точности определить из неоднородного волнового уравнения

1 д2 p 2

-V p =

c2 дт2

д 2Т-ду {ду j

Т] * AV1V]

(16)

где c0 - скорость распространения звука в данной среде. Заметим, что тензор Грина Tij выражается через производные от компонент турбулентных скоростей, эффективному вычислению которых были посвящены предыдущие разделы работы. Решение уравнения (16) следует искать при граничных условиях твердой стенки, которые эквивалентны однородным граничным условиям Неймана.

Известно, что решение нестационарного уравнения (16) в безграничной среде может быть построено

в явном виде p( х, t) =

1

52

Т.. (

--14л дх,дх,. J R

у, t--

R

Л

dy.

Если же среда имеет твердые границы, то 4лр( х, t) = JV (у,т) v j (у,т) dy -

-J

1

дх 1дх] у др(у, т)

S | х - у | доу

| х - у |

dSy -

-J

S

р( у,т)

1

| х-у |2 c0 | х-у \

др( у,т) дt

д(| х - у |)

до,.

dS„

т = t -

\ х - у|

1

0

0

+

c

0

о

c

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

Литература

1. Роуч П.Дж. Вычислительная гидродинамика. М., 1980. 616 с.

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

3. Пейре Р., Тейлор Т.Д. Вычислительные методы в задачах механики жидкости. Л., 1986. 352 с.

4. Белоцерковский О.М., Опарин А.М. Численный эксперимент в турбулентности. М., 2000. 223 с.

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

5. Moin P., Kim J. Numerical investigation of turbulent channel flow // J. of Fluid Mechanics. 1982. Vol. 118. P. 341377.

6. Kim J., Moin P., Moser R. Turbulence statistics in fully developed channel flow at low Reynolds number // J. of Fluid Mechanics. 1987. Vol. 177. P. 133-166.

7. Рождественский Б.Л., Симакин И.Н. Моделирование турбулентных течений в плоском канале // Журн. выч. мат. и мат. физики. 1985. Т. 25. С. 96-121.

8. Никитин Н.В. Прямое численное моделирование трехмерных турбулентных течений в трубах кругового сечения // Изв. АН СССР. Механика жидкости и газа. 1994. № 6. С. 14-26.

9. Голдстейн М.Е. Аэроакустика. М., 1981. 295 с.

23 октября 2008 г.