CC BY
152Математические заметки СВФУ Январь—март, 2017. Том 24, № 1
УДК 519.63
ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ДВУМЕРНОЙ
КОНВЕКЦИИ РЭЛЕЯ - БЕНАРА
В. В. Григорьев, П. Е. Захаров
Аннотация. Рассматривается конвекция Рэлея — Бенара (естественная конвекция). Это течение, которое образуется в вязкой среде при нагреве снизу и охлаждении сверху. В результате происходит образование вихрей (конвективных ячеек). Данный процесс описывается системой нелинейных дифференциальных уравнений в приближении Обербека — Буссинеска. В качестве управляющих параметров, определяющих режимы конвекции, выбраны число Рэлея и число Прандтля. Система решена методом конечных элементов с помощью вычислительного пакета РЕшС8. Получены численные результаты при различных числах Рэлея. Исследована интегральная характеристика (число Нуссельта) в зависимости от числа Рэлея.
Ключевые слова: естественная конвекция, приближение Обербека — Буссинеска, метод конечных элементов, конвекционные ячейки.
1. Введение
Конвекция, связанная с неоднородным нагревом, является одним из самых распространенных видов течений газа и жидкости в природе. Немалую роль она играет и в различных технических устройствах.
Конвекция в плоском горизонтальном слое жидкости, подогреваемом снизу и охлаждаемом сверху, или конвекция Рэлея — Бенара, является типом конвекции, который рассматривается чаще всего [1—3].
Конвекция Рэлея — Бенара играет важную роль в широком диапазоне явлений в геофизике, астрофизике, метеорологии, океанографии и инженерии, иными словами, представляет собой атмосферный феномен, который обычно наблюдается в широких пространственных и временных масштабах [4, 5]. Математическая модель этой конвекции в настоящее время успешно применяется, к примеру, в изучении мантии Земли (образование вулканов), при моделировании поверхности Солнца, в некоторых космических и атмосферных явлениях [6].
В данной работе проводится численное моделирование нелинейной модели двумерной конвекции Рэлея — Бенара в приближении Обербека — Буссинеска методом конечных элементов с помощью вычислительного пакета ЕЕшС8 [7, 8].
Работа выполнена при поддержке Мега-гранта Правительства Российской Федерации (№ 14.Y26.31.0013).
© 2017 Григорьев В. В., Захаров П. Е.
2. Математическая модель
Рассмотрим ограниченную область
П = {х е Б2 : х = (хьх2), 0 < х1 < Ь, 0 < х2 < Н},
заполненную вязкой несжимаемой жидкостью, где Ь — ширина, а Н — высота емкости (рис. 1).
Рис. 1. Схема области
Движение жидкости описывается уравнениями Навье — Стокса в приближении Обербека — Буссинеска, где плотность р предполагается линейно зависящей от температуры Т [9-11]:
1 (ди \
— ( + и Уи ) - Аи + Чр = ЯаТе2, (1)
Рг у дЬ )
V- и = 0, (2)
дТ
— + и • УТ = АТ, (3)
где и — поле скорости и р — давление. На границах Г1 и Г2 скорость равна нулю:
и = 0, хеГ1 иГ2, о<кге.
Температура равняется 1 на границе Г1 и 0 на границе Г2. Начальное распределение температуры равно 0.5,
Т = 1, х е Г1, 0 <1<1е,
Т = 0, х е Г2, 0 <1<1е, т = 0.5, хе п, 0<г<1е-
По боковым границам Г3 и Г4 могут стоять условия непротекания или периодические граничные условия.
В безразмерном уравнении остаются только два управляющих параметра определяющие режимы конвекции: Рг — число Прандтля и Иа — число Рэлея. По определению
Рг = 1
Ж
где V — кинематическая вязкость жидкости и Ж — коэффициент температуропроводности, и
д/ЗАТН3 Ка =--—,
где д — ускорение свободного падения, в — коэффициент теплового расширения жидкости, ДТ — разность температур между верхней и нижней границами.
3. Аппроксимация
Для аппроксимации по времени используем метод конечных разностей. Определим равномерную сетку по времени
шт = {Г : Г = тп, п = 0,1,..., Ж, тЖ = ¿е>.
Для краткости используем обозначение Тп = Т (£", х). Для решения подобных задач можно использовать неявную схему или схемы расщепления. Для решения этой задачи используется схема расщепления: сначала вычисляется скорость и давление и затем с помощью найденной скорости вычисляется температура:
1 Л.П+1 _ \
— -+ ип+1 • Уип+1 - Аи™+1 + Урп+1 - Иа Тпе2 = 0, (4)
Рг у т у
V • ип+1 = 0. (5)
ТП+1 _ ТП
-+ ип+1 • УТп+1 - АТп+1 = 0. (6)
т
Проведем аппроксимацию по пространству задачи (1)—(3) с использованием метода конечных элементов. Запишем эти задачи в вариационной формулировке. Для пробных и тестовых функций для уравнения скорости введем соболев-ское пространство W(Q), состоящее из векторных функций и таких, что и2, и) и | gгadи|2 имеют конечный интеграл в О:
W = {и е н1(О): и = 0, х е г1 и г2}, W = {V е н1(О): V = 0, х е г1 и г2},
теперь определим пространства функций давления:
д = н1(о), (3 = н1(о),
аналогично определим пространства функций температуры:
V = {Т е н1(О): Т = 0, х е гь Т = 1, х е г2},
V = {г е н1(О) : г = 0, х е г1 и г2}.
Умножим уравнения на тестовые функции и проинтегрируем по области с использованием формулы Грина. Получаем задачу: найти функции ип+1 € W, рп+1 € Q, удовлетворяющие уравнению
1 Г иП+1 _ ип 1 г
— /--+ — / (ип+1 •gradu'г+1)Wx
Рг у т Рг ]
о о
+ 1 gгad ип+1 • grad V ¿х — J рп+1 div V ¿х о о
+ У ип+1д ¿х — И^у Т^ • е2 ¿х = 0 Vv € W, Уд € (7) оо
Аналогично получаем задачу для температуры: найти функцию Тп+1 € V удовлетворяющую уравнению
/Т п+1_т п г
-гсЬс+ / gradTrl+1 • gradг¿х
о о
+ У ип+1 • gгad Тп+1г ¿х = 0 Уг € V. (8) о
Для области О строим сетку О^ и для этой области определим конечномерные пространства W^ с W, W^ с W, С (ь С Vh С V и ^ с V. Для устойчивости решения выберем элементы Тау1ог-Нос^ [8,12], где в качестве базисов используются полиномы Лагранжа второй и первой степеней соответственно. Для пространства Vh используем элементы Лагранжа первой степени.
4. Численное исследование
Задача была решена с помощью вычислительного пакета ЕЕшС8 [13]. Были проведены несколько расчетов с разными значениями управляющих параметров Рг и Иа. Значения брались в рамках стационарного режима конвекции. Нелинейная часть в скорости решалась итерационным методом Ньютона. Также для проверки реализации математической модели была взята статья [14], в которой авторы провели натурный эксперимент по визуализации конвективных ячеек. В эксперименте была взята дистиллированная вода, а число Рэлея авторы брали в диапазоне от 1500 до 4000. Аспектное отношение равно 16. Мы преследовали цели проверки соответствия количества конвективных ячеек при одном и том же значении числа Рэлея. Из-за того, что результаты были взяты после достижения предельного стационарного режима, число Прандтля было взята равным единице. Число Рэлея для сравнения было выбрано равным 2500.
Рис. 2. Результаты натурного эксперимента
Рис. 3. Линии тока численного эксперимента
На рис. 2 и 3 отображаются результаты натурного и численного экспериментов. Количество конвективных ячеек совпадает. Результаты сравнения удовлетворительные, поэтому приступаем к основной задаче. Для решения задачи в данной статье Ь = 3, Н = 1. Управляющий параметр Рг во всех расчетах равен 0.7, а параметр Иа изменяли в диапазоне от 106 до 108. Расчеты проводились на сетке с 12508 узлами и с шагом по времени 0.0001.
Как видно из рис. 4-6, шлейфы образуются в разные моменты времени. Приходим к выводу, что с ростом числа Рэлея процессы идут быстрее и увеличивается количество шлейфов в начале процесса. Конвективные ячейки успешно образовались при всех рассматриваемых числах Рэлея и все достигли предельных стационарных режимов. Для каждого из случаев первые возникшие шлейфы соединялись в несколько больших шлейфов. Путем такого соединения процессы достигали своих стационарных режимов, образуя конечные конвективные ячейки. Заметим, что в случае Иа = 107 присутствуют только две конвективные ячейки вместо четырех в остальных случаях. Объясняется это тем, что аспектное отношение равно трем, поэтому количество конвективных ячеек будет равно или двум, или четырем.
В подобного рода задачах представляет интерес оценка эффективности потока. Наиболее естественной и принятой мерой, которая квалифицирует усиление вертикального потока тепла за счет конвекции, является число Нус-сельта Ми.
Интегрируя по ж2 е (0, Н) усредненный поток тепла по горизонтальному пространству иТ — VT в вертикальном направлении, получим математическое определение числа Нуссельта [9]:
н
о
6 ТО 561
Рис. 4. Температура и линии тока при Яа = 106, а) Ь = 0.0209, Ъ) Ь = 0.0230, о) Ь = 0.0240, а) Ь = 0.0250, е) Ь = 0.0270, £) Ь = 0.03, g) Ь = 0.14
Рис. 5. Температура и линии тока при Яа = 107, а) 4 = 0.0042, Ь) 4 = 0.0055, с) 4 = 0.0065, а) 4 = 0.0075, е) 4 = 0.0095, £) 4 = 0.0125, g) 4 = 0.033
Рис. 6. Температура и линии тока при Яа = 108, а) Ь = 0.0008, Ъ) Ь = 0.0012, о) Ь = 0.0018, а) Ь = 0.0022, е) Ь = 0.003, £) Ь = 0.0035, g) Ь = 0.0125
Рис. 7. Влияние числа Рэлея на установление среднего числа Нуссельта во времени
На рис. 7 сравниваются результаты расчетов по установлению во времени среднего числа Нуссельта для различных чисел Рэлея. Как видно из рис. 7, интенсивность конвективного переноса тепла увеличивается с ростом числа Рэлея. При Иа = 108 наблюдается колебание, которое показывает, что число Рэлея приближается к критическому значению. При дальнейшем росте числа Рэлея будет расти частота колебаний. Если число Рэлея стало больше чем его критическое значение, то режим конвекции сменится со стационарного на турбулентный. Само критическое значение числа Рэлея зависит от аспектного отношения: чем оно меньше, тем меньше критическое значение числа Рэлея, и наоборот, чем больше аспектное отношение, тем больше критическое значение числа Рэлея [1].
5. Заключение
В работе методом конечных элементов численно исследована конвекция Рэлея — Бенара в безразмерных переменных.
Показано, что с увеличением числа Рэлея процесс достигает своего предельного стационарного режима быстрее, также увеличивается количество начальных шлейфов, которые, соединяясь между собой, образуют большие шлейфы, тем самым образуя конвективные ячейки.
Число Нуссельта принимает колебательный характер при числе Рэлея, равном 108, и с дальнейшим ростом будет расти частота колебаний.
ЛИТЕРАТУРА
1. Гетлинг А. В. Конвекция Рэлея — Бенара. Структуры и динамика. М.: Эдиториал УРСС, 1999.
2. Колмычков В. В., Мажорова О. С., Попов Ю. П., Федосеев Е. Э. О численном моделировании конвекции Рэлея — Бенара / Препринт. Ин-т прикладной математики им. М. В. Келдыша РАН, 2007. 3. Гетлинг А. В. Формирование пространственных структур конвекции Рэлея — Бенара. // Успехи физ. наук. 1991. Т. 161, вып. 9. С. 1-80.
4. Закинян Р. Г., Сухов С. А., Ларченко И. Н. Математическое моделирование тепловой конвекции // Современные проблемы науки и образования. 2011. Т. 6. Режим доступа: www.science-education.ru/100-5016
5. Chilla F., Schumacher J. New perspectives in turbulent rayleigh-benard convection // Europ. Phys. J.(E). 2012. V. 35, N 7. P. 1-25.
6. Spiegelman M. Myths and methods in modeling // Columbia Univ. Course Lecture Notes. Available online at http://www.ldeo.columbia.edu/- mspieg/mmm/course.pdf, accessed. 6:2006, 2004.
7. Самарский А. А., Вабищевич П. Н. Вычислительная теплопередача. М.: Книжный дом ЛИБРОКОМ, 2009.
8. Logg A., Mardal K.-A., Wells G. Automated Solution of Differential Equations by the Finite Element Method. Berlin: Springer Science & Business Media, 2012. (The FEniCS book; V. 84).
9. Otto F., Seis C. Rayleigh-Benard convection: improved bounds on the nusselt number //J. Math. Phys. 2011. V. 52, V. 8. P. 083702.
10. Sakievich P. J., Peet Y. T., Adrian R. J. Large-scale thermal motions of turbulent rayleigh-benard convection in a wide aspect-ratio cylindrical domain // Intern. J. Heat Fluid Flow. 2016. V. 61. P. 183-196.
11. Vynnycky M., Masuda Y. Rayleigh-Benard convection at high rayleigh number and infinite prandtl number: Asymptotics and numerics // Physics of Fluids. 2013. V. 25, N 11. P. 113602.
12. Вычислительные технологии. Профессиональный уровень / М. Ю. Антонов, Н. М. Афанасьева, В. С. Борисов и др. 2014.
13. Alnjs M., Blechta J., Hake J., Johansson A., Kehlet B., Logg A., Richardson C., Ring J., Rognes M. E., and Wells G. N. The fenics project version 1.5 // Arch. Numer. Software. 2015. V. 3. P. 9-23.
14. Lir J. T. and Lin T. F. Visualization of roll patterns in Rayleigh-Benard convection of air in a rectangular shallow cavity // Intern. J. Heat Mass Transfer. 2001. V. 44, N 15. P. 2889-2902.
Статья поступила 28 ноября 2016 г.
Григорьев Василий Васильевич, Захаров Петр Егорович Северо-Восточный федеральный университет им. М. К. Аммосова, Институт математики и информатики, ул. Кулаковского 48, Якутск 677891 d_alighieri@rambler.ru, zapetch@gmail.com
Математические заметки СВФУ Январь—март, 2017. Том 24, № 1
UDC 519.63
NUMERICAL MODELING OF THE TWO-DIMENSIONAL RAYLEIGH-BENARD CONVECTION V. V. Grigoriev and P. E. Zakharov
Abstract: We study Rayleigh-Benard convection which is a type of natural convection occurring in a plane horizontal layer of viscous fluid heated from below and cooled from above, in which the fluid develops a regular pattern of convection cells known as Benard cells. The process of rotation is described by a system of nonlinear differential Oberbeck-Boussinesq equations. As convection parameters, the Rayleigh number and the Prandtl number are taken. The system is solved using the finite element method by FEniCS. We obtain numerical results for varying Rayleigh numbers and study the dependence of the Nusselt number on the Rayleigh number.
Keywords: natural convection, Oberbeck-Boussinesq approximation, finite element method, convection cells.
REFERENCES
1. Getling A. V. Rayleigh-Benard Convection: Structures and Dynamics, World Sci., Singapore (1998). (Adv. Ser. Nonlinear Dyn.; V. 11).
2. Kolmychkov V. V., Mazhorova O. S., Popov Yu. P., Fedoseev E. E. "On numerical modeling of the Rayleigh-Benard convection," preprintKeldysh Inst. Appl. Math. (2007).
3. Getling A. V. "Formation of spatial structures in Rayleigh-Benard convection," Sov. Phys. Usp., 34, No. 9, 737-776 (1991).
4. Zakinyan R. G., Sukhov S. A., and Larchenko I. N. "Mathematical modeling of the heat convection [in Russian]," Contemporary Problems in Science and Education, 6,
5. Chilla F. and Schumacher J. "New perspectives in turbulent Rayleigh-Benard convection," Eur. Phys. J. (E), 35, No. 7, 1-25 (2012).
6. Spiegelman M. Myths and Methods in Modeling, Columbia Univ. Course Lect. Notes (available online at http://www.ldeo.columbia.edu/~mspieg/mmm/course.pdf.
7. Samarskii A. A. and Vabishchevich P. N. Numerical Heat Transfer [in Russian], Librokom, Moscow (2009).
8. Logg A., Mardal K.-A., and Wells G. Automated Solution of Differential Equations by the Finite Element Method, Springer, Berlin (2012). (The FEniCS book; V. 84).
9. Otto F. and Seis C. "Rayleigh-Benard convection: Improved bounds on the Nusselt number," J. Math. Phys., 52, No.8, 083702 (2011).
10. Sakievich P. J., Peet Y. T., and Adrian R. J. "Large-scale thermal motions of turbulent Rayleigh-Benard convection in a wide aspect-ratio cylindrical domain," Int. J. Heat Fluid Flow, 61, 183-196 (2016).
11. Vynnycky M. and Masuda Y. "Rayleigh-Benard convection at high Rayleigh number and infinite Rrandtl number: Asymptotics and numerics," Phys. Fluids, 25, No. 11, 113602 (2013).
12. Computational Technologies. Profesional Level [in Russian] (P. N. Vabishchevich, ed.), URSS, Moscow (2014).
© 2017 V. V. Grigoriev and P. E. Zakharov
13. Alnjs M., Blechta J., Hake J., Johansson A., Kehlet B., Logg A., Richardson C., Ring J., Rognes M. E., and Wells G. N. "The FENiCS project version 1.5," Arch. Numer. Software, 3, 9-23 (2015).
14. Lir J. T. and Lin T. F. "Visualization of roll patterns in Rayleigh-Benard convection of air in a rectangular shallow cavity," Int. J. Heat Mass Transfer, 44, No. 15, 2889-2902 (2001).
Submitted November 28, 2016
Vasiliy V. Grigoriev and Petr E. Zakharov M. K. Ammosov North-Eastern Federal University, Institute of Mathematics and Informatics, 48, Kulakovsky St., Yakutsk 677000, Russia d_alighieri@rambler.ru, zapetch'Sgmail .com
