Модель турбулентности k-ƒ в задаче пограничного слоя на вращающихся телах Текст научной статьи по специальности «Физика»

УДК 533.6

МОДЕЛЬ ТУРБУЛЕНТНОСТИ к,-е В ЗАДАЧЕ ПОГРАНИЧНОГО СЛОЯ НА ВРАЩАЮЩИХСЯ ТЕЛАХ

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

Самарский государственный аэрокосмический университет им. академика С.П. Королева (национальный исследовательский университет)

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

В статье использована модель турбулентности к _ е для решения задач пограничного слоя на вращающихся телах. Модель реализована в системе МЛТЬЛБ и показала хорошее соответствие с результатами экспериментов Парра и Лутандера-Ридберга.

Ключевые слова: к — е модели турбулентности, пограничный слой, вращающиеся тела.

1. ПОСТАНОВКА ЗАДАЧИ

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

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

д¥х д¥у

—- + —— + —- = 0, дх Я дх ду

У_ VI дЯ + у ду= ии +1 Т х дх Я дх у ду дх р ду

У дк + уудЯ+У ду = 1 дт

.(1)

дх Я дх у ду р ду

где Ух,у,Уг - компоненты вектора скорости жидкости, и (х) - скорость внешнего (невязкого) обтекания тела, Я (х) - радиус тела враще-

ния, т = /

дК

ту = /

дУ_

касательные

Рис. 1. Схема обтекания жидкостью вращающегося осесимметричного тела [1]

мой жидкости уравнения модели к _е могут быть представлены в виде:

дК тг дК д Ух — + V — = — х дх у ду ду

