автореферат диссертации по информатике, вычислительной технике и управлению, 05.13.18, диссертация на тему:Вычислительные алгоритмы и комплекс программ для численного моделирования течений неньютоновских жидкостей в кольцевом канале

На правах рукописи

Гаврилов Андрей Анатольевич

ВЫЧИСЛИТЕЛЬНЫЕ АЛГОРИТМЫ И КОМПЛЕКС ПРОГРАММ ДЛЯ ЧИСЛЕННОГО МОДЕЛИРОВАНИЯ ТЕЧЕНИЙ НЕНЬЮТОНОВСКИХ ЖИДКОСТЕЙ В КОЛЬЦЕВОМ КАНАЛЕ

05.13.18 - Математическое моделирование, численные методы и комплексы программ

АВТОРЕФЕРАТ

диссертации на соискание ученой степени кандидата физико-математических наук

1 5 ЯНЗ 2015

Новосибирск — 2014

005557374

005557374

Работа выполнена в Федеральном государственном бюджетном учреждении науки Институте теплофизики им. С.С. Кутателадзе Сибирского отделения Российской академии наук, г. Новосибирск

Научный руководитель: кандидат технических наук,

Дектерёв Александр Анатольевич

Официальные оппоненты: Старченко Александр Васильевич,

доктор физико-математических наук, профессор, Национальный исследовательский Томский государственный университет, декан механико-математического факультета

Яковенко Сергей Николаевич кандидат физико-математических наук, доцент, Институт теоретической и прикладной механики им. С.А. Христиановича СО РАН, г. Новосибирск, старший научный сотрудник

Ведущая организация: Федеральное государственное бюджетное

учреждение науки Институт вычислительного моделирования СО РАН, г. Красноярск

Защита состоится 20 февраля 2015 г. в 10:00 часов на заседании диссертационного совета ДМ003.046.01 на базе Федерального государственного бюджетного учреждения науки Института вычислительных технологий Сибирского отделения Российской академии наук по адресу: 630090, г. Новосибирск, пр. Академика Лаврентьева, 6, конференц-зал ИВТ СО РАН.

С диссертацией можно ознакомиться в библиотеке и на сайте Федерального государственного бюджетного учреждения науки Института вычислительных технологий Сибирского отделения Российской академии наук

http://www.ict.nsc.ru/ruyStructure/disCouncil/gavrilov2014

Автореферат разослан « 9 » декабря 2014 г.

Ученый секретарь диссертационного совета, к.ф.-м.н. А.С. Лебедев

Общая характеристика работы

Актуальность темы

Задача о течении жидкости в канале между двумя цилиндрическими трубами является классической задачей гидродинамики, имеющей широкое практическое применение. На практике важно знать характеристики течений в широком диапазоне параметров (геометрия канала, свойства жидкости и режим течения), а также их зависимость от этих параметров. Задачу существенно усложняет неньютоновская реология рабочей среды, вязкость которой нелинейно зависит от скорости деформации течения. Универсальным инструментом исследования течений жидкости является численный эксперимент, который позволяет получить информацию о структуре течения при относительно небольших временных, трудовых и материальных затратах.

Отсутствие систематических экспериментальных данных относительно турбулентных течений неньютоновских жидкостей инициировало появление работ по численному моделированию турбулентных течений неньютоновской жидкости. Однако сложность состоит в том, что в настоящее время отсутствует замкнутая, хорошо обоснованная и экспериментально подтвержденная модель турбулентности для неньютоновских сред, в частности, для вязкопла-стических и псевдопластических жидкостей.

Несмотря на многолетнюю историю изучения данных течений, имеющийся материал все еще не может в полной мере дать необходимую информацию обо всех параметрах течения. В связи с этим необходимо иметь инструментарий, с помощью которого с хорошей точностью и в широком диапазоне параметров можно было бы предсказывать характеристики течений рассматриваемого класса.

Цель работы заключается в разработке численного алгоритма на основе методов вычислительной гидродинамики и программного обеспечения для моделирования установившихся течений неньютоновских сред в кольцевом канале в широком диапазоне изменения параметров течения.

Задачи, решенные в ходе достижения поставленной цели:

1. Сформулированы математические модели для полностью развитых ламинарных и турбулентных режимов течений нелинейно-вязкопластичных сред в кольцевом канале.

2. Разработан численный алгоритм решения уравнений гидродинамики установившихся течений неньютоновских сред на основе метода конечного объема и методики расщепления.

3. Создан программный комплекс для моделирования течений неньютоновских сред в кольцевом канале в широком диапазоне изменения параметров течения.

4. Проведена верификация разработанных алгоритмов и программного обеспечения посредством сопоставления результатов численных расчетов с аналитическими решениями, данными экспериментов и результатами, полученными другими авторами.

Научная новизна изложенных в диссертационной работе результатов заключается в следующем:

1. Предложена и численно обоснована оригинальная процедура определения градиента давления для установившихся течений с использованием интегральных балансных соотношений.

2. С использованием осреднения по Рейнольдсу обобщенных уравнений Навье-Стокса сформулирована новая модель осредненной молекулярной вязкости для турбулентных течений неньютоновских жидкостей.

