Конечномерная инвариантная аппроксимация и периодические режимы течения Блазиуса Текст научной статьи по специальности «Математика»

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

Том 13, № 4, 2008

Конечномерная инвариантная аппроксимация и периодические режимы течения Блазиуса

Т. Г. ДАРМАЕВ

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

e-mail: [email protected]

In this article the method of the finite-dimensional invariant projection of the Navier—Stokes equations is applied for the parallel Blasius flow of viscous incompressible fluid on a semi-infinite flat plate. A numerical study of the periodical regimes was presented.

Введение

Известно, что в области ламинарно-турбулентного перехода в пограничном слое быстро нарастают возмущения. Первая попытка оценить влияние нелинейности на судьбу возмущений была предпринята Л.Д. Ландау [1]. С помощью качественных рассуждений он показал, что нелинейность может как стабилизировать нарастающие возмущения, создавая новый устойчивый режим течения, так и вызвать рост возмущений, устойчивых в линейном приближении. Дж. Стюарт и Дж. Ватсон [2, 3] количественно нашли уравнение Ландау для слабонеустойчивых возмущений плоскопараллельных течений в виде асимптотических рядов. Численные расчеты для течения Пуазейля выполнили В. Рейнольдс и М. Поттер [4]. Существенный вклад в развитие нелинейной теории был сделан В.В. Струминским [5]. Он получил решение уравнений Навье—Стокса для слабонелинейных возмущений в виде сходящихся рядов. Просуммировав бесконечные ряды для амплитуд возмущений, он показал, что нарастающие возмущения стабилизируются, а затухающие в линейном приближении затухают и при нелинейном рассмотрении. Решения Струминского сходились на некотором расстоянии от линейной нейтральной кривой. Но вышеуказанные слабонелинейные теории гидродинамической устойчивости [1-5] и прямое численное интегрирование уравнений Навье—Стокса [6] хорошо описывают начальную стадию ламинарно-турбулентного перехода. Слабонелинейная теория не применима для описания последующих стадий, поскольку нужно рассматривать члены высших порядков по амплитуде. Прямое численное интегрирование уравнений Навье—Стокса неэкономично и не гарантирует правильного описания асимптотического поведения решения. Б.Ю. Скобелев [7] разработал метод инвариантной конечномерной проекции уравнений Навье—Стокса, который позволяет получать нелинейные решения в окрестности нейтральной кривой в виде сходящихся рядов. Существенным достоинством метода инвариантной проекции является то, что он гарантирует правильное описание асимптотического поведения решений (т.е. при t ^ то), а также то, что начально-краевая задача для возмущений ламинарного течения сводится к конечномер-

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

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

В данной работе метод инвариантной конечномерной проекции уравнений Навье— Стокса применяется для плоскопараллельного течения Блазиуса вязкой несжимаемой жидкости над плоской полубесконечной пластиной.

1. Эволюционные уравнения для трехмерных

моногармонических периодических возмущений

Выберем систему координат так, чтобы ось х была направлена вдоль пластины, г-координата — поперек пластины, у-координата — перпендикулярно пластине, начало координат совпадает с передней кромкой пластины.

Рассмотрим уравнение Навье—Стокса для завихренности:

дш

— = (ш ■ У)и - (и -У)ш + vДш, (1)

и уравнение неразрывности:

(-и) - IX1 + Ж + £ = "• (2)

где V — коэффициент кинематической вязкости, х = (х,у,г) — координаты, и = («1,^2,^3) — скорость, ш = (^1,^2,^3) — завихренность, ш = Ух и или

дщ ди2 ди1 ди3 ди2 дщ д д д

= ^---, ^2 = ^---, ^3 = ^---X— , ш -V = Ш1— + + .

ду дг дг дх дх ду дх ду дг

Предположим, что уравнения (1) и (2) имеют стационарное решение и(0), ш(0). Положим и = и(0)(у)+и, ш = ш(0)(у)+шй. Решение и удовлетворяет естественным граничным условиям на пластине. В дальнейшем будут рассматриваться уравнения для возмущений й, ш, и поэтому знак тильды будем опускать. Подставляя и, и ш в (1), получаем уравнение для возмущений:

^ = (ш(0) ■ У)и + (ш ■ У)и(0) - (и(0) ■ У)ш-

-(и ■ У)ш(0) + (ш ■ У)и - (и ■ У)ш + VДш. (3)

Для плоскопараллельных течений скорость и10) = и (у) — профиль стационарно-

(0) (0) п (0) йи (0) (0) „ л,

го решения, и2 = и3 = и, а завихренность и3 =--—, и1 = и2 = и. Учитывая

йу

соотношения для ш, а также уравнение (2), и преобразовывая (3), получим следующую систему уравнений:

д и д д\ д сРи ди2

дг дх ) 2 Су2 дх д д

— [(ш ■ У)иэ - (и ■ - тт" [(ш ■ У)щ - (и ■ У)^]

дх д

дх

д

дг дх

д^ о^

дх2 дх2 1

д2 д

—- + и---VД ш2 = (ш ■ У)и2 — (и ■ У)^2 —

д2и2 дш2 2 + 2

Си ди2 Су дх

дх2 + дх2

из = -

дхду д 2и2

дх дШ2

дхду дх

(4)

Положим давление р = р(0) (у) + р, тогда из уравнения для скорости, опуская знак тильды для р, имеем

д д\ ттди1 Си др

— - VД щ + и—--и^— = ----(и ■ У)щ,

дЬ ) дх Су дх

д 4 \ ттди2 _ - vд ) щ2 + и_

д Л , ттди3

дг -д)из + и1х

др

-Т,--(и ■ У)и2,

ду

др

- дх - (и ' ^

(5)

Рассмотрим трехмерные моногармонические возмущения, периодические по продольной координате х с волновым числом а и по поперечной координате х с волновым числом в: и)т(у, вкт) = ~й,кт(у)егвкт, 3 = 1, 2, 3; -то < к,т < то, вкт = ках + твх, г — мнимая единица.

Положим:

и = Е икт(у,вкт), 3 = 1, 2, 3;

к,т „кт

р = £ ркт(у,9кт)

,т

и и0° = 0 с учетом граничных условий.

Подставляя и^ в уравнение неразрывности (2) и приравнивая члены при е кт, полу-

чаем:

ди кт ди кт ди кт

кад±- + тв(1 - 6л°)+ (1 - 6л°)ди2

дв,

т

дв

т

ду

0,

где 6 — символ Кронекера, Л = |к| + |ш|.

Таким образом, уравнение неразрывности удовлетворяется введением функций тока:

дф

т

ду

каи кт + три кт + 6Ло(и?° + си°°),

т

/1 X \дф „ кт

(1 - 6Л0)^- = -и2 ,

дв

т

где с — произвольная постоянная.

Из последних двух уравнений системы (4) и (6) находим:

(1 - ¿А0

„ кт

^кт

дфкт ка—т:--+ 1шршкт

дУ

+ $А0Щ

00

< и

кт

-¿(1 - 6ао)Ф

кт

(7)

и

кт

(1 - 5.

А0

^кт

дфкт

тв ----¿каш1кт

дУ

+ 5аоиз

00

где 11т = (ка)2 + (тв)2, ш = (1 - 5ао) Е ш2кт(У, вкт).

к,т

Складывая в (5) уравнения для щт, умноженное на ка, для щт — на тв (к, т = 0),

для и"0 — на 5а0 и для щ0 — на с5а0 , и учитывая (7), получим

дд — - VД + каи —— дЬ двкт

дфкт , ди дфкт

- ка——

-с5ао [(ик1т1 -у^2]00 - 7к2

др

кт

т дв,

кт

дУ дУ двкт

[ка(ик1т1 ■ У)ик2тз + тв(ик1т1 ■ У)ик2т2] кт - 5ао [(икт ■ У)ик2т2] 00 . (8)

- \ка

Из второго уравнения системы (5) с учетом (7) получаем

дфкт дркт

— - иД + каП д

дЬ двкт/ дв к т

дУ

+ [(ик1т1 ■ У)и22т2]

т

Используя соотношения (7), продифференцируем первое уравнение системы (5) по у,

второе — по в к т, умножим на и сложим:

д ^* д \ - ^кт дф

— - VД + каи , дЬ двкт

Афкт - ка

кт д I,

^т^ [(и к 1т1 -У)и22т2] ^ -

-5

А0

д

д (и к 1т1 ■ У)и к2т2 дУ

¿У2 двкт двкт

00

- с5ао

д

(ик1т1 ■ У)ик

ду 3

00

где Д

д2

+ 1кт

д

- — [ка(ик1т1 ■ У)ик2т2 + тв (ик1т1 ■ У)ик2т2] д2

т

(9)

дУ2 ' 'ктдв2т'

Обозначим:

%к 1 к2 = (к\к2а2 + тгт2в2)

дфк 1

т1

дУ

+ гав(тгк2 - к!т2)ш.

1 т1

Якгкъ = (к\к2а2 + тг т2в2) ^

1 т1

+ гав(тгк2 - к!"^)^1т1,

¿У

в2

¿ф'

т

Д = ^- , Ект = каиО0 + тви30,Укте = (ка + ств+ ¿(тв - ска)^,

где ш2 = Е т(У,Ь)в^0^т, фкт = фкт(У^увкт, пао = 1 - 5ао, причем Ш°0 = 0 и Хкф1 =

к,т

4 1к2 .

2

Выражая икт через фкт и и^, с помощью (6) и (7) получаем из (9):

IА + ^А - ^) - ^А

ф

кт = ЬктП\10ПЛ20

А

к1к2 Лк2т2

фк2т2 - Ф

кхтх

дфк2

т2

1т1

—г-

• пл1о пл2о

^22т2

1

^Ук1т1

%к2к ду^к1к2 + Ак1к2 ~ду^к2к

дф^1 д ^к — Ффк1т1

ду ду

ду д 2 ^к2к

ду2

кт

кт

-г

. / д2ф

кт

ду2

-Ккт — Ф

кт д Нкт

ду2 , гг1кктг}ктКкт - г5Хо ^^ ду / ¡к

X

2 т2

X

т2

Ук

к2т2С

д^к1к2 . У д^к2т2С \ 0 Ук1т1 д ^к2т2С дф 1 1 д^к2т2С

+ Ак1к2-^- -2ф -

1 т1

ду

ду

ду2

ду ду

00

Используя (7), из второго уравнения системы (4) находим и

кт

д л , г г д

— - VА + каи——

дг дбкт

кт и2 =

дик2т2

и

к1т1 д и2 2

+ (к2аик1т1 + Ш2вик31т1)

ду

к1т1 \ д и2

ди22т2

дв

2 т2

(10)

-(и к 1т1 • Ч)ик2т2 - тр

2 т2

(и дикт

1у дв,

т

т

Обозначая

Ук1к2 =

1 - 6

Л10

1 т1

оп , \ д2Фк1т1 -Г, , 2 а2\ дик1т1

ар(к2тг - кгт2)———--г(кгк2а + т1т2р )

ду2

ду

(11)

и учитывая, что

получаем далее

дик2т2

и

к1т1и и'2 2

ду

-г(1 - 6л20)и

1 т1

дф'

2 т2

ду

д_ дг

~~ — VА + гкаи ) икт

+ Фк2т2 ук2 - гпл10

гпл2

1

дфу к2т2

к1т1 дф

%

ду

2 Ук1к2и22т2 - ф

+ ар (кг т2 - к2тг)ф к1т1 фк2т2 +

т

1 т1 д и2

1 т1

дик 2т2\

дУ )

+

д (и ~

+фктд (каиО0 - три°°) - гикт(каи°° + три3°) - тр^фк ду (у

Учитывая и°° = 0, получаем

и

00

дф

00

ду

- си

00 з •

С помощью (5) находим уравнение для и3

00

д

д2

--V-

дг ду2

и

00

Пл1°Пл2о

^22т2

, к1т1 I од2фк2т2 дик2тЛ

ф I т2р ду2 - г ду )

(12)

(13)

7'к 1т1

рдфк2т2 к Л к2т2

к1к2 \ т2р—ду--гк2аи2 2

00

1

1

1

Система уравнений (10), (12)-(14) полностью определяет величины —кт, о)и°° и «3°, при этом ставятся следующие граничные условия:

Л " " Л "

^т(0) = = = Туфкты =0; (15)

= = 0, м?°(0) = и°°(то) = 0.

Таким образом, нужно решить систему трех уравнений для —кт, ))^т и и°°, а и°° находится по формуле (13).

Систему уравнений (10), (12), (14) можно записать в виде системы эволюционных уравнений для векторов vkm = (Д—кт, ))^т, ^0«°°) в виде