(

v + -

'k

дК

ду

+s

( ят/ У (ям V

дУх_ v дУ у

ду ду

напряжения между слоями жидкости. Вязкость / определяется как сумма молекулярной / и турбулентной / вязкостей.

Модель k — s добавляет в систему еще два уравнения - кинетической энергии турбулентности к и диссипации турбулентности Е [3]. В случае тонкого пограничного слоя несжимае-

Куркин Евгений Игоревич, аспирант кафедры аэрогидродинамики. E-mail: [email protected] Шахов Валентин Гаврилович, кандидат технических наук, профессор, заведующий кафедрой аэрогидродинамики. E-mail: [email protected]

ТТ дЕ тг дЕ д Vx — + V — = — х дх у ду ду

дК

v ду у f

- Е;

v + -

дЕ

V

W

ду

(2)

Е

+c , —s

s к »

f ЙТУ \2 f PtV V

V

v ду

л л

где V = —, sm = — = c

дУ1

\ ду у К2

-c

s2

К

р

р

Е

параметры

c^ = 0.09, cei = 1.44, се2 = 1.92, ^ = 1, = 1.3.

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

у = 0: Ух = Уу = 0, У2 = Яш;

к(х,п) = Ц2К(ху), фп) = ЦзЕ(x,У),

К =s,

ал = К

У = Уе : Ух = и (X) У = 0,

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

Модель к -а имеет вычислительные особенности на границе у = 0. Поэтому граничные условия для нее задаются на внешней границе и на границе, расположенной на некотором расстоянии у0 от поверхности тела. Учитывая это, граничные условия для уравнений к -а модели турбулентности можно записать в виде:

+ mfs+s+m [у2 + ОУ ]-а-Ъщык =

= х(икх- & )

(%) + тЛ+с а а+ [у2+О 2 р

у = у0 : Е0 =(ат )

СБ

(Зу„Л2 (ЗУ,^2

ду

У0

ду

-Са2 к - (3т2 - 1)иа = х(иах - Л ) •

Граничные условия системы примут вид: П = 0: Л = 0, и = 0, g = 1;

К 0 = Л Е0

(ат }

СБ

П = П0: а = (ат)с8 [у2 + Ор ],

к=

у = уе:

йКе

Ее йЕе

йх и(х)' йх

= -са2

Е2

■(А

хСм

(5)

КЦ (х )'

где (ат ^са = М- - турбулентная вязкость, оп-р

ределяемая в пристенной области у е [0, у0 ] по какому либо альтернативному закону турбулентности, к примеру по модернизированным формулам модели Себеси-Смита [4].

2. МАТРИЧНЫЙ ВИД УРАВНЕНИЙ

ПОГРАНИЧНОГО СЛОЯ

дк

ц = г!е: и = 1, g = 0, х--+а+2т2к = 0,

дх

да

иь ь / Л

дх + Са2к + (3т2-1)а = 0

а

дх Ь£ к

В этих уравнениях введены обозначения:

Ъ = 1 + а+, О =

шЯ

и '

т2 =

и'

т3 =

Я!

Я '

_ т2 +1 _Цх т = т3 +--2—, ^ех = , безразмерная тур-

Вводя функцию тока /(х,у), безразмерную +

переменную типа переменой Блазиуса булентная вязкость ат в пристенной области

вычисляется по алгебраическим соотношениям

(ат )с

зуется

(ат)

'СБ

V

а при Ц>Ц0 исполь-

соотношение

модели

к-а :

П=уЦ\х и представляя 1//=Я(х))ух(/(х)/(х,ц), У =ш1(х,г/), систему (1) приводим к виду [4]:

и=П у=^ р=gv, (3)

(ЪУ)п + тЛУ + т2 (1 - и2) + тзО2g2 =

= х(иих- Ж) (4)

(ЪР) + т1Лр - 2mзug = х их - Л) •

где Л, g - безразмерные функции. Нижние индексы п и х обозначают дифференцирование ную переменную, характеризующую коэффици-

+ „ к2

ат = СМ Кех~ .

а

Разделение пограничного слоя на пристенную область, в которой турбулентная вязкость описывается алгебраической моделью, и область, где реализуется к - а модель турбулентности, удобно проводить, вводя дополнитель-

по соответствующим переменным.

ент сопротивления трения на стенке. Такой

Уравнения турбулентности можно предста- подход, к примеру, предложен Себеси [3]. Гра-

вить в безразмерном виде, записывая, :

ницу п0 выберем исходя из постоянства вели-

2

чины У+ = 0 г = 50 , где ит = w0U (х). Для определения Пважно значение w на стенке твердого тела (Wo), остальные значения не рассматриваются, а производная w по толщине пограничного слоя считается равной 0 = 0) . Величина w0 может быть добавлена в список

переменных модели в виде w0 =

4 Vo2 + Q2 Po2

Re

1/4

•, а

граница n0 в таком случае рассчитывается как

По = Уо

1

w,

1^/Rex

Система уравнений в частных производных решается методом сеток с использованием конечно-разностной схемы "прямоугольник". Получившиеся нелинейные уравнения линеаризуются по Ньютону [5]. Следующая итерация переменных /, и, V, g, р, к, 8,8, q, w находится в

виде: Л+ = Л + 5 л ■

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

А5 = г. (7)

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

A=

A) со " «о" "0

А с 01 Г

B A C , 8= 8 , r= rj

в« AJ-I J 8J-i rJ-1

в, Aj _ . «j _ _ rJ _

Векторы переменных

^ = \_5fj 5 5 5gj 5р] 5 58] 5е} 5qj 5wj ],

] = 0,...,3.

Блоки А],В],С] - матрицы размером 10 х 10, блоки г - матрицы 10 х 1. Матрицы разрежены. Значения их ненулевых элементов приведены в Приложении.

Для начала работы итерационного метода приближений требуется задания начальных профилей используемых переменных. Начальные профили Л, и, V, g, р могут быть определены как начальные профили для ламинарного течения [4]. Профили же к, 8, 8, q можно задать следующим образом:

к = ф 2 + Q2 p2 ,

S = с,

'iVRe

S+)

V2 + Q2 p2

a

s = кп q = s

