Non-Newtonian fluid flow in a liddriven cavity at low Reynolds numbers.pdf Задача о течении вязкой несжимаемой жидкости в квадратной каверне с верхней стенкой, движущейся с постоянной скоростью в горизонтальном направлении, является известным и широко используемым тестом для верификации и оценки эффективности численных методов, применяемых для решения уравнений Сто^а. В связи с широким использованием такого варианта тестирования, имеется большое количество накопленных данных, позволяющих определить правильность новых результатов, тем самым гарантируя их достоверность. Основное количество результатов имеется для течений при больших числах Рейнольдса [1-6]. В работах [7, 8] рассмотрены случаи малых чисел Рейнольдса Re = 1 и 0. Рассмотрение в этих работах проведено для ньютоновской жидкости. В то же время актуальным представляется проведение исследований данного течения при малых числах Рейнольдса, но с учетом неньютоновского поведения жидкости. Основные уравнения Основными уравнениями для описания двумерного течения неньютоновской жидкости при малых числах Рейнольдса являются уравнения Стокса. Математическая постановка задачи записывается в безразмерных переменных. В качестве характерного размера выбрана длина стенки каверны, в качестве характерной скорости - скорость движения верхней крышки каверны. В приближении ползущего течения и при отсутствии массовых сил уравнения движения имеют вид = 0, i, j = 1,2, daj dx, (1) где a j = - p5 j + т j - компоненты полного тензора напряжений, p - давление, 5j - символ Кронекера, тj - компоненты тензора вязких напряжений, Xj - декартовы координаты. В качестве реологической модели, описывающей неньютоновское поведение жидкости, выбран степенной закон, также известный как модель Оствальда де-Виля, характеризующий зависимость вязкости от скорости сдвига: Tj = 2nej, (2) где n = Y"-1 - коэффициент эффективной вязкости, п - показатель нелинейности, 1/2 Y = (2eijeji) - интенсивность скоростей деформаций, etj = (dui /dXj + duj /dxi)/2 - компоненты тензора скоростей деформаций, ui - компоненты вектора скорости. Жидкости, поведение которых описывается данной моделью, можно подразделить на три типа: псевдопластичные (п < 1), ньютоновские (п = 1) и дилатантные (W > 1) . На рис. 1 представлена область решения. Верхняя стенка каверны движется в горизонтальном направлении с постоянной скоростью. Рис. 1. Область решения Систему (1) необходимо дополнить уравнением неразрывности (3) - = 0 дх, и граничными условиями, которые состоят в задании компонент скорости на стенках каверны: на верхней подвижной стенке пх = 1, п2 = 0 и на остальных стенках (4) (5) u = 0. Гранично-интегральная формулировка задачи и метод решения Представим уравнение (1) в виде N dai. _ь_ (6) дх, где aj = - p5j + 2eif - линейная часть полного тензора напряжений; - нелинейная векторная функция, которую будем dj ij дх, д * i = &- [2(1 -n* j j рассматривать как плотность источников, распределенных по области течения Q . Тогда, в соответствии с положениями непрямого метода граничных элементов [9], можно записать щ (x) = jGl} (x,|)ф; (|)dГ(|) +J Gl} (x,z)Tj (z)dQ(z), (7) Г Q где ф j (|) - плотность фиктивных источников, распределенных по границе области течения Г . Функции G, являются фундаментальными решениями уравнений Стокса и определяются формулой [10] 1 U , 1 >'iyj G«(*- 5) = -ds„ ta- -у-1, |8) -I 5:,- Ы - H--~ 1 1 r r 2 1 где у, = x, -1,, r = (yy )2 . Если на границе области течения Г заданы значения скорости, то уравнения (7) позволяют получить значения неизвестных граничных сил ф j (|) (| е Г). Это возможно сделать при известной функции Т, (z) (z е Q). Так как эта функция заранее неизвестна, возникает необходимость организации итерационного процесса. Для численного решения уравнений (7) используются постоянные элементы и постоянные ячейки. Граница области течения Г разбивается на N элементов. Функции фj (|) считаются постоянными на каждом элементе. Область течения разбивается на N2 квадратных ячеек. Функции Т, (z) считаются постоянными внутри ячейки. Тогда уравнения (7) в дискретной форме приобретают вид N N N u (xp ) = И AGf +Х kL AGf”, (9) q=1 k=1 -=1 где AG,Pq = J G,, (xp,dГ(|), AGpk- = J G, (xp,z)dQ(z), xp - середина A^ AQkL элемента p (узел). Для вычисления 2N неизвестных ф^ используются 2N уравнений (9), соответствующие N элементам, на которых заданы ut (xp). Коэффициенты получаемой системы линейных алгебраических уравнений AGijpq в случае постоянных элементов можно вычислить аналитически. Технология вычисления изложена в [11]. Для вычисления интегралов по области AGijpqm используются стандартные квадратурные формулы Гаусса без выделения особенностей. Особенности в этих интегралах имеют вид ln(1 / r). Следовательно, при интегрировании по области эти интегралы существуют в обычном смысле. Такой подход значительно упрощает алгоритм решения. При проведении расчетов использовалась квадратурная формула с 64 узлами. Для решения системы нелинейных алгебраических уравнений (9) относительно ф^ применялся метод простой итерации. На первой итерации использовались значения Тk- определенные по ньютоновскому полю течения. Для решения соответствующей системы линейных алгебраических уравнений использовался метод Гаусса. Далее использовались значения *^ , рассчитанные в соответствии с полем течения, полученным на предыдущей итерации. Функции *в центре ячейки (k, m) вычислялись конечно-разностным способом с использованием рассчитанных (t;Nn ) в вершинах ячеек (узлах сетки). Значения (t;Nn ) полностью определяются производными (dui / Xj ) , значения которых находились с использованием центральных разностей во внутренних узлах и односторонних разностей в приграничных узлах в соответствии со значениями uikm , вычисленными в узлах сетки. Для улучшения сходимости итерационного процесса в случае сильной нелинейности использовался метод релаксации с динамическим подбором коэффициента релаксации. Кроме того, процесс вычисления оптимизировался путем распараллеливания алгоритма. Анализ полученных результатов Вышеописанный метод расчета позволил провести вычисления в диапазоне изменения параметра нелинейности п от 0.2 до 1.2. Результаты расчетов представлены на рис. 2 - 8. Исследование аппроксимационной сходимости (рис. 2) показало, что приемлемая точность расчетов обеспечивается при разбиении границы Г на 160 и более элементов. Представленные в работе расчеты были проведены при N = 256 . 0.2 0 -0.4 0 0.4 0.8 щ -0.4 -0.2 0 0.2 0.4 *i Рис. 2. Аппроксимационная сходимость решения для п = 0.5 : а - профили составляющей скорости u1 (х1 = 0), б - профили составляющей скорости u2 (х2 = 0.5) *2 0.8 0.6 0.4 Изменение параметра n приводит к смещению центра основной вихревой зоны. Причем уменьшение п по сравнению с ньютоновским случаем (п = 1.0) сопровождается смещением центра к верхней крышке, а увеличение - к нижнему основанию. Соответствующие качественные и количественные результаты показаны на рис. 3 и 4. -0.5 -0.3 -0.1 0.1 0.3 *1 -0.5 -0.3 -0.1 0.1 0.3 *1 Рис. 3. Изолинии функции тока Ф ( u1 = ; u2 = -dl/dq ) с шагом АФ = 0.2 для раз личных значений показателя нелинейности: и = 0.5 (а), и = 0.7 (б), и = 1.0 (в), и = 1.2 (г) Рис. 4. Зависимость координаты центральной точки основного вихря от величины n 0.84 0.80 0.76 0.2 0.4 0.6 0. В работах, касающихся исследования данного течения, стандартно приводятся профили компонент скорости u1 и u2 при х1 = 0, х2 = 0.5 соответственно. Для рассматриваемого в настоящей работе нелинейного случая и малых чисел Рейнольдса указанные данные содержатся на рис. 5 и 6. х2 0.8 0.6 0.4 - a - Н-1-1- П-0.2 -е-^^я=о.5 ■ -в-п-0.7 - -/2=1.0 ООО и-1.2 . 1 1 1 1 0 0.4 0.8 щ -0.4 0.8 mi 0.2 -0.4 Рис. 5. Профили составляющей скорости u1 вдоль линии х1 = 0 для 0.2 < п < 1.2 (а) и сравнение с известными данными для ньютоновской жидкости (б) 0.4 х\ Рис. 6. Профили составляющей скорости u2 вдоль линии х2 = 0.5 для 0.2 < п < 1.2 (а) и сравнение с известными данными для ньютоновской жидкости (б) Распределение величины эффективной вязкости n по области течения соответствует распределению интенсивности скоростей деформаций у и параметру нелинейности п (псевдопластичному и дилатантному реологическому поведению). Это соответствие показано на рис. 7, 8. Зоны высоких значений эффективной вязкости наблюдаются в областях низких скоростей сдвига для псевдопла-стичных жидкостей (центр и нижние углы каверны) и в областях высоких скоростей сдвига для дилатантных жидкостей (верхние углы каверны). Размеры таких зон увеличиваются при уменьшении п (п < 1). Рис. 7. Поле интенсивности скоростей деформаций у для различных значений величины и : и = 0.2 (а), и = 0.5 (б), и = 0.7 (в), и = 1.0 (г), и = 1.2 (д) Рис. 8. Поле эффективной вязкости п внутри области течения для различных значений показателя нелинейности: и = 0.2 (а), и = 0.5 (б), и = 0.7 (в), и = 1.2 (г) (область течения, закрашенная черным цветом, соответствует значениям п>55 для псевдопластичной жидкости и п> 1.8 для дилатантной) Заключение Непрямым методом граничных элементов проведено исследование течения неньютоновской жидкости в квадратной каверне при малых числах Рейнольдса. Значение показателя нелинейности и варьировалось в диапазоне от 0.2 до 1.2. Получены профили компонент вектора скорости в характерных сечениях каверны. Для ньютоновской жидкости проведено сравнение полученных данных с результатами других исследователей. Показано, что уменьшение величины и приводит к смещению центра, вокруг которого вращается жидкость, к верхней крышке каверны. Представлены поля распределения эффективной вязкости и интенсивности скоростей деформаций по области течения. Полученные результаты могут быть использованы для тестирования численных методов решения уравнений Стокса при степенной зависимости вязкости от скорости сдвига.

Скачать электронную версию публикации