Нахождение стационарных решений задачи о конвекции Рэлея - Бенара - Кармана Текст научной статьи по специальности «Физика»

УДК 532.5 Б01 10.18522/0321-3005-2015-4-44-48

НАХОЖДЕНИЕ СТАЦИОНАРНЫХ РЕШЕНИЙ ЗАДАЧИ О КОНВЕКЦИИ РЭЛЕЯ - БЕНАРА - КАРМАНА*

© 2015 г. Т.Ф. Долгих, М.Ю. Жуков, Е.В. Ширяева

Долгих Татьяна Федоровна - аспирант, кафедра вычислительной математики и математической физики, Институт математики, механики и компьютерных наук им. И.И. Воровича Южного федерального университета, ул. Мильчакова, 8а, г. Ростов-на-Дону, 344090, e-mail: dolgikh@sfedu. ru

Dolgikh Tatiana Fedorovna - Post-Graduate Student, Department of Calculus Mathematics and Mathematical Physics, Vorovich Institute of Mathematics, Mechanics and Computer Science of the Southern Federal University, Milchakov St., 8a, Rostov-on-Don, 344090, Russia, e-mail: [email protected]

Жуков Михаил Юрьевич - доктор физико-математических наук, профессор, заведующий кафедрой вычислительной математики и математической физики, Институт математики, механики и компьютерных наук им. И.И. Воровича Южного федерального университета, ул. Мильчакова, 8а, г. Ростов-на-Дону, 344090, e-mail: [email protected]

Ширяева Елена Владимировна - кандидат физико-математических наук, доцент, кафедра вычислительной математики и математической физики, Институт математики, механики и компьютерных наук им. И.И. Воровича Южного федерального университета, ул. Мильчакова, 8а, 344090, г. Ростов-на-Дону, e-mail: [email protected]

Zhukov Mikhail Yurievich - Doctor of Physical and Mathematical Science, Professor, Head of Department of Calculus Mathematics and Mathematical Physics, Vorovich Institute of Mathematics, Mechanics and Computer Science of the Southern Federal University, Milchakov St., 8a, Rostov-on-Don, 344090, Russia, e-mail: zhuk@math sfedu.ru

Shiryaeva Elena Vladimirovna - Candidate of Physical and Mathematical Science, Associate Professor, Department of Calculus Mathematics and Mathematical Physics, Vorovich Institute of Mathematics, Mechanics and Computer Science of the Southern Federal University, Milchakov St., 8a, Rostov-on-Don, 344090, Russia, e-mail: [email protected]

Численно исследована задача Рэлея - Бенара - Кармана о стационарных режимах конвекции для цилиндров с отношениями высоты к радиусу 0,25 < 2к!Я < 4. На основе метода конечных элементов в сочетании с методом Ньютона реализован алгоритм прямого расчета стационарных режимов и определения критического параметра потери устойчивости.

Ключевые слова: конвекция Рэлея - Бенара - Кармана, стационарные режимы течения жидкости, потеря устойчивости.

The problem on the Rayleigh-Benard-Karman convection in cylinder is studied. Qualitative change of stationary modes offluid flow is investigated. To solve the problem we use the finite element method and the Newton method. The results of calculations for different relations height of the cylinder to its radius is presented.

Keywords: Rayleigh - Benard - Karman convection, stationary regimes offluid flow, instability.

Для определения критических параметров потери устойчивости основного режима течения в задаче о конвекции Рэлея - Бенара - Кармана в работе [1] в случае, когда отношение высоты цилиндра к радиусу равно единице, использован алгоритм численного решения нестационарной задачи в сочетании с движением по параметру (числу Рэлея). Обычно для этих целей традиционно используется решение спектральной задачи теории гидродинамической устойчивости [2]. Однако в данном случае такой подход сопряжен со значительными трудностями. Явные соотношения для основного режима течения отсутствуют, и его приходится определять численно, и лишь затем решать задачу устойчивости, линеаризованную на приближенном

решении исходной задачи. Ниже показано, что гораздо эффективнее строить численное решение стационарной задачи при помощи метода конечных элементов в сочетании с методом Ньютона. Дополнительное преимущество такого метода заключается также в том, что при движении по параметру легко следить за качественными изменениями структуры течения.

Основные уравнения