-Ьктукт + Nкт(^7), -то < I < то,

ЛГ

или в виде одного уравнения для вектора V = |укт(у)егбкт } , —то < к, т < то,

= —¿V V + N (V), (16)

где Ь и N — замкнутые неограниченные операторы в гильбертовом пространстве Н.

2. Двумерная инвариантная проекция для течения Блазиуса

Согласно теории, изложенной в [7], спектр р оператора (—Ь) допускает разбиение: р(^) = р^) и р2(^), р1^) П р2(^) = 0, где р1^) — ограниченная часть спектра, и существует проекционный оператор Р, такой, что пространство Н представимо в виде прямой суммы ортогональных подпространств:

Н = Р,н © (I — Р,)Н.

Оператор (—¿1) = — Р, Ь действует в подпространстве Р, Н, его спектр ограничен и равен р1(^). Спектр неограниченного оператора (—¿2) = —(I — Р/равен р2(V). Ограниченная часть спектра р1^) состоит из п пар простых изолированных собственных значений (Ак, Ак), где Ак — комплексно-сопряженное к Ак. Обозначим Р(к) проектор на собственное пространство оператора (—Ь), отвечающее паре собственных значений

(Ак, Ак):

Р(к) V = (V,—к )н ^к + (V,—к )н <Ак = Ук, (^кк) = 1,

где ^к — собственная функция оператора (—Ь), —к — собственная функция сопряженного оператора (—Ь*), (., .)я — скалярное произведение в комплексификации вещественного пространства Н.

Система уравнений (16) принимает вид