3. Разработан двухшаговый алгоритм построения расчетной сетки для моделирования установившегося турбулентного течения в кольцевом канале.

4. На основе предложенных математических моделей и численных алгоритмов создан оригинальный инструментарий вычислительного эксперимента для моделирования течений неньютоновских сред в кольцевом канале. Программный комплекс позволяет проводить исследования пользователями, не обладающими специальными знаниями из области вычислительной гидродинамики.

Практическая значимость. Разработанные модели и их программная реализация позволяют проводить детальное исследование течения неньютоновской среды в кольцевом зазоре. Программный комплекс, созданный автором, используется для исследования характера течения в кольцевом зазоре и построения базы данных, представляющей зависимость силовых факторов от геометрии задачи и режимных параметров.

Метод исследования. Для решения прямых задач гидродинамики используются современные методы математического моделирования и численные методы.

Обоснованность и достоверность основных результатов, полученных в диссертации, основывается на строгом математическом описании используемых численных алгоритмов, детальных методических расчетах тестовых задач, сопоставлении результатов численных расчетов с данными экспериментов и результатами, полученными другими авторами.

На защиту выносятся следующие результаты, соответствующие трём пунктам (3, 4, 5) паспорта специальности 05.13.18 — «Математическое моделирование, численные методы и комплексы программ» по физико-математическим наукам:

1. Численный алгоритм решения уравнений гидродинамики установившихся течений неньютоновских сред в канале произвольного сечения на основе метода конечного объема и методики расщепления.

2. Двухшаговый алгоритм построения расчетной сетки для моделирования установившегося турбулентного течения в кольцевом канале.

3. Модель осредненной по Рейнольдсу молекулярной вязкости для турбулентных течений неньютоновских сред, подчиняющихся реологической модели Гершеля-Балкли.

4. Программный комплекс для моделирования течений неньютоновских сред в кольцевом канале в широком диапазоне изменения параметров течения.

5. Результаты тестовых расчетов с использованием разработанного программного обеспечения.

6. Результаты численного моделирования винтовых режимов течения неньютоновских сред, в частности, данные о формировании переходных режимов с разделением течения на область турбулентного движения, примыкающей к вращающемуся цилиндру, и область ламинарного течения.

Представление работы. Основные результаты диссертации докладывались на следующих научных мероприятиях: Международная конференция «Sixth International Symposium On Turbulence, Heat and Mass Transfer» (Рим, 2009); Всероссийский семинар по теплофизике и теплоэнергетике (Красноярск, 2009); Международная конференция «Современные проблемы прикладной математики и механики: теория, эксперимент и практика» (Новосибирск, 2011); Международная конференция «Математические и информационные технологии» MIT-2011 (Врнячка Баня, Сербия, 2011); Всероссийский семинар «Фундаментальные основы МЭМС- и нанотехнологий» (Новосибирск, 2011); 4-я Всероссийская конференция «Фундаментальные основы МЭМС- и нанотехнологий» (Новосибирск, 2012).

Публикации. По материалам диссертации опубликовано 11 работ, в том числе 4 статьи в рецензируемых журналах, рекомендованных ВАК для представления основных результатов диссертации. Список работ приведен в конце автореферата.

Личный вклад автора. Личный вклад автора заключается в разработке математической модели и численного метода, создании специализированной расчетной программы и выполнении ряда тестовых расчетов. Автор участвовал в постановке задачи и анализе полученных результатов.

Структура и объем работы. Диссертация состоит из введения, четырех глав, заключения, списка литературы из 102 наименований, 6 таблиц и 38 рисунков. Общий объем диссертации составляет 145 страниц.

Содержание работы

Во Введении обосновывается актуальность темы исследования, формулируется цель и ставятся основные задачи работы. Раскрываются научная новизна и практическая значимость работы. Формулируются основные положения, выносимые на защиту.

В главе 1 рассматриваются вопросы гидродинамического описания установившегося течения несжимаемой среды в кольцевом канале. Приводится обзор научной литературы по изучаемой проблеме. В параграфе 1.1 формулируется математическая постановка задачи о течении в кольцевом канале на основе уравнений Навье-Стокса для обобщенной ньютоновской среды. В этой постановке ламинарный режим течения описывается уравнениями Навье-Стокса:

V-(pv) = 0 (1)

V-(/?vv) = -Vp + V- т, (2)

где р - плотность жидкости, V - вектор скорости, р - давление. Для обобщенной ньютоновской модели тензор вязких напряжений т определяется следующим образом:

т = 21&,

где ¡г - эффективная или кажущаяся молекулярная вязкость, Б - тензор скоростей деформации. Кажущаяся вязкость является нелинейной функцией второго инварианта тензора скоростей деформации У:

у = 4= (3)

В общем случае рабочая среда является неньютоновской вязкопластической жидкостью, поведение которой описывается реологией Гершеля-Балкли:

т = КГ+т0,т>т0 П

■ ^ т = \— тт

7=0 т<т0 \2

при этом кажущаяся вязкость вводится выражением

М(у) = {КГ+Ч)/У. (4)