Уравнения, описывающие конвекцию Рэлея -Бенара - Кармана в вязкой несжимаемой жидкости в цилиндре, вдоль оси которого действует сила тяжести, с предположением вращательной симмет-

*Работа выполнена при финансовой поддержке базовой части технического задания 213.01-11/2014-1 Минобрнауки РФ, Южный федеральный университет, и в рамках проектной части государственного задания в сфере научной деятельности (задание № 1.1398.2014/К).

рии и стационарности, в цилиндрических координатах и безразмерных переменных имеют вид [1-5]

V7 V 1 MV1

v-Vv--k„ +— k, =

w 1 Ra

= -vp--rot rot v ч--- вk^

Re PrRe

(1)

-v-k +v-V6> = -

1

PrRe

Ав ,

n л 1 5 5 52 divv = 0 , A =--r--1--- ,

r dr dr dz

^ 5 d 15 dw

V = k).--hkz—, divv =--гмн--

dr dz r dr dz

л л л (1 ЗА dJL

rot Ar,A.,A7 =--^--*-

r Ф z Л Л. Л

l г оф dz

k.+

dAr dAz dz dr

(

S гАф dr

d4 d</>

\

На основаниях цилиндра, которые вращаются с одинаковой по величине О (П = 1), но противоположной по направлению угловой скоростью, заданы условия прилипания и температура, а на боковых стенках цилиндра — условия прилипания и отсутствие потока тепла

и = м/ = 0 , у = Ог , в = 0 , г = к ; (2)

и = м> = 0 , V = -От , О = О , г =-к ; и=м>г =0 , V = 0 , вг =0 , г = 0 ; и = м> = 0 , у = 0 , вг =0 , г =К . Здесь у= и. у. м' - скорость: р - давление: в — отклонение от равновесной температуры Т0 г ; к . кф, к. орты цилиндрической системы координат; Я, Ъ — радиус и полувысота цилиндра.

В качестве равновесной температуры Т0 г выбрано распределение

Т z =

Тп h-z

2h

Т =1 -'о 1 •

(3)

К =м„,

где g, ~ ускорение силы тяжести; О,, Н ,. Т , -величины угловой скорости, радиуса цилиндра и разности температур; 8, — коэффициент температурного расширения; /3, - температуропроводность; //, — кинематическая вязкость; А, - полувысота цилиндра; р, — плотность жидкости.

Метод решения

Для отыскания стационарных решений нелинейную задачу (1)-(4) приводим к операторному виду, который обычно используется для решения задачи Стокса методом конечных элементов

Р и,у,м>,р,е =/1+/2+/3=0 , (5)

где /[ = | —г V-Уу -У + у211 -т>У с1гс1г ,

в

/2 = |r —Р div v — Vp • V + w© drdz

Л =11 ¿i-^WV-UU + VV

1 Ra

--rV<9-V0 +--rßW \drdz .

где Т0 — температура нижнего основания цилиндра.

Такой выбор равновесной температуры соответствует решению задачи V - 0 при О - О.

Т° -к = Т0, Т° к = 0. Температура в области имеет вид

Т х,г =Т° г +6 х,г . (4)

Числа Рейнольдса Re, Прандтля Рг и Рэлея Яа, а также размерные и безразмерные величины связаны соотношениями

/4 д,[м, 5,

Re PrRe'

Здесь V = U,V,W , Р, 0 - пробные функции. Вычисляя производную Фреше, получим

DF u,v,w,p,0 Su,öv,Sw,Sp,Se = (6)

= DIl +DI2 +DI3,

DIX = J -r <5v-Vv + v-Vi>v -V +

D

+2vSVU — Suv + uSv V drdz ,

DI2 = JV —Pdivc>v-VSp-V + Sw® drdz ,

1 Ra

--rVSe-W& +-—rSeiV \drdz .

Re PrRe

Для нахождения решения применяется метод Ньютона

ОР Ук,рк,вк д\,8р,36 =Р Ук,рк,ек ,

к+1 к У,р,в = \,р, в - д\,др,дв .

Обратим внимание, что при решении задачи (5) не использовался стабилизирующий член и не осуществлялось преобразование члена Р div V в VI' ■ у при помощи интегрирования по частям.

D

D

D+ = Ü<r<RO<z<±h

Результаты вычислительных экспериментов

Численное решение получено с помощью пакета конечных элементов FreeFem++ [6]. Зеркальная симметрия и,у,м>,р,9 г,г = и,у,-м>,р,-в г,-г решения задачи (1)-(3) контролировалась с использованием коэффициента симметрии функции тока ц/

£= \4>с1гс1г / \4drdz

/ о_

Обеспечение возможной симметрии относительно г — () при вычислениях осуществлялось выбором симметричной квадратной сетки. При расчетах с относительной погрешностью, не превышающей 1(Г9, задавалось 1800 кусочно-квадратичных элементов. Вычисления проводились с некоторым шагом по Я а. а все остальные параметры фиксированы ( Яс = 90. Рг = 1).

Для визуализации решения использовалась функция тока цг = г , где

ду-А-ч* = 5 .

На рис. 1 приведены изолинии функции тока для различных значений числа Рэлея.

Хорошо видно, что на интервале 8540 <11а< 8640 происходит изменение структуры течения - режим с двумя вихрями теряет устойчивость и переходит в режим с четырьмя вихрями. Обратим внимание на то, что потеря устойчивости происходит без потери зеркальной симметрии. В частности, это означает, что если в качестве биффуркаци-онной кривой выбирать зависимость коэффициента симметрии от числа Рэлея: .V = .V Я а , то такой переход между режимами не будет наблюдаться.

На рис. 2 показаны изолинии скоростей их,2 ,

УХ, 2 , ^ Х, 2 , функции

тока у/ х,г , температуры в х, г , азимутальной компоненты вихря скорости ^ х, 2 и давления р х, г для значений числа Рэлея 11а = 17000 в случае 2Ь/ Н - 0.25. Вплоть до значений Яа «18000 основной режим не теряет устойчивости, что находится в хорошем соответствии с асимптотическими результатами работы [4].

Рис. 1. Изолинии функции тока для различных значений Яа, = 4

Рис. 2. Изолинии характеристик течения Ra = 8640, 2h/R = 0,25

Асимптотическая модель

Точность и достоверность полученных результатов позволяет контролировать асимптотическая модель, построение которой в основном следует работам [4, 7]. Для поиска вращательно-сим-метричных режимов (о / од = 0) в случае, когда высота цилиндра много меньше радиуса (г »¡г), используем асимптотическую модель. Ищем решение задачи (1)-(4) в виде

г = R + hx, z = hy, р = hp, t = ht, (7)

u=hu, w = hw, v = v, Ju = hJu0, S = hS0, p = h-lp0, Mo = Oil), S0 = O( 1), J30=0( 1),

T = T<0\z) + he, |L = r = h-1y0, y0=O(l), 2 h

2 h 2 ° где T(0) (z) - «равновесный» профиль температуры; у — градиент температуры.

Если подставить (7) в (1)-(4) и пренебречь членами более малого порядка по h, чем сохраненные, то получим, опуская знак «~», уравнения

Вычисляя a (О), получим

и,-— = -рх +дА0И , К

(8)

= Qry + 2QR^T

(-1)*

■ежкх sin яку .

(10)

к=1 л к

Ряд (10) суммируется. Действительно, вводя комплексную переменную ^ = х + гу и функцию

к=i лк

дифференцируя

получим сумму геометрической прогрессии с показателем (/ = ех'. и после суммирования и интегрирования с учетом условия на бесконечности будем иметь

к=1 1 + е 1

F(g) =--\а(\ + еж().

л

Вычисляя мнимую часть функции F(£), окончательно получим

ein -п-,,

(11)

2QR ежх sin жу v0=Çlry--arctg-

w, =-py+^\w + pae, ux+wy= 0, et-y0w= SoAo0, vt=MoA0v, A0=()„+()^ и краевые условия

u = w = 0, 0 = 0, y = ±l, (9)

v = ±Qr = ±a{R + hx), y = ± 1,

u = w = 0, v = 0 , x = -R/h = -R0,

0Х = О, x = -R/h = -R0,

и = wx=0, v = 0 , 0X = 0, x = 0.

В асимптотической модели (8), (9) задача для определения скорости v отщепляется и может быть проинтегрирована независимо от остальных уравнений.

Ищем стационарное решение задачи для

v = v0 (х, у) в виде v0 = Qгу + (х) sin яку. Для

к=1

определения ак(х) имеем, заменяя R0 на +оо при /г —>0, ак"(х)-я1к1ак(х) = 0 , ак(^) = ак(-оо) = 0 . Тогда ак(х) = ак(0)елЬ .

Величины ак (0) определяются краевым

условием для v0 (0, у) = 0 0 = QRy + 'Ур,.. (0) sin яку .

к \ + ежх cos я у

Профиль азимутальной скорости, рассчитанный для полной и асимптотической моделей, отличается от профиля, заданного формулой (11), не более чем на 3 %, начиная со значений параметров R/h> 3.

Литература

1. Bordja L., Tuckerman L.S., Witkowski L.M., Navarro M.C., Barkley D., Bessaih R. Influence of counter-rotating von Karman flow on cylindrical Rayleigh-Benard convection // Phys. Rev E., 2010. Vol. 81. P. 036322.

2. Гершуни Г.З., Жуховицкий Е.М. Конвективная устойчивость несжимаемой жидкости. М., 1972. 392 с.

3. Жуков М.Ю., Ширяева Е.В. Расчет стационарных режимов конвекции Рэлея - Бенара - Кармана // Современные проблемы МСС : тр. XVI Междунар. конф. Т. 1. Ростов н/Д., 2Q12. С. 99-103.

4. Молодчаная А.В., Ширяева Е.В. Асимптотическая модель конвекции Бенара - Кармана для тонких цилиндров // Современные проблемы МСС : тр. XV Междунар. конф. Т. 2. Ростов н/Д., 2Q11. С. 176-180.

5. Долгих Т. Ф. Вычисление бифуркационных кривых для стационарной задачи конвекции Рэлея - Бенара - Кармана // Современные проблемы МСС : тр. XVll Междунар. конф. Т. 1. Ростов н/Д., 2Q14. С. 166-170.

6. Жуков М.Ю., Ширяева Е.В. Использование пакета конечных элементов FreeFem++ для задач гидродинамики, электрофореза и биологии. Ростов н/Д., 2QQ8. 256 с.

7. Владимиров В.А., Дениссенко П.В., Жуков М.Ю., Юдович В.И. Течение вращающейся жидкости в зазоре между двумя параллельными пластинами // Мат. моделирование. 2QQ1. Т. 13, № 2. С. 27-38.

v

2

k=1

References

1. Bordja L., Tuckerman L.S., Witkowski L.M., Navarro M.C., Barkley D., Bessaih R. Influence of counter-rotating von Karman flow on cylindrical Rayleigh-Benard convection. Phys. Rev E, 2010, vol. 81, p. 036322.

2. Gershuni G.Z., Zhukhovitskii E.M. Konvektivnaya us-toichivost' neszhimaemoi zhidkosti [Convective stability of incompressible fluid]. Moscow, 1972, 392 p.

3. Zhukov M.Yu., Shiryaeva E.V. [Calculation of steady-state regimes Rayleigh-Benard-Karman convection]. Sovremennye problemy MSS [Modern problems of MCC] : materials of XVI Intern. Conf. Vol. 1. Rostov-on-Don, 2012, pp. 99-103.

4. Molodchanaya A.V., Shiryaeva E.V. [Asymptotic model Benard convection-pocket for thin cylinders]. Sovremennye

problemy MSS [Modern problems of MCC] : materials of XV Intern. Conf. Vol. 2. Rostov-on-Don, 2011, pp. 176-180.

5. Dolgikh T.F. [The calculation of the bifurcation curves for stationary problem of Rayleigh-Benard-Karman convection]. So-vremennye problemy MSS [Modern problems of MCC] : materials of XVII Intern. Conf. Vol. 1. Rostov-on-Don, 2014, pp. 166-170.

6. Zhukov M.Yu., Shiryaeva E.V. Ispol'zovanie paketa ko-nechnykh elementov FreeFem++ dlya zadach gidrodinamiki, elektroforeza i biologii [Using the finite element package Free-Fem ++ for hydrodynamics problems, electrophoresis and biology]. Rostov-on-Don, 2008, 256 p.

7. Vladimirov V.A., Denissenko P.V., Zhukov M.Yu., Yu-dovich V.I. Techenie vrashchayushcheisya zhidkosti v zazore mezhdu dvumya parallel'nymi plastinami [Vortex flows in the gap between two parallel plates]. Mat. modelirovanie, 2001, vol. 13, no 2, pp. 27-38.

Поступила в редакцию_25 сентября 2015 г.