МЕХАНИКА ЖИДКОСТИ И ГАЗА № 5 • 2009
УДК 532.529.6: 532.6: 519.6
© 2009 г. Л. Б. ДИРЕКТОР, И. Л. МАЙКОВ ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ДИНАМИКИ КАПЛИ ВЯЗКОЙ ЖИДКОСТИ
Предложена численная модель, разработанная на основе метода функции уровня, позволившая описать как нелинейные колебания одиночной капли жидкости, так и процессы дробления и слияния капель. Уравнения Навье—Стокса, записанные в переменных "скорость—давление" на прямоугольной равномерной сетке в цилиндрической системе координат, решаются при помощи метода расщепления по физическим процессам. Получены неос-циллирующие решения для двухфазных сред с характерным отношением плотностей фаз меньше 10-3 и Яе > 1000. Возможности разработанного подхода демонстрируются на примере решения задачи о падении капли из капилляра (отрыв от капилляра, образование "шарика Плато", полет капли, соударение с плоской поверхностью, осцилляции капли на поверхности и растекание). Сравнение результатов моделирования с известными расчетными моделями и экспериментальными данными показало удовлетворительное согласие как по фазам колебаний, так и форме капель.
Ключевые слова: математическое моделирование, вязкая жидкость, функция уровня, шарик Плато, свободная поверхность, нелинейные колебания.
Сложность описания динамики капли вязкой жидкости связана с необходимостью решения системы уравнений Навье—Стокса с учетом граничных условий на свободной поверхности, меняющей свое положение. Эволюция свободной поверхности во многом определяет основные характеристики гидродинамических процессов.
Проблеме разработки численных методов решения задач со свободными границами посвящено большое число публикаций. Следует отметить, что большинство методов разработано для решения задач динамики одиночной вязкой капли, и обобщение этих методов на ансамбль капель проблематично. Наиболее перспективно использование методов сквозного счета. Двухфазная система моделируется однофазной средой с резко меняющимися параметрами на границе раздела областей в предположении, что изменение физических свойств при переходе из одной области к другой происходит внутри этой границы [1—7]. Поверхностные силы на границе раздела областей рассматриваются как объемные и вводятся в уравнение движения. Для определения границы раздела используется уравнение переноса скалярной функции, которая претерпевает разрыв на границе областей [3, 4]. Другой подход заключается в использовании вспомогательной функции — функции уровня [6, 7], имеющей противоположные по знаку значения в разных областях, а нулевой уровень функции соответствует положению границы раздела. Уравнение переноса решается для вспомогательной функции, а значения исходных функций (плотности, вязкости) определяются по функции уровня. Использование функции уровня позволяет описывать динамику областей с изменяющейся топологией.
Построенная с использованием метода адаптивных сеток модель, представленная в [8], позволяет рассчитывать процесс колебаний одиночной капли и за счет аналитического преобразования расчетной области дает возможность с высокой точностью определять форму капли. Однако для описания ансамбля капель и для расчета колебаний большой амплитуды она непригодна из-за нарушения однозначности преобразования расчетной области.
Представленная в настоящей работе модель описывает как процессы нелинейных колебаний одиночной капли жидкости (колебания в процессе свободного падения, колебания на горизонтальной поверхности), так и процессы дробления и слияния нескольких капель.
Возможности разработанного подхода демонстрируются на примере решения задачи о падении капли из капилляра (отрыв капли от капилляра, полет, соударение с плоской поверхностью и растекание капли).
1. Математическая модель. Система из двух несмешивающихся и несжимаемых жидкостей описывается уравнениями Навье — Стокса и неразрывности
р,^ + PXUV)U, = -Vpt +V • T, + Pig (1.1)
dt
V- U, = 0 (1.2)
где i — номер области (среды) (i = 1, 2); pi — плотность; U — вектор скорости; pi — давление; t — время; = (VU, + VU,r) — тензор вязких напряжений; ^ — динамическая вязкость; g — вектор ускорения свободного падения; верхний индекс " T' означает операцию транспонирования. Индекс "1" относится к более плотной среде.
На границе раздела сред выполняется динамическое условие
(ti -T2)n = (pi -p2 + ok)n (L3)
где n — вектор нормали, направленный из первой среды во вторую; к — коэффициент кривизны поверхности; ст — коэффициент поверхностного натяжения.
Представим систему из двух сред как однофазную с резко меняющимися параметрами (плотность и вязкость) на границе раздела областей. Предположим, что плотность и динамическая вязкость каждой среды постоянны.
Введем функцию уровня ф (x, t) (фиг. 1, а) ф(х, t) > 0, ф е Q1; Ф (x, t) = 0, феГ; ф(х, t) < 0,ф е Q2, где x — радиус-вектор точки пространства Qi иГиQ2.
Такое представление двухфазной среды позволяет вводить любое количество областей Qi и Q^. Геометрические характеристики границы раздела Г (в трехмерном случае — поверхности) однозначно определяются по функции уровня.
Вектор нормали к поверхности n = Уф/|Уф|, ф = 0, кривизна границы раздела k = V-Уф/ |Уф|, ф = 0.
Введем скорость (в дальнейшем будем называть ее скоростью жидкости) U = U1, ф> 0; U = U2, ф< 0, причем скорость — непрерывная функция на границе раздела Г областей Qi и Q2 (Ui = U2 при ф = 0).
Граница раздела движется со скоростью частиц жидкости, и уравнение переноса функции уровня записывается в виде
^ + U-Уф = 0 (1.4)
dt
Поверхностная капиллярная сила может быть сведена к объемной силе. Объемная капиллярная сила действует только в узком переходном слое и в пределе толщины переходного слоя, стремящейся к нулю, превращается в динамическое условие на границе раздела (1.3). Выражение для объемной силы, согласно [9], имеет вид
F = -a kb (ф) Уф где 8 (ф) — дельта-функция.
П2(ф < 0) Г(Ф = 0)
г
тЯ
гт
СП б
0
кЯ г
Фиг. 1. Расчетные области (а) и (б): 0.2 — несмешивающиеся жидкости; Г граница раздела сред; ф (х, г) — функция уровня; п — вектор нормали
а
Уравнения (1.1), (1.2) можно представить как
р— + р( и V )и = -Ур + V т + р% + Р (1.5)
дг
V- и = 0 (1.6)
Р = Р2 + (Р1 -р2)#(ф), Ц = Ц 2 + (Ц -Ц 2)#(ф), где Н( ф) - функция Хевисайда: Н(Ф) = 0, Ф < 0, Н(Ф) = 1, ф > 0.
Определим безразмерные переменные следующим образом:
х* = х, и* = и, г* = —, р* = -р--, р* = £, ц* = , е = £
Ь — Ь Р1— Р1 g
где Ь, и — характерные длина и скорость.
Уравнения (1.5), (1.6) в безразмерном виде можно записать как
— + (и*У)и* = -Ур* + -^У-т* + —е + (1.7)
дг* р* Яе Бг Wep*
V и * = 0 (1.8)
где e — единичный вектор в направлении вектора g; т * — безразмерный тензор вязких напряжений; Р* — безразмерная объемная сила; Яе — число Рейнольдса; Бг — число Фруда; We — число Вебера и
р * = Х + (1 -X )#(ф), ц * = п + (1 -П )Н( Ф), , Г1= —
Р1
2. Численный метод. Пусть в момент времени г*п известны поля скоростей и *п и
давления р* . Тогда для расчета неизвестных функций в момент времени г * используется схема расщепления [10], состоящая из трех шагов. Запишем уравнение (1.7) в виде
ЯТТ* 1
ди- = -1 Ур* + г (2.1)
дг* р*
1 1 1 р*
г = -(и* У)и* + — Ут* + — е + —— Яе Бг Wep*
На первом шаге решения уравнения (2.1) предполагается, что перенос количества движения осуществляется только за счет конвекции и диффузии
и** - и*" = ^п
^п + 1
где А1п +1 = I *" +1 -1 *", и ** — промежуточное поле скорости, которое не удовлетворяет уравнению неразрывности.
На втором шаге по найденному промежуточному полю скорости и ** с учетом
и*п +1
* рассчитывается поле давления у| 1 _ур*п +1|=—]_ у. и** (2.2)
Ф*(Фп) ) д*
п + 1
Для решения уравнения Пуассона (2.2) на каждом шаге по времени обычно используются либо итерационные, либо прямые методы. Здесь используется метод сопряженных градиентов [11].
На третьем шаге предполагается, что перенос количества движения осуществляется только за счет градиента давления (конвекция и диффузия отсутствуют)
и*п +1 = и**-м
г 1 ^
^ ур*п + 1 ^р* )
(2.3)
Уравнение Пуассона (2.2) получается применением оператора дивергенции к обеим частям равенства (2.3) и с учетом уравнения неразрывности (1.8)
V-и*п +1 = 0.
Для устойчивости численной процедуры проведем дополнительное сглаживание (размазывание) границы раздела сред (в результате граница приобретает конечную толщину) с помощью модельной функции Хевисайда [5]
ф) =
0, ф < -£
2(+ф+п„„о),
1, ф > е
где е — малая положительная величина.
При численном решении задачи используется шахматная (сдвинутая) сетка. Рассмотрим конвективные члены уравнений (1.4) и (1.7). В общем случае конвективный перенос по оси х1 для произвольной переменной Е, = {и, v, ф} можно записать в виде
р5 _ р ^ +1/2.3 - ^ -1/2.3
3 л
Ах1
где Р, 3(и) — функция, зависящая от компоненты скорости и и аппроксимированная на сетке для переменной £. В качестве базовой схемы используем схему против потока первого порядка точности. Используя интерполяционную формулу Ньютона [12], можно получить конечные разности различных порядков точности и определить поправки к базовой схеме (ЕМО — схемы) [7]. Аналогично определяются конвективные потоки в направлении оси х2.
II
II
I
I
Фиг. 2. Форма капли (I — расчет по разработанной модели, II — расчетные данные [13]) в процессе ее колебаний на горизонтальной поверхности (время, мс): ^ = 0, 0.6, 1.7, 2.8, 11.4, 14.8, 17.6, 21.6 (а-з)
При аппроксимации вязкостных членов (второй член в правой части уравнения (1.7)) и уравнения для кривизны поверхности используется центрально-разностная схема второго порядка точности.
3. Результаты. Все расчеты проводились с использованием прямоугольной равномерной сетки с шахматным расположением узлов в цилиндрической системе координат с симметрией по углу. Расчетная область представляла из себя цилиндр радиусом кЯ и высотой тЯ, где Я — радиус капли; к, т — целые числа, определяющие область вычислений; координаты х1 и х2 соответствуют координатам г и z в цилиндрической системе координат (фиг. 1, б).
Для проверки адекватности разработанной численной модели рассматривался процесс колебаний капли, возникающих
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.