Здесь реология среды задается тремя коэффициентами: показателем степени п (п< 1), показателем консистенции к и пределом текучести вязкопластической жидкости 25. Ниже определенного предельного значения напряжений среда ведет себя как твердое тело (гс^), выше этого предела - как несжимаемая вязкая жидкость

Результатом решения задачи являются силовые факторы, обеспечивающие заданную кинематику жидкости. Основная искомая величина - это продольный градиент давления (перепад давления на единицу длины канала), соответствующий заданному расходу. Сопротивление среды вращению цилиндра формирует суммарный момент вязких сил, действующих на цилиндр. При наличии эксцентриситета суммарный вектор нормальных напряжений на стенке отличен от нуля. Его величина и направление также представляют практический интерес.

Для жидкостей с предельным напряжением вводится монотонная регуляризация кажущейся вязкости, что позволяет рассматривать вязкопластиче-скую среду как нелинейную вязкую жидкость. Для жидкости Гершеля-Балкли (4) в работе используется следующее выражение:

^ = [КГ + г0 (1 - ехрС-т г/ Г)]/г, (5)

которое при стремлении скорости деформации к нулю позволяет ограничить эффективную вязкость. Для выбора величины коэффициента регуляризации т выполнены методические расчеты течения жидкости Бингама.

В параграфе 1.2 изложен подход для моделирования турбулентных режимов течения на основе осреднения по Рейнольдсу. Приведено описание используемой дифференциальной модели турбулентности.

Широкое распространение в инженерной практике получили различные полуэмпирические модели, связанные с тем или иным способом замыкания осредненных по Рейнольдсу уравнений Навье-Стокса (ЯАЫБ). Согласно подходу Рейнольдса, любые мгновенные значения гидродинамических парамет-

ров потока представляются в виде суммы осредненной величины (по времени) и ее пульсационной составляющей: <р = (ф) + ф', где угловые скобки обозначают операцию осреднения по времени, штрих относится к пульсацион-ным величинам, а символ «крышечка» - к мгновенным величинам. В рамках модели вихревой вязкости уравнения Рейнольдса имеют вид:

^+V • (рШ) = V • [(//+//, )2в], от

где и - средняя скорость, Р - осредненное давление, в - осредненный тензор скоростей деформации, ¡1 - среднее значение коэффициента молекулярной вязкости, - турбулентная вязкость.

В работе используются две двухпараметрические дифференциальные модели турбулентности: к-е модель Abe-Kondoh-Nagano и к-со ББТ модель Ментера. Модель Abe-Kondoh-Nagano является модификацией двухпарамет-рической диссипатнвной к-е модели турбулентности на случай пристеночных течений. Модель Ментера, являющаяся суперпозицией моделей к-Е и к-со, хорошо зарекомендовала себя при моделировании пристеночных течений, в том числе течений с положительным градиентом давления. Для моделирования ограниченных стенками течений большую роль играет описание пристеночного слоя течений. В параграфе приводится обзор методов пристеночных функций.

Параграф 1.3 посвящен вопросам моделирования турбулентных режимов неньютоновских сред. Приводится обзор по существующим моделям турбулентности для неньютоновских сред. В турбулентных течениях молекулярная вязкость неньютоновской среды зависит от флуктуирующей скорости деформации. Зависимость молекулярной вязкости от скорости деформации течения приводит к дополнительным сложностям при формулировании осредненных по Рейнольдсу уравнений Навье-Стокса. При описании турбулентного течения необходимо оперировать осредненной величиной молекулярной вязкости, которая является функцией турбулентных характеристик потока.

На основе метода осреднения по Рейнольдсу предложена новая модель осредненной молекулярной вязкости для течений жидкостей Гершеля-Балкли. Эффективная молекулярная вязкость среды зависит от значения скорости деформации на малых масштабах, а следовательно, и от флуктуаций поля скорости. В выражении для определения осредненной скорости диссипации пренебрегаем тройной корреляцией (корреляция пульсаций коэффициента молекулярной вязкости с пульсациями скорости деформации), в результате имеем соотношение

ре = /г( 2«).

Выражение для средней скорости деформации представляется суммой двух слагаемых (3):

Первое слагаемое в правой части рассчитывается по градиентам средней скорости, а второе определяет среднюю пульсационную скорость деформации. Используя формулу для скорости диссипации, получаем выражение для средней скорости деформации:

Г =2 5,+0*0///. (6)

В этом же приближении можно считать, что среднее значение коэффициента молекулярной вязкости снова связано со средним значением скорости деформации соотношением (4).

Таким образом, значение коэффициента средней молекулярной вязкости можно найти, решив систему нелинейных уравнений (4) и (6). Для нахождения распределения коэффициентов средней молекулярной вязкости предлагается использовать двухшаговый итерационный алгоритм. На первом шаге по текущим значениям средней скорости, скорости диссипации и коэффициенту молекулярной вязкости рассчитывается значение локальной средней скорости сдвига (6). Затем по найденной средней скорости сдвига определяются значения коэффициента средней эффективной молекулярной вязкости (4), которые используются при решении уравнений гидродинамики на новом итерационном слое.

Глава 2 посвящена описанию численной методики. Численный метод основывается на методе контрольного объема для сетки произвольного типа с выпуклыми многогранными ячейками. Даны формулы для аппроксимации нестационарного, конвективного и диффузионного членов обобщенного дифференциального уравнения переноса. Представлен способ нахождения градиента сеточного поля. Описаны используемые методы решения системы линейных алгебраических уравнений, получающихся после дискретизации дифференциальных уравнений.

В параграфе 2.3 рассматриваются особенности дискретизации и численного решения уравнений Навье-Стокса для несжимаемых сред. Одним из основных подходов к построению численных алгоритмов для уравнений На-вье-Стокса является подход с использованием методики расщепления. В параграфе приведено построение схемы расщепления на основе нестационарных уравнений гидродинамики.

Запишем нестационарные уравнения Навье-Стокса (1) и (2), дискрети-зированные по времени:

= (7) У-(/7у") = 0, (8)

где п - номер временного слоя, г- временной шаг. В правой части динамического уравнения оператор Н объединяет линеаризованное конвективное и диффузионное слагаемые:

Н (V""1) V" = - V • (/7У"-1 у") + V • Т ( V").

Полученная схема обобщена на случай стационарной постановки при помощи введения псевдо-временного шага г, согласованного с процедурой

нижней релаксации для уравнения количества движения. Методика расщепления для стационарной задачи реализуется следующей последовательностью шагов:

1. Определение промежуточной скорости V* , которая является решением динамического уравнения (7) с давлением с предыдущего итерационного слоя:

где к - номер итерационного слоя. Величина массового потока на гранях контрольных объемов и эффективная вязкость берутся с предыдущей итерации.

2. Решение эллиптического уравнения Пуассона для поправки давления с использованием промежуточной скорости и выполнение последующей коррекции давления:

3. Определение поправки скорости и коррекция скорости на новом А>ом итерационном слое:

<5У = --У(<ЗД. у<=УЧ5У Р

Скорректированная скорость удовлетворяет уравнению неразрывности (8).

Цикл повторяется до момента достижения сходимости. Выбор параметра г осуществляется на основе рассмотрения дискретизированного уравнения переноса количества движения (7). Соответствующее выражение для определения локального псевдо-временного шага имеет вид (БГМРЬЕС процедура):

pJ а

т = ---

аР 1 -а'

где а - параметр релаксации, аР - диагональный коэффициент, J - объем контрольного ячейки Р.

В рассматриваемых расчетных задачах задана величина массового расхода Схс, рабочей среды, протекающей через поперечное сечение канала. Для выполнения глобального баланса массы в такой постановке используется процедура коррекции расхода. Процедура коррекции массового потока обеспечивает заданный расход через границы расчетной области с одновременным выполнением граничных условий д\/дп = 0. Расход через некоторую поверхность Г рассчитывается по формуле:

I

Если скорость на поверхности на данной итерации не удовлетворяет заданному расходу, вводится аддитивная поправка скорости, обеспечивающая заданный расход:

= |/?(у + Дг/п)-п<^5' = |/7у-П£/5' + Дм| рс/Б,

А«=7— pv-rnff

J/Ji/5

z

Эта поправка добавляется к вектору скорости на поверхности, в результате чего расход через границу становится равным заданному. При решении уравнения количества движения на поверхностях с заданным расходом ставятся граничные условия Дирихле.

Поскольку рассматриваются полностью установившиеся течения, то расчетная сетка может иметь один слой контрольных объемов в аксиальном направлении (вдоль оси z). Пример расчетной сетки показан на Рисунке 1.1. Для определения аксиального градиента давления предложена оригинальная процедура с использованием интегральных балансных соотношений. Алгоритм следует традиционной формулировке периодических вдоль течения граничных условий и предполагает выделение из поля давления периодической составляющей. Но в отличие от классического итерационного варианта, в новом алгоритме линейная составляющая давления определяется явным способом из интегрального баланса сил.

Поле давления разбивается на две составляющие р-р{+р2. Первая из них является периодической и не зависит от координаты z, так что Pi = РЛХ,у), где плоскость х,у - это плоскость поперечного сечения канала. Вторая составляющая поля давления зависит только от г и ее градиент в аксиальном направлении постоянен: dp/dz = др2 /dz = const Эта константа неизвестна (она является целевой функцией всего алгоритма) и должна быть получена в процессе решения задачи. Из условия баланса интегральных сил, действующих на объем жидкости в установившемся течении вдоль оси z, легко получить величину градиента давления вдоль оси z:

где £ - поверхность твердых стенок канала, так что интегрирование выполняется по этой поверхности, п - единичный вектор нормали к поверхности стенки. Индекс г соответствует аксиальной составляющей вектора.

Рисунок 1 - Пример сетки для расчета течения в кольцевом канале (1), разбиение расчетной сетки на слои (2)

Явный способ определения продольного градиента давления не требует введения дополнительного релаксационного множителя. Фактически итерационное изменение градиента давления регулируется коэффициентом релаксации для уравнения количества движения (7).

В параграфе 2.4 сформулирован алгоритм постановки граничных условий на стенке с помощью метода пристеночных функций для используемых моделей турбулентности. Метод пристенных, или пристеночных, функций подразумевает моделирование пристеночного слоя, т.е. введение модели, описывающей поведение течения вблизи стенки. Особенности пристеночного турбулентного слоя, которые не разрешаются на грубой сетке, вносятся в модель посредством специальных граничных условий. Дано детальное описание используемого варианта метода для к-е и к-(0 моделей. Метод применяет бесшовное объединение стандартных пристеночных функций и асимптотики для вязкого подслоя на основе метода с весами. Этот метод позволяет использовать единый алгоритм для сеток с произвольным разрешением турбулентного пограничного слоя.

Результаты тестирование программы для моделирования ламинарных и турбулентных режимов течения ньютоновской среды представлены в главе 3. Расчетная программа тестировалась на решении большого числа задач данного класса, описываемых либо аналитическими решениями, либо существующими экспериментальными данными. В параграфе 3.1 показано сравнение по профилю скорости расчетных данных для установившегося ламинарного течения в концентрическом канале с аналитическим решением. Показано сравнение расчетов с экспериментальными данными для вариантов течения в кольцевом канале с эксцентриситетом и вращением внутреннего цилиндра. Для всех задач получено отличное совпадение с эталонными данными.

Параграф 3.2 посвящен развитым турбулентным режимам течения ньютоновских сред. Приведены результаты моделирования простого сдвигового течения с использованием алгоритма пристеночных функций на различных по детализации сетках.

Из результатов моделирования развитых турбулентных течений в трубах сформулировано следующее требование к разностной сетке в пристеночной области (условие разрешимости градиентной области). Первый узел пристеночной зоны должен лежать в ламинарном подслое и иметь у+ ~1. Сетка должна быть равномерной в диапазоне 1<>>+<30 для разрешения градиентов турбулентных величин. Методические расчеты показывают, что безразмерный шаг пристеночной сетки Ду+ должен быть не менее 1.

Разработан алгоритм определения параметров экономичной расчетной сетки, позволяющий избегать попадания первого пристеночного узла сетки в буферный подслой (5<у+<30). Для построения сетки, позволяющей разрешить пристеночный слой, вводятся три сеточных слоя ячеек (Рисунок 1.2). Две области, непосредственно примыкающие к стенкам, имеют заданную толщину и равномерное распределение узлов по радиальному направлению.

Третья область сетки - внутренняя область, расположенная между пристеночными слоями с заданным сгущением узлов к границам.

Заданное количество узлов в пристеночной зоне, безразмерное расстояние от стенки до первого узла у+=1 и значение касательного трения на стенке однозначно задают толщину пристеночных слоев сетки. Отношение толщин примыкающих контрольных объемов пристеночной области и внутренней области сетки не должно превышать заданного числа, а количество узлов внутренней области сетки не должно превышать значения Л-п, < 100. Если для данного трения на стенке не удается построить сетку, удовлетворяющую этим условиям, то строится «однозонная» сетка без выделения пристеночных слоев с равномерным сгущением узлов к стенке. Построение такой грубой сетки происходит при больших числах Рейнольдса в тех случаях, когда для разрешения всех подобластей пограничного слоя требуется большое количество узлов.

Процедура построения окончательной сетки состоит из двух этапов. На первом этапе трение на стенках оценивается на основе имеющихся эмпирических корреляций для турбулентного течения в кольцевом канале с учетом эксцентриситета. На втором этапе по полученным в предварительном расчете величинам вязкого трения на стенке анализируется, какая в данном случае должна быть расчетная сетка, и выполняется построение окончательной расчетной сетки.

Для тестирования процедуры построения сетки были проведены две серии расчетов. В первой серии расчетов моделировалось напорное течение без вращения внутренней трубы при разных числах Рейнольдса. Расчетные значения коэффициентов сопротивления и данные по корреляции Jones and Leung приведены на Рисунке 2.1.

0 001

л+<1 Л >30

0.12 — F

0.08 —

ф расчет - — корреляция

"1 I I ! I . .ill I I ; ill:

Л+<1 yf > 30

.Г

1 "11111- I I I rill

1.0Е+5 Re

1 ОЕ+6 1.0Е+4 1.0Е+5 1.0Е+6 1 0Е+7

Re,,

Рисунок 2 - Коэффициент сопротивления для напорного течения без вращения (1) и для напорного течения с вращением внутреннего цилиндра (2)

На Рисунке 2.1 вертикальная полоса в районе Яе~2х105 обозначает область переключения с трехзонной сетки (<2х105) на однозонную сетку (>2x10 ). В области переключения алгоритма наблюдается незначительный

скачок на распределении коэффициента сопротивления. В целом отклонение от корреляции не превосходит 6%, а в среднем составляет 4%. На трехзонной сетке сходимость решения достигалась за ~5000 итераций, на однозонной за -400 итераций.

Вторая серия расчетов была выполнена для напорного течения при Яе=20х103 с вращением внутреннего цилиндра при различных вращательных числах Рейнольдса 11ет. Зависимость интегрального результата расчета - коэффициента сопротивления/- от числа Кеш показана на Рисунке 2.2. Там же на Рисунке 2.2 область переключения алгоритма отмечена вертикальной полосой. Слева от полосы при 11еш <1.1х105 работает алгоритм построения трехзонной сетки, справа от полосы при Яеы >1.1x105 строится однозонная сетка. Полученные распределения интегральных величин гладкие и не содержат разрывов или скачков в области переключения алгоритма. На трехзонной сетке сходимость решения достигалась за ~30000 итераций, на однозонной за ~10000 итераций.

В главе 4 рассматриваются вопросы моделирования течений неньютоновских сред, в общем случае описываемых реологией Гершеля-Балкли, для ламинарных и турбулентных режимов течения.

Тестирование представленного алгоритма было проведено на большом числе различных задач в широком диапазоне изменения параметров течения и геометрии кольцевого канала. Расчетные данные сопоставлялись с известными аналитическими решениями, экспериментальными данными и расчетами других авторов. Во всех случаях численные результаты хорошо согласовывались с тестовыми данными. В параграфе 4.1 приведены три примера такого сопоставления. Сопоставление расчетов и аналитических решений для течений степенной и бингамовской жидкостей в концентрическом кольцевом канале проводилось как по интегральным данным, в частности по продольному градиенту давления, так и по распределениям аксиальной скорости.

На Рисунке 3.1 аналитические решения (сплошные линии) сопоставлены с численными решениями для течений степенной среды в кольцевом канале с отношениями радиусов 0 = 0.5 и Яе=67. Использовалась достаточно грубая равномерная сетка (40х 120 узлов). С уменьшением показателя степени п градиент скорости в пристеночной области увеличивается. На Рисунке 3.2 численное решение сопоставлено с аналитическим решением для течения бингамовской жидкости в концентрическом канале с отношением радиусов 0 = 0.5 для трех значений числа Бингама: Вп =5, 50 и 250. Численное решение получено на той же сетке, что и для степенной среды. Согласование численного решения с аналитическим решением также везде хорошее. Однако стоит отметить, что точность расчетов при больших значениях числа Бингама сильно зависит от детализации расчетной сетки в области пограничного слоя, который может быть очень тонким. Для значения числа Бингама равного 104 толщина погранслоя составляет около 0,03 % от ширины кольцевого зазора. Для разрешения таких узких пристеночных слоев необходимо использовать неравномерную сетку со сгущением сеточных линий к стенке.

Рисунок 3 - Профили скорости для течения степенной (1) и бингамовской (2) жидкостей в концентрическом кольцевом канале в зависимости от безразмерного радиуса. Сплошные кривые - аналитическое решение, символы - расчет

Важным тестом является установившееся ламинарное течение жидкости Гершеля-Балкли в концентрическом кольцевом канале с вращением внутренней трубы. Численные результаты сопоставлены с экспериментальными данными. Расчеты выполнены для двух режимов со следующими критериями, задающими характер течения: число Рейнольдса Re = 0.527, число Тейлора Та = 2.58 и Re = 0.12 и Та = 0.0014. В обоих случаях наблюдается хорошее согласование расчетных и экспериментальных данных по профилям аксиальной и тангенциальной компонент скорости в выбранных сечениях.

Параграф 4.2 посвящен вопросам моделирования турбулентных режимов неньютоновских сред. Тестовые расчеты выполнены для установившихся течений степенной жидкости и жидкости Гершеля-Балкли в прямой круглой трубе. Данные, полученные с помощью предлагаемого численного алгоритма, сопоставлялись с результатами прямого численного моделирования (DNS).

Для всех режимов получено хорошее предсказание распределений сдвиговых турбулентных напряжений (Рисунок 4). С увеличением предельного напряжения или уменьшением показателя степени п наблюдается уменьшение интенсивности турбулентного переноса импульса. Расположение максимума напряжений немного смещается в сторону от стенки.

На Рисунке 5 сопоставлены данные распределений нормированной средней молекулярной вязкости fl*=/u/jum полученные с использованием предлагаемой RANS модели и прямого численного моделирования. Полученные распределения хорошо согласуются с данными DNS во всей области течения. Для всех данных наблюдается монотонный рост вязкости при отдалении от стенки. Непосредственно от стенки до значения у+- -10 средняя вязкость слабо меняется. В области максимальной генерации турбулентности (у+ -10) наблюдается резкий рост средней вязкости.

0.8 —I

т+

0.6 —

0.4 —

0.2 —

0

0.8 —

О 1 т+

if --- 6 Ii 'S

8 I

t

0.6 —

0.4

0.2 —

3'

й/7

J»

i.

TTTTTl-1 I I III

Лттт

0

** r?Tillll|-1 I I 11 lll|-1 I I 111111

1000

1 10 У+ 100 1000 1 10 У+ 100

1) 2) Рисунок 4 - Распределения сдвиговых турбулентных напряжений: 1) степенная среда (DNS: 1 - п=1, 2 - п=0.75, 3 - n=0.69, RANS: 4 - п=1, 5 -п=0.75, 6 -п=0.69), 2) жидкость Гершеля-Балкли (DNS: 1 - 11=1,^=0.0; 2 - n=0.6,^=0.0, 3 - n=0.6,Zo=0.1 rw, RANS: 4 - n= 1,5 - n=0.6, to=0.0, 6 - n=0.6,io=0.1 rw)

