Решение задачи Орра-Зоммерфельда псевдоспектральным методом по полиномам Чебышева Текст научной статьи по специальности «Физика»

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

Том 13, № 3, 2008

Решение задачи Орра—Зоммерфельда псевдоспектральным

методом по полиномам Чебышева

Т. Г. ДАРМАЕВ

Бурятский государственный университет, Улан-Удэ, Россия

e-mail: dtg@bsu.ru

In this article a numerical method, relying on the Chebyshev polynomials is used to calculate the spectrum of eigenvalues and define curves of constant increment of the Orr—Sommerfeld problem with high accuracy.

В линейной теории гидродинамической устойчивости исследование развития малых возмущений в плоскопараллельных течениях вязкой жидкости сводится к известной задаче Орра—Зоммерфельда [1]:

D4v — 2a2D2v + a4v = iaRe[(U — c)(D2v — a2v) — vD2U ], (1)

где D = d/dy, a — волновое число, Re — число Рейнольдса, c = cr + ici, v = vr + ivi —

искомые комплексные собственное число и функция, i — мнимая единица.

Для течения Пуазейля в трубе профиль скорости основного потока

U(У) = 1 — у2, — 1 < У < 1,

а v(y) удовлетворяет граничным условиям:

v(—1) = Dv(—1) = v(1) = Dv(1) = 0. (2)

Для течения Блазиуса над плоской полубесконечной пластиной U(у) = f(у), где f (у) находится численно из уравнения [2]:

f + 1 ff" = 0, f' ^ 1 f" ^ 0, при У ^ (3)

со следующими граничными условиями:

v(0) = Dv(0) = v(to) = Dv(<x>) = 0. (4)

Для течения Блазиуса задача на бесконечном интервале (0, то) сводится к анализу задачи на конечном интервале (0,L) при больших L, а затем к линейной замене

у = 1 —n/L — на интервале (0,1). Симметрично отображая профиль на интервал (0, — 1),

получаем аналогичную течению Пуазейля задачу.

Исследованию задачи Орра—Зоммерфельда посвящено множество работ [2]. Главными целями численного решения при заданном основном течении являются: поиск

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

кривой нейтральной устойчивости (ci = 0), кривой постоянной скорости роста

(aci = const), вычисление собственных значений и функций для данной пары положительных значений a и Re.

Основные методы решения принадлежат двум классам: 1) конечно-разностные методы [2-5]; 2) спектральные методы, использующие разложение по полной системе ортогональных функций [6, 7].

Конечно-разностные методы являются эффективными при нахождении нейтральной кривой, наименее устойчивого собственного значения и собственных функций. Применение методов спектрального разложения позволяет отыскать спектр собственных значений, но при этом собственная функция вычисляется хуже. В работе [6] с помощью разложения по полиномам Чебышева исследованы первые собственные значения для течения Пуазейля, а в [7] с помощью метода Галеркина по полиномам Лежандра исследуется спектр течения Блазиуса.

В данной работе применяется метод коллокации или псевдоспектральный метод по полиномам Чебышева для исследования устойчивости плоскопараллельных течений, в частности течений Пуазейля и Блазиуса.

Решение задачи (1), (2) ищем в виде следующего разложения:

N

v(y) = znTn{y), — 1 < y < 1, (5)

n=0

где Tn(y) — полиномы Чебышева n-го порядка 1-го рода:

Tn(cos ив) = cos ив, n > 0.

Заданный профиль скорости основного течения также разлагаем в ряд:

N

U(y) = unTn(y), -1 < y <1 (6)

n=0

и находим коэффициенты ип.

Подставляя (5), (6) в (1), используя свойства полиномов Чебышева и приравнивая коэффициенты при различных Тп, получаем уравнения для искомых коэффициентов гп\

1N

[p3(p2 — 4)2 — 3n2p5 + 3n4p3 — pn2(n2 — 4)2]zp — (2a2 — ia Re c) x

24

p=n+4 p=n( mod 2)

N

x P(P2 — n2)zp + (a4 — ia Re c)dnzn—

p=n+2 p=n( mod 2)

1 N

— ^ia Re Y zp Y P(P2 — n2)Un-m — a2 ^ zpu,n—p—

p=2 m=p( mod 2) |p|<N

|m|<p—2 |n—p|<N

|n—m|<N

Y Zp Y m(m2 — (n — p)2)un = 0, (7)

|p|<N m=|n—p|+2