3. ПРИМЕРЫ РАСЧЕТА И СРАВНЕНИЕ С ЭКСПЕРИМЕНТАЛЬНЫМИ ДАННЫМИ

Реализация предложенной модели проведена в системе МЛТЬЛБ в модернизированной программе УеЛеЬ

Результаты расчета для тела Парра [1] при О = 1,--.,4 и = 3 -Ю5...9 -105 представлены на рис. 2, 3. Сравнение рассчитанных профилей пограничного слоя с данными эксперимента [1] показывает хорошее соответствие.

Сравнение толщины вытеснения импульса

Уе

0х = | и (1 - и~)с1у на рис. 4 и коэффициента мо-

Рис. 2. Сравнение безразмерных профилей меридиональной и и окружной g компонент скорости пограничного слоя при различных скоростях вращения О, Яе = 3 • 105, х / Ят = 3,1. Сплошная линия - расчеты, пунктир - эксперимент Парра [1]

Рис. 3. Сравнение безразмерных профилей меридиональной и и окружной g компонент скорости пограничного слоя при различных числах Рейнольдса, О = 1, х / Ят = 3,1. Сплошная линия - расчеты, пунктир - эксперимент Парра [1]

Рис. 4. Сравнение толщины вытеснения импульса в меридиональном направлении, Яе = 3 • 105, х / Ят = 3,1. Сплошная линия - расчеты, пунктир - эксперимент Парра [1]

мента силы трения в окружном направлении М

С =

^ Л /Г

М 1 гл2г,5 на рис. 5 также показывают хо-

2 Р Ят

Хорошее соответствие результатов эксперимента Лутандера-Ридберга [2] по исследованию положения точки отрыва потока на вращающемся в осевом направлении шаре при Яе = 1,6 • 10 с результатами моделирования (рис. 6) подтвер-

рошее соответствие расчета и результатов экс- ждает возможность использования описанной

перимента.

модели и ее программной реализации в системе

Рис. 5. Сравнение коэффициента момента трения в окружном направлении, Яе = 3 -105 Сплошная линия - расчеты, пунктир - эксперимент Парра [1]

136 3

134 132 130 128

ч 1 1 1 1 1

1 у \ V \ \

1 ^^^^^ 1 ^^ 1 1 1 V

0 0,25 0,5 0,75 1

а

1,25 1,5

Рис. 6. Сравнение положения точки отрыва потока на вращающемся шаре: сплошная линия - расчет, пунктир - эксперимент Лутандера-Ридберга [2]

Элементы матриц А0, С0 и г0, определим из

МЛТЬЛБ для нахождения точек отрыва потока в задачах вращающихся осесимметричных тел в условий на поверхности тела вращения: осевом потоке.

А1-1 = Д2.2 = А 3,4 = _д5,2 = _ А 6,4 = _д7,6 =

4. ВЫВОДЫ

Модель турбулентности к представлена в форме, предназначенной для описания осесим-метричного пограничного слоя на вращающихся телах, и реализована в системе МЛТЬЛБ. Сравнение результатов расчета по этой модели с данными экспериментов Парра и Лутандера-Ридбер-га показывает хорошее соответствие, что позволяет использовать модель турбулентности для решения задач пограничного слоя на вращающихся телах, в том числе до точки отрыва потока.

ПРИЛОЖЕНИЕ. КОЭФФИЦИЕНТЫ ЛИНЕАРИЗОВАННЫХ УРАВНЕНИЙ ТУРБУЛЕНТНОСТИ

Коэффициенты зависят от значения координаты , определяющей положение точки сетки на толщине пограничного слоя.

= _Д8,7 = _А9,8 = _А10,9 = 1

До3 = А503 = /2, А43 = 2^, А*5 = 2О2р

0

А^5 = _4Яе

р5,2 _ р6,4 _ р7,6 _ р8,7 _ р9,8 _ р10,9 _ 1 С0 _ С0 _ С0 _ С0 _ С0 _ С0 _ 1,

C3,3 = С^5 = _УК/2.

