Пограничный слой на вращающихся осесимметричных телах при их осевом обтекании Текст научной статьи по специальности «Физика»

УДК 533.6 Вестник СПбГУ. Сер. 10, 2008, вып. 4

Е. И. Куркин, В. Г. Шахов

ПОГРАНИЧНЫЙ СЛОЙ НА ВРАЩАЮЩИХСЯ ОСЕСИММЕТРИЧНЫХ ТЕЛАХ ПРИ ИХ ОСЕВОМ ОБТЕКАНИИ

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

Матричная форма уравнений положена в основу программы УегЬе1, разработанной в среде МА^^АВ. Программа позволяет рассчитывать профили скоростей и напряжений трения пограничного слоя при ламинарном и турбулентном режимах потока, а также делать выводы о коэффициентах сопротивления тела и находить точку отрыва потока. Для описания турбулентного режима используется модифицированная модель Себиси-Смита.

С помощью программы УегЪе1 исследовано обтекание шара и полутела вращения. Результаты сравнены с ранее известными.

1. Постановка задачи. Описание пограничного слоя проводится в криволинейной системе координат: х - меридиональная координата, измеряемая вдоль направляющей тела вращения, у - вертикальная координата, измеряемая по вертикали к поверхности тела, г - окружная координата (рис. 1).

Куркин Евгений Игоревич - студент (магистрант) факультета летательных аппаратов Самарского государственного аэрокосмического университета им. акад. С. П. Королева, лаборант-исследователь кафедры аэрогидродинамики Самарского государственного аэрокосмического университета им. акад. С. П. Королева. Автор 12 работ. Научное направление: численные методы аэрогидродинамики, теория пограничного слоя. E-mail: [email protected].

Шахов Валентин Гаврилович - профессор, кандидат технических наук, заведующий кафедрой аэрогидродинамики Самарского государственного аэрокосмического университета им. акад. С. П. Королева. Автор 170 работ. Научные направления: аэродинамика летательных аппаратов, теория пограничного слоя и турбулентности, численные и экспериментальные методы механики жидкости и газа. E-mail: [email protected].

(§ Е. И. Куркин, В. Г. Шахов, 2008

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

<914 14 дН дУу

дх К дх ду

дУх

дх

Ух.

дУI

дх

V;2 дя

Я дх 414 дН Я дх

0,

дУх = ду т. дУг у ду

І/— + 1 дТхУ дх р ду ’

1 дтгу

(1)

р ду

где Ух, Уу, Уг - компоненты вектора скорости жидкости; и(х) - скорость внешнего (невязкого) обтекания тела; Н(х) - радиус тела вращения; тху = №^4 тгу = ~~

касательные напряжения между слоями жидкости. Вязкость определяется согласно модифицированной модели Себиси-Смита как сумма молекулярной ц и турбулентной ^ вязкостей.

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

при у = 0: Ух = Уу =0, Уг = Яш,

при у = уе: Ух = и(х), Ух = 0,

(2)

где уе - толщина пограничного слоя.

2. Матричный вид уравнений пограничного слоя. Введение функции тока ф{х,у) : Ух = Уу = ~обращает уравнение неразрывности в тождество. Остав-

шиеся уравнения системы (1) приводятся к безразмерному виду переменной типа переменной Блазиуса г/ = • Течение жидкости в пограничном слое может быть описано

безразмерными функциями /(ж, гу), д(х,г]), такими, что ф = К\п^си/, [Ух(х,г]) = и/п\, Уг = шКд и их производными

и /п, ^ /пп, Р дп' (3)

В результате перехода к безразмерному виду система (1)-(2) примет вид

(Ьу)п + ші/у + Ш2 (1 - и2) + шз9?д2 = х (иих - у/х) , (Ьр)п + ш1/р - 2шзид = х (идх - р/х)

(4)

с граничными условиями

при п = 0: / = 0, и = 0, д =1, при п = Пв'- и =1, д = 0

(5)

и коэффициентами

Ь = 1 + Є

+ т,

Ші = Ш3 +

то2 + 1 2 ’

и'

то2 = —ж,

Я

то3 = —ж,

П- —

и ■

Коэффициент турбулентной вязкости єт в случае ламинарного пограничного слоя равен нулю. При расчете турбулентного пограничного слоя коэффициент єт считается отдельно для внутренней (пристенной) є+ті и внешней єто областей пограничного

слоя [2], причем внешняя область пограничного слоя е+ = е^о начинается с момента превышения е^г над ет,о. Используя [3], формулам модели Себиси-Смита можно придать вид