|n—p|<N m+n=p( mod 2)

гр ир ^|р|и|рЬ N — р — 4 2 ^р

р = n(mod2) обозначает, что р — п делится нацело на 2. Используя следующие свойства полиномов Чебышева:

^р = 1, р > 0,

^-1п2,

Тга(±1) = (±1)П, Тга(±1) = (±1)П

можно привести граничные условия (2) при у = —1 и у = 1 соответственно к виду

N

Е

п=С п=С( тоё 2)

N

0,

N

Е

п=С п=С( тоё 2)

N

п2гп = 0,

Е

пгП

Е

п2г„

0.

(8)

(9)

П=1

п=1( тоё 2)

п=1 п=1( тоё 2)

Для течений с симметричным профилем размер системы уравнений (7) сокращается вдвое, так как можно рассматривать только четные гП, п = 0, 2,..., N = 2М. При этом автоматически выполняются граничные условия (9).

Для нахождения М +1 неизвестных коэффициентов г2т, т = 0,1,...,М = N/2 из М + 3 уравнений применяется тау-метод [8], который заключается в отбрасывании двух последних уравнений при т = М — 1,М так, что с двумя граничными условиями (8) получаются М + 1 уравнение для М + 1 неизвестных коэффициентов г2т, т = 0,1,..., М — 1. В результате из (7) получаем алгебраическую задачу на собственные значения общего вида:

Аг = сВг, (10)

где г — вектор неизвестных длины М +1;

А, В — известные матрицы размерности (М +1) х (М +1) следующего вида:

1 1 1 ... 1 \ 0 0 0 ... 0

0 4 9 Ю 0 0 0 ... 0

А= 0 0 а11 ... а1,М-1 , В = 0 0 &11 •••Ь1,М-1

0 0 ЙМ-1,1 ••ЙМ-1,М-1 / 0 0 ЬМ-1,1 ••ЬМ-1,М-1 /

Первые две строки матриц А и В получены из граничных условий (8), а элементы Ягл и — из уравнений (7). Так как матрица В из-за нулевых первых двух строк

особенная, то возникают трудности с ее обращением. Чтобы избежать этого, обнуляя внедиагональные элементы и нормируя диагональные элементы первых двух строк матрицы А и применяя это преобразование одновременно к матрице В, получаем:

1 0 0 ... 0 \ 0 0 0 ... 0

0 1 0 ... 0 0 0 0 ... 0

А= 0 0 Й11 •••Й1,М-1 , В = 0 0 &11 -1 1,М 1 1-0

0 0 ЙМ-1,1 ••ЙМ-1,М - 1 0 0 ЬМ-1,1 о- - 1 ,М 1

Собственные значения системы

(А — сВ )- = 0,

(11)

П

0

л и ~ иМ^—1 Э и I ||M—1 где A =|| aik ||i и B = || bik ||i — Уже невырожденная матрица размерности

(M — 1) х (M — 1), те же, что у исходной системы (7). Применяя далее QZ-алгоритм [9] к полученной системе (11), получаем собственные значения. Для нахождения собственных функций из вычисленных z обратным преобразованием находим z, затем используем разложение (5).

На основе алгоритма была разработана универсальная программа на языке FORTRAN Power Station 4.0 для расчета устойчивости плоскопараллельных течений с двойной точностью.

В данной работе расчеты проведены для течений Пуазейля и Блазиуса. В табл. 1 даны первые 10 собственных чисел течения Пуазейля, вычисленных при a =1, Re = 10 000, M = 200.

Первое собственное число с наибольшей мнимой частью определяет устойчивость течения. В табл. 2 приведены его значения, вычисленные для течения Блазиуса при a = 0.16869, Re = 500 и разных M и L.

Расчеты позволяют сделать следующие выводы. Выбором L при фиксированном M или N, или увеличением M при фиксированном L можно достичь повышения точности. Спектральный метод целесообразно применять для вычисления одновременно несколь-

Т аблица1

(2.374486880407458E—001, 3.791166309875633E-003)

(1.900557202258878E—001,—1.828207620297516E—001) (3.491230585037001E—001,—1.245165788603773E—00l) (3.685068183504831E—001,—2.388279976521671E—00l) (4.748377439650907E—001,—2.086553689039426E—001) (5.128807977577818E—001,—2.868928204096470E—001) (5.843306096842552E—001,—2.668386487364604E—001) (6.208853976536459E—001,—3.359016412279827E—001) (6.623888858297891E—001,—2.797281175566295E—001) (7.194333618198993E—001,—2.617623955442141E—001)