3

2.5 —

A A A. n=0.75 DNS □ □ 0(1=0.69 DNS

- n=0.75 RANS

-------n=0.69RANS

1.5 —

100

RANS DNS T0=0 A

io=0.1iw - - - • О

+ j t0=0.5Tw--

1

0 0.2 0.4 06 0.8 1 0 0.2 0.4 0.6 0.8 1

1-Г/Я Л-rlR

1) 2) Рисунок 5 - Распределения нормированных средних молекулярных вязкостей вдоль радиуса трубы: 1) для степенной среды, 2) для жидкости Гершеля-Балкли. Символы - DNS, линии - RANS модель

В параграфе 4.3 описываются результаты численного моделирования винтовых режимов течения неньютоновских сред с доминированием вращательного течения, индуцированного вращением внутреннего цилиндра. Данная задача интересна тем, что условия развития турбулентности в течении существенно меняются по сечению канала.

Особенности указанных течений рассмотрены на примере течений жидкости Бингама с параметрами А\= 0.005 Па с т0=0.3 Па в коаксиальном канале с

отношением радиусов 0.3 при Яе=200. Для течений сред с предельным напряжением область ламинарно-турбулентного перехода сдвигается в область больших вращательных чисел Яеш. В области ламинарно-турбулентного перехода наблюдается существенное снижение перепада давления вдоль канала (Рисунок 6). Уменьшение сопротивления объясняется качественной перестройкой течения, связанной с формированием области турбулентного течения вблизи вращающегося цилиндра. Для оценки характера течения в целом в работе вводится интегральный критерий степени турбулизации потока /г- Это число равно доли объема течения, в котором турбулентная вязкость больше осредненной молекулярной вязкости.