^ = —¿(к)ук + N(кЧ к =1, 2,..., п, ЛГ (17)

— = —+ N2 V,

(Г

где V = у1 + у2 + ... + Уп + г, г = (I — Р,)v; Ь(к) = Р(к)Ь, N(к) = Р(k)N, N2 = (I — Р,)N.

Введем в подпространствах Р(к) Н полярные системы координат и будем искать решение (17) в виде

Ук = 2Ке(т'ке^р) = 2Ие(тке^р) + Ук(гь гп, Оь 9П), г = X(ть ...,Гп,вг, ...,9п).

(18)

Динамическая система

От = (Пк + Ь(к))гк, ^ = (ак + с(к))гк, к = 1, 2,..., п, аъ аъ

определяет поведение траекторий уравнения (18) на 2п-мерном инвариантном многообразии, где Лк = Пк+г.к, а Ь(к), с(к) — некоторые функции, зависящие только от координат

г к, Ок.

Предельное многообразие определяется функциями д*(г, О), Ь(к)(т, О), с(к)(т, О) в виде рядов по степеням Тк'

те те ,,. те

д* = £ д^, ь(к) = Е Ьк , с(к) = Е сРт3, |в|=2 |в|=1 |в| = 1

п

где дв = Е ^ + Х9, т9 = т1 ■ т22 ■... ■ тП , | = 81 + ... + зп. Функции дв определяются к=1

следующей рекуррентной системой линейных уравнений'

.112 + ^ = — 2Яе £ [(6™ + гС£Рк] - £ с£1) & + ,

1 к=1 к+р=в 1 (19) те • '

N = £ N9Т9, 5к =(81 ,...,8к - 1,...,8п).

|9|=2

В соответствии с теорией, изложенной в [7], рассмотрим собственные векторы оператора (—Ьи):

Л + Ь„ V = 0.

В покомпонентной записи это уравнение принимает следующий вид:

{ а2 \ . (РТТ * {а2 \2 ,

(Л + гкаи) ^— — ^ Фкт — гка/З^^ — » (¿уа — т!^ ^ = 0,

а2 аТ

(Л + гкаи)ик2т — » ( — — си1т = —т^—'фкт, (20)

аУ

Ли00 — »¿^ = 0.

Легко видеть, что первое уравнение — это уравнение Орра—Зоммерфельда для трехмерных возмущений (Л = —гкас и без ограничения общности можно положить к =1, т = ±1).

Для двумерной инвариантной проекции рассмотрим первое собственное число задачи Орра—Зоммерфельда А = п + Тогда оператор (—Lv) имеет однопараметрические семейства собственных функций вида

ip(x, y, z) = ^Аф (x, y, z) , u2 (x, y, z), 0 j ,

ф (x, y, z) = f (y) [ei(ax+ez) + , (21)

(x,y,z) = fi (y) [ei(ax+ez) - Tei(ax-ez)] ,

где f (y) — собственная функция задачи Орра—Зоммерфельда для течения Блазиуса, а fl(y) — решение второго уравнения в (20) при k = m — 1.

Значение параметра 0 < т < 1 определяется постановкой задачи. При т =1 мы имеем дело с возмущением в виде стоячей волны в направлении z-координаты:

ф (x, y, z) = f (y)eiax cos fiz.

При т = 0 имеем бегущую трехмерную волну, при промежуточных значениях т — смесь двух предыдущих возмущений.

Сопряженный собственный вектор у*km = (ф*кт, u2^km, 8\ou3°0) удовлетворяет системе

Г / й2 \ й ( \2 йи

(А - гкаЩ - 72кпг) ф*кт - 2гкаи~ф*кт - V — - у2кт) ф*кт = -ш^—u*2km,

(А - гкаи) и*кт - V - т2т) <т = 0,

Аи300 - v^--u*300 = 0. х 3 йу2 3

Очевидно, что собственному вектору системы (21) соответствует сопряженный собственный вектор вида

р* (х, у, г) = (ф * (х, у, г), 0, 0).

Здесь ф * (х,у,г) = / * (у) (ег(аж+в2:) + т1ег(ах-вг^ , где / * (у) — решение следующего (сопряженного с уравнением Орра—Зоммерфельда) уравнения (т1 можно положить равным нулю):

(А • тт\( й2 А ** О- йий/* ( й2 а2\2 0 а 2 2 + в2 (А - гаи){йу-2 - й) / - - ^ - й) / =0> й =а +

В случае двумерной инвариантной проекции предельное многообразие определяется функциями

тете те

9 * = Е 9я г3, Ь(г) = £ Ь2пГ2п, с(г) = Е С2п,г2п,

|Я|=2 п=1 п=1

93 = Е 98кегкГ, 0 * = в - ах. к=-в

Из (19) получаем следующую рекуррентную систему уравнений для дак: гк [(а - аи) Дк + аБ2и] д8к - ^ Дкдзк = = -6к1(6в_1 + гсв-1)Д1/ - 6к,-1(Ьв-1 - гсв-1)Д1 ¡--гк ^ сдДкдРк - га [1дФВД)дрз - 3Вддгд22др^ , (22)

д+р=в д+р=в 1+]=к

С

где И,е — число Рейнольдса, О = —, Дк = О2 - (ка)2.

Су

Функции д3к удовлетворяют следующим граничным условиям: д3к (0) = Од^к (0) = д3к (то) = Бд3к (то) = 0.

Будем искать решения (18) в классе 2п-периодических функций от в к, удовлетворяющих условиям ортогональности вида

2п

У е_гГ (д„ р*)в Св* = 0, (23)

0

где — собственный вектор сопряженного оператора (-Ь*), а скалярное произведение (■, ■)в в данном случае имеет вид

2п/а 2п/в те

= / Сх / Сх / (хГх-2- + х».х22> + х»^) Су,

0 0 0

х1 = ф, х2 = Ш2, х3 = и30. Так как gmk — вектор вида

gs = (дфтк,шщк,бл0и30)) (в > 2),

то условие ортогональности (23) принимает следующий вид:

М ^ - а2) ф-с*=°-

Нормируем М * и М условием

те

М СС^ - ^ = '■

0

Тогда коэффициенты Ь5, с5 определяются формулой

Ь3_1 + гс5_1 = М *

("Л, - С-2 - а2) ф.

« А " ^ I Су2 " ) ф1

д+Р=в

Су,

где через ("811)1 обозначена правая часть уравнения (10) для ф]1.

те

3. Численные расчеты

Амплитуда периодических режимов определяется из уравнения

N

+ 5] ь2пА2п

0, N ^ то,

П=1

где ^ = аа — линейный коэффициент нарастания.

Рекуррентная система (22) для нахождения коэффициентов при N = 2, 3,4, 5 решалась методом ортогональной прогонки [8] на неравномерной разностной сетке, сгущающейся в пограничном слое.

В данной работе проводилось численное исследование поперечных возмущений (в = 0). В результате вычислений выявилась следующая картина (рис. 1). При некотором значении волнового числа а = а* от нижней ветви линейной нейтральной кривой течения Блазиуса (рис. 1 — сплошная линия) ответвляются устойчивый и неустойчивый режимы. Неустойчивый режим соответствует верхней части, а устойчивый — нижней части амплитудной поверхности. На рис. 2 представлен срез амплитудной поверхности при а = 0.206906. При некоторых числах Рейнольдса происходит слияние этих режимов в точках тангенциальной бифуркации [9], соответствующих точкам складки в теории катастроф [10]. С некоторого а передняя складка амплитудной поверхности из нефизической области отрицательных значений квадратов амплитуд выходит в область положительных значений. Эти точки передних складок, полученные численными расчетами при в = 0, отображены на рис. 1 пунктиром. С увеличением а амплитудная поверхность периодических решений отрывается от линейной нейтральной кривой, и при некотором значении а периодические решения исчезают.

Рис. 1. Линейная и нелинейная нейтральные кривые

Рис. 2. Срез амплитудной поверхности при а = 0.206906

В работе [11] были проведены эксперименты при малых значениях волнового числа а и выяснено, что линейная теория достаточно хорошо описывает развитие таких возмущений. Результаты данной работы показывают, что требуются тщательные исследования в области значений а > а* для обнаружения новых нелинейных эффектов в экспериментах.

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

[1] Ландау Л.Д. К проблеме турбулентности // Докл. АН СССР. 1944. Т. 44. С. 339-342.

[2] Stuart J.T. On the non-linear mechanics of wave disturbances in stable and unstable parallel flows. Part 1. The basic behavior in plane Pouseuille flow //J. Fluid Mech. 1960. Vol. 9, pt. 3. P. 353-370.

[3] Watson J. On the non-linear mechanics of wave disturbances in stable and unstable parallel flows. Part 2. The development of a solution for plane Pouseuille flow and for plane Couette flow //J. Fluid Mech. 1960. Vol. 9, pt. 3. P. 371-389.

[4] Reynolds W.C., Potter M.C. Finite-amplitude instability of parallel shear flows // J. Fluid Mech. 1967. Vol. 27, pt. 3. P. 465-492.

[5] Струминский В.В. К нелинейной теории развития аэродинамических возмущений // Докл. АН СССР. 1963. Т. 153, № 3. С. 547-550.

[6] Rist U., Fasel H. Direct numerical simulation of controlled transition in a flat-plate boundary layer // J. Fluid Mech. 1995. Vol. 298. P. 211-248.

[7] Скобелев Б.Ю. Конечномерная инвариантная аппроксимация уравнений Навье— Стокса и автоколебательные режимы течения Пуазейля // Приклад. мат. и механика. 1990. Т. 54, вып. 3. С. 416-429.

[8] Годунов С.К. О численном решении краевых задач для систем линейных обыкновенных дифференциальных уравнений // Успехи мат. наук. 1961. Вып. 3(99). С. 171-174.

[9] Йосс Ж., Джозеф Д. Элементарная теория устойчивости и бифуркаций. М.: Мир, 1983. 301 с.

[10] Постон Т., Стюарт И. Теория катастроф и ее приложения. М.: Мир, I980. 607 с.

[11] Бойко А.В., Грек Г.Р., Сбоев Д.С. Спектральный анализ локализованных возмущений в пограничном слое при докритических числах Рейнольдса. Новосибирск, 2002. (Препр./ИТПМ СО РАН; № 1-2002).

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