г; = __О2р2, Г05 =(Г4)о, Г06 =(г5)0.

Уравнения импульсов задаются одинаково при всех ] = 1,..., J, поэтому первые 4 строки матриц А^, В^, С] и г}. будут одинаковы во всех областях пограничного слоя:

АУ0 = А21 = 1, А2'2 = _И] /2, А3;1 =(,3);,

А32 =(,5);, А3;3 = (,1);, А- = (,7);,

А41 = (в ), ,А4 2 = (в),, А- = (в),,

АГ = (в)у А-5 = (Д);, ВУ0 = В21 = -1, в;2 = -И]/2, В31 =(,4)в- =(,6)

ВУ 3 =()у В34 =(5, )у В- =(в4 )у

ВУ 2 =(вб ) у, В4 3 = (вю )В4 4 = (Д )

ВТ = (в ); . г; = 0, г; = (/ );, гЗ = (/);,

г;4 = (г3);, первые четыре строки матриц С; -нулевые.

(Г1)(г5); , а также (51); (5*);,

(Д ); ,---,(Д0); соответствуют коэффициентам линеаризации уравнений импульсов в меридиональном и окружном направлениях [4].

Задание модели турбулентности в различных областях пограничного слоя приводит к определению нижних шести строк матриц А;, В;, С;, Г; в каждой из областей отдельно:

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

у = ^.Л - 1:

А52 = А" = = А» = =А;,; =А9*=А;05=-1,

А» = А-5 = : А;- = :а99 = -Йу+1 /2, ,

с;-3 = с6,> = с';' = с? = '--и,+1/2 , г5 =(Г.),,

с« = с" = с;6 = сГ: =с9,*=с;0,'=1,

г;= (/5 )у г; = (Гйк , »у. г; = ()у, коэффициен-

ты матриц нижних шести строк Внулевые,

(к );-1 = к;-1 - к; + ¿Л-!/; , Ы; -1 = У-8у + Ь;Ч;

1]-1/2 .

- граница между различными моделями турбулентности ; = ;0. На этом слое матрица В;0 вычисляется по правилам пристенной области, матрица С 0 вычисляется по правилам для области пограничного слоя с к - 8 моделью турбулентности, а 5-10 строки матриц А;0 и г;0 имеют вид:

А7,2 = А8,4 = А9,6 = А10,8 = -1

А;0 = А;0 = А;0 = А;0 = 1,

А7,3 = А8,5 = А9,7 = А10,9 = -и /2

А;0 = А;0 = А;0 = А;0 = П;+1 1 2,

А53 = П А55 = П А56 = П А5'8 = П

;0 А./0 _ и2> А;0 ~ и3-> А;0 ~ и4->

А6,3 = П А6,5 = П А6,6 = П А6,8 = П

А;0 _ и5> А;0 _ иб> А;0 _ и7 > А;0 _ ^8 ,

У0

(ГВ1 ) ;0 , Г;0 (ГВ2 ) ;0 , Гу0 (Г4 )

Гу0 (Г5 )у0 , Г;0 (Гйк )у0 , г10 (Гйе )у0 ,

д д где П1 = — (<)С8 , П = 8 ф (+т)с8 ,

А = -2Яехсрк;0,Д = (8+)^ ,

П5 = 2Яех ^^0, П6 = 2Яех смк]0О2Р;0,

А = 2Яех с^ 0 [^ + О2 р20 ], Д = -28

(ГП1);0 = хс^к^0 - (<)С5 8 у0

(ГП2 );0 = 820 - Кех сА20 [^0 + ОР20 ] .

У 0 >

- область пограничного слоя с к -е моделью турбулентности, у = ;0 +1, •••, 3 -1:

А72 = А8,4 = А9,6 = А10,8 = -1,

АГ = АГ = АГ = АГ = ^+1 /2;

А5;к =(а2к-1)у ,А6к = (%к-1)у,

В6к =(^2к )у ,В5к = (а2к )у, к = 1_9 ;