1000 г

д □

100

А □

X-

Л ЛАЛ А л / /'

□ □ □ □□ ЛД

-X— X х □ л / <

х □ V '

р/ л' .

1

'т

0.8

Рисунок 6 - Зависимость коэффициента сопротивления / и степени турбулизации потока/г от числа Яещ

0.005

А

\

\ х

¿-е оо ^о е-о-* ^ ^

----7-1-

X

Яей,103

- 20

30

-Х- 40

----60

-О- 80 - 100

ее

3

0.2

0.4 0.6 г/ДИ

1)

Рисунок 7 - Распределения в кольцевом зазоре: 1) аксиальной скорости нормированной на среднерасходную скорость, 2) кинетической энергии турбулентности нормированной на скорость вращения цилиндра. Левая граница соответствует стенке вращающегося цилиндра

На Рисунке 7 представлены профили аксиальной скорости в кольцевом канале для разных вращательных чисел Яеш. В ламинарном режиме течения наблюдается течение с движущимся квазитвердым ядром, занимающим большую часть канала, и узкими пристеночными слоями. Характер течения кардинально меняется при развитии турбулентности вблизи вращающегося внутреннего цилиндра, что наблюдается при числах Рейнольдса Яеш>20х103. Вблизи цилиндра эффективная вязкость достаточно мала для развития не-устойчивостей и переходу к турбулентности. Развитие турбулентных пульсаций приводит к интенсификации обмена количеством движения между пристеночным слоем и квазитвердым ядром потока и увеличивает размеры слоя с малой вязкостью. Происходит разделение течения в канале на область турбулентного движения, примыкающей к вращающемуся цилиндру, и область ламинарного течения. С увеличением числа Рейнольдса Яею размер кольца, занятого турбулентным течением, и аксиальный поток через турбулентное кольцо увеличиваются (Рисунок 7). Расчеты показывают, что при некотором значении 11ет (-40x103) весь аксиальный поток сосредотачивается в области турбулентного течения, а аксиальное движение в квазитвердой части практически отсутствует. При Кеш~100х103 весь канал занят турбулентным течением с относительно симметричным профилем аксиальной скорости.

