CC BY
16
ISSN 1992-6502 (P ri nt)_
2016. Т. 20, № 2 (72). С. 74-80
Ъыьмт QjrAQnQj
ISSN 2225-2789 (Online) http://journal.ugatu.ac.ru
УДК 532.5.013:536.24:519.6
Численное моделирование конвективных течений
аномально термовязкой жидкости
в. с. Кулешов 1, к. в. Моисеев2
1kuleshovvs@gmail.com, 2constgo@mail.ru Институт механики им. Р. Р. Мавлютова УНЦ РАН, Уфа Поступила в редакцию 04.04.2016
Аннотация. Представлены результаты численного исследования тепловой конвекции аномально термовязкой жидкости в квадратной полости. Рассмотрена зависимость в виде гауссовой кривой (колоко-лообразной функции). Математическая модель основана на уравнениях динамики сплошной среды и записана в приближение Обербека-Буссинеска. Вычислительный код реализован по методу контрольного объема и алгоритма SIMPLE с применением многопроцессорных технологий. Изучено влияние параметров аномалии вязкости на режимы конвективных течений, вычислены интегральные коэффициенты теплоотдачи в плоской ячейке, подогреваемой снизу, с вертикальными теплоизолированными границами. Выявлены зоны параметров задачи, при которых аномалия вязкости не влияет на картину течений с сохранением интегральных коэффициентов теплоотдачи.
Ключевые слова: тепломассоперенос, свободная конвекция, аномально вязкая жидкость, термовязкая жидкость, вязкий барьер.
Во многих процессах, происходящих на Земле, в том числе обусловленных деятельностью человека, встречается естественная (свободная) конвекция, способная оказать значительное влияние на их развитие. К таким процессам относятся океанические и атмосферные течения или процессы, протекающие при добыче высоковязкой нефти тепловыми методами, или в аппаратах химической промышленности.
В связи с истощением запасов традиционных месторождений углеводородного сырья доля трудноизвлекаемых запасов (аномально вязкой, высокосернистой нефти и др.) в последние годы постоянно увеличивается, поэтому моделирование технологических режимов добычи, транспортировки, хранения и переработки сред с аномалией вязкости по температуре и сложной реологией самой жидкости представляется актуальной задачей.
Особенности течений аномально термовязких жидкостей (далее аномально термовязкой будем называть жидкость с немонотонной зависимостью вязкости от температуры) впервые
Работа поддержана грантом РФФИ № 14-08-97060-р_поволжье_а и Программой фундаментальных исследований ОЭММПУ РАН № 13 «Вихри и волны в сложных средах».
изучены в плоском [1] и в цилиндрическом каналах [2, 3]. В данных работах представлены результаты моделирования течения жидкости с немонотонной температурной зависимостью вязкости и влияние перераспределения поля температур на картину течения и расходные характеристики.
В работе [4] рассматривалась задача о тер-мо-гравитационной конвекции жидкости в плоской ячейке с немонотонной зависимостью вязкости от температуры и влияние квадратичной зависимости вязкости от температуры с параметром отношения вязкостей на глобальную картину течения и теплообмен на изотермических стенках. Обнаружено существование периодических и квазипериодических течений.
В настоящей работе численно исследуется тепловая конвекция аномально термовязкой жидкости в квадратной полости с вертикальными адиабатическими и горизонтальными изотермическими стенками. В качестве модельной рассмотрена жидкость с немонотонной зависимостью вязкости от температуры, описываемая гауссовой кривой, определяющая общий характер изменения вязкости жидкой серы и высоковязкой нефти с учетом двух параметров аномалии: область заполнения и отношение вязкостей [1].
ПОСТАНОВКА ЗАДАЧИ И ОСНОВНЫЕ УРАВНЕНИЯ
В качестве моделируемой системы рассматривается двумерное конвективное течение несжимаемой ньютоновской жидкости с немонотонной зависимостью вязкости от температуры (аномальная жидкость) в замкнутой полости Q = [о, Z]x[0, Z]e R2 (рис. 1). Будем считать, что две вертикальные стенки адиабатические, на одной из горизонтальных стенок поддерживается постоянная температура Th (hot - горячая стенка), а на другой Tc (cold - холодная стенка), при этом Th > Tc. На всех стенках области для продольной и поперечной составляющих скорости жидкости задаются условия прилипания.
Система дифференциальных уравнений тепловой конвекции, записанная в безразмерном виде в приближении Обербека-Буссинеска [5], состоит из уравнения неразрывности, уравнений импульсов и уравнения теплопереноса:
^ = 0,
дхг
d(v • )
(1)
/
+
Pr v Ra
dv;
~dt
-
2
dx,.
dp = —— +
dx.
dx,,
dv, dv, —L +
Kdxk dx.
k
«0+«^ = (pr. Ra)
dt dx
г J, - (
2 •
■ + 8г 20,
d 20
V
dx, dx.
i J
ц(0) = 1 + Ae
(2)
(3)
(4)
где V/ - /-ая составляющая вектора скорости жидкости; X/ - пространственная координата; ^ -время; р - отклонение от гидростатического давления; 9 - безразмерная температура, опре-
Т — Т 1
деляемая как: 0 = -—Тт-, То =1 (Тя + Тс);
Тн — Тс 2
ц - динамическая вязкость жидкости; Рг = — и
а
„ gP(Tн — Тс )Ьг
Ra = ——-—--безразмерные числа
аУо
Прандтля и Рэлея соответственно [5, 6]; У0 = ц0/р0 - кинематическая вязкость жидкости при температуре То; цо, ро - значения динамической вязкости и плотности жидкости при температуре То; g - ускорение свободного падения; р - коэффициент теплового расширения; Ь - линейный размер области; А и В - параметры аномалии вязкости жидкости. Параметр А = Цтах / Цшт - 1 показывает отношение макси-
Рис. 1. Схема расчетной области
мального и минимального значений вязкости в заданном диапазоне изменения температуры [Тс:Тн], а параметр В характеризует степень наполненности заданного температурного интервала.
Переход от размерных к безразмерным величинам в системе дифференциальных уравнений (1-4) осуществляется при помощи следующих характерных параметров: линейный размер
области Ь; скорость V* = )12; время
= Ь / vо; давление р* = р0V2.
Следует отметить, что физические переменные и теплофизические константы жидкости отсчитываются от заданной температуры То. Также полагается, что теплофизические переменные, за исключением вязкости, не зависят от температуры.
Рассматриваемая система (1-4) в начальный момент времени находится в невозмущенном состоянии, то есть жидкость в полости неподвижна vi | = 0, определяется при средней
температуре То: 0| = 0о = 0 . На границах расчетной области задаются следующие параметры: скорость на стенках равна нулю, то есть V. = 0, на адиабатических границах зада-
I ХГЧЛУПииТ!
ются условия Неймана —
dx
= 0, а на изотер-
x=0,1
мических поддерживается постоянная темпера-
1 1
тура 0 =0~= - и 0| = 0Г =— .
\y=0 н 2 у=0 C 2
Полученная модель решается численно с использованием метода контрольного объема и алгоритма SIMPLE [7], модифицированного до второго порядка аппроксимации нестационарных членов [8]:
d
2
3ф"+1 - 4ф" + ф"-1 dt 2Ы '
где индекс n + 1 обозначает новые значения физической переменной ф = (u, v, 0), n - в текущий и n - 1 - в предыдущий момент времени.
Система дифференциальных уравнений (1-3), записанная в дискретном виде, сводится к системе линейных алгебраических уравнений, которая в матричном представлении имеет пятидиагональный вид, поэтому для эффективного ее решения применяется свободная высокопроизводительная библиотека решателей Hypre (High performance preconditio -ners) [9]. При реализации вычислительного кода использовался метод AMG (algebraic multigrid method) [10].
Разработанный вычислительный код был протестирован на таких задачах, как термогравитационная конвекция в квадратной полости несжимаемой жидкости с постоянной вязкостью, подогреваемой сбоку [11] и снизу [12].
В табл. 1-3 приведены результаты тестовых вычислительных экспериментов при фиксированном числе Прандтля Pr = 0,71 и для чисел Рэлея Ra = 104; 105; 106. В этих таблицах показаны характеристики течений, полученные на разных пространственных сетках и при шаге по времени At = 10-4, в последнем столбце приводятся для сравнения результаты работы [12], полученные также с помощью метода контрольного объема на сетке 256x256 узлов.
Результаты тестовых расчетов показали, что разработанный вычислительный продукт позволяет получить достоверные решения для поставленной задачи.
Таблица 1 Результаты вычислений для Ra = 104
Характе- Результаты вычислений Работа[12]
ристика 50x50 100x100 200x200 256x256
NuH 2,1577 2,1580 2,1580 2,1581
NuC 2,1577 2,1580 2,1580 2,1580
umax 0,25110 0,25200 0,2520 0,25228
vmax 0,26235 0,26338 0,26360 0,26369
Таблица 2 Результаты вычислений для Ra = 105
Характе- Результаты вычислений Работа [12]
ристика 50x50 100x100 200x200 256x256
NUh 3,9027 3,9078 3,9098 3,9103
Nuc 3,9027 3,9078 3,9098 3,9103
umax 0,33225 0,34086 0,34343 0,34434
vmax 0,36514 0,37287 0,37497 0,37569
Таблица 3 Результаты вычислений для Ra = 106
Характе Результаты вычислений Работа[12]
ристика 50x50 100x100 200x200 256x256
NUh 6,5208 6,3951 6,3373 6,3092
Nuc 6,5208 6,3951 6,3375 6,3092
umax 0,32152 0,35016 0,36425 0,37088
vmax 0,37172 0,39274 0,401978 0,40600
РЕЗУЛЬТАТЫ МОДЕЛИРОВАНИЯ
В рамках представленного исследования авторы поставили перед собой вопрос, как влияют параметры колоколообразной зависимости вязкости от температуры (4) на режимы течения жидкости в квадратной полости. Для получения ответа был исследован диапазон изменения параметров задачи (1-4): безразмерные числа Прандтля полагались равными Рг = 0,1; 1 ;10, числа Рэлея - Ra = 103; 104; 105; 106, параметр аномалии А изменялся в диапазоне 1 < А < 103 второй параметр аномалии В - в диапазоне 10 2 < В < 104. Расчеты проводились на суперкомпьютере Уфимского государственного авиационного технического университета [13] на сетке размером 100*100 узлов и шагом интегрирования по времени А( = 10-4.
В результате вычислительных экспериментов установлены критические значения параметров аномалии (В > 100), при которых отчетливо наблюдается влияние так называемого «вязкого барьера» (зоны максимальной вязкости) на картину течений в исследуемой области. Так, при В = 1000 с ростом величины А происходит формирование третичных вихрей, расположенных примерно на высоте у ~ 0,2 и у ~ 0,8 (рис. 2). При этом наблюдается увеличение этих вихрей и уменьшение вторичных вихрей (расположенных в левом верхнем и правом нижнем углах), об этом свидетельствует поднятие фронтальной поверхности температуры, иначе говоря, линии средней температуры при 00 = 0.
На рис. 3 показана зависимость высоты поднятия линии фронта на левой адиабатической стенке полости от величины А, где сплошная линия соответствует аппроксимации значений высоты, полученных в результате численного эксперимента.
Для изучения теплопереноса вычислим средний коэффициент теплообмена по всей заданной области число Нуссельта Ки, которое представляет собой отношение полного потока тепла к потоку тепла за счет теплопроводности жидкости, и определяемое как
Рис. 2. Распределения температуры (верхняя часть) и линии тока (нижняя часть) в квадратной ячейке:
а - постоянная вязкость; б - А = 10; В = 103; Рг = 0,1; Ка = 106; в - А = 300; В = 103; Рг = 0,1; Ка = 106
_ 1 1
Ми = Ц
о о
(Рг- Ra )2 V -0 —
Ш су
dxdy .
Следует отметить, что при изменении параметра А от 1 до 300 образовавшиеся дополнительные третичные вихри не приводят к увеличению или уменьшению интенсивности теплоотдачи (максимальное изменение - 2,5%) в сравнении с постоянной вязкостью (число Нуссельта Ки = 5,82).
На рис. 4 представлено распределение продольной скорости в центральном вертикальном сечении х = о,5 и центральном горизонтальном сечении у = о,5. Видно, что с ростом параметра А от 1 до 10 для поперечной и продольной компонент скорости наблюдается увеличение максимального значения скорости, однако при дальнейшем увеличении числа А отмечается отклонение горизонтальной компоненты скорости. Для вертикальной компоненты характерно распределение с резким скачком значения скорости вблизи адиабатических стенок с линейным перепадом в центре ячейке.
С дальнейшим увеличением параметра аномалии А система теряет устойчивость, и стационарный режим сменяется периодическим, о чем свидетельствует изменение среднего числа Нуссельта
Ки в момент времени от 900 до 1000 (рис. 5). Следует отметить, что несмотря на то, что разница достигает трехкратного значения, среднее значение числа Нуссельта, осредненное на заданном интервале, составляет 5,02, что на 13,7% процентов меньше в сравнении с постоянной вязкостью.
Рис. 3. Высота фронтальной поверхности температуры
Рис. 4. Эпюры продольной составляющей вектора скорости (а) и поперечной (б)
Для данного режима характерно динамическое изменение третичных вихрей, что приводит к перестройке поля температуры в центральной области и к формированию обширной зоны вязкого барьера, как это представлено на рис. 6. В первоначальный момент времени (рис. 6, а) наблюдается разрыв зоны «вязкого барьера» на две подобласти. С течением времени происходит слияние «вязкого барьера» в одну область (рис. 6, б). Далее эта зона вращается в глобальном течении (рис. 6, в), после чего происходит образование разрывной зоны «вязкого барьера» и процесс повторяется периодически.
Рис. 5. Изменение среднего числа Нуссельта Nu со временем
ЗАКЛЮЧЕНИЕ
В данной работе проведено численное исследование естественной конвекции аномально термовязкой жидкости в плоской ячейке с вертикальными адиабатическими и горизонтальными изотермическими стенками с колоколооб-разной зависимостью вязкости от температуры. Результаты получены с помощью метода контрольного объема и алгоритма SIMPLE.
Изучено влияние параметров аномалии вязкости на режимы конвективных течений и интегральные коэффициенты теплоотдачи в плоской ячейке, подогреваемой снизу, с вертикальными теплоизолированными границами.
Выявлены зоны параметров задачи, при которых гауссовская или колоколообразная зависимость вязкости от температуры не влияет на картину течений с сохранением интегральных коэффициентов теплоотдачи в отличие от жидкости с постоянной вязкостью.
Также в ряде случаев учет зависимости вязкости от температуры приводит к существенному изменению характера процессов тепломассообмена: обнаружены значения параметров аномалии A, при которых происходит перестройка течения и формирование вторичных и третичных вихрей, которые, в свою очередь, приводят к перераспределению поля температуры и смене стационарного режима на периодический с соответствующими колебаниями тепловых потоков и вихревых структур.
Полученные результаты окажутся полезными при проектировании аппаратов в нефтехи-
Рис. 6. Распределения температуры (верхняя часть), вязкости (средняя часть) и линии тока (нижняя часть) в квадратной ячейке при А = 1о3, В = 1о3, Рг = о.1, Яа = 1о6 в моменты времени:
а - г = 992; б - г = 997; в - г = 998
мической промышленности и разработке новых технологий для добычи, транспортировки, хранении и переработки аномально вязкой нефти со сложной реологией.
Результаты работы обсуждались на Российской научно-технической конференции «Мавлю-товские чтения», посвященной 90-летию со дня рождения член-корреспондента РАН, доктора технических наук, профессора Р. Р. Мавлютова и были рекомендованы к публикации в Вестнике УГАТУ.
СПИСОК ЛИТЕРАТУРЫ
1. Урманчеев С. Ф., Киреев В. Н. Установившееся течение жидкости с температурной аномалией вязкости // Доклады академии наук. 2004. Т. 396, № 2. С. 204-207. [ Urmancheev S.F., Kireev V.N. Steady flow of a fluid with an anomalous temperature dependence of viscosity // Doklady akademii nauk. 2004. T. 396, № 2. S. 204-207. ]
2. Хизбуллина С. Ф. Численное исследование течения жидкости с немонотонной зависимостью вязкости от температуры // Вестник Башкирского университета. 2006. № 2. С. 22-25. [ Khizbullina S.F. Numerical research of liquid flow with non-monotonous dependence of viscosity on
temperature // Vestnik Bashkirskogo universiteta. 2006. № 2. S. 22-25. ]
3. Хизбуллина С. Ф., Урманчеев С. Ф., Киреев В. Н.
Математическое моделирование течения аномально термовязкой жидкости в цилиндрическом канале // Труды Четвертой Российской национальной конференции по теплообмену: в 8 т. Т. 2. Вынужденная конвекция однофазной жидкости. М.: Издательский дом МЭИ. 2006. С. 145-148. [ Khizbullina S.F., Urmancheev S.F., Kireev V.N. Mathematical modeling of flow with anomalous thermoviscous fluid in cylindrical channel // Trudy Chetvertoi Rossijskoi natsional'noi konferentsii po teploobmenu: v 8 tomakh. T. 2. Vynuzhdennaya konvektsiya odnofaznoi zhidkosti. - M.: Izdatel'skii dom MEI. 2006. S. 145-148. ]
4. Ильясов А. М., Моисеев К. В. Урманчеев С. Ф. Численное моделирование термоконвекции жидкости с квадратичной зависимостью вязкости от температуры // Сиб. журн. индустр. матем. 2005. Т. 8, № 4. С. 51-59. [ Il'yasov A.M., Moiseev K.V. Urmancheev S.F. Numerical simulation of fluid's thermal convection with square dependence of viscosity on temperature // Sib. zhurn. industr. matem. 2005. T. 8, № 4. S. 51-59. ]
5. Гершуни Г. З., Жуховицкий Е. М. Конвективная устойчивость несжимаемой жидкости. М.: Наука, 1972. 392 с. [ Gershuni G.Z., Zhukhovitskii E.M. Convective stability of incompressible fluids. M.: Nauka, 1972. 392 s. ]
6. Гетлинг А. В. Конвекция Рэлея-Бенара. Структуры и динамика. М.: Эдиториал УРСС. 1998. 248 с. [ Getling A.V. Rayleigh-Benard convection. Structure and Dynamics. M.: Editorial URSS. 1998. 248 s. ]
7. Патанкар С. Численные методы решения задач теплообмена и динамики жидкости. М.: Энергофизмат. 1984. 152 с. [ Patankar S. Numerical heat transfer and fluid flow. Hemisphere Publishing Corporation, New York, 1980. ]
8. Андерсон Д., Таннехилл Дж., Плетчер Р. Вычислительная гидромеханика и теплообмен: в 2-х т. Т. 1. / пер. с англ. М.: Мир. 1990. 384 с. [ D. Anderson, Tannehill J., R. Pletcher Computational Fluid Mechanics and Heat Transfer. Hemisphere Publishing Corporation, New York, 1984. ]
9. Библиотека решателей (High performance preconditioners) [Электронный ресурс]. Режим доступа: http://computation.llnl.gov/project/linear_solvers/ (дата обращения: 01.02.2016). [ Biblioteka reshatelei (High performance preconditioners) [Elektronnyi resurs]. Rezhim dostupa: http://computation.llnl.gov/project/linear_solvers/. ]
10. Falgout Robert D. An introduction to algebraic Multigrid // Computing in Science and Engg. 2006. V. 8. № 6. P. 24-33.
11. Le Quere P., Alziary De Roquefort T. Computation of natural convection in two-dimensional cavities with Chebyshev polynomials // J. Comput. Phys. 1985. V. 57. № 2. P. 210-228.
12. N. Ouertatani et al., Numerical simulation of two-dimensional Rayleigh - Benard convection in an enclosure // C. R. Mecanique. 2008. V. 336, № 5. P. 464-470.
13. Вычислительный кластер УГАТУ [Электронный ресурс]. Режим доступа: http://www.ufa-rb.ru/supercomputer/ (дата обращения: 05.02.2016). [ Vychislitel'nyi klaster UGATU [Elektronnyi resurs]. Rezhim dostupa: http://www.ufa-rb.ru/supercomputer/ ]
ОБ АВТОРАХ
КУЛЕШОВ Василий Сергеевич, асп. лаб. механики многофазных систем. Дипл. матем., системный программист (УГАТУ, 2014). Готовит дисс. о конвективных течениях аномально термовязкой жидкости. Иссл. в обл. механики неоднородных и реологических сложных сред в полях внешних сил.
МОИСЕЕВ Константин Валерьевич, с.н.с. лаб. механики многофазных систем, доц. каф. механики и конструирования машин УГНТУ. Дипл. математик (Баш. гос. ун-т, 2004). К.ф.-м.н. по механике жидкости газа и плазмы (Тюм. гос. ун-т, 2009). Иссл. в обл. механики неоднородных и реологических сложных сред в полях внешних сил.
METADATA
Title: Numerical simulation of convection anomalous
thermoviscous flow. Authors: V.S. Kuleshov, C.V. Moiseev
Affiliation: Mavlutov's Institute of Mechanics of Ufa Branch, RAS
Email: kuleshovvs@gmail.com, constgo@mail.ru Language: Russian.
Source: Vestnik UGATU (scientific journal of Ufa State Aviation Technical University), vol. 20, no. 2 (72), pp. 74-80, 2016. ISSN 2225-2789 (Online), ISSN 1992-6502 (Print). Abstract: The simulation results for thermal convection of anomalous thermoviscous liquid in square cavity are presented. The liquid is called anomalous thermoviscous if the viscosity depends on temperature nonmonotonically. The gauss curve dependence (bell function) is considered for convenience of numerical study. The mathematical model is based on the classical fluid mechanics equations written in Oberbeck - Boussinesq approximation. The SIMPLE algorithm developed using multiprocessor technologies is used for calculations. The influence of viscosity anomalous on convective flow regime is studied. The integral heat transfer coefficients in a heated from below 2D-cavity with thermo insulated vertical walls. The range of parameters, for which at constant integral heat transfer coefficients the viscous abnormality does not affect the fluid flow are identified. Key words: heat transfer, mass transfer, natural convection, anomalous viscosity of fluid, thermoviscous fluid, viscous barrier.
About authors:
KULESHOV, Vasily Sergeyevich, Postgrad. (PhD) Student Dept. "Mechanics of multiphase systems" Mavlutov's Institute of Mechanics of Ufa Branch, RAS, Dipl. Mathematician and system programmer (UGATU, 2014). MOISEEV, Konstantin Valer'evich, senior researcher Dept. "Mechanics of multiphase systems" Mavlutov's Institute of Mechanics of Ufa Branch, RAS, Cand. Of Phys. And math Sci. (Tum. State Univ., 2009), Dipl. Mathematician (Bash. SU, 2004).