с7,3 = с85 = с97 = с1/,9 = -ь}+1 /2, с72 = с84 = с9;6 = с;07 = 1,г5 =()у,

г; = (гБ)у, =(/4)у, г;=(Г5)у, г9=(Г1&)

г0 =(Ге )у ,

где коэффициенты уравнения к : т„ х

а

,) = тт1_+ ^ 1

1 )у = 2 5у + 2 5у-1/2 д/ 1дх) {'

Ч тП п хп п д \д/

а2). ^^1 +—1/2— I — I 2/у 2 у-1 2 у-1/2 д/ |дх) -1,

аз)у = -тПк" - Хх-(|к

2 I дх

у-1/2

«4 = -т1кП-1 -

х(дк_ 2 I дх

у-1/2

«).=1 га, (а.).=1 ет

а9

1 ( дР 1 1

(дР 1

др ),-1,

0

5

(aii) j = hj15

-1 " ^и2

д К ) 1 ( д P

'j 11 1 д к ) j 2 ( д к yj

+ —

12 j

1 (д0 )п х

д (дк

2 (дк )j 2 j-1/2 дк (дху;

(«12 V = - h- S"n-1

д b

д к ) j-1 2 1 д к ) j-1

1 I дР + —

2 1-1 2 1 дк ) j-1 2

д ( д к

• j-1/2

д к I дх

j-1

(а,) j = „ Г(Ь 2 )П + ffп + (f

j-1/2

(«14 ) = - h j-' (b2 )П-1 + j + —1 (f )

(a 1 5) = h -1s" l1^) +1 V 1 5'1 j j1 де ). 2

эр

де

Э0

де

к v=-h- s- 11 ^ее

дЬ^Т +1

де) j- 2

ср у -(д£ у

,де) i-1 (де) J

(a7), = (a), = (кт), = (a8), =0,

- )j=-f (f) П .,„ - 2 (3 - п - 1)еп.

- )=- т (ff ) П .„,- 1 - ое-1

1 ( дР1

- )j=1 (17) J ■ -)j=2 (f

j-1

(-9>J = 2№)J. - )j = 2$)^

(i')j=hJ1 qn №).+1

-)=-hJ 1 q-^it) П-+2

д P1 д к

'CP |

v дк ) t _ 1 ( дк

д 01)

дк j

( д01

j-1

) = hjlqn ^

