CC BY
52
2015
ВЕСТНИК ТОМСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА Математика и механика
№ 6(38)
УДК 532.516.5
DOI 10.17223/19988621/38/11
М.А. Пономарева, М.П. Филина, В.А. Якутенок
ТЕЧЕНИЕ НЕНЬЮТОНОВСКОЙ ЖИДКОСТИ В КВАДРАТНОЙ КАВЕРНЕ ПРИ МАЛЫХ ЧИСЛАХ РЕЙНОЛЬДСА1
Рассмотрен вопрос о распределении кинематических и динамических характеристик течения неньютоновской жидкости в квадратной каверне. В качестве реологической модели использовался степенной закон. Численное решение получено в приближении ползущего течения непрямым методом граничных элементов. Исследования проведены в диапазоне изменения показателя нелинейности n от 0.2 до 1.2. Приведены профили компонент вектора скорости в характерных сечениях каверны. Результаты для ньютоновского случая совпадают с данными других авторов. Показано, что уменьшение величины n приводит к смещению центра, вокруг которого вращается жидкость, к верхней крышке каверны. Представлены поля распределения эффективной вязкости и интенсивности скоростей деформаций по области течения.
Ключевые слова: неньютоновская жидкость, течение в каверне, непрямой метод граничных элементов.
Задача о течении вязкой несжимаемой жидкости в квадратной каверне с верхней стенкой, движущейся с постоянной скоростью в горизонтальном направлении, является известным и широко используемым тестом для верификации и оценки эффективности численных методов, применяемых для решения уравнений Стокса. В связи с широким использованием такого варианта тестирования, имеется большое количество накопленных данных, позволяющих определить правильность новых результатов, тем самым гарантируя их достоверность. Основное количество результатов имеется для течений при больших числах Рейнольдса [1-6]. В работах [7, 8] рассмотрены случаи малых чисел Рейнольдса Re = 1 и 0. Рассмотрение в этих работах проведено для ньютоновской жидкости. В то же время актуальным представляется проведение исследований данного течения при малых числах Рейнольдса, но с учетом неньютоновского поведения жидкости.
Основные уравнения
Основными уравнениями для описания двумерного течения неньютоновской жидкости при малых числах Рейнольдса являются уравнения Стокса. Математическая постановка задачи записывается в безразмерных переменных. В качестве характерного размера выбрана длина стенки каверны, в качестве характерной скорости - скорость движения верхней крышки каверны. В приближении ползущего течения и при отсутствии массовых сил уравнения движения имеют вид
Исследование выполнено при финансовой поддержке Гранта Президента РФ (МК-3687.2014.1) и РФФИ в рамках научного проекта № 14-08-31579 мол_а.
Течение неньютоновсной жидкости в квадратной каверне при малых числах Рейнольдса 91
до-
дх,-
= 0, i, j = 1,2,
(1)
где о- = -p5j + тj - компоненты полного тензора напряжений, p - давление, 5-- символ Кронекера, т- - компоненты тензора вязких напряжений, х- - декартовы координаты.
В качестве реологической модели, описывающей неньютоновское поведение жидкости, выбран степенной закон, также известный как модель Оствальда де-Виля, характеризующий зависимость вязкости от скорости сдвига:
т- = 2Пеу, (2)
где n = Yn 1 - коэффициент эффективной вязкости, п - показатель нелинейности, 1/2
Y = (2ёуё^) - интенсивность скоростей деформаций, ё- = (dui /дх- + ди- /dxt)/2
- компоненты тензора скоростей деформаций, ut - компоненты вектора скорости.
Жидкости, поведение которых описывается данной моделью, можно подразделить на три типа: псевдопластичные (п < 1), ньютоновские (п = 1) и дилатантные
(п >1).
На рис. 1 представлена область решения. Верхняя стенка каверны движется в горизонтальном направлении с постоянной скоростью.
Систему (1) необходимо дополнить уравнением неразрывности
^ = 0 (3)
дХ;
и граничными условиями, которые состоят в задании компонент скорости на стенках каверны: на верхней подвижной стенке
и1 = 1, и2 = 0 (4)
и на остальных стенках
U = 0. (5)
Гранично-интегральная формулировка задачи и метод решения
Представим уравнение (1) в виде
до-
ij = ш дХ : ^
(6)
где oN = -pby + 2ё- - линейная часть полного тензора напряжений;
д г п дт f
Шi =----1 2(1 -п)ёу 1 =-- - нелинейная векторная функция, которую будем
дх у J Л дх-
рассматривать как плотность источников, распределенных по области течения Q .
92
М.А. Пономарева, М.П. Филина, В.А. Якутенок
Тогда, в соответствии с положениями непрямого метода граничных элементов [9], можно записать
U (x) = jGj (x,£)фj (|)dШ +j Gj (x,z)Tj (z)dQ(z), (7)
Г Q
где ф j (§) - плотность фиктивных источников, распределенных по границе области течения Г . Функции Gj являются фундаментальными решениями уравнений Стокса и определяются формулой [10]
G»<x->=- ^), (8)
где y = xi ~‘%i, r = (yy )2 . Если на границе области течения Г заданы значения скорости, то уравнения (7) позволяют получить значения неизвестных граничных сил ф j (§) (| е Г). Это возможно сделать при известной функции Т(z) (z е Q).
Так как эта функция заранее неизвестна, возникает необходимость организации итерационного процесса.
Для численного решения уравнений (7) используются постоянные элементы и постоянные ячейки. Граница области течения Г разбивается на N элементов. Функции фj (|) считаются постоянными на каждом элементе. Область течения
разбивается на N2 квадратных ячеек. Функции Т(z) считаются постоянными внутри ячейки. Тогда уравнения (7) в дискретной форме приобретают вид
N N N
u (xp ) = И Gq +^Т У AGpkm, (9)
q=1 k=1 m=1
где AGjq = j Gj (xp,dГ(§), AG?km = j Gj (xp,z)dQ(z), xp - середина
AFq AQkm
элемента p (узел).
Для вычисления 2N неизвестных ф-j используются 2N уравнений (9), соответствующие N элементам, на которых заданы у (xp). Коэффициенты получаемой системы линейных алгебраических уравнений AGjpq в случае постоянных
элементов можно вычислить аналитически. Технология вычисления изложена в
[11].
Для вычисления интегралов по области AGpqm используются стандартные
квадратурные формулы Гаусса без выделения особенностей. Особенности в этих интегралах имеют вид ln(1/ г). Следовательно, при интегрировании по области эти интегралы существуют в обычном смысле. Такой подход значительно упрощает алгоритм решения. При проведении расчетов использовалась квадратурная формула с 64 узлами. Для решения системы нелинейных алгебраических уравнений (9) относительно ф^^ применялся метод простой итерации. На первой итерации использовались значения Тkm определенные по ньютоновскому полю тече-
Течение неньютоновской жидкости в квадратной каверне при малых числах Рейнольдса 93
ния. Для решения соответствующей системы линейных алгебраических уравнений использовался метод Гаусса. Далее использовались значения ^^ , рассчитанные в соответствии с полем течения, полученным на предыдущей итерации. Функции ^Jkn в центре ячейки (k, m) вычислялись конечно-разностным способом с использованием рассчитанных (tNn ) в вершинах ячеек (узлах сетки).
Значения (tNn ) полностью определяются производными (dut / Xj ) , значения
которых находились с использованием центральных разностей во внутренних узлах и односторонних разностей в приграничных узлах в соответствии со значениями и1™, вычисленными в узлах сетки. Для улучшения сходимости итерационного процесса в случае сильной нелинейности использовался метод релаксации с динамическим подбором коэффициента релаксации. Кроме того, процесс вычисления оптимизировался путем распараллеливания алгоритма.
Анализ полученных результатов
Вышеописанный метод расчета позволил провести вычисления в диапазоне изменения параметра нелинейности n от 0.2 до 1.2. Результаты расчетов представлены на рис. 2 - 8.
Исследование аппроксимационной сходимости (рис. 2) показало, что приемлемая точность расчетов обеспечивается при разбиении границы Г на 160 и более элементов. Представленные в работе расчеты были проведены при N = 256 .
*2
0.8
0.6
0.4
0.2
0
-0.4 0 0.4 0.8 щ -0.4 -0.2 0 0.2 0.4 *i
Рис. 2. Аппроксимационная сходимость решения для n = 0.5 : а - профили составляющей скорости u1 (х1 = 0), б - профили составляющей скорости u2 (х2 = 0.5)
Изменение параметра n приводит к смещению центра основной вихревой зоны. Причем уменьшение n по сравнению с ньютоновским случаем (n = 1.0) сопровождается смещением центра к верхней крышке, а увеличение - к нижнему основанию. Соответствующие качественные и количественные результаты показаны на рис. 3 и 4.
94
М.А. Пономарева, М.П. Филина, В.А. Якутенок
-0.5 -0.3 -0.1 0.1 0.3 *1 -0.5 -0.3 -0.1 0.1 0.3 *1
Рис. 3. Изолинии функции тока Ф ( u1 = дФ/дх2; u2 = -дФ/дх1 ) с шагом АФ = 0.2 для различных значений показателя нелинейности: n = 0.5 (а), n = 0.7 (б), n = 1.0 (в), n = 1.2 (г)
Рис. 4. Зависимость координаты х2 центральной точки основного вихря от величины n
В работах, касающихся исследования данного течения, стандартно приводятся профили компонент скорости u1 и u2 при х1 = 0, х2 = 0.5
соответственно. Для рассматриваемого в настоящей работе нелинейного случая и малых чисел Рейнольдса указанные данные содержатся на рис. 5 и 6.
Течение неньютоновской жидкости в квадратной каверне при малых числах Рейнольдса 95
Рис. 5. Профили составляющей скорости u1 вдоль линии х1 = 0 для 0.2 < n < 1.2 (а) и сравнение с известными данными для ньютоновской жидкости (б)
-0.4 -0.2 0 0.2 0.4 х\ -0.4 -0.2 0
0.4 х\
Рис. 6. Профили составляющей скорости u2 вдоль линии х2 = 0.5 для 0.2 < n < 1.2 (а) и сравнение с известными данными для ньютоновской жидкости (б)
Распределение величины эффективной вязкости п по области течения соответствует распределению интенсивности скоростей деформаций у и параметру нелинейности n (псевдопластичному и дилатантному реологическому поведению). Это соответствие показано на рис. 7, 8. Зоны высоких значений эффективной вязкости наблюдаются в областях низких скоростей сдвига для псевдопластичных жидкостей (центр и нижние углы каверны) и в областях высоких скоростей сдвига для дилатантных жидкостей (верхние углы каверны). Размеры таких зон увеличиваются при уменьшении n (n < 1).
96
М.А. Пономарева, М.П. Филина, В.А. Якутенок
Рис. 7. Поле интенсивности скоростей деформаций у для различных значений величины n : n = 0.2 (a), n = 0.5 (б), n = 0.7 (в), n = 1.0 (г), n = 1.2 (d)
Течение неньютоновсной жидкости в квадратной каверне при малых числах Рейнольдса 97
Рис. 8. Поле эффективной вязкости п внутри области течения для различных значений показателя нелинейности: n = 0.2 (а), n = 0.5 (б), n = 0.7 (в), n = 1.2 (г) (область течения,
закрашенная черным цветом, соответствует значениям п>5.5 для псевдопластичной жидкости и п> 1.8 для дилатантной)
Заключение
Непрямым методом граничных элементов проведено исследование течения неньютоновской жидкости в квадратной каверне при малых числах Рейнольдса. Значение показателя нелинейности n варьировалось в диапазоне от 0.2 до 1.2. Получены профили компонент вектора скорости в характерных сечениях каверны. Для ньютоновской жидкости проведено сравнение полученных данных с результатами других исследователей. Показано, что уменьшение величины n приводит к смещению центра, вокруг которого вращается жидкость, к верхней крышке каверны. Представлены поля распределения эффективной вязкости и интенсивности скоростей деформаций по области течения. Полученные результаты могут быть использованы для тестирования численных методов решения уравнений Стокса при степенной зависимости вязкости от скорости сдвига.
98
М.А. Пономарева, М.П. Филина, В.А. Якутенок
ЛИТЕРАТУРА
1. Елизарова Т.Г., Милюкова О.М. Численное моделирование течения вязкой несжимаемой жидкости в кубической каверне // Журнал вычислительной математики и математической физики. 2003. Т. 43. № 3. С. 453-466.
2. Haque S., Lashgari I., Giannetti F. Stability of fluids with shear-dependent viscosity in the lid-driven cavity // J. Non-Newtonian Fluid Mechanics. 2012. No. 173. P. 49-61.
3. Anderson P.D., Galaktinov O.S., Peters G.W. Mixing of non-Newtonian fluids in time-periodic cavity flows // J. Non-Newtonian Fluid Mechanics. 2000. No. 93. P. 265-286.
4. Lashckarbolok M., Jabbari E. Collocated Discrete Least Squares (CDLS) meshless method for the simulation of power-law fluid flows // Scientia Iranica. 2013. No. 20(2). P. 322-328.
5. Nejat A., Sharbatdar M., Jalali A. A Newton-Krylov finite volume algorithm for power-law non-Newtonian fluid flows using pseudo-compressibility method // 20th AIAA Computational Fluid Dynamics Conference AiAA. 2011. С. 1-19.
6. Li Q., Hong N., Shi B. Simulation of power-law fluid flows in two-dimensional square cavity using multi-relaxation-time lattice boltzmann method // Commun. Comput. Phys. 2014. V. 15. No. 1. P. 256-284.
7. Liu M.B., Xie W.P., Liu G.R. Modeling incompressible flows using a finite particle method // Applied Mathematical Modeling. 2005. V. 29. P. 1252-1270.
8. Shamekhi A., Aliabadi A. Non-Newtonian lid-driven cavity flow simulation by mesh free method // ICCES. 2009. V. 11. No. 3. С. 67-72.
9. Бреббия К., ТеллесЖ., Вроубел Л. Методы граничных элементов. М.: Мир, 1987.
10. Ладыженская О.А. Математические вопросы динамики вязкой несжимаемой жидкости.
М.: Наука, 1970.
11. Ponomareva M.A., Filina M.P., Yakutenok V.A. The indirect boundary element method for the two-dimensional pressure- and gravity-driven free surface Stokes flow // WIT Transactions on Modelling and Simulation. 2014. V. 57. P. 289-304. DOI: 10.2495/BE370241.
Статья поступила 15.11.2015 г.
Ponomareva M.A., Filina M.P., Yakutenok V.A. NON-NEWTONIAN FLUID FLOW IN A LID-DRIVEN CAVITY AT LOW REYNOLDS NUMBERS
DOI 10.17223/19988621/38/11
In this paper, the question about the distribution of kinematic and dynamic properties of a non-Newtonian fluid flow in a lid-driven square cavity is considered. The power-law model is used as the rheological model. The numerical solution is received using the indirect boundary element method in the creeping flow approximation. The study is performed in the range of the power-law index from 0.2 to 1.2. The velocity component profiles at the mid-span of the cavity are obtained. For the case of the Newtonian fluid, a comparison with known results showed a good agreement. It is shown that the position of the main vortex shifts towards the upper moving lid as the power-law index decreases. The fields of effective viscosity and deformation rate intensity inside the flow domain are presented.
Keywords: non-Newtonian fluid, flow in lid-driven cavity, Indirect Boundary Element Method.
PONOMAREVA Maria Andreevna (Candidate of Physics and Mathematics, Tomsk State University, Tomsk, Russian Federation)
E-mail: pma@ftf.tsu.ru
FILINA Maria Petrovna (Tomsk State University, Tomsk, Russian Federation)
E-mail: filina.mari@mail.ru
YAKUTENOK Vladimir Albertovich (Doctor of Physics and Mathematics, Prof., Tomsk State University, Tomsk, Russian Federation)
E-mail: yva@ftf.tsu.ru
Течение неньютоновской жидкости в квадратной каверне при малых числах Рейнольдса 99
REFERENCES
1. Elizarova T.G., Milyukova O.M. Chislennoe modelirovanie techeniya vyazkoy neszhi-maemoy zhidkosti v kubicheskoy kaverne. Zhurnal vychislitel'noy matematiki i matematicheskoy fiziki, 2003, vol. 43, no. 3, pp. 453-466. (in Russian)
2. Haque S., Lashgari I., Giannetti F. Stability of fluids with shear-dependent viscosity in the lid-driven cavity. J. Non-Newtonian Fluid Mechanics, 2012, no. 173, pp. 49-61.
3. Anderson P.D., Galaktinov O.S., Peters G.W. Mixing of non-Newtonian fluids in time-periodic cavity flows. J. Non-Newtonian Fluid Mechanics, 2000, no. 93, pp. 265-286.
4. Lashckarbolok M., Jabbari E. Collocated Discrete Least Squares (CDLS) meshless method for the simulation of power-law fluid flows. Scientia Iranica, 2013, no. 20(2), pp. 322-328.
5. Nejat A., Sharbatdar M., Jalali A. A Newton-Krylov finite volume algorithm for power-law non-Newtonian fluid flows using pseudo-compressibility method. 20th AIAA Computational Fluid Dynamics Conference AIAA, 2011, pp. 1-19.
6. Li Q., Hong N., Shi B. Simulation of power-law fluid flows in two-dimensional square cavity using multi-relaxation-time lattice boltzmann method. Commun. Comput. Phys., 2014, vol. 15, no. 1, pp. 256-284.
7. Liu M.B., Xie W.P., Liu G.R. Modeling incompressible flows using a finite particle method. Applied Mathematical Modeling, 2005, vol. 29, pp. 1252-1270.
8. Shamekhi A., Aliabadi A. Non-Newtonian lid-driven cavity flow simulation by mesh free method. ICCES, 2009, vol. 11, no. 3, pp. 67-72.
9. Brebbiya K., Telles Zh., Vroubel L. Metody granichnykh elementov. Moskow, Mir Publ., 1987. (in Russian)
10. Ladyzhenskaya O.A. Matematicheskie voprosy dinamiki vyazkoy neszhimaemoy zhidkosti. Moskow, Nauka Publ., 1970. (in Russian)
11. Ponomareva M.A., Filina M.P., Yakutenok V.A. The indirect boundary element method for the two-dimensional pressure- and gravity-driven free surface Stokes flow. WIT Transactions on Modelling and Simulation, 2014, vol. 57, pp. 289-304. DOI: 10.2495/BE370241.