Таблица2

M L C

50 20 (3.669001421988583E—001, 9.008614672213490E—003)

50 25 (3.674198471510465E—001, 8.884444385117093E—003)

50 30 (3.675159113309415E—001, 8.861197869179240E—003)

50 30 (3.675372666980753E—001, 8.856108912264829E—003)

50 40 (3.675370055823886E—001, 8.856033114205948E—003)

50 50 (3.675378470516833E—001, 8.854982808473398E—003)

50 40 (3.675401124612824E—001, 8.859160141251153E—003)

50 60 (3.675431967080640E—001, 8.856010388764098E—003)

100 30 (3.675159317493594E—001, 8.861206560442054E—003)

100 40 (3.675369826998388E—001, 8.856105188630574E—003)

100 50 (3.675377039644471E—001, 8.855951933329087E—003)

100 60 (3.675377492667720E—001, 8.855973804581585E—003)

200 40 (3.675346324654210E—001, 8.855022110732356E—003)

200 50 (3.675354342513038E—001, 8.858435217178896E—003)

Рис. 1. Спектр собственных значений течения Рис. 2. Спектр собственных значений течения Пуазейля: а = 1, Яе = 10 000, М = 100 Пуазейля: а = 1, Яе = 10 000, М = 200

Рис. 3. Спектр собственных значений течения Рис. 4. Собственная функция течения Пуазей-Блазиуса: а = 0.179, Яе = 580, М = 200 ля (А-мода)

Рис. 5. Собственная функция течения Пуазей- Рис. 6. Собственная функция течения Пуазейля (Р-мода) ля (Я-мода)

AL

0.7 J---------1--------'--------1--------'---------1--------1--------1---------'--------т—

0.6 0.7 0.8 0.9 1.0

Re/10000

Рис. 7. Кривые постоянных скоростей роста

ких собственных значений и функций. Вычисленные спектры при разных M = 100 и 200 даны для течения Пуазейля (рис. 1 и 2) и Блазиуса (рис. 3). Видно, что при M = 200 достигается приемлемая точность. Различаются пристеночные моды [10], которые иногда называют модами Эйри, или А-модами, для которых cr ^ 1, с = O[(aRe)-1/3], центральные моды, иногда называемые модами Пекериса, или P-модами, для которых cr ^ 0, c = O[(aRe)-1/2] при aRe ^ то, и имеют место моды Шенстеда, или S-моды, для которых cr ^ 2/3, c ~ —i[n(n + 2)]2a/4Re при aRe/n ^ 0. На рис. 4-6 приведены собственные функции (сплошные линии — реальные части, пунктир — мнимые части) для мод A, P, S соответственно. На рис. 7 изображены полученные спектральным методом кривые по значениям с.

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

[1] Линь Цзяо-цзяо. Теория гидродинамической устойчивости. М.: ИЛ, 1958.

[2] Гольдштик М.А., Штерн В.Н. Гидродинамическая устойчивость и турбулентность. Новосибирск: Наука, 1977. 421 с.

[3] Thomas L.N. The stability of plane Poisseuille flow // Phys. Rev. 1953. Vol. 90. P. 780-783.

[4] Mack L.M. A numerical study of the temporal eigenvalue spectrum of the Blasius boundary layer // J. Fluid Mech. 1976. Vol. 73, N 3. P. 497-520.

[5] Дармаев Т.Г. Анализ эффективности различных методов численного решения задачи Орра—Зоммерфельда // Математика, ее приложения и математическое образование: матер. междунар. конф. Ч. 1. Улан-Удэ, 2002. С. 156-160.

[6] Orszag S.A. Accurate solution of the Orr-Sommerfeld stability equation // J. Fluid Mech. 1971. Vol. 50, N 4. P. 689-703.

[7] Melenk J.M., Kirchner N.P., Schwab C. Spectral Galerkin discredization for hydrodynamic stability problems // Computing. 2000. Vol. 65. P. 97-118.

[8] Lanczos C. Applied Analysis. Prentice-Hall, 1952.

[9] Moler C.B., Stewart G.W. An algorithm for generalized matrix eigenproblems // SIAM J. Numer. Anal. 1973. Vol. 10. P. 241-256.

[10] Дразин Ф. Введение в теорию гидродинамической устойчивости: пер. с англ. М.: Физ-матлит, 2005. 288 с.

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