1 (C01

Се Jj 2 ( де j

2 (,-n -1)

Xn с ( де )

1 2 "j-1/2 Се ( Cx ) j'

)j = -[ (b2 S )П - (b2 S )П-1 ] hj 1 -

- [ PU/2 - 0П-1/2 - 2 - 2 ("к )П-1/2 + -П ( fs )П-1/2 ] +

+х

j-1/2

j-1/2

' j-1/2

' j-1/2

b2 =

+ С к2 h1

= 1 = 1 +^Rex-, р = . Re — [v2 + Q2p2], ak a е " х е [ ]

0 = е, = ^Reх , = ^Rex "к2

дк a.

е де a.

10 = 1, CP = 2 c Re к- [v2 + Q 2 2 ] де дк ц х е [ ]'

CP П - к V 2 гл2 2 "1 т- = c^ Re х — [V2 + Q2 p2 J ,

де е 2

CP о о к2 CP 2 _ к2 _ 2

2c-Re , ф=2c^Re х TQp

уравнение е:

(-16 )j=-h - q -1|

се

1 (с 01

--1

- -1

.1(3-; - j xn»" 1/21е(|е)" ,

' j -1

(-17 )j = h-' (b, )n+^f; + f )

(-18 )■ = - h-1 (b,)" 1 + -nfn1 + f^^

vmz j j V ^ j-1 2 Jj-1 2 ( дх у

(-7 )j =(-)J =(-13 )J =(-14 )J = 0 ,

(E )J = -[(b,q )" - (b,q )"-1J h J- 1 -[ (P )"-1/ (01)"-1/2 + -n (fq )"-1/2 - (3-n - 1)(«е)П

+ x

дх J " qi-1/2 ^ ' j-1/2

j-1/2

(-)= -Lq; + X^Lq" Aff

2 qi + 2 qi-1/2 if (дх j

е+ c к2

b, = 1 + ^ = 1 Re x — a„ ст„ е '

(-2) = -Lq" 1 + ——q"" 1/2 )"

yri)j 2 4j-1 2 4j-1/2 if (с— ) ^

P1 = СеlС^ Re —к [ v2 + Q2 p2 J , 01 = cs2 к,

О=2с £ О=-с £1

д8 82 к' дк" с£2 к2,

дР дР

— = 0, — = 2с81с „Яе хку,

> ^ 81 „ х >

д£ ду

дР 2

-дР- = гс^Яехк°2р ,

^ = с81с„ Яе х [у2 + О2 р21 дк 8 „ х [ ]

- внешняя граница пограничного слоя, у = 3. В3 в таком случае может быть вычислена по правилам для основной части пограничного слоя, матрицы же А3 и г3 примут вид:

АГ =(«2к-1)у, А6-к =(^2к-1)у, к = \,...,9 ;

Л7'2 = Л8, = 1,Af = Ei, Л9/ = E2, J = E3, АГ = E4, rj =(г,, rj =(Ге)j , rj = 0, J = 0 , rj =( rEi )j, J =( Ге 2) j ,

где E = 2< + x" — \ — dk ldx ,

(ГЕ 2 )j = -

de

x" \ -I + ce

д x

J)

k"

(3 m 2 -1)

Вычисление производных проводится по

схеме (/] = А /_-2 + Лг/_-1 + А/_ , где в случае использования первого порядка А1 = 0,

A =•

-1

A3 =

•, в случае второго

порядка A1 =

(x" - 2 Хп-1 )(xn - 2 Xn )

A2 =_

(xn-1 - xn-2 )(xn-1 - xn )

A = 2 xn - xn-1 - xn - 2

3 (xn - xn- 2 )(xn - xn-1 )'

причем производные вида

Äff

df [d x

= A3.

СПИСОК ЛИТЕРАТУРЫ

E3 = Ce2

(kn)

2 ' 2

E = 1,

E4 =(3m2n -1) + xT^-l^^

^ ' йс? I Pbr

deldx

+ ce

2en

fe )j = -

x" \ — | + en, + 2m7nkn

laxj j 2 J

1. Pörr О. Untersuchungen der dreidimensionalen Grenzschicht an rotierenden Drehkörpern bei axialer Anströmung//Applied Mechanics. 1963. Vol. 32. N. 6. P. 393-413.

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

3. Cebeci T. Analys of Turbulent Flows, Elsevier, 2004. 376 p.

4. Куркин Е.И., Шахов В.Г. Пограничный слой на вращающихся осесимметричных телах при их осевом обтекании // Вестник Санкт-Петербургского университета. 2008. Сер. 10. Вып. 4. С. 38-49.

5. Себиси Т., Брэдшоу П. Конвективный теплообмен, М.: Мир, 1987. 590 с.

K-8 MODEL OF TURBULENCE IN THE PROBLEM OF THE BOUNDARY LAYER AT ROTATING BODIES

© 2012 E.I. Kurkin, V.G. Shakhov

Samara State Aerospace University named after Academician S.P. Korolyov (National Research University)

1

n n-1

xn xn-1

n

J

It is suggested k — e turbulence model for solving the boundary layer on rotating bodies. The model is implemented in MATLAB and showed good agreement with experimental results of Parr and Lutander-Rydberg.

Keywords: k —e turbulence model, boundary layer, rotating body.

Evgeni Kurkin, Postgraduate Student at the Aero-Hydrodynamics Department. E-mail: [email protected] Valentin Shakhov, Candidate of Technics, Professor, Head at the Aero -Hydrodynamics Department. E-mail: [email protected]