Содержание к диссертации
Введение
1 Моделирование течений жидкости со свободной поверхностью 16
1.1 Численные методы решения задач о течении невязкой жидкости со свободной поверхностью 16
1.2 Численные методы решения задач о течении вязкой жидкости со свободной поверхностью 20
2. Колебания капли невязкой жидкости под действием сил поверхностного натяжения 25
2.1 Введение 25
2.2 Постановка задачи 29
2.3 Методы решения 32
2.3.1 Метод граничных элементов 32
Плоский случай 32
Осесимметричный случай 36
2.3.2 Метод конечных разностей 45
2.3.3 Эволюционные алгоритмы расчета движения свободной поверхности 49
2.4 Результаты исследования 55
3 Моделирование вязкого течения со свободной поверхностью внутри вращающегося горизонтального цилиндра 74
3.1 Введение 74
3.2 Постановка задачи 79
3.3 Метод расчета 83
3.4 Результаты 89
Заключение 115
Литература 118
- Численные методы решения задач о течении вязкой жидкости со свободной поверхностью
- Постановка задачи
- Эволюционные алгоритмы расчета движения свободной поверхности
- Метод расчета
Введение к работе
Диссертационная работа посвящена решению фундаментальных и прикладных задач о течениях идеальной и вязкой несжимаемых жидкостей со свободной поверхностью в плоской и осесимметричной постановках. В качестве инструмента исследований применяется метод граничных элементов.
Класс задач о течении жидкости со свободной поверхностью представляет большой интерес для самых различных областей науки и техники, будь то химическая технология, гидрометеорология, охрана окружающей среды или аэрокосмические исследования. Наличие свободной поверхности в области течения является характерной особенностью таких процессов, как распыление аэрозолей, нанесение покрытий, литье и многих других. Данное обстоятельство обуславливает неослабевающий интерес исследователей как к развитию методов решения, так и к более тщательному изучению конкретных задач, имеющих теоретическое и практическое значение.
В связи с нерегулярностью границ областей, характерной для исследуемых процессов, при количественном анализе трудно рассчитывать на получение аналитических результатов, и, как правило, решения большинства практических задач приходится искать с использованием численных методов. Многообразие задач, связанных с исследованием течений жидкости со свободной поверхностью, породило значительное количество численных методик, учитывающих особенности той или иной рассматриваемой проблемы. Зачастую при создании алгоритмов и их численной реализации исследователи используют специальные приемы, связанные с отслеживанием эволюции свободной границы и выполнением граничных условий на ней. Также остро стоит вопрос о повышении точности расчетов в окрестности линий трехфазного контакта. Выбор того или иного подхода к решению определяется спецификой задачи и особенностями методов, имеющихся в арсенале исследователя. Эффективность метода
5 напрямую влияет на возможность проведения точного численного эксперимента в заданном диапазоне изменения определяющих параметров.
В настоящей работе исследуются две задачи о течении несжимаемой жидкости со свободной поверхностью: 1) смоделирован процесс колебаний капли идеальной жидкости под действием сил поверхностного натяжения в невесомости; 2) проведено исследование течения высоковязкой жидкости в частично заполненном горизонтальном цилиндре, вращающемся вокруг собственной оси с постоянной скоростью.
Детальное изучение поведения жидкостей в капельном состоянии интересует исследователей на протяжении вот уже более полутора сотен лет. В первую очередь, этот класс задач интересен в связи со своей непосредственной практической значимостью, так как жидкости в капельном состоянии встречаются во многих природных и технологических процессах. Возрождение интереса к изучению колебаний капель в безгравитационной среде в последние два десятилетия связано с применением космических технологий, включающих бесконтейнерное производство материалов [1]. В геологии и кристаллографии исследование процессов разрыва капель и газожидких пузырьков в минералах позволяет судить об истории возникновения минералов. Возможность предсказания поведения капель важна в таких явлениях как распыление аэрозолей [2], нанесение покрытий методом напыления, гетерогенное горение, технологии бесконтейнерного производства материалов в космосе, клеточное деление в биологических системах, взаимодействие поверхностей радаров с влагой дождевых облаков [3], а также непрямое измерение реологических параметров жидкостей.
Существенную роль в изучении поведения жидкостей в капельном состоянии в том или ином процессе играют, конечно, натурные эксперименты и эксперименты, проводимые на модельных установках. Недавно было предложено оптическое измерение частоты колебаний как неагрессивный метод оценки динамики поверхностного натяжения и, возможно, прочих физических констант, связанных с поверхностной
6 реологией [4-6]. Получение аналитических результатов для большинства рассматриваемых процессов затруднено сложностью используемого математического аппарата и нерегулярностью границ областей, поэтому важнейшим инструментом изучения здесь являются численные методы. Применительно к гидродинамике использование численного эксперимента часто оказывается наиболее информативным и дешевым инструментом исследования, позволяющим получать необходимые сведения о процессе в широком диапазоне изменения входных параметров.
В окружающей среде с нулевой гравитацией масса жидкости самопроизвольно уменьшает площадь своей поверхности, преобразуясь в сферическую каплю, которая может колебаться в различных модах. Природа колебаний определяется начальными условиями, а амплитуда движения определяется энергией, сообщенной капле в процессе формирования.
Малые колебания невязких сферических капель, помещенных в другую невязкую окружающую жидкость, были впервые изучены Рэлеем [7] и Ламбом [8] в контексте линейной теории и много позднее их исследовали Tsamopoulos и Brown [9] в контексте слабо нелинейной теории. В настоящее время известно множество работ, посвященных анализу указанного процесса [10-37]. Однако в большинстве работ исследователи ограничиваются рассмотрением малых деформаций и, соответственно, колебаний капель. Предыдущие попытки моделирования процесса колебаний капель показали, что для проведения длительных по времени расчетов необходима высокая точность, более того, также необходимо применение какой-либо сглаживающей процедуры (или внесение другой значительной поправки в расчет), с тем, чтобы предотвратить рост вычислительной неустойчивости, причина появления которой неизвестна. Другим важным фактором является точность вычисления кривизны свободной поверхности, которая является основным параметром исследуемого процесса. Кроме того, значительное влияние как на устойчивость, так и на получаемые результаты оказывает и
7 выбранный исследователем алгоритм расчета эволюции поверхности во времени.
Таким образом, в рассматриваемой задаче, представляет интерес не только изучение течения, возникающего при колебаниях капель и эволюции свободной поверхности, но также рассмотрение вычислительных особенностей, возникающих при численном моделировании процесса и, в конечном итоге, разработка эффективного вычислительного алгоритма расчета конкретных течений.
Другой задачей, исследуемой в рамках диссертации, является математическое моделирование двумерного вязкого течения со свободной поверхностью внутри вращающегося горизонтального цилиндра.
Исследование вязкого течения со свободной поверхностью в горизонтальном вращающемся цилиндре представляет существенный прикладной интерес для процессов химической технологии. Характер движения жидкости определяет эффективность реализации таких технологических процессов, как смешение, центробежное литье, нанесение покрытий и других.
Математическое моделирование процесса перемешивания полимерных композиций в объемном смесителе осложняется пространственным характером процесса перемешивания, наличием свободной поверхности, а также изменением реологических свойств среды во времени. Именно поэтому до настоящего времени нет работ, использующих полную математическую постановку задачи исследования процесса смешения. Вместе с тем, актуальной является задача разработки инженерных методик расчета гидродинамических процессов, реализуемых в технологии переработки высоконаполненных полимерных композиций, позволяющих прогнозировать режимы перемешивания с целью правильной организации технологического процесса изготовления изделий с заданными эксплуатационными характеристиками. Процесс смешения с использованием объемных смесителей является одним из этапов переработки полимеров
8 методом свободного литья. Стоит отметить, что число работ в этом направлении ограничено.
Данное обстоятельство, по-видимому, связано с необходимостью рассматривать течения при наличии свободных поверхностей, меняющихся со временем, что существенно усложняет их математическое описание и, как следствие, возможность получения количественных зависимостей. Первые попытки математического моделирования гидродинамических процессов переработки полимеров методом свободного литья с учетом эволюции свободных поверхностей представлены в работах [38-39].
В настоящем исследовании рассматривается установившееся течение вязкой несжимаемой жидкости в частично заполненном вращающемся горизонтальном цилиндре в приближении ползущего движения. В этом случае основным безразмерным критерием является величина, равная отношению числа Рейнольдса к числу Фруда и характеризующая соотношение вязких и гравитационных сил в потоке жидкости. В работе [40] показано, что в зависимости от значения этого критерия можно выделить два режима течения. При первом режиме, когда значение определяющего критерия мало, жидкость покрывает внутреннюю поверхность цилиндра слоем постоянной толщины, т.е. осуществляется течение, близкое к квазитвердому движению. В случае второго режима, который реализуется при достаточно больших значениях основного критерия, практически вся жидкость находится в нижней части цилиндра, и течение разделяется на зону циркуляционного движения внизу и пленочного течения в остальной части потока.
В существующих работах, использующих приближенные методы, в основном рассматривается случай, когда форма свободной поверхности незначительно отличается от окружности [41-46]. В [47-48] рассматривается течение, при котором форма свободной поверхности близка к горизонтальной, и реализуется течение с циркуляционной зоной. Эти работы выполнены при значительных упрощениях исходных дифференциальных
9 уравнений, сводящих их к обыкновенным. Имеются также попытки использования метода конечных разностей для решения данной задачи [49]. Однако в большинстве работ представлены лишь кинематические характеристики течения. Вопрос об исследовании характера течения, формируемого на начальном этапе процесса смешения, ставился лишь в работе [49].
Более полный обзор численных и экспериментальных работ, посвященных исследованию рассматриваемого процесса, приводится в третье главе настоящей диссертации.
Таким образом, в рассматриваемой задаче актуальными вопросами являются построение устойчивых вычислительных методик расчета режимов течения, позволяющих также определять необходимые кинематические и динамические характеристики внутри области течения.
Для решения плоских задач гидродинамики течений со свободной поверхностью можно использовать различные конечно-разностные методы [39]. При рассмотрении течений в областях сложной геометрии, с большими деформациями свободной границы, при наличии нескольких, значительно различающихся характерных размеров, возникают вычислительные трудности, связанные с построением конечно-разностных сеток, которые, кроме того, должны адаптироваться к изменениям свободной поверхности. В связи с этим актуальность приобретают поиски подходов, позволяющих упростить и унифицировать алгоритмическую процедуру численного решения. Одним из путей преодоления этих трудностей является переход к системе граничных интегральных уравнений, эквивалентной исходной системе дифференциальных уравнений. В этом направлении большие успехи достигнуты в механике деформируемого твердого тела [50], опираясь на которые далее строится в общем виде непрямой вариант метода граничных элементов для задачи о плоском течении вязкой жидкости со свободной поверхностью.
В работе большое внимание уделяется построению вычислительных процедур, включающих численную реализацию гранично-интегральных уравнений, в том числе особенностям построения процедуры в каждой конкретной задаче. К задачам настоящего исследования также относится разработка эффективных численных алгоритмов решения поставленных задач на основе метода граничных элементов.
В связи с вышеизложенным, целью диссертационной работы является
Разработка устойчивых вычислительных методик расчета плоских и осесимметричных течений несжимаемой жидкости со свободной поверхностью на основе непрямого метода граничных элементов.
Исследование процесса колебаний капли невязкой несжимаемой жидкости под действием сил поверхностного натяжения в невесомости.
Моделирование течения высоковязкой жидкости в частично заполненном вращающемся горизонтальном цилиндре, проведение параметрических исследований процесса.
Научная новизна работы заключается в следующем: > Разработан и протестирован алгоритм реализации метода граничных элементов для осесимметричного случая; апробированы различные численные алгоритмы для расчета эволюции свободной поверхности, проведен их сравнительный анализ.
Проведенное исследование колебаний капель идеальной жидкости позволило выявить влияние начальной деформации формы свободной поверхности на характеристики колебательного процесса.
Сформулирована и реализована математическая постановка, позволяющая в рамках модели ползущего движения исследовать двумерное вязкое течение со свободной поверхностью внутри вращающегося горизонтального цилиндра. Разработана универсальная вычислительная методика на основе метода граничных элементов, использование которой
11 позволило выявить установившиеся режимы течения высоковязкой жидкости внутри вращающегося горизонтального цилиндра; получить в широком диапазоне изменения определяющих параметров распределения кинематических и динамических характеристик внутри рассматриваемой области; исследовать эволюцию свободной границы.
Практическая значимость.
Практическая ценность работ, посвященных исследованию процесса колебания капель, обусловлена широкими возможностями использования полученных результатов [51] применительно к технологии спекания в порошковой металлургии, метеорологии, теории двухфазных потоков.
Исследуемое в работе течение вязкой жидкости в частично заполненном вращающемся горизонтальном цилиндре рассматривается как модель гидродинамического процесса, реализуемого в объемном (гравитационном) смесителе на стадии перемешивания. На основе метода граничных элементов разработана методика и создана программа расчета, проведены параметрические исследования. Выявлены режимы течения с образованием циркуляционной зоны и с формированием сдвигового течения по всей поверхности. Кроме того, выявлен режим течения, при котором часть поверхности цилиндра не смачивается жидкостью. Определены критические значения основных параметров (соотношения гравитационных и вязких сил в потоке W и степени заполнения цилиндра жидкостью X), разделяющие разные режимы течения. Впервые для данной задачи получены поля скоростей, распределения второго инварианта тензора напряжений, сдвиговых и сжимающих напряжений в области, занятой жидкостью, в широком диапазоне изменения определяющих параметров. Распределение сдвиговых и сжимающих напряжений позволяет оценить силовое воздействие на агломераты порошкообразных компонентов, а второго инварианта тензора напряжений - выделить зоны с максимальными
12 значениями диссипативной функции и наибольшей интенсивностью смешения.
В качестве характеристики эффективности перемешивания используется величина объемной мощности потока, равная произведению скорости сдвига и напряжения сдвига. Получены зависимости объемной мощности от определяющих параметров, установлен экстремальный характер этих зависимостей, что впервые показано расчетным путем и также подтверждено экспериментально. Максимуму удельной мощности отвечает оптимальная скорость вращения ротора смесителя. Определены диапазоны изменения определяющих параметров, для которых при определении характеристик течения можно использовать приближенное решение [48].
Исследование эволюции свободной границы на начальном этапе вращения цилиндра, из положения горизонтальной свободной границы к установившейся позволило определить четыре варианта поведения свободной поверхности.
В результате физического моделирования подтверждено наличие режимов течения, выявленных расчетным путем [52].
Разработанные методики и пакеты прикладных программ проходили апробацию в федеральном центре двойных технологий «Союз». Получено три акта внедрения программ расчета на предприятия ФГУП «ФЦДТ «Союз».
Достоверность и обоснованность полученных результатов обеспечивается корректностью математических постановок задач, внутренними проверками разработанных методик и программ расчета (проверка аппроксимационной сходимости, выполнение законов сохранения и пр.), а также согласованием с известными экспериментальными данными, а также с результатами, полученными другими авторами.
Апробация работы. Материалы диссертации по мере их получения докладывались и обсуждались на международной конференции «Моделирование процессов в синергетических системах» (Улан-Удэ, 2002);
13 на VIII Всероссийской научно-технической конференции «Механика летательных аппаратов и новые материалы» (Томск, 2003г); на IX Всероссийской научной конференции студентов-физиков и молодых ученых (Красноярск, 2003); на Всероссийской научной конференции молодых ученых (НГТУ, Новосибирск, 2003); на X всероссийской научно-технической конференции «Физика и химия высокоэнергетических систем» (Томск, 2004); на научной сессии молодых ученых научно-образовательного центра «Физика и химия высокоэнергетических систем» (Томск, 2004); на IV и V Всероссийских научных конференциях «Фундаментальные и прикладные проблемы современной механики» (Томск, 2004, 2006); на I, II и III Всероссийских конференциях молодых ученых «Физика и химия высокоэнергетических систем». (Томск, 2005, 2006, 2007); на VI Международной конференции «Лаврентьевские чтения по математике, механике и физике» (ИГиЛ СО РАН, Новосибирск, 2005); на Международной школе-конференции молодых ученых «Физика и химия наноматериалов» (Томск, 2005); на VI Всероссийском съезде по теоретической и прикладной механике (Ниж. Новгород, 2006); на XXII научной конференции стран СНГ «Дисперсные системы» (Украина, Одесса, 2006); XLV Международная научная студенческая конференция «Студент и научно-технический прогресс» (НГУ, Новосибирск, 2007); на XIV Международном симпозиуме «Оптика атмосферы и океана. Физика атмосферы» (Бурятия, 2007) и на Всероссийской конференции «Проблемы механики сплошных сред и физики взрыва» (ИГиЛ СО РАН, Новосибирск, 2007).
Публикации. Основные результаты диссертации представлены в трудах вышеперечисленных конференций, а также в журналах «Теоретические основы химической технологии», «Theoretical Foundations of Chemical Engineering», «Вычислительные технологии», «Известия вузов. Физика», «Оптика атмосферы и океана». Всего по материалам диссертации
14 опубликовано 27 работ [53-79], из них 21 в соавторстве с профессорами Г.Р. Шрагером и В.А. Якутенком.
Краткое содержание работы по главам
Во введении раскрывается актуальность и практическая значимость математического моделирования течений несжимаемой жидкости со свободными границами. Сформулированы цель и основные задачи исследований.
В первой главе приведен обзор современных численных методов расчета течений несжимаемой идеальной и вязкой жидкости со свободными границами.
В второй главе в плоской и осесимметричной постановках исследуется процесс нелинейных колебаний капли идеальной жидкости под действием сил поверхностного натяжения в отсутствии гравитации. Разработана методика использования метода граничных элементов применительно к моделированию осесимметричных течений невязкой несжимаемой жидкости со свободными границами. Приведен сравнительный анализ результатов использования конечно-разностного и гранично-элементного подходов к решению, в осесимметричном случае протестированы различные алгоритмы отслеживания эволюции свободной границы капли. Полученные результаты сравниваются с известными данными.
В третьей главе в рамках численного моделирования с использованием модели ползущего течения, исследовано поведение объема вязкой жидкости, заполняющего часть горизонтального цилиндра, вращающегося с постоянной угловой скоростью. Изучен характер эволюции свободной поверхности жидкости на начальном этапе вращения из начального положения горизонтальной свободной границы. Выявлена возможность реализации четырех различных вариантов трансформации свободной поверхности. Показано, что установившееся течение
15 характеризуется двумя режимами движения: при первом режиме, который реализуется при быстром вращении цилиндра, жидкость в цилиндре образует малоподвижный пристеночный квазикольцевой слой; второй режим, возникающий при относительно медленном вращении цилиндра, характеризуется наличием зон циркуляционного течения в нижней части цилиндра и тонкого пристеночного слоя в верхней его части. Наличие указанных режимов течения подтверждается экспериментальными данными. На основе метода граничных элементов разработана методика расчета кинематических и динамических характеристик течения, оценки величины удельной мощности перемешивания.
В заключении подведены основные итоги проведенных исследований, намечены пути дальнейшего развития работы.
Основные положения, выносимые на защиту:
К* Алгоритм реализации метода граничных элементов для расчета осесимметричных и плоских течений идеальной несжимаемой жидкости со свободной границей;
Сравнительный анализ численных методик для расчета эволюции свободной поверхности;
Результаты численного анализа нелинейных колебаний капель идеальной несжимаемой жидкости под действием сил поверхностного натяжения;
Алгоритм решения задачи о течении вязкой жидкости, частично заполняющей вращающийся горизонтальный цилиндр;
Режимы течения, полученные при исследовании стационарных и квазистационарных гидродинамических процессов во вращающемся цилиндре;
Результаты расчетов кинематических, динамических характеристик и зависимости величины удельной мощности потока от определяющих параметров течения во вращающемся цилиндре.
Численные методы решения задач о течении вязкой жидкости со свободной поверхностью
Первым [107] и по настоящее время наиболее распространенным численным методом, применяемым для численного решения задач гидродинамики вязкой жидкости со свободными границами, является метод конечных разностей (МКР). В [107] описан оригинальный метод расчета плоских течений вязкой жидкости со свободной границей (ныне известный как МАС-метод), использующий явную разностную аппроксимацию уравнений Навье-Стокса на прямоугольной эйлеровой сетке. Специальный выбор сетки позволяет выделить элементарные ячейки, для которых тождественно выполняется разностный аналог уравнения неразрывности. В центре ячейки рассчитывается давление, а на ее сторонах - скорости. Для определения положения свободной границы с течением времени на эйлеровой сетке располагаются частицы-маркеры, траектории движения которых следуют течению жидкости. В данном методе свободная поверхность описывается ступенчатой функцией, что приводит к потере ориентации границы, и, как следствие, невозможности точной аппроксимации граничных условий на ней. Однако при исследовании течений жидкостей, обладающих незначительной вязкостью, вклад вязких сил в граничное условие на свободной поверхности невелик, и МАС-метод позволяет получать достаточно точные результаты.
Дальнейшие многочисленные модификации МАС-метода [108- 116] связаны с применением более точных схем аппроксимации конвективных членов, более точном выполнении граничного условия для давления на свободной поверхности. Модификация МАС-метода для расчета высоковязких течений со свободными границами предложена в [ПО]. В [117-120] представлено применение данной модификации при решении таких задач, как заполнение осесимметричных каналов, растекание столбика жидкости по горизонтальной плоскости. Основным достоинством МАС-метода является выполнение разностного аналога условия неразрывности за счет использования разнесенной сетки для расчета давления и скоростей. В то же время, это несет дополнительные трудности при расчете данных характеристик на свободной границе. К недостаткам относится также трудности построения разностной сетки для областей со сложной геометрией.
В работах [121 - 124] предложен численный метод расчета течений вязкой жидкости со свободной поверхностью, который позволяет аппроксимировать естественные граничные условия на свободной поверхности с сохранением её ориентации. С помощью данного метода стало возможным удовлетворить условия на границах раздела [125, 126], учесть капиллярные эффекты [127], а также неньютоновское поведение жидкости [128 - 131]. Большинство задач рассмотрено в осесимметричной постановке, обобщение метода на пространственный случай выполнено в [132 - 133].
Наряду с методами эйлерова типа, рассмотренными выше, существуют методы лагранжева типа [134], которые базируются на лагранжевых координатах, а расчетные ячейки перемещаются вместе с жидкостью. Основное преимущество такого подхода состоит в способности точно отслеживать движение свободной границы или поверхности раздела, однако эти методы пригодны лишь для расчета течений, в которых ячейки не претерпевают сильные искажения. Наряду с перечисленными методами существуют смешанные эйлерово-лагранжевые методы - метод ALE [135], метод VOF [136]. Метод объема жидкости (метод VOF) основан на рассмотрении частичного объема жидкости в ячейках расчетной сетки. Его особенностью также является введение некоторой функции F, значение которой принимается за единицу, в ячейках, полностью заполненных жидкостью и нулю, где жидкость отсутствует. Таким образом, ячейки, содержащие свободную поверхность, имеют 0 F 1. Направление быстрейшего изменения функции F дает направление нормали к свободной поверхности. В настоящее время развитие МКР для решения уравнений Навье-Стокса идет в основном по пути повышения порядка аппроксимации [137- 140].
Помимо конечно-разностных методов последнее время всё большее распространение в задачах моделирования течений вязких жидкостей со свободными границами получают методы конечных и граничных элементов. Широкое их распространение в значительной степени объясняется гибкостью в отношении геометрии исследуемой области и краевых условий задачи, а также развитием вычислительной техники.
Одна из первых работ, использующих МКЭ, опубликована в 1965 году [141]. В ней рассмотрено течение жидкости со свободной поверхностью между двумя параллельными пластинами. В работе делается вывод о применимости МКЭ для численного исследования широкого класса задач, связанных с течениями нелинейно-вязких жидкостей со свободными границами. Дальнейшее развитие метода привело к решению целого ряда практически важных задач: исследованы задачи о зарождении волн на поверхности и водосливе [142-144], задачи струйного течения с учетом упругих и капиллярных эффектов [145, 146], процессы заполнения каналов сложной формы [147]. В [43] исследовано течение вязкой жидкости во вращающемся горизонтальном цилиндре; в [148, 149] - растекание столба жидкости под действием силы тяжести. В работах [150, 151] для моделирования течений расплавов полимеров с учетов неизотермичности и нелинейно-вязкопластичных свойств разработан вариант МКЭ, основанный на представлении нелинейной задачи в виде ряда сходящихся линейных задач о течении ньютоновской жидкости с соответствующей эффективной вязкостью. Подход, основанный на применении метода Галеркина, реализован в цикле работ [152-154]. Эти исследования имеют большую методическую ценность, так как в них рассмотрено влияние на эффективность численных процедур вида конечного элемента, способов построения конечно-элементных сеток (что, вообще говоря, выделяется в отдельную проблему при изучении областей, имеющих сильно отличающиеся масштабы), методов решения получающихся систем нелинейных уравнений. Большое внимание уделено постановке граничных условий на твердой стенке и линии трехфазного контакта. И в настоящее время МКЭ не только не утратил своей актуальности, но продолжает широко использоваться при решении задач о течении вязкой жидкости со свободными границами [155 - 158].
В настоящее время большинство работ по МГЭ относится к задачам линейной теории упругости, что объясняется, прежде всего, тем, что сама история развития линейной теории упругости тесно связана с теорией граничных интегральных уравнений. Как уже отмечалось выше, МГЭ в гидродинамике используется в основном при исследовании потенциальных течений идеальной жидкости и для решения уравнений Стокса, описывающих ползущее течение вязкой жидкости. Например, течение жидкости при заполнении полостей различной конфигурации, растекание жидкости на плоскости [159- 160]. В [40, 161, 162] МГЭ используется для исследования медленных течений вязкой жидкости со свободной поверхностью. Рассмотрены задачи об определении стационарной формы свободной поверхности жидкости, находящейся внутри вращающегося горизонтального цилиндра, о заполнении вертикального канала, об истечении из клиновидных и прямоугольных емкостей. Данный подход основан на применении фундаментальных решений линеаризованных уравнений Навье-Стокса и использует аналитические способы интегрирования последних. Обобщение данного подхода на трехмерный случай приведено в [163]. Задачи о пространственной экструзии неньютоновских жидкостей рассматриваются в работах [164 - 166]. В [167, 168] с использованием МГЭ исследуется заполнение плоского канала жидкостью.
Постановка задачи
Рассматривается процесс нелинейных колебаний капли идеальной несжимаемой жидкости под действием сил поверхностного натяжения в отсутствии силы тяжести. Тем самым, в плоском приближении моделируется процесс колебаний бесконечного цилиндра идеальной жидкости, который в начальный момент времени покоится и имеет в поперечном сечении форму эллипса с полуосями a, b (рис. 2.1). В осесимметричном случае моделируется процесс колебаний капли, которая в начальный момент времени покоится и имеет форму вытянутого эллипсоида вращения с полуосями а, Ь, с = а (рис. 2.2).
Когда гравитация незначительна, капля жидкости стремится принять сферическую форму, с тем чтобы уменьшить поверхностную энергию. Такая капля начинает колебаться, и в отсутствие окружающей среды единственным источником затухания колебаний являются слабовязкие эффекты, генерируемые в тонком слое вдоль поверхности [10].
Пионерскими работами в изучении колебаний капель являются труды Рэлея ([7], 1879) и Ламба ([8], 1932), которые позволили получить линеаризованное решение для малых отклонений капли от сферической формы. В [7] с помощью полиномов Лежандра Pn, п = 2,3,4... были выражены основные моды колебаний и вычислены соответствующие частоты. В работе [8] был выяснен эффект малой вязкости, и позднее в [11] решено уравнение, описывающее эффект произвольной вязкости для бесконечно малых колебаний жидкости. Параллельно PI позднее аналитические исследования были дополнены малоамплитудной теорией и были рассчитаны эффекты вязкости и поверхностной реологии, как описано, например, в [5, 21]. Теоретические исследования дополнялись лабораторными наблюдениями, проводимыми в условиях микрогравитации или под влиянием акустической радиации [6, 13].
Детальное численное исследование бесконечно малых колебаний капли с произвольной вязкостью было проведено в [12]. Малые и умеренные колебания невязких капель рассматривались в [9] в контексте слабо нелинейной теории, расширив исследование [7] для начальных возмущений гармоник п = 2,3,4. Было показано, в числе прочего, что нелинейность "слабая", и что частоты первых четырех мод уменьшаются с возрастанием амплитуды. Количественно это подтверждено для второй моды, в том числе и экспериментально. Численно данная задача решалась и другими авторами (Foote, 1973 [14]; Alonso, 1974 [15] - МАС-методом, Benner, 1983 [17] - МКЭ), однако до 1988 гг. не существовало исследования колебаний большой амплитуды, впервые изученных в [10].
Для изучения колебаний большой амплитуды несколькими авторами выполнено численное моделирование в осесимметричном случае. Ранние исследования моделировали деформацию вязких капель, используя МАС-метод [14-16]. Durr и Siekmann [23] смоделировали колебания невязких капель умеренной амплитуды, используя представление о гармоническом потенциале простого слоя. Lundgren и Mansour [10] смоделировали колебания большой амплитуды, используя обобщенный метод вихря, развитый в работах Baker а и соавторов [24 - 26], и рассчитали эффекты малой вязкости, линеаризовав основные уравнения внутри пограничных слоев вдоль свободной поверхности. Их результаты показали, что сферическая капля может разрушиться, когда сообщенная ей энергия достаточно велика. Pelekasis и др. [21] провели всесторонний анализ трех альтернативных формулировок: эйлеровой формулировки в обозначениях канонических переменных, гранично-интегральной формулировки, основанной на третьей формуле Грина и формулировки, основанной на обобщенном методе вихря.
Другими авторами изучены двумерные и осесимметричные колебания невязких и вязких капель с использованием конечно-элементных и гранично-элементных методик [5, 18, 27, 28]. Mashayek и др. [29] рассчитали осесимметричные колебания вязких капель при наличии внутренней циркуляции, образуемой приложенной единичной тангенциальной скоростью вдоль поверхности капли. Apfel и соавторы исследовали осесимметричные колебания вязких капель, взяв в рассмотрение эффекты реологии поверхности раздела фаз [6, 30, 31]. Basaran и соавторы изучили процесс образования [32, 33] и формы колебаний подвешенных осесимметричных капель [34, 35]. Basaran и Patzek [37] применили конечно-элементный метод для моделирования трехмерных колебаний невязких капель, но, к сожалению, их работа доступна лишь в форме тезисов. Насколько известно, в упомянутой работе авторами не развит какой либо специальный метод моделирования пространственных течений невязкой жидкости.
В работе [19] представлен численный метод моделирования трехмерных колебаний невязкой капли, базирующийся на обобщенном методе вихря, раработанном Вакег ом и соавторами [24 - 26], и использовавшийся для изучения плоских и осесимметричных течений. Представленный в работе алгоритм успешно применим для исследования малых и умеренных деформаций капель.
Таким образом, в настоящее время накоплено большое число работ, посвященных анализу указанного процесса. Однако до сих пор исследователям не удалось смоделировать распад капли под действием поверхностного натяжения (т.е. собственной деформации), ни даже определить, какие значения деформаций являются критическими для капель, приводящими их к распаду. В большинстве работ (например, [10], [18], [19]) авторы ограничиваются рассмотрением малых колебаний капель в плоском [18], осесимметричном [10], [19] и трехмерном случаях [19]. Это связано с явлением численной неустойчивости применяемых схем. Как отмечается практически всеми исследователями в данной области, получение точных результатов возможно лишь при малой деформации поверхности. Поэтому в каждой из указанных работ авторам приходится прибегать к различного рода ухищрениям с целью улучшения устойчивости алгоритма, таким как перераспределение узлов вдоль свободной поверхности на каждом шаге по времени [18], [19], "сглаживание" функции потенциала вдоль свободной поверхности через определенное число шагов по времени с помощью сплайн-интерполяции [18] или с использованием полиномов Лежандра [19, 20], наложение дополнительных условий на значения потенциала в точках свободной поверхности [10] и другим. Естественно, возникает вопрос, каким образом подобные вмешательства в численную схему искажают физическую сущность процесса. Наиболее распространенной "корректировкой" алгоритма является внесение поправки в значение потенциала на свободной границе [10] или в динамическое граничное условие. Это эквивалентно добавлению некоторого диффузионного слагаемого в последнее уравнение, что, несомненно, стабилизирует схему, приводя, однако к некоторому затуханию колебаний. В [18] было показано, что указанная добавка в динамическое граничное условие может рассматриваться как учет вязкости исследуемой жидкости, тем самым, снимая вопрос о физическом соответствии и правомерности указанных действий.
Кроме того, значительное влияние как на устойчивость, так и на результат оказывает и выбранный исследователем алгоритм эволюции поверхности во времени. Этот вопрос впервые был детально рассмотрен в [171], где было показано в частности, что использование метода Рунге-Кутта третьего и четвертого порядков, являющихся условно устойчивыми по своим свойствам, приводит к слабой диссипации схемы, вследствие чего колебания идеальной капли, смоделированные с их использованием, окажутся слабо затухающими.
Эволюционные алгоритмы расчета движения свободной поверхности
Для наблюдения изменения профиля свободной поверхности вдоль свободной границы располагается конечное число частиц, движение которых осуществляется в соответствии с эйлеровым (плоская и осесимметричная задачи) или лагранжевым (осесимметричная задача) представлением кинематического условия. В настоящей работе были исследованы следующие виды разностного представления кинематического условия на свободной границе: в лагранжевой форме со вторым порядком, в эйлеровой -с первым порядком, схема «предиктор-корректор» второго порядка и схема Рунге-Кутта четвертого порядка. Для реализации любого из представленных алгоритмов прежде всего необходимо решить краевую задачу (2.1), (2.4), (2.6), т.е. определить значения компонент вектора скорости на свободной поверхности. Ниже представлены основные этапы используемых алгоритмов.
После решения исходной задачи (2.1), (2.4), (2.6) непрямым вариантом МГЭ или методом конечных разностей становятся известными нормальные компоненты скорости дер/дп во всех точках свободной поверхности. Касательные компоненты dep/ds вычисляются с помощью процедуры численного дифференцирования [177] 5ф/а5 = ар,ах= ap/fr (268) dt дерdXj дер"ай" п,- Зф „ п2,5s 2 2 dr dt" дердх2 дер"ай" "п2 дер ds ] где х номер узла, % рассматривается как непрерывная вещественная функция. Для вычисления производных по % используются приближенные пятиточечные формулы центральных и односторонних разностей. Далее могут быть вычислены компоненты скорости v]5 v2 на свободной поверхности (2.69) (2.70) Лагранжевы производные второго порядка Рассмотрим теперь производные второго порядка, входящие в (2.65)-(2.67): = di = avL + VSv,+v2Svi = 9L + Vev,+v2SvLi(271) dt" dt dt 3XJ dx2 dx\ 5xj dx2 ==ІЇ2= + = 3 (2.72) dt2 dt dt дх\ dx2 дх2 dx{ dx2 d2cp dv, dv? die —L + V- - + dt " dt dt T = y\ (2.73) dtz " dt - dt где ф[ = Зф/St. Очевидно, что ф!, также как и ф, удовлетворяет уравнению Лапласа. Тогда для ф1 имеем следующую краевую задачу: V !=0, dф Фі - (vf + vi = —- - + к. dt v l 2J 2 которая также решается методом граничных элементов, описанным в части 2.3.1. В результате решения определяются значения на свободной поверхности 5фі/5п; dq \/ds вычисляются численным дифференцированием аналогично dcp/8s - (2.68); дщ/дх\, дц \/дх2 вычисляются по формулам, аналогичным (2.69), (2.70). Выражения для dvi/dx]} dv2/5xj, входящие в (2.71), (2.72), получаются с помощью процедуры численного дифференцирования по % следующих выражений [177]: д\] _ 1 D дх.\ dv2 dxY D 5vj dx] dv2 dx2 d% dx dx d% dw\ dx2 dv2 dxj Л (2.74) 3x 3v ? dv где D = (dXi/dx)2 + (dx2/dX)2, a —2 5x2 определяются из уравнения неразрывности и условия отсутствия вихря dv2 _ 3vj Зх2 dvj _ Sv2 (2.75) cbq 3x2 5xj Производная dK/dt в (2.73) вычисляется разностями назад на каждом шаге по времени. Таким образом, становятся известными все величины, необходимые для вычисления нового положения свободной границы и нового значения потенциала на ней со вторым порядком точности по At. Движение частиц свободной поверхности с использованием кинематического условия в эйлеровой форме
Схема предиктор-корректор (метод Эйлера-Коши) Кинематическое условие на свободной границе также используется в эйлеровой форме (2.76), а новая форма свободной поверхности определяется в два этапа. На первом этапе вычисляется «грубое приближение» по методу Эйлера fjn+1=fin+At- E j\ где Ф" - правая часть (2.77), вычисленная на шаге п or(t,f,n(t,e))=vr?- fin А0 На втором этапе происходит уточнение f (t, 0) по формуле П , жП 1 (2.81) Здесь Ф" = o(t, f" j, а Ф" - правая часть (2.77), при вычислении которой используется приближенное решение, полученное на первом этапе: ФР Ф + ДІ,? 1). Аналогичный алгоритм используется и для вычисления потенциала ф"+ для нового положения свободной границы. Схема Рунге-Кутта четвертого порядка Во всех использованных схемах шаг по времени At ограничен сверху условием Куранта At ASmin/Vmax, (2.84) где ASmin - длина наименьшего элемента в случае использования МГЭ/ минимальный шаг сетки в случае МКР, a Vmax - наибольшая скорость, достигаемая на каком-либо элементе свободной поверхности. На первоначальном этапе расчета бралось значение At = 0.001. Как показало исследование, при достаточно малых начальных деформациях (к 1.5) значение At оставалось постоянным в ходе расчетов At = 0.001. При расчетах больших начальных деформаций, начиная с к = 1.5, выбор шага по времени осуществляется автоматически по (2.84). 2.4 Результаты исследования
В плоском случае методом конечных разностей исходная задача решалась с использованием как метода Гаусса-Зейделя, так и алгоритма продольно-поперечной прогонки. Величины периодов колебаний и последовательности форм свободной поверхности, полученные с применением различных методов, совпадают в диапазоне начальных деформаций 1.005 к 2.0. Времена расчета одного периода колебаний с использованием различных методов отличаются незначительно, несмотря на результаты, представленные в таблице 1. Это объясняется тем, что при использовании метода Гаусса-Зейделя для расчета потенциала из уравнения
С использованием метода граничных элементов в плоском случае был проведен цикл расчетов с различными значениями начального отношения осей к в диапазоне 1.002 к 2.0. При этом граница разбивалась на N = 96 постоянных элементов. В расчетах, проведенных cN 96, отличий ни в величине периода, ни в наблюдаемых формах свободной поверхности не отмечалось. Эволюция свободной поверхности в процессе колебаний иллюстрируется рис. 2.3.
Метод расчета
Исследованию течения внутри цилиндра посвящено значительное количество работ [39, 41 - 48, 52, 64, 155, 161, 183-191]. Кроме важности для практических целей (нанесение покрытий на внутреннюю поверхность труб, центробежное литье, использование цилиндра в качестве объемного смесителя) внимание к изучению данного процесса обусловлено многообразием форм его протекания, включающим очевидные предельные случаи при со -»0 и со - оо, где со - частота вращения цилиндра. В первом случае объем жидкости располагается практически целиком в нижней части цилиндра, а на стенках вне объема образуется тонкая пленка. Возникающее при этом возвратное течение с циркуляционной зоной хорошо описывается в рамках приближенного решения, полученного в [47 - 48]. Эта работа выполнена при значительных упрощениях исходных дифференциальных уравнений, сводящих их к обыкновенным. Достоверность результатов [47 -48] подтверждается их сравнением с данными [64], полученными численным методом для полных уравнений Стокса с точным выполнением граничных условий на свободной поверхности. Имеются также попытки использования метода конечных разностей для решения данной задачи [49]. Однако в большинстве работ представлены лишь кинематические характеристики течения.
Второй предельный случай ( у- оо) характерен тем, что жидкость почти равномерно распределяется вдоль твердой стенки. В существующих работах, использующих приближенные методы, в основном рассматривается случай, когда форма свободной поверхности незначительно отличается от окружности [41 -46]. Решение данной задачи, полученное в приближении теории смазки [183] неоднократно уточнялось, проверялось экспериментально и численно и обобщалось на случай степенной жидкости [43, 46, 184 - 189]. Кроме этого, в [186] в приближении теории смазки получено условие перехода от первого вида течения ко второму. Это условие в безразмерной форме для ньютоновской жидкости можно сформулировать в 1/9 г виде: XW =-41, где А, - безразмерный коэффициент заполнения цилиндра жидкостью, W - критическое значение безразмерного параметра W, W = pgR/дхо (R - радиус цилиндра, р - плотность жидкости, g ускорение силы тяжести, ju - коэффициент динамической вязкости). При W W происходит потеря устойчивости слоя жидкости на поднимающейся стенке, что проявляется в скачкообразном изменении решения для толщины слоя. На самом деле в этот момент на свободной поверхности появляется наплыв, в окрестности которого зарождается зона циркуляционного течения так, как это описано в работах [64], [40]. Результаты численного исследования, полученные в настоящей работе, позволяют уточнить условие существования пленочного течения: AWJ л/3, что согласуется с экспериментальными результатами, полученными в [40].
Помимо вышеуказанных двух предельных типов установившегося течения возможен и другой характер поведения свободной поверхности. Например, в [46] утверждается, что существует четыре различных конфигурации установившихся форм свободной поверхности: две с полным смачиванием твердой стенки и две с частичным. В [184] проведено экспериментальное исследование осевой неустойчивости свободной поверхности жидкости в цилиндре, вращающемся с большими угловыми скоростями. Выявлены различные трехмерные периодические структуры. Отмечается возможность капельного поведения жидкости, похожего на эффект «дождевания», наблюдаемый при центробежном литье [192]. Изучению трехмерной структуры в полностью заполненном вращающемся наклонном цилиндре посвящена работа [190]. С помощью разработанного псевдоспектрального метода получены картины течения для чисел Рейнольдса порядка 103. Обсуждение возможных форм свободной границы на внутренней поверхности цилиндра приводится в [189]. На основе решения в приближении теории смазки в [183] получены формы с резким скачком толщины слоя, такие же как, например, в работах [46, 185]. В [185] сделана попытка получения стационарных форм свободной границы с использованием пакета прикладных программ NEKTON (Fluent, Inc.). Решение удалось получить только при Са 6, где Са — (vR/u/a (с коэффициент поверхностного натяжения). В случае больших капиллярных чисел не удавалось получить сходящиеся решения, ввиду сильных деформаций расчетной сетки. Численное моделирование процесса представлено также в работе [155], где на основе метода VOF исследовано течение в миксерах, применяемых в пищевой промышленности. Рассмотрено горизонтальное и вертикальное расположение цилиндра. При этом исследование ограничено случаем, когда основная часть перемешиваемой жидкости находится на дне цилиндра. Таким образом, существует необходимость детального описания поведения свободной поверхности при переходе от начального горизонтального положения к некоторой стационарной форме, получения этих стационарных форм и соответствующих кинематических и динамических характеристик течения.
Для понимания механизма перемешивания, качественного представления степени влияния определяющих факторов используют различные упрощения математической формулировки задачи. Одним из таких упрощений является предположение о том, что смеситель представляет собой бесконечно длинный горизонтальный цилиндр, частично заполненный высоковязкой жидкостью и вращающийся вокруг собственной оси. Данное упрощение является основным в настоящей работе. Исследование рассматриваемого гидродинамического процесса приводит к необходимости решения задачи из класса плоских задач гидродинамики со свободной поверхностью. Для решения подобных задач можно использовать различные конечно-разностные методы [39]. При рассмотрении течений в областях сложной геометрии с большими деформациями свободной поверхности, при наличии нескольких значительно отличающихся характерных размеров, возникают вычислительные трудности, связанные с построением расчетных сеток, которые, кроме того, должны адаптироваться к изменениям свободной поверхности. В связи с этим актуальность приобретают поиски подходов, позволяющих упростить и унифицировать алгоритмическую процедуру численного решения. Одним из путей преодоления этих трудностей является переход к системе граничных интегральных уравнений, эквивалентных исходной системе дифференциальных уравнений. Именно такой подход используется в данной работе.
Для определения стационарной формы свободной поверхности используется математическая постановка, с начальным распределением жидкости слоем равномерной толщины на боковой поверхности цилиндра. С целью исследования характера течения при переходе от начального горизонтального положения к некоторой стационарной форме начальная форма области, занятой жидкостью задается в виде сегмента в нижней части цилиндра, при этом свободная поверхность горизонтальна. Математическая постановка задачи формулируется в приближении ползущего движения. Для численного решения задачи используется метод граничных элементов.
Разработана методика и создана программа расчета течения. Проведены параметрические исследования кинематических и динамических характеристик процесса в зависимости от значений определяющих параметров. Выявлены различные режимы течения и получены критические значения параметров, разделяющие эти режимы. Представлена зависимость удельной мощности потока от основных параметров, характеризующая интенсивность перемешивания в рассматриваемом течении. Проведено сравнение результатов расчета с данными численного и физического моделирования, с аналитическим решением задачи в приближении пленочного течения. Определены диапазоны изменения определяющих параметров, для которых при определении характеристик течения можно использовать приближенное решение [48]. Полученные результаты подтверждаются экспериментальными данными [52, 191] и являются основой для выбора наиболее эффективных режимов процесса перемешивания.
Похожие диссертации на Численное моделирование течений жидкости со свободной поверхностью методом граничных элементов