Є+* =0, 16^

і і у

1 — ехр ( — —

\/Лех\/ь2 + П2р27(Г7, є+ 0 = 0,0168у^ё^|??е - /е|

77*г

гДе А = ^ (ае®)4 #; N = (1 - И, 8р+)2; р+ = (11еж) 4 т2Я; Б = \Л>о + ЇЇ2РІ; 11е;

~г; 7*г = 1 — ехр

-Сіг (Х-ХІГ)Ц^

,+ )2. п+ _

С4г = 8,33-1(Г4^(11еж)

о Т “ Р0>

0,66

х/ ; ^о, р0 - значения

функций V и р на границе тела вращения.

Система уравнений в частных производных (3)-(5) решается методом сеток с использованием конечно-разностной схемы «прямоугольник» [2, 4]. Получившиеся нелинейные уравнения разрешаются методом последовательных приближений Ньютона. Следующую итерацию представим линейно через предыдущую:

/*+1 = / + 6/*, = и* + 5п\, V+1 = Vі* + 6уІ; , д*+1 = д* + 6д*, р!+1 = р* + 5р\.

Итоговую систему линейных алгебраических уравнений удобно рассматривать в матричном виде:

А6 = г. (6)

Матрица коэффициентов системы (6) имеет блочную трехдиагональную структуру:

А

Ао Со 6о го

В1 А1 С1 61 г1

Bj аЬ Су ,6 = 6; , г = гj

В ч 1 А ч 1 С— 6.1 -1 ГJ-1

BJ AJ 6J ГJ

. (7)

Блоки Aj, Bj, О, 6у имеют вид г ь:

А

о =

6/ 6и 6V 6gj 6Pj \Т,з = 0.. з; го = [ 0 0 0 (г4)1 (г5)1 ]Т

(r1)j (г2 )з Ы j (г*Ь {r5)j ]Т,І = 1..3, ГJ = [ (rl)J (г2)J (r3)J 0 0 ]Т

1 0 0 0 0 1 ~Ну 2 0 0 0

0 1 0 0 0 (йз)Ь (sъ)j (sl)j (s7)j 0

0 0 0 1 0 ,А = (вз ^ (Ръ)ь ш3 вь (в1); , 3 = 1..3-1,

0 - і —Л-1 1 2 0 0 0 -1 -к, 2 0 0

0 0 0 -1 — Ні 2 J 0 0 0 -1 — Ну 2 .

1 — hJ 2 0 0 0 0 0 0 0 0

(s3)J (s7)J 0 0 0 0 0 0

AJ = (вз^ (в5)J (в9 ) J (в7^ (в1) J , С = 0 0 0 — Ну 2 0 0 1 - О 3

0 1 0 0 0 0 1 0 0

0 0 0 1 0 0 0 0 1 — Ну 2 -1

2

У

г

У

Б,

1 Ну 2 0 0 0

(«4) 3 («в)з («2)3 («8)3 0

)4 (в4 (вв), (в10)з (в2),

0 -1 — Ну 2 0 0

0 0 0 -1 — Ну 2

1..7.

(8)

Коэффициенты матриц (8) равны

( \ _ і і іі , гп1 + ап 4 и.п х

(ві). - ^ 6^+ - із~-^

,П

/ \ т1 + ап і ап п-1

(5з)7 = -------о------У3 + "У V1/2’

/ \ _ і - її® . т1 + ап лі —/і. т-

^2)3 — ~)1з Ьз-И о

Хп т— 1

2

■* 3-І 2 І-1/2 ’

2

+ а„ («4)7' =--------^-------«

і „.і—1

(«б), = - (т^ + ап) и,, («в), = - (т^ + а„) и,—і, (вг), = т^ід,, (в8), = т^П2^— і;

“пгп-1 ч _ и-іг, г , К + <*

(Рі)з - «3 ьз +---2---(Р2>з-~Ыз Ь3-1 +

_____п_гі ^п лп-1

2 ■'і-1 2 і-1/27

ті + а.

ҐО \ _ ""і І ^і і , ^І і—і

(&), - ----------2-----^ + ~2~рз-і/2>

в), =

ті + ап 2

І „і— 1

■^•-1 + -у^-1/2

2ті і

(ДО, = + ^7-1,2’ ш, = -

2ті і

і „і—1

г **' п— х

9] — 1 ^ ^'-1/2

(/Зг)7 = -

2ті і

і і— 1 М3 “ “^3-1/2’

(/%),• = -

2ті і

і і—1

(во), =0, (віо), =0;

(г1), = /,—і -1, + ^зUj—1|2, (г4),—і = и,— і -и, + З^ — і^і (г5), — і = д, —1 -д, + hjpj — 1|2, (г2), = ^3 —1|2 - 1 (Ь)^3 - Ц — 1^' — і) + (ті + аі) (/^3 —1|2 - (ті + аі) (м2)^ —1|2 +

+ ті^і (д2)^ 1/2 + аі (3172/З—1|2 - З-т^' —1

(гз)з = ^Г—1|2 - ^ 1 {Ь3Р3 - ЬЗ—іРЗ—і) + (ті + аі) (/р)3 —1|2-(2ті + аі) ^Уз—іП-

— ап и

■і—1 д3 — ді—1 и3 + /і—1 Р — Рі—1 Я

з —1|2дЗ —1|2 дЗ —1|2 иЗ—1|2 + і, —1|2РЗ—1|2 pj — 1|2/j — 1|2

і-1 иі

і-1 і

і-1

где

і-1 31|2

і-1

Lj—1|2

т і—1 1j — 1|2

-Ь

— 1 + а

з— 1|2 + аі 1

і-1 ( 2)і-1

(/,і|2 - (и )з—1|2

1(Ьзvj-Ьз—і,0+ті (М—|+т2 1 - (У)3—1|2 + тз^2 (д2)3—1|2}

1

і-1 31|2

1

(/Р)і——1|2 - (ид)—|2

1

і — 1

мі— і1|2 = {Ь— 1 (ЬзРз-Ьз—іРз—і)+ті (/р)з—1|2-2тзидз—1|2}

Матричное уравнение (6) позволяет последовательно находить приближения величин следующего слоя сетки по их значениям на предыдущем слое. Такой подход не применим к первому слою, поэтому зададимся в качестве его начального приближения аппроксимацией многочленами начального профиля скоростей и напряжений, удовлетворяющими граничным условиям (5):

а

41

п

п

2

2

II

т

2

4 \Пе

3-U-

2 \?7е

3^_

2%

2

- Ё J_

2 Г7е

1-1^

Пе

gj — 1 - Uj ,

Pj =

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

Sfo — 0, 5щ — О,

— О,

Suj — О, Sgj — О.

3. Программа Уегіеі. Предложенная модель описания пограничного слоя реализована в среде МА^^АВ 7.0 программой УегЪе1 [5]. Результат программы - безразмерные функции f, и, V, д, р, физически интерпретируемые в дальнейшем отдельными процедурами авторов. Все переменные МА^^АВ имеют двойную точность. Алгоритм работы УегЪе1 представлен на рис. 2.

Begin (Vertel)

v_ input ввод модели из файла или формы

X

zero_profile задание начального приближения первого слоя многочленами (11)

I

Заполнение сетки решением for ind.n=l:num.x

1) расчет слоя layer_n

2) проверка точки отрыва по касат. напряжению: v3 — 0, если да-остановка цикла

С ЕпсКУсПс!)

Рис. 2. Алгоритм программы УегЬе1.

Процедура у_три1 производит постановку задачи: формирует сетку переменных, вычисляет коэффициенты. Начальное приближение задается процедурой 2его_ргоЯ1е. Циклическое вычисление последующего слоя методом последовательных приближений 1ауег_п (рис. 3) заполняет сетку х, п решением.

Трехдиагональный вид матрицы коэффициентов (7) позволяет использовать частный метод решения системы линейных алгебраических уравнений - метод матричной прогонки [2, 5]. Он дает выигрыш машинного времени по сравнению с решателем м\м системы МА^АВ в 3 раза, а по сравнению с решателем МА^АВ для разряженных матриц в 2,3 раза.

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

Изменение функций /, и, V, д, р при турбулентном течении существенно более сильное у стенки тела, чем вблизи внешней границы пограничного слоя. В этом случае используется неравномерная сетка переменной п, шаг которой нарастает как геометрическая прогрессия. Результаты вычислений показали, что наиболее удобным знаменателем прогрессии является значение 1.05, а начальный шаг у стенки тела равен 0.075.

4. Тестирование программы Уег1е1. Оно проводилось путем сравнения результатов этой программы с известными точными решениями уравнений Навье-Стокса.

2

2

2

1

v

j

Рис. 3. Алгоритм процедуры 1ayer_n.

Переход между слоями сетки по переменной х проверялся расчетом осесимметричного обтекания вращающегося цилиндра [4]. Достоверность расчета передней критической точки затупленных тел подтверждается сравнением параметров обтекания передней критической точки вращающегося шара единичного радиуса и результатов Тиффорда и Чжу Шень-до [1] осевого обдува вращающегося диска. Замена а = ^ и

С = У\/= ^7л/1 + ^ обеспечивает единый вид результатов обоих задач. Сравнение профилей скоростей в критической точке при расчете пограничного слоя с равномерным распределением 300 точек с шагом 0,01 по п показано на рис. 4.

Приведенные коэффициенты сопротивления осевому потоку равны

С = сги -^= \/ а2 + и2 = 2^о \ 1 + ^21 3^4 -

- / V 1 J

Сравним коэффициенты сопротивления вращающегося диска осевому потоку СуеПе\ с результатами Тиффорда и Чжу Шень-до С [ 1]:

17-1 0 0,1 0,25 0,5 1 1,5 2 4 6 оо

С 1,03 1,05 1,14 1,38 1,83 2,17 2,33 2,52 2,57 2,61

Д/е пеі 1,021 1,037 1,124 1,373 1,871 2,171 2,333 2,539 2,585 2,624

01234567 01234567

Рис. 4- Профили осевой (а) и окружной (б) проекции скорости при обтекании диска. Штриховая линия - результаты Тиффорда и Чжу Шень-до [1], сплошная - Уег1е1.

Эти коэффициенты также характеризуют высокую точность результатов расчета программой УегЪе1.

5. Исследование обтекания полутела вращения. Полутело вращения образовано наложением пространственного источника и однородного потока (рис. 5).

Параметрическое задание геометрии полутела вращения и скорости его невязкого обтекания дано в [1]. Рассчитав в программе УегЪе1 безразмерные величины касательных напряжений V и р, определим сопротивление трения полутела вращения набегающему потоку, а также момент сопротивления его вращению.

Коэффициент сопротивления тела вращения вычислим по формуле

Ал/й [ха [ХР с Ха 1г = „2 р2 \j-VoRdx.

Зависимость коэффициента сопротивления от относительной скорости вращения тела при разных его длинах, полученная программой УегЪе1 и определенная Шлих-тингом [1] с помощью приближенных методов, представлена на рис. 6, где число

Рис. 5. Схема обтекания полутела вращения.

Рейнольдса Ие =

и ос Ягг

Ут = шЯт. Результаты программы удовлетворяют теоре-л/Йё^ха гг от Тт^у а значения коэффициента

тически полученному ВИДУ зависимости \/КеЬха £г ОТ ТТ

и °Р

сопротивления хорошо соответствуют ранее известным [1]. Коэффициент момента трения вычисляется так же:

См = —

47Гу^ Ґ

ияът Л

-Я ро сіх.

Сравнение результатов со значениями Шлихтинга [1] представлено на рис. 7.

Рис. 6. Коэффициент сопротивления при ламинарном режиме обтекания полутела вращения, полученный программой УегЬе1 (а) и приближенно Шлихтингом [1] (б).

Рис. 7. Коэффициент момента трения при ламинарном режиме обтекания полутела вращения, полученный УегЬе1 (1) и приближенно Шлихтингом [1] (2).

6. Исследование обтекания шара. Проведем сравнение ламинарного и турбулентного пограничного слоев при обтекании вращающегося шара (рис. 8).

V

Число Рейнольдса Ие = 2Д"^/ос = 1.3 • 106 допускает существование как ламинарного, так и турбулентного пограничного слоев. Профиль скорости турбулентного режима течения (рис. 9) быстрее нарастает в пристенной области. К тому же такой профиль медленнее достигает величины невязкого обтекания, что показывает большую толщину турбулентного пограничного слоя. На ряде сеток используемая турбулентная модель вносит неустойчивость в решение. Сглаживание вычисления коэффициента полиномом пятой степени стабилизирует решение и позволяет получать гладкие профили. Хотя точка отрыва потока заметна и на профиле скорости, механизм ее появления удобнее изучать, основываясь на значении возникающих в жидкости напряжений, стремящихся к нулю на границе стенки и сосредотачивающихся во внутренней области пограничного слоя (рис. 10).

Рис. 8. Схема обтекания шара.

Рис. 9. Профили безразмерной скорости пограничного слоя в осевом и (а) и окружном д (б) направлениях (И,е =1.3 • 106, Ут= 1).

Сплошная линия - турбулентный режим, штриховая - ламинарный.

Напряжения, возникающие в жидкости на границе тела, определяют действующие на тело силы, а их равенство нулю является индикатором отрыва потока. Достоверность получения графиков поверхностных напряжений подтверждается их сравнением с результатами Шлихтинга [1] (рис. 11). Характер течения сильно влияет на поверхностные напряжения. При неизменном числе Рейнольдса касательные напряжения при турбулентном режиме вдвое больше, чем при ламинарном. Возрастание числа Ие на порядок увеличивает возникающие напряжение еще вдвое (рис. 12).

Отрыв потока при ламинарном режиме происходит при в8 ~ 106° (рис. 13). При турбулентном режиме положение точки отрыва существенно смещается вниз по потоку. Увеличение угловой скорости вращения шара вызывает малое смещение точки отрыва потока вперед как при ламинарном, так и при турбулентном режиме течения. Полученные результаты хорошо согласуются с замерами Лутандера-Ридберга [1].

Рис. 10. Пристенная часть профилей безразмерного напряжения пограничного слоя в осевом V (а) и окружном р (б) направлениях (И,е = 1.3 • 106, Ут/П^ = 1).

Сплошная линия - турбулентный режим, штриховая - ламинарный.

Рис. 11. Сравнение продольного (а) и окружного (б) касательного напряжения на поверхности вращающегося шара при ламинарном режиме.

Сплошная линия - решение Уег1е1, штриховая - результаты Шлихтинга [1].

9, град. 0, град.

Рис. 12. Касательное напряжение на поверхности вращающегося шара в осевом (а) и окружном (б) направлениях при его турбулентном обтекании.

05, град.

Рис. 13. Положение точки отрыва потока.

Увеличение скорости вращения осесимметричного тела изменяет не только положение точки отрыва потока, но и влияет на характер профилей скорости его продольного обтекания. При скорости вращения поверхности тела, не превышающей скорости его осевого обтекания, продольная скорость монотонно возрастает до значения на внешней границе (рис. 14). Быстрое вращение осесимметричного тела приводит к «вентиляторному» эффекту - превышению скорости течения внутренних слоев жидкости пограничного слоя над скоростью невязкого обтекания тела. На рис. 14 представлены продольные проекции скорости на направление внешнего потока при турбулентном и ламинарном режимах. Пограничный слой при турбулентном режиме имеет большую толщину, и «вентиляторный» эффект быстрее затухает при смещении вниз по потоку.

I-------------1-------------1------------1-------------, I__________I__________I___________I__________I__________I___________L_

0 2 4 0 2 4 6

Рис. 14. Профили продольной проекции скорости при «вентиляторном» эффекте в ламинарном (а) и турбулентном (б) режимах (Re = 1.3 • 106, Vm/U^ = 10).

Summary

Kurkin E. I., Shakhov V. G. Boundary layer of axial flow around revolving axisymmetric bodies.

A system of non-linear parabolic Prandtl boundary layer equations is brought to a finite difference form and solved using successive approximation method. A “Vertel” code, based on the matrix form of equations, allowed to calculate velocity profiles and friction stress values for both laminar and turbulent flows. The turbulent mode is described using a modified Cebici-Smith model. The author gives a benchmark analysis of Vertel calculations for a revolving sphere and a half-body of revolution in comparison with earlier results.

Key words: laminar flow, turbulent flow, axial flow, boundary layer, Prandtl equations, Cebici-Smith model.

Литература

1. Дорфман Л. А. Гидродинамическое сопротивление и теплоотдача вращающихся тел. - М.: Физ-матлит, 1960. - 260 с.

2. Себиси Т., Брэдшоу П. Конвективный теплообмен | Пер. с англ.; Под ред. У. Г. Пирумова. - М.: Мир, 1987. - Б90 с.

3. Шахов В. Г. О гипотезе турбулентности в пространственных пограничных слоях II Труды Куйбышевск. авиац. ин-та. - 1971. - Вып. ЗБ. - С. 102-109.

4. Куркин Е. И., Шахов В. Г. Турбулентный пограничный слой при осевом обтекании вращающегося осесимметричного тела 11 Сб. трудов XIII Всерос. науч.-техн. семинара «Управление движением и навигация летательных аппаратов». - Самара: Изд-во Самарск. аэрокосм. ун-та, 2007. - Ч. II. -С. 32-38.

Б. Куркин Е. И., Шахов В. Г. Исследование пограничного слоя при осевом обтекании вращающегося осесимметричного тела в среде MATLAB 11 Труды Всерос. науч. конференции «Проектирование научных и инженерных приложений в среде MATLAB». - СПб.: Изд-во С.-Петерб. ун-та, 2007. -C. 29Б-307.

Статья рекомендована к печати проф. Е. И. Веремеем.

Статья принята к печати 29 апреля 2008 г.