Известия Тульского государственного университета Естественные науки. 2015. Вып. 4. С. 5-13 = Математика
УДК 519.632.4
О спектральной задаче для оператора Орра—Зоммерфельда *
С.Д. Алгазин, Г.Х. Соловьёв
Аннотация. Методом вычислительного эксперимента исследуется задача о распределении собственных значений оператора Орра-Зоммерфельда.
Ключевые слова: численный алгоритм без насыщения, уравнение Орра-Зоммерфельда.
Введение
Запишем полную систему уравнений движения вязкой несжимаемой жидкости, состоящую из уравнения неразрывности и трёх уравнений Навье-Стокса:
аиг.
йт и = 0; р—— = —Угр + дАиг аЛ
(1)
Неизвестными являются три компоненты вектора скорости и = {и, и, и:} и давление р. Обезразмерим уравнения, выбрав в качестве размерно-независимых параметров характерную длину, скорость и плотность жидкости Ь, и и р. Тогда в системе останется единственный безразмерный параметр — число Рейнольдса К = Ьир/д. Раскрывая производные в (1), получаем следующую систему уравнений:
ди ди ди о дх ду дх '
ди ди ди ди др 1 (д2и д2и д2и дЛ дх ду дх дх К \дх2 ду2 дх2
ди ди ди ди др 1 / д2и д2и д2и \
дЛ дх ду дх ду К \ дх2 ду2 дх2)
ди ди ди ди др 1 (д2и д2и д2и
---1-и-----\-и— +--—г- + —г- + —г-
(2)
дЛ + и дх + У ду + и дх дх + К \ дх2 + ду2
дх2
* Работа выполнена при финансовой поддержке РФФИ (проекты № 15-01-01739, № 16-01-00189).
Будем рассматривать малые возмущения некоторого заданного установившегося плоскопараллельного течения
и = и0 (г), и = 0, ш = 0, р = р0,
вектор скорости которого параллелен оси х и зависит только от координаты г, а давление постоянно. Возмущённое движение имеет вид
и = и0 (г) + и' (х, у,г,Ь), и = и' (х, у, г, Ь), ш = ш' (х, у,г,Ь), р = р0 + р (х, у, г, Ь).
Здесь штрихом обозначены возмущённые величины, причём \и'\ , \и'\ , \ш'\ << тах |ио|, \р'\ << ро. Устойчивость будем изучать в линейном приближении. Для этого подставим эти выражения в (2), линеаризуем и вычтем те же уравнения для невозмущённого течения. В результате получим следующую систему линейных уравнений относительно возмущенных величин:
ди + ди + дш о дх ду дг '
ди' ди' дио
-Т7Т + ио — + ш —-
дЬ дх дг
др' 1 дх К
д2и' д2и' д2и'
дх2 ду2
дг2
ди'
ди'
др' 1 дЬ и° дх ду К
д2и' д2и' д2и'
дш' дш' дЬ дх
= др' + 1
= - ~дх + К
дх2
д 2ш' дх2
+
ду2
+
+ д2ш' + ду2
дг2 д2ш'
дг2
(3)
Система (3), дополненная граничными и начальными условиями, определяет поведение малых возмущений течения.
Заметим, во-первых, что хотя физическая постановка задачи предполагает, что возмущения действительны, в силу линейности уравнений можно рассматривать и комплексные решения. Такие решения предполагают, что физический смысл имеет их действительная и мнимая части, а комплексная форма используется исключительно из-за математического удобства.
Во-вторых, полученная линейная система имеет коэффициенты, зависящие только от г и, следовательно, имеет решения, экспоненциальные по х, у и Ь. Это свойство используется в основном методе исследования устойчивости — методе нормальных (или собственных) мод, который заключается в следующем. Рассматриваются не произвольные возмущения с начальными условиями, а возмущения специального вида:
( и' (х, у, г, Ь) и' (х, у, г, Ь) ш' (х, у, г, Ь) \ р (х, у, г, Ь)
)
( и (г) Я(г)
ш (г) \ Р (г)
)
е
(4)
Здесь а и в — заданные вещественные числа. Такие решения, называемые модами, являются бегущими волнами: их вещественные (и мнимые) части имеют вид
Re f' (x, y, z, t) = Re ( f(z) ei(ax+ßy-^ =
= ^ Re/ (z) cos (ax + ßy — Reut) — Imf (z) cos (ax + ßy — Reut) j е1тш* = = f(z) cos (ax + ßy + ф (z) — Reut) е1шш1.
В каждом слое z = const движение имеет вид волны, перемещающейся без деформации в направлении, заданном a и ß Одновременно происходит усиление или затухание волны, зависящее от знака Imu Величины a и ß называются волновыми числами; вектор {a,ß} в плоскости xy называется волновым вектором. Направление волнового вектора задает направление движения волны, а его длина определяет длину волныА:
2п
у а2 + в2 =
л
Величина и называется частотой волны. Она должна находиться из решения системы (3) после подстановки туда (4), как будет показано ниже. Таким образом, и = и (а, в) Если 1ти (а, в) > 0 для каких-нибудь а и в, то течение неустойчиво. Если же 1ти (а, в) ^ 0 для всех а, в £ К, то можно говорить лишь об устойчивости возмущения вида (4). Однако часто из этого следует и устойчивость возмущений произвольного вида.
Итак, поведение возмущений в виде бегущей волны (4) часто определяют устойчивость течения по отношению к произвольным возмущениям, поэтому сначала остановимся на исследовании таких возмущений. Возмущения произвольного вида будут рассмотрены позже. Подставим (4) в (3), введя для удобства вместо и величинус = и/а, называемую фазовой скоростью:
. - дш
гаи + гви + —— = 0, дг
. , ^дио . ^ 1 ( 2^ д2и\
га (и0 — с) и + ш —— = -гар +--—а и — в и + ,
дг К \ дг2
ia (uo — c) и = —ißp + R (—a2и — ß2v + ,
■ , s - dp Ii д2 w\
ia (u0 — c) w = — -—b— —a w — ß w + 1 dz R \
(5)
dz2)'
Граничные условия для этой системы зависят от рассматриваемого основного течения. Если, например, это течение между двумя твёрдыми стенками г = г\,г2, то это условия прилипания:
и = и = и = 0, х = х\,х2.
Тривиальное решение (5) с однородными граничными условиями — нулевое. Очевидно, нас интересуют решения, отличные от нуля; они существуют лишь при некоторых определённых значениях с. Таким образом, получаем задачу на собственные значения с, которые могут быть выражены в виде
^ (а, в, с, К) = 0. (6)
Целью исследования уравнения (6) является определение такого числа Рейнольдса Ксг, что течение устойчиво при К ^ Ксг и неустойчиво при К > Ксг. Такое исследование упрощается благодаря теореме Сквайера, который показал, что достаточно исследовать только плоские возмущения (т.е. возмущения с в = 0), поскольку они являются наиболее неустойчивыми.
Теорема Сквайера. Для вычисления критического числа Рейнольдса Ксг течения вязкой жидкости достаточно рассматривать только плоские возмущения.
Уравнения Орра—Зоммерфельда и Рэлея. Таким образом, дальше будем рассматривать плоские возмущения. Система (5) принимает вид
гаи+ Ц =0,
ia (uo - c) u + w = -iap+ R (-a2u + ia (uo - c) w = -§ + R (-a2w + 0) .
Преобразуем её к одному уравнению относительно w. Выразим u из первого уравнения
- i dw
u = adz (7)
и подставим во второе:
^ i ( ) dw i ^du0 1 / dw 1 d 3w P a Uo C dz aW dz R \ dz a2 dz3
Подставляя p в третье уравнение, получаем 1 (cftw 2 d 2w 4 Л d (, . dw ^du0\ 2, . ^
— 'dz2 + aw) = dz [(uo-c) dz-wu -dz)-a (uo-c) w
Это уравнение называется уравнением Орра-Зоммерфельда. Часто его записывают в операторной форме:
1 /г, 2 2\2 ~ =( ) / гч 2 2\ ~ d 2u0
(.D2 — a2) w = (u0 — c) (D2 — a2) w--w,
iRa dz2
где D = d/dz.
Граничное условие прилипания на твёрдых стенках в силу (7) формулируется так:
— дши
ш = — = 0, г = г\, г2-дг
С этими граничными условиями уравнение Орра-Зоммерфельда определяет задачу на собственные значения с.
В случае невязкой жидкости К — с и уравнение Орра-Зоммерфельда вырождается в уравнение, называемое уравнением Рэлея:
д (, . дш ^ ди0\ 2 , . ^
- и— с) дг — ш — а (ио — с) ш = 0,
или
д2и
(ио — с) (Б2 — а2) ш —др°ши = 0-
Поскольку оно имеет не 4-й, а 2-й порядок, то ставится лишь по одному граничному условию на жёстких стенках — условию непротекания: ши = 0,
г = гь г2-
Итак, мы получили два уравнения, описывающие поведение возмущений: в вязкой жидкости это уравнение Орра-Зоммерфельда, в невязкой — уравнение Рэлея. Каждое из них вместе с граничными условиями определяет задачу на собственные значения с. Если существует а £ К, такое, что 1тс (а) > 0, то течение неустойчиво. В противном случае течение устойчиво, при этом различают два типа устойчивости: если 1тс (а) < 0, то устойчивость асимптотическая (возмущение экспоненциально затухает), если же 1тс(а) = 0, то такая устойчивость называется нейтральной — амплитуда возмущения не растёт и не затухает. Изучению уравнения Орра-Зоммерфельда посвящены работы [1, 2, 5]. Численные методы можно найти в [1, 3, 4].
1. Численный алгоритм без насыщения для уравнения Орра-Зоммерфельда
Уравнение для собственных функций. Пусть ¡а = ^Х? — а2. Тогда
¡1 • ф + гаК (и" — та) Ф — КМаф = 0, Ф (±1) = Ф' (±1) = 0. (8)
Здесь а — вещественный параметр; И — число Рейнольдса; и=1 — у2 — невозмущенный профиль; Л — спектральный параметр (неустойчивому течению соответствуют Л с положительной мнимой частью). Запишем уравнение (8) в виде
ф1У (х) = р (х) ф" (х) + д (х) ф,
где
p (x) = 2a2 + XR + iaRU (x), —q (x) = a4 + iaRU" + ia3RU + R\. Пусть K (x,y) — функция Грина оператора с граничными условия-
ми
ф (±1) = ф' (±1) = 0.
Тогда
1 1
ф (x) = У A (x, у)ф (y) dy + xj B (x, y) ф (y) dy, -1 -1
A (x, y) = — (a4 + ia3RU (y)) K (x, y) + 2 ia RU' (y) Ky (x, y) + + (2a2 + iaRU (y)) Kyy (x, y); B(x, y) = R [Kyy(x, y) — a2K(x, \
Причем
K(xy) = l 24 (1+ x)2(1 — y)2(1 + 2y — 2x — xy), y>x, K (x,y) 1 24 (1 — x)2(1 + y)2(1 + 2x — 2y — xy), y<x,
(1+x)
2
Kyy(x,y)H ^ —^ т
I [—1 — 2y — xy] , y>x,
K (x, y) = i (1+x)§(1 (x — 2y+xy), y>x,
[ (1-x)8(1+y) (x — 2y — xy), y < x. 2. Дискретизация
Имеем
ф (x) ^ ¿] Inj (x) ф, Inj (x) = (x — (xj)' xj = cos2—= ^ 2'"''
1 1 = Е/A fey)«* (y)dyфj + x WB(^n (y) dy j + л,
j -1 j -1 1 1 aij = j A (xi, y) Inj (y) dy, bij = j B (xi, y) Inj (y) dy, 11
— приближенные значения собственных функций в узлах, ц — собственные значения дискретной задачи, х = (Х1,Х2, ■■■,Хп)',
Х = Ах + М^Х, (I - А) Х - В = 0. Получаем алгебраическую задачу на собственные значения
-1 - (I - A)"1 B
X = 0.
3. Уравнение Орра-Зоммерфельда. Формулы для программирования
Входные параметры: а — вещественное, И — число Рейнольдса, т размер матрицы (п=2т),
, « + », чётный случай,
« — », нечётный случай.
Конечномерная задача имеет вид: и = Аи + ХБь, где А и В — матрицы размера т х т. Имеем
2т-1
А- =1
a- j —
m
k=0
*3-m E '( i1 ± (-1)^ COS j Imk X)
0j = -nZ- = -7ZZ- ' j = 1, 2, 3,..., m;
(2j-1)n = (j-1) n 2n 4m
х- = cos 0-, i = 1, 2, 3,..., m; причем «+» для четных собственных функций и «-» для нечетных, штрих у знака суммы означает, что слагаемое при k=0 берётся с коэффициентом 1/2. Формулы для Imk (х-) имеют вид
Imk (Xi) = [2а2Ilk (Xi) - aAhk (х-)] + i [aRhk (x-) + 2aRI3k (x-) - a3RI5k (x-)]
Поэтому следует написать две подпрограммы для вычисления ReA и ImA. Далее
1 2m-1
Bij = m 2 '(1 ± (-1)k) cos k0j ■ Ink (Xi) ,
k=0
Ink (xi) = i [a3Rhk (xi) - aRIik (x-)] .
Таким образом, матрица B — чисто мнимая (сравни с ReA ).
Формулы для вычисления интегралов Ijk,j = 1, 2, 3, 4, 5, имеют вид
Ilk (x) = -1 (1 - x2) [Iko (x) + (x + 2) Iki (x)] +
+ 4 (1 + x)2 [-Jko (x) + (2 - x) Jki (x)] ,
I2k (x) = + 4 (1 - X2) \Jk2 (X) + (X + 2) Ik3 (x)] -
- ^(1+ x)2 \-Jk2 (x) + (2 - x) J k3 (x)] + Ilk,
I3k (x) = - 4 (1 - x2) \xIkl (x) - 2I k2 (x) - (2 + x) I k3 (x)] -- 4 (1 + x)2 [xJki (x) - 2Jk2 (x) + (2 - x) Jk3 (x)] ,
I4k(x) = 24 (1 - x2)x
x [(1 + 2x) I\ (x) + 3xIki (x) - 3Ik2 (x) - (2 + x) Ik3 (x)] + + 24 (1 + x)2 [(1 - 2x) Jko (x) + 3xJki (x) - 3Jk2 (x) + (2 - x) Jk3 (x)] ,
I5k (x) = I4k (x) -
-1 (1 - x2) \(1 + 2x) Ik2 (x) + 3xIk3 (x) - 3Ik4 (x) - (2 + x) Ik5 (x)] -- 24 (1 + x)2 [(1 - 2x) Jk2 (x) + 3xJk3 (x) - 3Jk4 (x) + (2 - x) Jk5 (x) .
Список литературы
1. Скороходов С.Л. Численный анализ спектра задачи Орра-Зоммерфельда // Ж. вычисл. матем. и матем. физ. 2007. Т. 47. № 10. С. 1672-1691.
2. Шкаликов А.А. Спектральные портреты оператора Орра-Зоммерфельда при больших числах Рейнольдса // Современная математика. Фундаментальные направления. 2003. Т. 3. С. 89-112.
3. Алгазин С.Д. Численные алгоритмы классической математической физики. М.: Диалог-МИФИ, 2010. 240 с.
4. Бабенко К.И. Основы численного анализа. М.: Наука, 1986. 744 с.
5. Курочкин С.В. Метод выявления неустойчивости и поиска неустойчивых собственных значений в задаче Орра-Зоммерфельда // Журнал вычисл. матем. и матем. физ. 2001. Т. 41. № 1. С. 86-94.
Алгазин Сергей Дмитриевич ([email protected]), д.ф.-м.н., в.н.с., Институт проблем механики им. А.Ю. Ишлинского РАН, Москва.
Соловьёв Гариф Хусаинович ([email protected]), к.т.н, доцент, Московский государственный машиностроительный университет, Москва.
About the spectral problem for the Orr—Somerfield S.D. Algazin, G.X. Solov'ev
Abstract. The method of computing experiment investigates a problem about a distribution of eigenvalues of an operator of Orr-Somerfield.
Keywords: numerical algorithm without saturation, equation of Orr-Somerfield.
Algazin Sergey ([email protected]), doctor of physical and mathematical sciences, leading scientific researcher, Institute for Problems in Mechanics of RAS, Moscow.
Solov'ev Garif ([email protected]), candidate of technical sciences, associate professor, Moscow State Machine-Building University, Moscow.
Поступила 24-10-2014