В заключении сформулированы основные результаты данной работы:

1. На основе метода контрольного объема и методики расщепления разработан численный алгоритм решения уравнений гидродинамики неньютоновских сред. Предложена оригинальная процедура определения градиента давления для задач с периодическими граничными условиями с использованием интегральных балансных соотношений.

2. Сформулирована новая модель осредненной по Рейнольдсу молекулярной вязкости для турбулентных течений неньютоновских сред, подчиняющихся модели Гершеля-Балкли. В тестовых расчетах получено хорошее предсказание распределений осредненной вязкости как в пристеночном слое, так и в ядре течений. Предложенная модель в рамках метода осреднения по Рейнольдсу позволяет правильно описать уменьшение турбулентного переноса импульса с уменьшением показателя степени и увеличением предельного напряжения.

3. Разработан двухшаговый алгоритм построения конечно-объемной сетки для расчета полностью развитого турбулентного течения в кольцевом канале. Построение сетки, позволяющей уменьшить ошибки, возникающие при использовании автоматических пристеночных функций, основывается на результатах предварительного расчета.

4. Проведено большое количество тестовых расчетов для ламинарных и турбулентных режимов течения заданного класса жидкостей. Тестовые расчеты продемонстрировали хорошую точность и устойчивость вычислительного алгоритма в широком диапазоне параметров течения и характеристик рабочей среды.

