CC BY
104
Вычислительные технологии
Том 17, № 1, 2012
Численный алгоритм для моделирования установившихся ламинарных течений неньютоновских жидкостей в кольцевом зазоре с эксцентриситетом
A.A. Гаврилов1, А. В. Минаков2, А. А. Дектерев1, В. Я. Рудяк3 1 Институт, теплофизики СО РАН, Новосибирск, Россия 2 Сибирский федеральный университет, Красноярск, Россия 3Новосибирский государственный архитектурно-строительный университет,
Россия
e-mail: gavandSyandex. ru, valery. rudyak@mail. ru
Представлен численный алгоритм расчета установившихся ламинарных течений неньютоновских вязкопластических жидкостей в кольцевых каналах с эксцентриситетом и вращением внутреннего цилиндра. Алгоритм основан на методе контрольного объема для структурированной сетки и БНИРЬЕ-процедуре для связи полей давления и скорости. Градиент давления вдоль течения определяется из интегрального баланса сил. Для моделирования эффективной вязкости жидкостей с предельным напряжением используется метод монотонной регуляризации. Алгоритм позволяет описывать данный класс течений в широком диапазоне изменений параметров кольцевого канала и режимов течения. Приведены детали реализации и результаты тестирования предложенной численной методики. Для ряда течений в кольцевом зазоре проведено сопоставление численных результатов с известными аналитическими решениями или экспериментальными данными. Во всех рассмотренных случаях получено хорошее согласование с тестовыми данными.
Ключевые слова: метод контрольного объема, Б1МРЬЕ-процедура, неньютоновская жидкость, ламинарное течение.
Введение
Течения в цилиндрических каналах широко распространены в различных технологических приложениях: в теплообменниках, подшипниках скольжения, центрифугах, в некоторых видах миксеров, при бурении скважин и т. д. Ламинарное течение в концентрическом канале для ньютоновской жидкости описывается известным точным решением [1]. На практике, однако, приходится иметь дело с течением в канале с эксцентриситетом, а внутренняя (или внешняя) труба может вращаться. Это существенно усложняет задачу описания течения. Кроме того, помимо обычной ньютоновской жидкости в таком канале могут иметь место течения различных неньютоновских жидкостей. Во всех таких случаях удовлетворительных аналитических решений задачи не существует. Но поскольку актуальность этой задачи с практической точки зрения чрезвычайно высока, то в последнее десятилетие для ее решения было развито несколько численных методов [2-4]. Алгоритм для описания течений ньютоновской жидкости, основанный на методе контрольного объема, предложен и авторами данной работы [5].
Первые данные о свойствах течений неньютоновских жидкостей в каналах с эксцентриситетом появились лишь в 1980-х годах. Наибольший успех связан с изучением течений так называемых степенных жидкостей (см., например, [6-8] и цитируемую там литературу), В работе [8] проведен сравнительный анализ нескольких формул, предложенных для нахождения падения давления в таких каналах. Оказалось, что все проанализированные соотношения, описывающие падение давления, — неудачны. Более успешной была корреляция, предложенная в [9]. Ее точность составляет 10-15%, Однако и она работает лишь при индексах жидкости п > 0.4 и в ограниченном диапазоне отношений радиусов внешней и внутренней труб. Самые полные и последовательные численные исследования течения степенной жидкости в эксцентричном канале проведены в [10], где результаты получены для различных степенных жидкостей, чисел Тейлора и Рейнольдса,
Для жидкостей Бингама вплоть до настоящего времени вообще отсутствовали данные о течении в эксцентричном канале, тем более с вращением. Некоторые экспериментальные результаты содержатся в работе [10] (см, также цитируемую там литературу), В [11] сделана попытка решения задачи о напорном течении с помощью некоторого вариационного метода. Почти такая же картина характерна для жидкостей Херше-ля —Балкли, здесь можно указать лишь работу [12], Сложность решения задачи для неньютоновских жидкостей связана главным образом с двумя обстоятельствами: зависимостью коэффициента эффективной вязкости от скорости сдвига и наличием предельного напряжения. Это не позволяет для решения задачи непосредственно использовать алгоритмы, разработанные для описания течений ньютоновских жидкостей, в частности, алгоритма, предложенного в [5], Необходима его модификация. Разработка такой модификации алгоритма [5] и является целью данной работы,
1. Математическая модель
Рассматривается изотермическое ламинарное полностью развитое стационарное течение несжимаемой неньютоновской жидкости с постоянной плотностью и заданным расходом Q, В этой постановке течение описывается следующими уравнениями гидродинамики:
где р — плотность жидкости, V — вектор скорости, р — давление, Т — тензор вязких напряжений, который для рассматриваемых неньютоновских жидкостей определяется выражением Т = /1еО, где эффективный коэффициент вязкоети жидкости /1е зависит от скорости сдвига 7:
Далее рассматривались следующие неньютоновские жидкости: жидкости со степенным реологическим законом (степенные жидкости), ^е = к^-1, жидкости Бингама, = (то + кА/)А/-1, жидкости Хершеля — Балкли, /1е = (т0 + кА/п)А/-1, и жидкости Роберт-сопа — Стиффа, ^е = (т0 + -1, Здее ь т0 — предел текучести вязко пластической
жидкости. Ниже определенного предельного значения напряжений среда ведет себя
пк
логических моделей.
V- рv = 0, V- р\\ = -Ур + V- Т
(1)
Расчетная область представляет собой кольцевой канал, образованный двумя гладкими прямыми трубами круглого сечения. Радиус внутренней трубы R1; внешней R2, Ось внутренней трубы может быть смещена относительно оси внешней трубы на величину AR, определяющую эксцентриситет. Далее используется безразмерное значение эксцентриситета: e = AR/(R2 — R1), Внутренняя труба может вращаться вокруг своей оси с постоянной угловой скорость ш. На поверхностях цилиндров выполняются обычные условия прилипания. Течение жидкости в кольцевом канале в общем случае зависит от пяти безразмерных параметров: отношения радиусов в = R1/R2, эксцентриситета e, показателя степени n, числа Бннгама Вп = т0/(кГп), отношения вращательного числа Рейнольдса Rer к аксиальному Re, Здесь введена характерная скорость сдвига
Г
\
wRi \2 / U
R2 — R1 R2 — R1 а указанные числа Рейнольдса определяются выражениями
Re = Р^ № ~ Ri) ш=ЫЪ-Ъ)и
r кГn-1 ktn-1
где среднемассовая скорость U = Q/(pS) (Q — расход жидкости), S — площадь поперечного сечения канала,
В приложениях первые три параметра обычно изменяются в следующих пределах: 0.1 < в < 0.99, 0 < e < 10 < Вп < 18000, В отсутствии вращения внутреннего цилин-в
Пуазейля, Основным интегральным результатом моделирования является продольный градиент давления dp/dz, безразмерный перепад давления вводится числом
faE,e = 2{Rl~R2)2
U (kfn-1
2. Алгоритм и его численная реализация
Численный метод, развиваемый в данной работе, основан на методе контрольного объема, Детали его реализации для описания течений ньютоновских жидкостей подробно описаны в [5]. Здесь отмечаются только наиболее принципиальные моменты и специфические черты алгоритма, обусловленные неньютоновской реологией жидкости и граничными условиями для установившегося течения.
Основная сложность при моделировании течений вязкопластической среды связана с существованием сингулярности в тех областях, где скорость деформации равна нулю. Для вязкопластических жидкостей с ненулевым предельным напряжением т0 в области течения |Т| < то (область твердого тела) скорость сдвига стремится к нулю, в силу чего эффективная вязкость жидкости неограниченно возрастает. Чтобы преодолеть данную трудность, используют различные регуляризации исходных реологических моделей. Простейшая из них состоит в искусственном ограничении эффективной вяз-
| Т| < то
величиной тах = шах(^е,^етах). Однако точность получающегося алгоритма в общем случае существенно зависит от выбора величины ^е тах. Заниженное значение ^е тах приводит к искусственному уменьшению области канала, занятой твердым телом, С другой
стороны, выбор заведомо завышенного значения ^е тах существенно снижает скорость сходимости задачи. Подобрать оптимальное значение ^е тах для всех рассматриваемых задач не удается. По этой причине данный метод при создании универсального алгоритма неприемлем.
Второй подход, предложенный Папанастасно [13], заключается в аппроксимации выражения для эффективной вязкости некоторым сглаженным выражением. Так, например, для жидкости Бингама вместо выражения /1е = (т0 + кА/п)А/-1 вводится следующее:
которое при стремлении y к нулю позволяет ограничить эффективную вязкость. Экспоненциальный множитель m должен быть достаточно большим (чаще всего его значение колеблется от 100 до 3000),
Конечно, последний подход также имеет свои недостатки, среди которых можно отметить проблему определения границ жестких зон в потоке и плохую обусловленность матриц диекретизированных систем. Последнее при больших параметрах регуляризации удается преодолеть путем построения подходящих предобуелавливателей [14], Вместе с тем метод регуляризации отличается простотой в реализации и его достаточно успешно можно применять при решении практических задач. При разработке алгоритма, описываемого в данной статье, используется именно этот метод. Ниже будет представлен результат методических расчетов по изучению влияния величины коэффициента регуляризации.
Поскольку рассматривается установившееся течение в канале с постоянной геометрией поперечного сечения при заданном расходе, то для его моделирования можно использовать метод периодических граничных условий [15], Предложенный алгоритм следует традиционной формулировке периодических граничных условий и состоит в выделении из поля давления периодической составляющей. Обычно падение давления между входным и выходным сечениями определяется с помощью метода расчета поправки градиента давления, основанного на условии постоянства заданного расхода, В нашем алгоритме периодическая составляющая градиента давления определяется из интегрального баланса сил.
Для моделирования течения в цилиндрическом зазоре применяется одноблочная структурированная сетка, построенная в ортогональной биполярной системе координат с использованием процедуры перераспределения узлов в азимутальном направлении. Это позволяет проводить практически однородное сгущение сетки к границам канала. Сетка имеет один слой контрольных объемов в направлении течения (вдоль оси z). Использование одного слоя ячеек означает, что производные всех полевых величин по z
давления, равны нулю. Для восстановления давления последнее разбивается на две составляющие: p = pi + p2, Первая pi вычисляется с использованием итерационного SIMPLE-подобного алгоритма и имеет нулевой градиент по z, т.е. pi = pi(x,y), где x, y — плоскость поперечного сечения капала, Вторая составляющая поля давления p2
z
постоянен: dp/dz = dp2/dz = const. Эта константа определяется из условия баланса
z
(2)
где Ь — длина расчетного участка канала, Е — поверхность стенок канала, в силу чего интегрирование выполняется по этой поверхности, п — единичный вектор нормали к поверхности стенки. Индекс г указывает на то, что в подынтегральное выражение входит аксиальная составляющая вектора, стоящего в скобках.
Определение продольного градиента давления явным способом только на основе сил трения на стенках (3) имеет недостаток. Абсолютные ошибки значений продольной и поперечной компонент трения на стенке имеют одинаковый порядок. Отношение характерных величин продольной и поперечной компонент трения на стенке пропорционально отношению чисел Рейпольдса Т/Т ~ Не/Ке,., Для течений с преимущественно вращательным движением, заданным вращением внутреннего цилиндра, Нег >> Ее, относительная ошибка определения продольного трения может достигать больших значений ~ £ (Кег/Ее), Соответственно относительная ошибка расчетного значения продольного градиента давления будет высока, В численном алгоритме большая погрешность определения градиента давления приводит к появлению массового дисбаланса, что выражается в неравенстве между заданным и расчетным расходами через центральное поперечное сечение канала. Для решения этой проблемы в работе используется подход, идеологически близкий к процедуре, применяемой при расчете периодических течений для восстановления градиента давления [15]. Чтобы массовый поток вдоль оси г сохранялся, в градиент давления, рассчитанный из условия баланса сил (3), вводится поправка
др др2
я- = "я- + д"> дг дг
пропорциональная массовому дисбалансу 8в = а(Яс — Я)/(£ь$)> гДе Я ~ заданный, Яс — расчетный массовый поток через попере чное сечение; ¿ь = Ь/и— характерное
время, рассчитываемое по характерной скорости и^ = ^(иЯ\)2 + [/2; а — множитель, регулирующий величину поправки. Введенная поправка при сходимости итерационного процесса и уменьшении дисбаланса массового потока стремится к нулю. Множитель а определяется на основе решения тестовой задачи течения ньютоновской жидкости в зазоре между коаксиальными цилиндрами с вращением внутреннего цилиндра со следующими параметрами: Не= 0.25, Кег = 250, Выполнены расчеты с различными значе-
а
ошибки определения коэффициента трения от величины а: для а = 0.5 £z = 0.01%, для а = 0.1 = 0.15%, для а = 0.01 = 1%, для а = 0.0 = 2.5%. Ир и а = 1.0 наблюдались неустойчивость численного алгоритма и отсутствие сходимости расчета. Для сходящихся расчетов количество итераций до сходимости не зависит от значений аа
коэффициента трения Исходя из результатов проведенного теста все дальнейшие расчеты выполнялись с коэффициентом а = 0.1, обеспечивающим хорошую точность (< 1 %) и устойчивость итерационного процесса. Таким образом, окончательно выражение для продольного градиента давления имеет вид
Применение алгоритма [5] для описания течений неньютоновской жидкости потребовало его модификации. Модифицированный итерационный процесс вычислений полей
скорости и давления на каждой итерации состоит из четырех этапов. На первом, решая линеаризованное второе уравнение системы (1), вычисляется промежуточное значение скорости V*: V • (р^у*) — V • Т* = —Ур\ где V1, р1 — значения скорости и давления на предыдущей итерации. Полученная после дискретизации линеаризованного уравнения система линейных алгебраических уравнений решается итерационным способом с применением метода неполной факторизации Б-1Ы1, в котором факторизованы только диагональные члены [16],
Эффективная вязкость в тензоре касательных напряжений определяется по полю скорости с предыдущей итерации. Используя градиент поля скорости с предыдущей 1-й итерации, рассчитывается локальная скорость сдвига 7, по которой затем вычисляется значение эффективной вязкости с применением процедуры регуляризации для заданной реологии.
Расчет эффективной вязкости на стенке ведется итерационным способом по реологическим формулам. Считая величину трения на стенке известной и используя значение вязкости па стенке с предыдущего итерационного слоя определяем скорость сдвига вблизи стенки = Тщ вде Тщ — модуль вектора касательных напряжений на стенке. Затем по реологическим соотношениям с учетом аппроксимации вида (2) рассчитывается текущее значение молекулярной вязкости на стенке.
На втором этапе алгоритма рассчитывается поправка давления р' и определяется "поперечная" составляющая поля давления на новой итерации для обеспечения выполнения условия неразрывности поля скорости в поперечном сечении канала У(т Ур') = У(ру*), где т — локальный псевдо-временной шаг. Затем по найденной поправке давления проводится коррекция скорости и давления в центрах контрольных объемов = V* — (т/р)Ур', р!+1 = р\ + р'. Для определения давления на стенках применяется процедура интерполяции из расчетной области.
После нахождения поля скорости в центрах контрольных объемов, удовлетворяющего уравнению неразрывности, выполняется процедура коррекции расхода [5] для определения распределения скорости на входной и выходной границах. Процедура состоит из следующих операций. Используя экстраполяцию из внутренней области, реализуем "мягкое" граничное условие для скоростей на входной и выходной границах ду/дп = 0, Полученное распределение скоростей может не удовлетворять заданному расходу, поэтому после экстраполяции значения скорости на границе корректируются с учетом заданного массового расхода. Скорректированные скорости на входной и выходной границах используются в качестве граничных для уравнения количества движения (условия типа условия Дирихле),
Наконец, на последнем этапе по полученному полю скорости рассчитываются касательные напряжения и по формуле (4) восстанавливается продольный градиент давления, В качестве критерия сходимости процесса глобальных итераций задается условие малости средних величин относительных приращений компонент скорости,
В силу того что в уравнения движения и в определяющие соотношения для течения несжимаемой жидкости входит только градиент давления, знание величины абсолютного давления для расчета не требуется, и для моделирования течения используется поле относительного давления. Для численного метода необходимо задать положение точки пулевого относительного давления г0 (расположение этой точки произвольно). Принимается, что значение абсолютного давления в этом узле известно, а величина относительного давления всегда должна быть равна нулю, С помощью найденного градиента и заданного значения абсолютного давления в точке нулевого относительного
давления поле абсолютного давления переечитываетея по формуле
(др\
¿Ы • (г - -го) • 3. Тестирование численного алгоритма
Для определения влияния величины коэффициента регуляризации т, входящего в выражение (2), выполнялись расчеты течения бингамовской жидкости в концентрическом канале со следующими параметрами: отношение радиусов в = 0.5, Ее = 500, Вп = 50. Использовалась сетка из 60 х 120 узлов со сгущением к стенке. Сходимость считалась полученной при достижении относительного итерационного отклонения продольной компоненты скорости значения < 10-7, Расчетные данные сравнивались с эталонным решением, в качестве которого использовался численный расчет с т = 2000. В табл. 1 приведены интегральные результаты методических расчетов: значения числа /аКе, число итераций г до сходимости, относительное отклонение А числа faRe от эталонного и отклонение расчетного профиля скорости '(г) от эталоппого '0(г)
А' = /|„ - 'оИгД^г.
Я.1 ' Яг
Видно, что при малых значениях т (т < 100) возможны большие ошибки определения как интегральных параметров, так и пространственных профилей скорости. С т
т
значение можно считать компромиссным вариантом, который обеспечивает требуемую точность расчетов при умеренных вычислительных затратах.
Тестирование представленного алгоритма было проведено на большом числе различных задач в широких диапазонах изменения параметров течения и геометрии кольцевого канала. Расчетные данные сопоставлялись с известными аналитическими решениями, экспериментальными данными и расчетами других авторов. Во всех случаях полученные численные результаты хорошо согласовывались с тестовыми. Ниже приведены три примера такого сопоставления.
3.1. Течение степенной жидкости в концентрическом кольцевом канале
Степенные жидкости являются достаточно типичными неньютоновскими жидкостями. Данной реологической модели удовлетворяют течения ряда растворов полимеров, некоторых лаков и красок, многих суспензий твердых частиц, некоторых видов буровых
Таблица 1. Результаты методических расчетов
т /аЫе А, % Ат, % г
10 265.3 1.2 2.4 700
100 267.3 0.48 0.68 1400
1000 268.4 0.07 0.1 1800
2000 268.6 0 0 2000
растворов и т.д. В случае ламинарного установившегося течения степенной жидкости в концентрическом кольцевом канале известно аналитическое решение, предложенное в [17]. В частности, профиль аксиальной скорости имеет вид
г/В-2
2 кЬ
©
г/В-2
'в2 \т Г
га = —,
( в2\т Г п
(5)
где Др — перепад давления на длине канала Ь, а параметр в находится из уравнения
в / по \ т 1 / в2
/(?-<)«-/И)«-
© 4 7 в
Связь перепада давления Др с объемным расходом Я определяется выражением
На рис. 1, о решение (5) (сплошные линии) сопоставлено с расчетными данными. Использовалась достаточно грубая равномерная сетка (40 х 120 узлов). Во всех случаях наблюдается хорошее согласование расчетных данных с точным решением. С уменьшением показателя степени п градиент скорости в пристеночной области увеличивается, поэтому в общем случае для численного разрешения пристеночной области необходимо использовать неравномерные расчетные сетки со сгущением узлов к стенкам канала.
В табл. 2 приведено сопоставление расчетного и аналитически полученного перепада давления при течении степенной жидкости в кольцевом канале, выполненное для двух расчетных сеток, состоящих из 20 х 80 и 40 х 120 узлов. Видно, что с увеличением детализации сетки расчетное значение faRe приближается к аналитическому. Как уже
Ш 0.1
0.05
\¥ 0.08 0.06 0.04 0.02
0.2
0.4
0.6
0.8
0.2
0.4
0.6
0.8
Рис. 1. Профили скорости для течения степенной (а) и бингамовской (б) жидкостей в концентрическом кольцевом канале в зависимости от безразмерного радиуса; сплошные линии — решение (5), точки — расчет; в = 0.5, Ые = 67
Таблица 2. Сравнение расчетных и аналитически полученных значений /аЫе для течения степенной жидкости в концентрическом канале
п Аналитическое Расчетное решение
решение сетка 20 х 80 сетка 40 х 120
0.1 2.93 3.47 2.96
0.5 7.94 7.97 7.95
0.85 19.39 17.18 19.23
было отмечено, максимальное расхождение наблюдается для минимального из рассмотренных значений показателя степени, однако и в этом случае на сетке 40 х 120 ошибка не превышала 1 %.
3.2. Течение бингамовской жидкости в концентрическом кольцевом канале
Другой важный в приложениях пример неньютоновских жидкостей — бингамовекая жидкость. Для начала течения такой жидкости необходимо преодолеть предельное напряжение т0. После этого данная среда ведет себя как обычная жид кость с вязкостью к. Жидкостями Бингама являются различные пасты, буровые, цементные и глинистые растворы, сточные грязи, некоторые виды красок и т.д. Профиль скорости ламинарного установившегося напорного течения бингамовской жидкости в концентрическом кольцевом канале описывается следующим аналитическим решением [18]:
Ш (г):
Ар Я:
2кЬ
-То(/31-в)-1(/312-02)+/321п(|
-То 1 -
Я
- М-
2 I и2,
/31 < ^ < /32,
г Ё2
(6)
где
/---Т0 + л/т02 + 4/32
Р=у/%-Т0Р2, & = \ --
То
2 т0Ь
ЛрД2
а параметр в2 находится из уравнения
2/32 (132 - То) 1п ~ 1 + (То + ©)2 + 2Т0 (1 - /32) = 0.
Связь перепада давления Ар с объемным расходом Q выражается формулой
Я = ^^ ((1 - в4) - 2/32(/32 - То)(1 - в2) - |(1 + 03)То + 1(2/32 - То)3То) .
На рис, 1,6 решение (6) сопоставлено с численным решением для течения бинга-
в = 0.5
значений числа Бингама: Вп = 5, 50 и 250, Численное решение получено па равномерной расчетной сетке, состоящей из 40 х 120 узлов. Согласование решения (6) с численным решением также везде хорошее. Однако следует отметить, что точность расчетов при больших значениях числа Бингама сильно зависит от детализации расчетной сетки
г
Таблица 3. Сравнение расчетных и аналитически полученных значений /аЫе для течения бингамовской жидкости в концентрическом канале
Число Аналитическое Расчетное решение
Бингама решение сетка 20 х 80 сетка 40 х 120
5 38.4 38.8 38.6
50 153.4 160.7 155.6
250 595.9 652.7 610.3
в области пограничного слоя, который может быть очень тонким. Для Вп = 104 толщина погранслоя составляет около 0,03% от ширины кольцевого зазора. Ясно, что для разрешения таких узких погранслоев необходимо использовать неравномерную сетку с детализацией пристеночной области.
Напорное течение бингамовской жидкости имеет две характерные области (см, рис, 1, б), В центре канала расположена область твердого ядра потока, в которой т < т0 и скорость однородна. Во второй области, у стенок канала, имеет место течение вязкой жидкости (т > то). С ростом числа Бингама размеры квазитвердой области увеличиваются,
В табл. 3 приведено сравнение расчетного и аналитически полученных значений безразмерного перепада давления в кольцевом канале faRe для рассмотренного течения бингамовской жидкости. Для наибольшего из чисел Бингама максимальная ошибка не превышала 3 % на сетке 40 х 120, В целом точность решения с увеличением детализации сетки возрастает,
3.3. Течение жидкости Хершеля — Балкли в концентрическом кольцевом канале с вращением внутренней трубы
Жидкость Хершеля — Балкли представляет собой суперпозицию степенной и бингамовской жидкостей, в силу чего в течениях проявляются особенности, характерные и для степенной, и для бингамовской жидкостей. Наличие вращения внутренней трубы существенно усложняет течение в кольцевом канале, поэтому тестовым примером является установившееся ламинарное течение данной жидкости в концентрическом кольцевом канале с вращением внутренней трубы. Численные результаты сопоставлены с экспериментальными данными работы [19], Внешний диаметр канала равнялся 0,065 м, внутренний — 0,04 м, в = 0.62, В качестве рабочей жидкости использовался 0,2%-й раствор карбопола, реология которого соответствует жидкости Хершеля — Балкли с параметрами т0 = 32, п = 0.43, к = 12.09, На рис, 2 приведены результаты сопоставления расчета безразмерных профилей аксиальной IV = Ш/и и тангенциальной V = У/{шК\) скоростей с экспериментальными, В первом случае расход жидкости составлял 0,1529 кг/с, а внутренняя труба вращалась с угловой скоростью 13,8 рад/с. Это соответствует значениям числа Рейнольдса Ее = 0.527, числа Тейлора Та = 2.58, Во втором случае расход равнялся 0,1504 кг/с, внутренняя труба вращалась с угловой скоростью 2,78 рад/с, что отвечает Ее = 0.12, Та = 0.0014, Расчет получен на равномерной сетке, состоящей из 20 х 120 узлов. Число Тейлора определяется соотношением
1.5
Ш
0.5
1
V
0 0.2 0.4 0.6 0.3 б
а
Рис. 2. Профили нормированной аксиальной Та = 2.58 и нормированной тангенциальной Та = 0.0014 в зависимости от безразмерного эксперимент [19]
0 0.2 0.4 0.6 0.3 5
б
скорости (а) для в = 0.62, 11е = 0.527, е = 0, скорости (б) для в = 0.62, Ие = 0.12, е = 0, радиуса. Сплошные линии — расчет, точки —
Из рисунка видно, что в обоих случаях наблюдается хорошее согласование расчетных и экспериментальных данных.
Заключение
Представленный численный алгоритм позволяет моделировать течения для практически произвольной геометрии кольцевого канала в широком диапазоне изменения параметров течения в рамках ламинарного режима. Так, в расчетах соотношение радиусов изменялось в диапазоне 0.1 < в < 0.9, а эксцентриситеты были практически любые (0 < е < 1). Расход жидкости, ее плотность, вязкость, скорость вращения внутренней трубы также варьировались в очень широких пределах. С помощью алгоритма можно моделировать течения практически всех наиболее распространенных вязкоплаетичее-ких флюидов. Помимо представленных в статье тестов, были выполнены расчеты для весьма сложной реологической жидкости — жидкости Робертсона — Стифа [6].
В работе для жидкостей с предельным напряжением использовался метод регуляризации эффективного коэффициента вязкости (см. формулу (2)). Альтернативой такому подходу могут быть вариационные методы. Численные методы на основе вариационных неравенств требуют, однако, введения дополнительных итерационных шагов (для реализации неявной связи полей вязких напряжений и скоростей деформации) и особой аккуратности при дискретизации дифференциальных уравнений [20]. Реализация использованного метода регуляризации на практике значительно проще и пригодна для выполнения массовых расчетов.
Вместе с тем следует иметь в виду, что в реальных течениях часто наблюдается турбулентность. Представленный алгоритм легко адаптируется для численного решении осредненных по Рейнольдсу уравнений Навье —Стокса с использованием различных моделей турбулентности, что, однако, сопряжено с преодолением ряда специфических трудностей.
Список литературы
[1] Ландау Л.Д., Лифшиц Е.М. Гидродинамика. М.: Наука, 1986. 735 с.
[2] Ooms G., Kampman-Reinhartz В.Е. Influence of drill pipe rotation and eccentricity on pressure drop over borehole during drilling // Europ. J. of Mechanics. B. Fluids. 1996. Vol. 15, No. 5. P. 695-711.
[3] Ooms G., Kampman-Reinhartz B.E. Influence of drill pipe rotation and eccentricity on pressure drop over borehole with newtonian liquid during drilling // SPE Drill k, Completion. 2000. Vol. 15, No. 4. P. 249-253.
[4] Escudier M.P., Gouldson I.W., Oliveira P.J., Pinho F.T. Effects of inner cylinder rotation on laminar flow of a Newtonian fluid through an eccentric annulus // Intern. J. of Heat and Fluid Flow. 2000. Vol. 21. P. 92-103.
[5] Гаврилов А.А., Минаков А.В., Дектерев А.А., Рудяк В.Я. Численный алгоритм для моделирования ламинарных течений в кольцевом канале с эксцентриситетом // Сиб. журн. индустр. математики. 2010. Т. 13, № 4 (44). С. 46-61.
[61 Bird R.B., Hassager О. Dynamics of polymeric fluids // Fluid Mech. 1987. Vol. 1, ch. 4. P. 169-253.
[7] Uner D., Ozge C., Tosun I. An approximate solution for non-Newtonian flow in eccentric annuli // Industr. & Eng. Chem. Res. 1988. Vol. 27. P. 698-701.
[81 Sestak J., zltny R., Dostal M. Flow of purely viscous non-Newtonian fluids in straight non-circular ducts //A Review and Comparison of Procedures for Rapid Engineering Friction Factor Estimates. European Congress on Rheologv. Cambridge, 2000.
[9] Haciislamoglu XL. Cartalos U. Practical pressure loss predictions in realistic annular geometries // SPE Annual Technical Conf. New Orleans. 1994. Paper SPE. No. 28304.
[101 Escudier M.P., Oliveira P.J., Pinho F.T. Fully developed laminar flow of purely viscous non-Newtonian liquids through annuli, including the effects of eccentricity and inner-cylinder rotation // Intern. J. of Heat and Fluid Flow. 2002. Vol. 23. P. 52-73.
[Ill Wachs A. Numerical simulation of steady Bingham flow through an eccentric annular cross-section by distributed Lagrange multiplier/fictitious domain and augmented Lagrangian method //J. Non-Newtonian Fluid Mech. 2007. Vol. 142. P. 183-198.
[121 Hussein Q.E., Sharif M.A.R. Viscoplastic fluid flow in irregular eccentric annuli due to axial motion of the inner pipe // The Canadian J. of Chem. Eng. 1997. Vol. 75. P. 1038-1045.
[131 Papanastasiou Т.е. Flow of materials with yield // J. of Rheologv. 1987. Vol. 31. P. 385-404.
[141 Гриневич П.П., Ольшанский М.А. Итерационный метод решения регуляризованной задачи Бингама // Вычисл. методы и программирование. 2010. Т. 11. С. 78-87.
[151 Быстрое Ю.А., Исаев С.А., Кудрявцев И.А., Леонтьев А.И. Численное моделирование вихревой интенсификации теплообмена в пакетах труб. СПб.: Судостроение, 2005. 392 с.
[161 Barrett R., Berry XL. Chan Т. Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods. SIAM, 1994. 141 p.
[171 Fredrickson A.G., Bird R.B. Non-Newtonian flow in annuli // Ind. Eng. Chem. 1958. Vol. 50, No. 3. P. 347-352.
[181 LajR-D W.M. Slurry and suspension transport. Basic flow studies on Bingham plastic fluids // Ibid. 1957. Vol. 49, No. 1. P. 138-141.
[191 Escudier М.Р., Oliveira P.J., Pinho F.T., Smith S. Fully developed laminar flow of non-Newtonian liquids through annuli: comparison of numerical calculations with experiments // Experiments in Fluids. 2002. Vol. 33. P. 101-111.
[20] Муравлёва E.A. Разностные схемы для расчета течений вязко-пластической среды в канале // Матем. моделирование. 2008. Т. 20, № 12. С. 76-88.
Поступила в редакцию 24 апреля 2011 г., с доработки — 15 сентября 2011 г.