5. На основе предложенных математических моделей и численных алгоритмов создан оригинальный инструментарий вычислительного эксперимента для моделирования течений неньютоновских сред в кольцевом канале.

6. Выполнено численное исследование винтовых режимов течения неньютоновских сред с доминированием вращательного течения. Обнаружено существование переходных режимов с разделением течения на область турбулентного движения, примыкающей к вращающемуся цилиндру, и область ламинарного течения. Такие режимы течения характеризуются локализацией аксиального потока в турбулентной области.

Публикации в периодических изданиях, рекомендованных ВАК:

1. Гаврилов, A.A. Применение нового численного алгоритма решения уравнений Навье-Стокса для моделирования работы вискозиметра типа физического маятника / A.A. Гаврилов, A.B. Минаков, A.A. Дектерев, В .Я. Рудяк // Теплофизика и аэромеханика. - 2008.- Т. 15, № 2. - С. 353-365.

2. Гаврилов, A.A. Численный алгоритм для моделирования ламинарных течений в кольцевом канале с эксцентриситетом / A.A. Гаврилов, A.B. Минаков, A.A. Дектерев, В.Я. Рудяк // Сиб. журн. индустр. матем. - 2010. - Т. 13 № 4 -С. 3-14.

3. Гаврилов, A.A. Численный алгоритм для моделирования установившихся ламинарных течений неньютоновских жидкостей в кольцевом зазоре с эксцентриситетом / A.A. Гаврилов, A.B. Минаков, A.A. Дектерев, В.Я. Рудяк // Вычислительные технологии. - 2012. - Т. 17, № 1. - С. 44-56.

4. Гаврилов, A.A. Моделирование коэффициента молекулярной вязкости вяз-копластичных жидкостей в турбулентных течениях / A.A. Гаврилов, В.Я. Рудяк // Доклады АН ВШ РФ. - 2013. - № 2 (21 ). - С. 55-66.

Публикации в других изданиях:

5. Gavrilov, A.A. Numerical simulation of fully developed turbulent non-Newtonian flow in eccentric annulus / A.V. Minakov, A.A. Gavrilov, A.A. Dekter-ev, V. Ya. Rudyak // Proceedings of Sixth International Symposium On Turbulence, Heat and Mass Transfer. - New York: Begell House Inc., 2009. - P. 333-335.

6. Гаврилов, A.A. Математическая модель и численная методика моделирования развитого турбулентного течения неньютоновских вязкопластических жидкостей / A.A. Гаврилов, A.B. Минаков, A.A. Дектерев, В.Я. Рудяк // Тезисы докладов международной конференции «Современные проблемы прикладной математики и механики: теория, эксперимент и практика», посвященная 90-летию со дня рождения академика H.H. Яненко. - 2011. - С. 85-86.

7. Гаврилов, A.A. Численное моделирование установившихся ламинарных течений неньютоновских вязкопластических жидкостей в кольцевом зазоре / A.A. Гаврилов, A.B. Минаков, A.A. Дектерев, В.Я. Рудяк // Сборник тезисов международной конференции «Математические и информационные технологии» MIT-2011. -2011. - С.78-79.

8. Гаврилов, A.A. Моделирование турбулентного течения неньютоновской жидкости на основе двухпараметрической модели турбулентности / A.A. Гаврилов, A.B. Минаков, A.A. Дектерев, В.Я. Рудяк // Тезисы докладов 3-го Все-

российского семинара «Фундаментальные основы МЭМС- и нанотехноло-гий».— 2011. - Вып. 3. - С. 45-47.

9. Гаврилов, A.A. Современные возможности CFD кода SigmaFlow для решения теплофизических задач / A.A. Гаврилов, A.B. Минаков, A.A. Дектерев // Сборник научных статей. Современная наука: исследования, идеи, результаты, технологии. - Киев: «НПВК Триакон». - 2010. - №2(4). - С. 117-122.

10. Гаврилов, A.A. Моделирование задач гидродинамики, теплообмена и горения с использованием CFD программы SigmaFlow / A.A. Гаврилов, A.B. Минаков, A.A. Дектерев, Е.Б. Харламов // Сборник тезисов VI Всероссийского семинара по теплофизике и теплоэнергетике. - 2009. - С. 24.

11. Гаврилов, A.A. Прямое численное моделирование турбулентного течения степенной жидкости в круглой трубе / A.A. Гаврилов, В.Я. Рудяк, Е.В. Подрябинкин // 4 Всероссийский семинар «Фундаментальные основы МЭМС- и нанотехнологий». - 2012. - Вып. 4. - С. 113-116.

Подписано к печати 3 декабря 2014 г. Заказ № 29 Формат 60/84/16. Объем 1 уч.-изд. лист. Тираж 100 экз.

Отпечатано в ФГБУН Институте теплофизики СО РАН 630090, Новосибирск, пр. Академика Лаврентьева, 1