CC BY
14ВЕСТНИК ПЕРМСКОГО УНИВЕРСИТЕТА
2021
• ФИЗИКА •
Вып. 1
УДК 532.59; 519.24 PACS 47.20.Ma; 47.11.-j
Моделирование квазистационарного рельефа методом решеточных уравнений Больцмана
И. В. Володин1^, А.А. Алабужев2^
1 Пермский государственный национальный исследовательский университет
2 Институт механики сплошных сред УрО РАН t ivanwolodin@gmail.com
1 alabuzhev@mail.ru
Исследуется поведение двухслойной системы несмешивающихся жидкостей под действием горизонтальных линейных вибраций в поле тяжести. Задача решается численно методом решеточных уравнений Больцмана. Для описания двухфазной системы используется модель НС2, которая позволяет описать поверхность раздела двух жидкостей. Рассматриваются два типа граничных условий на боковых поверхностях: периодические - замкнутая (кольцевая) область для сравнения с аналитическими результатами и условие прилипания для сравнения с экспериментами. Показано наличие квазистационарного рельефа на поверхности раздела жидкостей. Изучено существование рельефа в зависимости от вязкостей жидкостей.
Ключевые слова: метод решеточных уравнений Больцмана; двухслойная система; квазистационарный рельеф; горизонтальные вибрации
Поступила в редакцию 17.12.2020; после рецензии 21.02.2021; принята к опубликованию 27.02.2021
1 Perm State University
2 Institute of Continuous Media Mechanics UB RAS t ivanwolodin@gmail.com
* alabuzhev@mail.ru
The dynamics of two-layer system of immiscible liquids under the action of horizontal linear vibrations in the field of gravity was investigated. The numerical simulation was carried out by the lattice Boltzmann method (LBM) with model D2Q9. For the first time LBM was used to achieve the appearance of frozen wave (quasi-stationary relief) at the interface of two fluids. There are two types of boundary conditions for the sidewalls: a periodic condition for comparison with analytical results and no-slip condition for comparison with experiments. Various computational domains were considered. Both cases with the same viscosities of both phases and different viscosity ratios were studied. HCZ model was used to describe two-phase system and the interface of two liquids. The presence of a frozen wave on the interface of liquids was found. The dependence of liquids viscosity on the relief was studied. The obtained critical wave number coincides well with the theoretically predicted value for liquids with the equal viscosity and vanishing viscosity. The results of numerical calculations show a weak viscosity effect for a more viscous lower liquid. However, the destabilizing effect of viscosity is more significant for a more viscous upper liquid.
© Володин И. В., Алабужев А. А., 2021
Frozen wave simulation by the lattice Boltzmann method
I. V. Volodin1^, A. A. Alabuzhev2i
распространяется на условиях лицензии
Creative Commons Attribution 4.0 International (CC BY 4.0).
Keywords: lattice Boltzmann method; two-layer system; frozen wave; quasi-stationary relief; horizontal vibrations
Received 17.12.2020; revised21.02.2021; accepted27.02.2021
doi: 10.17072/1994-3598-2021-1-59-68
1. Введение
Гидродинамические системы с поверхностью раздела встречаются как в природе, так и в технологических процессах [1-3]. Их нетривиальное поведение привлекает большое внимание исследователей. Важную роль в динамике таких систем играют неустойчивости поверхности раздела двух жидкостей (фаз) [4, 5]. Хорошо известно, что колебательное воздействие может оказывать как стабилизирующий, так и дестабилизирующий эффект на границы раздела жидкостей. Например, в работах [6-8] было показано, что вертикальные колебания стабилизируют неустойчивость Рэлея-Тейлора (более тяжелый слой, помещенный поверх более легкой жидкости, гравитационно неустойчив). Однако эти колебания также приводят к параметрическому возбуждению стабильно стратифицированной двухслойной системы (параметрическое возбуждение, рябь Фарадея) [9-11]. Устойчиво стратифицированные слои жидкостей при горизонтальных колебаниях по-разному ускоряются за счет различного значения плотностей и вязкостей. В этом случае на поверхности раздела жидкостей может возникнуть квазистационарный рельеф или "застывшая волна" (frozen wave) [6, 7, 12], при которой контактная граница между фазами не изменяется. Это связано с тем, что сдвиговый поток, создаваемый периодическим внешним воздействием, состоит из противоточных слоев, так как жидкости несжимаемы [13-15], что является частным случаем неустойчивости Кельвина— Гельмгольца [16].
В работах [16, 17] рассматривалась система двух бесконечных горизонтальных слоев несме-шивающихся жидкостей между параллельными жесткими пластинами под действием горизонтальной вибрации с малой амплитудой и высокой частотой в пределе исчезающей вязкости обеих жидкостях. Использовалось условие сохранения горизонтального объемного потока в слое конечной толщины. Анализ линейной устойчивости показывает, что критический вибрационный параметр b равен:
-,2 _ .
(Р1+Р2Г
2PIP2(PI-P2):
■{ak+ipi-p2)9-)Xgh{kh), (1)
где рг и р2 - плотность нижнеи и верхней жидкости соответственно, И - высота слоя, а -коэффициент поверхностного натяжения, д -ускорение свободного падения, к - волновое число возмущения. Таким образом, при отсутствии вязкости формула (1) описывает неустойчивость Кельвина—Гельмгольца, где коротковолновые и
длинноволновые возмущения подавляются стабилизирующими эффектами капиллярных и плавучих сил соответственно [13-15]. Случай конечных (и равных) кинематических вязкостей показал, что неустойчивость Кельвина— Гельмгольца лишь незначительно зависит от вязкости [17, 18]. Критерий существования наиболее опасной длинноволновой
неустойчивости [16] был распространен на случай произвольной амплитуды колебаний [18]. Было показано, что в предельном случае исчезающей вязкости результаты [18] хорошо согласуется с результатами линейной устойчивости [16]. Кроме того, существует и параметрическая неустойчивость в дополнение к неустойчивости Кельвина—Гельмгольца [19-21].
Экспериментально исследованы осредненные эффекты горизонтальной вибрации на границе раздела двух несмешивающихся жидкостей: фторинертная жидкость FC-40 р1 = 1,85 г/см3, v = 0,02 Ст и машинное масло р^/р2 =2.13, у2 = 1.87 Ст, у2/уг = 93.5, где у1 и у2 - коэффициенты кинематической вязкости нижней и верхней жидкостей соответственно [22]. Результаты этих экспериментов хорошо согласуются с теоретическими предсказаниями [16] при высоких частотах колебаний, т. е. критическое значение колебательного параметра монотонно
уменьшается к постоянному порогу. Однако следующие эксперименты ^С-40 и машинное масло, Р1/Р2 = 1.95, у2 = 17 Ст, у2/у^ = 850) [23] показали, что для конечных частот вибраций экспериментальное пороговое значение колебательного параметра ниже, чем предсказано теоретически [16]. В этом случае условия эксперимента не соответствуют высокочастотному пределу, используемому в теории [16, 17]. Хотя в нижней жидкости из-за ее низкой вязкости толщина поверхностного слоя на несколько порядков меньше, но толщина слоя 8 равна половине толщины слоя масла И (безразмерная частота колебаний ш = 2(й/сг)2 «30) [23]. Позже эффект большого контраста вязкости у2/у\ ~ (102 ^104) на линейном пороге устойчивости был исследован экспериментально и численно [13, 14]. Показано, что модель невязких жидкостей [16] последовательно недооценивает порог
неустойчивости для жидкостей равной вязкости, но в целом переоценивает порог для жидкостей неодинаковой вязкости. Кроме того, существуют четыре области значений параметра У2/У\, в которых имеет место качественно различная динамика, отражающаяся в немонотонной зависимости наиболее неустойчивого волнового
числа и критической амплитуды. Анализ длинноволновой устойчивости [15] показал также, что различное соотношение вязкостей оказывает существенное влияние на порог длинноволновой неустойчивости. Дестабилизирующий эффект вязкости более значителен на конечных частотах вибраций, чем при высокочастотных. В отличие от этого вязкость подавляет длинноволновые возмущения на низких частотах.
В настоящей работе изучалась замороженная волна численным методом решеточных уравнений Больцмана (Lattice Boltzmann Method, LBM). На сегодняшний день авторам не известны статьи, посвященные использованию данного метода для изучения влияния вибрации на динамику поверхности раздела.
2. Численный метод
Одним из универсальных методов численного моделирования аэро- и гидродинамических задач является метод решеточных уравнений Больцмана [24-26]. Идея основана на кинетическом уравнении Больцмана для одночастичной функции распределения f(r,p,t) в фазовом пространстве координат и моментов. Выбранный метод, в отличие от сеточных, не решает уравнение Навье—Стокса, а подходит к решению задачи моделирования движения жидкости с помощью кинетического уравнения Больцмана, являющегося основным в физической кинетике [27]. Классические уравнения гидродинамики в приближении сплошной среды могут быть выведены из кинетического уравнения Больцмана. В настоящее время LBM активно развивается и привлекает большое внимание исследователей. Кинетическое уравнение Больцмана описывает эволюцию функции плотности вероятности в фазовом пространстве. Макроскопические параметры среды, такие как плотность, скорость и энергия, выражаются через функцию распределения. Таким образом, метод газовой динамики используется для моделирования течений, рассматривая систему на мезоскопическом уровне. Основанием для использования методов газовой динамики в задачах непрерывной среды является тот факт, что характерное время макросистемы значительно больше времени длины свободного пробега. В безразмерном виде кинетическое уравнение Больцмана имеет следующий вид [27]:
f + W f + (F + G)Vpf = Q , (2)
где v - микроскопическая скорость жидкости, F -эффективная молекулярная сила взаимодействия, G - сила тяжести, П - оператор столкновения, Vr, Vp - оператор дифференцирования по координатам и моментам соответственно.
Межмолекулярную силу F в уравнении Больцмана (2) можно записать в относительно простой форме
F = -V^ + Fs, (3)
где Fs = KpVAp представляет силу, связанную с поверхностным натяжением (параметр к определяет силу поверхностного натяжения), а ц является функцией плотности:
ц(р) = bp2 RT х- ap\ (4)
Здесь a - параметр межмолекулярного парного потенциала (потенциал Леннарда—Джонса), b = 2xd3 (3m) 1, d - эффективный диаметр молекулы, m - масса молекулы, R - газовая постоянная, T - температура, х - характеризует увеличение вероятности столкновения из-за увеличения плотности жидкости. Отметим, что ц(р) связано с давлением соотношением
ц(р) = p - pRT, а давление удовлетворяет слдеующему уравнению состояния:
p = pRT (1 + bpx)- ap\ (5)
Общепринятым преимуществом LBM являются простота создания параллельного кода, обусловленная микроскопической кинетической идеей, простота программирования, так как вычислительный алгоритм включает в себя только простейшие арифметические операции. К достоинствам также можно отнести легкость задания граничных условий, так как не используются сеточные методы, а состояние системы описывается одной функцией. К недостаткам метода относятся исключительная сложность использования данного подхода в задачах с деформируемыми границами и необходимая малость числа Маха, хотя в последнее время предпринимаются попытки расширить алгоритм на случай произвольно больших чисел Маха, например, алгоритм Particles on Demand [28].
Для моделирования движения применяется двумерная девятискоростная модель D2Q9. Все физические параметры обезразмерены в единицах вязкости. Дискретизация уравнения (2) и его дальнейшее программирование стали возможными после введения приближения Бхатнагара— Гросса—Крука [29] для оператора столкновения.
3. Модель HCZ
Кратко обсудим основную идею расширения метода решеточных уравнений Больцмана на двухфазную систему с межфазной поверхностью раздела. Как было сказано ранее, все макроскопические переменные выражаются через
функцию распределения. Сама функция распределения требует нормировки и в нашем случае нормируется на интегральную плотность всей системы. В общепринятой формулировке эволюционное [24-26] уравнение для функции распределения, учитывающее единичные пространственные и временные шаги, имеет следующий вид:
f (r + Ar, t + At) - f (r, t) = Q.
(6)
Q = -
f - f?
(7)
feq =
P
(2жЯТ у
D/ 2
exp
(v - -
2RT
2
(8)
V f KV feq = - v_- feq .
RT
Макроскопические параметры жидкости в конкретной точке (узле) вычисляются по следующим формулам:
p(r, t) = £ f, - (r, t) = X f Cj,
c2 Ar2
V = -
At
X —
Следовательно, Дг = С;Д£, где скорость частицы ограничена конечным набором из девяти скоростей = 0 ^ 8 модели D2Q9, которая активно используется при решении двумерных задач гидродинамики. Компоненты скоростей для D2Q9 имеют следующий вид: с0(0,0), сх(1,0), с2(0,1), с3(-1,0), с4(0,-1), с5(1,1), с6(-1,1), с7(—1,—1) и с8(1, —1). Приближение Бхатнагара—Гросса—Крука [24-26, 29] используется для оператора столкновения:
где /¿еч - равновесная функция распределения, которая соответствующим образом предписывает зависимость от локальных макроскопических свойств среды, т - релаксационный параметр среды, время между столкновениями двух частиц. Обычно используется распределение Максвелла в качестве равновесной функции /¿еч с точностью до второго члена разложения [24-26]:
где и и р - макроскопические скорость и плотность, соответственно, О - размерность пространства. В нашем случае плоской области О = 2.
Производная Vр/ не может быть вычислена
напрямую, поскольку неизвестна зависимость функции распределения от микроскопической скорости. Учитывая, что /еч является ведущей частью распределения /, а Vр/щ имеет наиболее важный вклад в Vр/, обычно предполагают [2426]:
Введем число Рейнольдса Re = UL/v, где U и L - характерная скорость и характерная длина в макромасштабе, соответственно. Этот
безразмерный параметр характеризует отношение нелинейного и диссипативного членов в уравнении Навье—Стокса. Величина N = L/Дг -число ячеек решетки в направлении характерной длины. Для Дг выбирается единичное значение, что является обычной практикой в LBM, L = N, следовательно, Re = UN/v (назовем этот параметр решетчатым числом Рейнольдса).
Основная проблема при прямом моделировании многофазных потоков уравнением Больцмана (2) связана с вычислением межмолекулярной силы (3). Эта трудность становится весьма существенной для жидкостей, находящихся далеко от критической точки. В этой ситуации численные схемы очень неустойчивы по отношению к небольшим численным ошибкам при вычислении межмолекулярной силы. Для адекватного описания двухфазной системы с помощью LBM существует несколько подходов. Одним из таких является мультифазная модель, которая была предложена в работе [30] и названная по первым буквам фамилии ее авторов -HCZ (He, Chen, Zhang).
По аналогии с функцией распределения f и, учитывая (4)-(5), введем новую функцию
g = RTf + ^(р)Г(0);
(11)
где Г( и) - это функция макроскопической скорости - :
Г (-) =
1
2kRT
exp
( t \2 А (v - иу
2RT
(12)
(9)
Уравнение эволюции для функции g :
* = RTf + Г(0) MP!.
dt dt w dt
(13)
В результате, уравнение Больцмана (2) можно переписать в следующем виде:
df f- feq (v--)(F + G) —+vV f = —J— -^-feq . (10)
dt r x RT
Подставляя (4)-(5) и (8)-(12) в (13), окончательно получим
dg_ dt
g - g
X
+ (v - - )■
(Г(и +G )-(Г(-)-Г(0 ))(p))
i=0
;=0
T
где я89 =РКГ Г(и) + ^(р)Г(0).
Через g мы можем вычислить давление и скорость, используя соотношения:
р = | я^ , рКТи = | .
Однако добавление уравнения (14) не позволяет сформировать полную (замкнутую) систему уравнений для многофазного потока, потому что функция распределения я дает возможность определить только давление и скорость. Предстоит определить еще одну важную переменную в многофазном потоке - макроскопическую плотность Р. Отметим, что плотность не определена однозначно для данного давления. Для несжимаемых многофазных потоков плотности жидкостей вдали от границы раздела точно известны. Далее необходимо только отслеживать разные фазы и формировать четкую поверхность раздела. Для этих целей можно использовать любую индексную функцию ф, удовлетворяющую уравнению (2). Если индексная функция известна, то можно определить физические свойства многофазной среды, например, плотность и вязкость. Концепция индексной функции использовалась во многих численных алгоритмах расчета многофазных течений, от метода конечных объемов до метода фиксации уровня. Тем не менее, вопрос определения поверхности раздела остается предметом исследования традиционного подхода. В нашей формулировке включение межмолекулярной силы автоматически решает эту проблему.
Если жидкости (фазы) не смешиваются, то уравнения эволюции функций распределения [31] для обеих жидкостей будут иметь следующий вид:
С _сеЧ
/Кг + с0С + 1) = /Кг, 0 - п О,
еЯ
д1(г + с0С + 1) = д^г, О - 91 У1 +5/(г, О,
в предположении, что временной шаг ДС = 1. Слагаемые ^ (г, и (г, описывают энергию взаимодействия между фазами, т.е. силы поверхностного натяжения. Разложения функций распределения в ряд по скорости будут иметь следующую форму:
еЧ ( , 2 (с(а (сги)2 (ии)2
2 с?
2 с2
РП ( с\и (с1и)2 (ии)2
,сБ 3,
Ф=^ /;Жр) = Р -С52Р
¿=о
/8 1 ис1 с1д1+2с1Рс
Плотность, вязкость и время релаксации могут быть найдены из выражений:
Ф~Фд
У = Уд +
Р = Рд +
т = та + л л ф-ф1
ф- Ф1
ф- Фд
ф- Ф1
ф- Фд
О; -vg),
(Л -Pg),
О/ "ТА
где индекс I означает жидкую фазу (нижняя жидкость), индекс д - газообразную (верхняя жидкость). Слагаемое, отвечающее за добавление энергии за счет поверхностного натяжения, имеет следующий вид:
О = - 2^) (С- иЖКи) + + (1 - 27)(с - «) - ГКО)),
1 \(с-и)^еч
где р - гидродинамическое давление, с5 -скорость звука, - весовой коэффициент, зависящий от выбранной сетки. Макроскопические свойства поля физической величины могут быть вычислены через формулы:
др
¥* = ар-а,
где а - коэффициент поверхностного натяжения. 4. Тестирование модели HCZ
Для проверки численного метода было получено решение двух тестовых задач, результаты которых хорошо известны: неустойчивость Рэлея—Тейлора и затухание гравитационно-капиллярных волн [32]. Постановка задач предполагала наличие периодических граничных условий на боковых стенках и условия отсутствия скольжения на верхних стенках. В качестве первой задачи, которая была решена с использованием модели НС2, рассматривалась задача Рэлея-Тейлора -неустойчивость границы раздела жидкостей различной плотности в гравитационном поле, когда слой более тяжелой жидкости располагается над слоем более легкой жидкости. Начальное возмущение границы раздела будет расти с течением времени, так как объем более тяжелой жидкости над плоскостью раздела начнет «тонуть» в более легкой жидкости, а объем более легкой жидкости под плоскостью раздела начнет «всплывать» в более тяжелой жидкости.
Качественная эволюционная картина такой неустойчивости, численно смоделированная авторами настоящей статьи, представлена на рис. 1. Размер расчетной сетки равен 128x512, безразмерные плотности р1 = 0.035 и р2 = 0.125, время релаксации т1 =т2 = 0.6875, коэффициент поверхностного натяжения а = 0.01.
гравитационно-капиллярных волн на поверхности жидкость-газ имеет вид [33]:
Рис. 2. Начальное распределение в задаче о затухании гравитационно-капиллярной волны
а2 = + ,
Р
(15)
Рис. 1. Эволюция неустойчивости Рэлея— Тейлора в моменты времени 1=0, 10000, 20000
В случае, когда тяжелая жидкость находится под легкой, возможно существование гравитационно-капиллярной волны, если граница раздела находится вне равновесного состояния. В рассматриваемой задаче рг = 0.12, р1 = 0.04, размер сетки 1000*350. Начальное распределение второй численной задачи показано на рис. 2.
где И - толщина слоя, к - волновое число, д -ускорение свободного падения. График затухания амплитуды колебаний А(Ъ) в узле х = 250 показан на рис. 3. Видно, что амплитуда колебаний со временем затухает.
Полученные численные результаты хорошо совпадают с результатми трехмерного моделирования в работе [34], в которой использовался метод конечных разностей с большим количесвом узлов расчетнгой сетки.
5. Квазистационарный рельеф
Рассмотрено также движение
стратифицированной жидкости в поле горизонтальных линейных колебаний. В системе отсчета, связанной с сосудом, вибрационные и гравитационные ускорения могут быть записаны посредством эффективного ускорения силы тяжести:
9ег = 9 + aa)2Ysm(a)t),
где о) - частота гармонических колебаний, а -амплитуда колебаний, у - единичный вектор в направлении оси колебаний, перпендикулярный д. Задача в этой постановке позволяет вывести осредненные уравнения движения, разделив все процессы на пульсирующие (быстро меняющиеся) и медленные, если характерное время колебаний значительно меньше вязкого гидродинамического времени:
V
а--.
I
Для бесконечного горизонтального слоя в [17] показано, что в случае, когда волновое число к и глубина слоя И удовлетворяют соотношению кк >> 1, наиболее опасные возмущения кт можно найти из соотношения:
,2 _ (Р Р2)
кт =-
и
(16)
Рис. 3. Динамика затухания амплитудных колебаний гравитационно-капиллярной
При наличии поверхностного натяжения в системе дисперсионное соотношение для
Начальное распределение задачи в приведенной выше постановке совпадает с начальным распределением предыдущей задачи (см. рис. 2), а = 0.02, ш = 0.1. Для установления значения решеточного числа Рейнольдса определим скорость через вибрационные параметры: и = а/Т, Т = 2п/ш.
Квазистационарный рельеф полости показан на рис. 4, приближенная средняя область изображена на рис. 5.
Рис. 4. Квазистационарный рельеф на границе раздела двухслойной системы при горизонтальных линейных колебаниях. а = 0.02, ш = 0.1, рх = 0.12, р2 = 0.04, т1=т2 = 0.6875, Де « 5
Л/\Л/^ЛлЛл/Wwww\/^Л/^Л
Рис. 5. Квазистационарный рельеф в
середине области (см. рис. 4)
Подставляя параметры задачи в (6), получаем результат, хорошо согласующийся с аналитической оценкой: формула (6) - кт «0.19, результат численного моделирования (см. рис. 4) -кт «0.21. Наличие периодических граничных условий генерирует волну, которая начинает перемещаться по области, влияя на квазистационарный рельеф, который формируется внешним периодическим воздействием, поэтому рельеф выглядит несимметричным образом. Привнесение в систему диссипативных эффектов позволяет эффективным образом подавлять возникновение бегущей по области волны. Этого можно добиться путем увеличения вязкости, а также замены граничных условий боковых стенок на твердые (см. рис. 6).
а)
б)
Рис. 6. Квазистационарный рельеф в отсутствие бегущей волны а = 0.02, ш = 0.1, рх = 0.12, р2 = 0.04: а — т1 =т2 = 68, Де « 0.01; б — случай твердых границ, т1 — т, = 10, Де « 0.1
Исследовано также влияние вибраций на сосуд конечной длины с твердыми вертикальными боковыми стенками, на которых ставились граничные условия прилипания. Остальные параметры системы оставались прежними. Начальное распределение и результаты численных расчетов при частотах ш = 0.020,0.035,0.040 представлены на рис. 7. На рисунках видно возникновение квазистационарного рельефа, т.е. неустойчивость коротковолнового типа. Этот рельеф разрушается с дальнейшим увеличением частоты колебаний, что согласуется с аналитической теорией [16-18].
1<С 1® 200 КО ХО 360 аОО
б)
В)
г)
Рис. 7. Форма поверхности раздела при а = 0.02, р1 = 0.12, р2 = 0.04, т1 = т2 = 0.6875; а — начальное состояние; квазистационарнцый рельеф при: б — ш = 0.02, Де « 0.4; в — ш = 0.035, Де « 0.7; г — ш = 0.04, Де « 0.8
Для детального сравнения с
экспериментальными результатами [13, 14, 23] и численным моделированием [15, 34] были также рассмотрены системы с большим соотношением вязкостей жидкостей. Отношение вязкостей превышало 102, т. е. верхняя жидкость достаточно сильно вязкая [13, 14, 23]. Дополнительно был
исследован случай с более вязкой нижней жидкостью с аналогичным обратным соотношением. Результаты численного моделирования представлены на рис. 8. Последний случай не был аналитически исследован (см. [13, 14, 23]), но расчеты показывают слабое влияние вязкости при v1/v2 = 110. Возможно, дополнительную роль сыграла сила инерции, т.к. нижняя жидкость более тяжелая. Заметим, что дестабилизирующий эффект вязкости более значим для противоположного случая v2/vx »1.
200
150 >, 1ÜD
50
50 1GG 150 200 250 3DD 350 X
а)
200
150 >, 100
50
50 100 150 200 250 300 350
X
б)
Рис. 8. Поверхность раздела системы двух жидкостей с большим соотношением вязкости v± = 110 v2, т2 = 0.6875, Re « 0.7: a - ш = 0.035; б - ш = 0.04
6. Заключение
Представлены результаты численного исследования возникновения неустойчивости на поверхности раздела двух устойчиво стратифицированных слоев несмешивающихся вязких жидкостей, подверженных горизонтальным вибрациям. Впервые с помощью метода решеточных уравнений Больцмана LBM численно исследовано появление квазистационарного рельефа (замороженной волны) на границе раздела двух сред. Полученное критическое волновое число хорошо совпадает с теоретически предсказанным значением [13-17] для жидкостей с одинаковой вязкостью [15, 18] и исчезающе малой вязкостью [16, 17]. Кроме того, рассмотрено появление замороженной волны на границе раздела двух несмешивающихся жидкостей с большим коэффициентом вязкости аналогично экспериментам [13, 14, 22, 23]. Результаты численных расчетов показывают слабый эффект вязкости для более вязкой нижней жидкости. Однако дестабилизирующий эффект вязкости более значителен для более вязкой верхней жидкости, что согласуется с результатами экспериментов. В дальнейшем планируется
получение новых численных результатов путем изменения температуры, рассмотрение задач с термокапиллярным эффектом.
Работа выполнена при финансовой поддержке Правительства Пермского края (Программа поддержки научных школ Пермского края, грант № C-26/788).
Список литературы
1. Shyy M., Narayanan R. Fluid dynamics at interfaces. Cambridge: Cambridge University Press, 1999. 461 p.
2. Gualtieri C., Mihailovic D. T. Fluid mechanics of environmental interfaces. London: CRC Press, 2018. 500 p.
3. Тятюшкина Е. С., Козелков А. С., Куркин А. А., Курулин В. В., Ефремов В. Р., Уткин Д. А. Оценка численной диффузии метода конечных объемов при моделировании поверхностных волн // Вычислительные технологии. 2019. Т. 24. № 1. С. 106-119
4. Drazin P. G. Introduction to hydrodynamic stability. Cambridge: Cambridge University Press, 2002. 278 p.
5. Bostwick J. B., Steen P. H. Stability of constrained capillary surfaces // Annu. Rev. Fluid Mech. 2014. Vol. 47, P. 539-568.
6. Wolf G. H. The dynamic stabilization of the Ray-leigh-Taylor instability and the corresponding dynamic equilibrium // Z. Physik 1969. Vol. 227. P. 291-300.
7. Wolf G. H. Dynamic stabilization of the interchange instability of a liquid-gas interface // Phys. Rev. Lett. 1970. Vol. 24, P. 444-446.
8. Lapuerta V., Mancebo F. J., Vega J. M. Control of Rayleigh-Taylor instability by vertical vibration in large aspect ratio containers // Phys. Rev. E. 2001. Vol. 64, 016318
9. Kumar K., Tuckerman L. S. Parametric instability of the interface between two fluids // J. Fluid Mech. 1994. Vol. 279. P. 49-68.
10. Silber M., Skeldon A. C. Parametrically excited surface waves: two-frequency forcing, normal form symmetries, and pattern selection // Phys. Rev. E. 1999. Vol. 59. P. 5446-5456
11. Shklyaev S., Khenner M., Alabuzhev A. A. Enhanced stability of a dewetting thin liquid film in a single-frequency vibration field // Phys. Rev. E. 2008. Vol. 77, 036320
12. Bezdenezhnykh N. A., Briskman V. A., Lapin A. Y., et al The influence of high frequency tangential vibrations on the stability of the fluid interface in microgravity // Int. J. Microgravity Res. 1991. Vol. 4, P. 96-97.
13. Talib E., Jalikop S. V., Juell A. The influence of viscosity on the frozen wave instability: theory and experiment // J. Fluid Mech. 2007. Vol. 584. P. 45-68.
14. Talib E., Juell A. Instability of a viscous interface under horizontal oscillation // Phys. Fluids. 2007. Vol. 19, 092102
15. Lyubimov D. V., Khilko G. L., Ivantsov A. O., et al Viscosity effect on the longwave instability of a fluid interface subjected to horizontal vibrations // J. Fluid Mech. 2017. Vol. 814, P. 24-41.
16. Lyubimov D. V., Cherepanov A. A. Development of a steady relief at the interface of fluids in a vibrational field // Fluid Dyn. 1987. Vol. 21, P. 849854.
17. Любимов Д. В., Любимова Т. П., Черепанов A. A. Динамика поверхности раздела в вибрационных полях. М.: Физматлит, 2003. 216 с.
18. Khenner M. V., Lyubimov D. V., Belozerova T. S. et al Stability of plane-parallel vibrational flow in a two-layer system // Eur. J. Mech. B/Fluids 1999. Vol. 18, P. 1085-1101.
19. Kelly R.E. The stability of an unsteady Kelvin-Helmholtz flow //J. Fluid Mech. 1965. Vol. 22. P. 547-560.
20. Lyubimov D. V., Khenner M. V., Stolz M. M. Stability of a fluid interface under tangential vibrations // Fluid Dyn. 1998. Vol. 33, P. 318-323.
21. Shklyaev S., Alabuzhev A. A., Khenner M. Influence of a longitudinal and tilted vibration on stability and dewetting of a liquid film // Phys. Rev. E. 2009. Vol. 79, 051603
22. Ivanova A. A., Kozlov V. G., Evesque P. Interface Dynamics of Immiscible Fluids under Horizontal Vibration // Fluid Dyn. 2001. Vol. 36(3), P. 362368.
23. Ivanova A. A., Kozlov V. G., Tashkinov S.I. Interface Dynamics of Immiscible Fluids under Circularly Polarized Vibration (Experiment) // Fluid Dyn. 2001. Vol. 36(6), P. 871-879.
24. Succi S. The lattice Boltzmann equation for fluid dynamics and beyond. Oxford University Press, 2001.
25. Куперштох А. Л. Моделирование течений с границами раздела фаз жидкость-пар методом решёточных уравнений Больцмана // Вестник НГУ. Серия: Математика. 2005. Т. 5. Вып. 3. С. 29-42.
26. Kupershtokh A. L., Medvedev D. A., Gribanov I. I. Thermal lattice Boltzmann method for multiphase flows // Phys. Rev. E. Vol. 98, 023308.
27. Питаевский Л. П., Лифшиц Е. М. Теоретическая физика. Т. 10. Физическая кинетика. М.: Физматлит, 2007. 536 с.
28. Dorschner B., Bosch F., Karlin I. V. Particles on Demand for Kinetic Theory // Phys. Rev. Lett. 2018. Vol. 121, 130602
29. Bhatnagar P. L., Gross E. P., Krook M. A Model for collision processes in gases. I. Small Amplitude processes in charged and neutral one-component systems // Phys. Rev. 1954. P. 511525.
30. He X., Chen S., Zhang R. A lattice Boltzmann scheme for incompressible multiphase flow and its application in simulation of Rayleigh-Taylor instability // J. Comp. Phys. 1999. Vol. 152, P. 642663.
31. Huang H., Sukop M. C., Lu X. Multiphase lattice Boltzmann methods: theory and application. New Jersey: Wiley-Blackwell, 2015. 373 p.
32. Kull H. J. Theory of the Rayleigh-Taylor instability // Phys. Reports 1991. Vol. 206, P. 197-325.
33. Ландау Л.Д., Лифшиц Е.М. Теоретическая физика в 10 т. Т. 6. Гидродинамика. М.: Физмат-лит. 552 с.
34. Lee H.G., Kim J. Numerical simulation of the three-dimensional Rayleigh-Taylor instability // Computers & Mathematics with Applications. 2013. Vol. 66, No. 8. P. 1466-1474.
References
1. Shyy M., Narayanan R. Fluid dynamics at interfaces. Cambridge: Cambridge University Press, 1999, 461 p.
2. Gualtieri C., Mihailovic D.T. Fluid mechanics of environmental interfaces. London: CRC Press, 2018, 500 p.
3. Tyatyushkina E. S., Kozelkov A. S., Kurkin A. A., Kurulin V. V., Efremov V. R., Utkin D. A. Evaluation of numerical diffusion of the finite volume method when modelling surface waves. Computational Technologies, 2019, vol. 24, no. 1. pp. 106119. (In Russian)
4. Drazin P.G. Introduction to hydrodynamic stability. Cambridge: Cambridge University Press, 2002, 278 p.
5. Bostwick J. B., Steen P. H. Stability of constrained capillary surfaces. Annu. Rev. Fluid Mech, 2014, vol. 47, pp. 539-568.
6. Wolf G. H. The dynamic stabilization of the Ray-leigh-Taylor instability and the corresponding dynamic equilibrium. Z. Physik 1969, vol. 227, pp. 291-300.
7. Wolf G. H. Dynamic stabilization of the interchange instability of a liquid-gas interface. Phys. Rev. Lett. 1970. vol. 24, pp. 444-446.
8. Lapuerta V., Mancebo F. J., Vega J. M. Control of Rayleigh-Taylor instability by vertical vibration in large aspect ratio containers. Phys. Rev. E. 2001, vol. 64, 016318
9. Kumar K., Tuckerman L. S. Parametric instability of the interface between two fluids. J. Fluid Mech. 1994, vol. 279, pp. 49-68.
10. Silber M., Skeldon A. C. Parametrically excited surface waves: two-frequency forcing, normal form symmetries, and pattern selection. Phys. Rev. E. 1999, vol. 59, 5446-56
11. Shklyaev S., Khenner M., Alabuzhev A. A. Enhanced stability of a dewetting thin liquid film in a single-frequency vibration field. Phys. Rev. E. 2008, vol. 77, 036320
12. Bezdenezhnykh N. A., Briskman V. A., Lapin A. Y., et al The influence of high frequency tangential vibrations on the stability of the fluid interface in microgravity. Int. J. Microgravity Res. 1991, vol. 4, pp. 96-97.
13. Talib E., Jalikop S. V., Juell A. The influence of viscosity on the frozen wave instability: theory and
experiment. J. FluidMech. 2007, vol. 584, pp. 4568.
14. Talib E., Juell A. Instability of a viscous interface under horizontal oscillation. Phys. Fluids. 2007, vol. 19, 092102
15. Lyubimov D. V., Khilko G. L., Ivantsov A. O., et al Viscosity effect on the longwave instability of a fluid interface subjected to horizontal vibrations. J. Fluid Mech. 2017, vol. 814, pp. 24-41.
16. Lyubimov D. V., Cherepanov A. A. Development of a steady relief at the interface of fluids in a vibrational field. Fluid Dyn. 1987, vol. 21, pp. 849854.
17. Lyubimov D. V., Lyubimova T. P., Cherepanov A. A. Dinamika poverhnosti razdela v vi-bracionnyh polyah (Interface dynamics under vibrations). Moscow: Fizmatlit, 2003, 216 p. (In Russian)
18. Khenner M. V., Lyubimov D. V., Belozerova T. S. et al Stability of plane-parallel vibrational flow in a two-layer system. Eur. J. Mech. B/Fluids, 1999, vol. 18, pp. 1085-1101.
19. Kelly R.E. The stability of an unsteady Kelvin-Helmholtz flow. J. Fluid Mech. 1965, vol. 22, pp. 547-560.
20. Lyubimov D. V., Khenner M. V., Stolz M. M. Stability of a fluid interface under tangential vibrations. Fluid Dyn. 1998, vol. 33, pp. 318-323.
21. Shklyaev S., Alabuzhev A. A., Khenner M. Influence of a longitudinal and tilted vibration on stability and dewetting of a liquid film. Phys. Rev. E. 2009, vol. 79, 051603
22. Ivanova A. A., Kozlov V. G., Evesque P. Interface Dynamics of Immiscible Fluids under Horizontal Vibration. Fluid Dyn. 2001, vol. 36(3), pp. 362368.
23. Ivanova A. A., Kozlov V. G., Tashkinov S.I. Interface Dynamics of Immiscible Fluids under Circu-
larly Polarized Vibration (Experiment). Fluid Dyn. 2001, vol. 36(6), pp. 871-879.
24. Succi S. The Lattice Boltzmann Equation for Fluid Dynamics and Beyond. Oxford University Press, 2001, 304 p.
25. Kyneprnmox A. K. Modelirovanie techenij s granitsami razdela faz zhidkost'-par metodom reshyo-tochnyh uravnenij Bolzmana. Vestnik NGU. Series: Mathematics, 2005. vol. 5, no. 3. pp. 29-42.
26. Kupershtokh A. L., Medvedev D. A., Gribanov 1.1. Thermal lattice Boltzmann method for multiphase flows. Phys. Rev. E. vol. 98, 023308.
27. Pitaevskii L. P., Lifshitz E. M. Physical Kinetics: Volume 10 (Course of Theoretical Physics S). Butterworth-Heinemann, 452 p.
28. Dorschner B., Bosch F., Karlin I. V. Particles on Demand for Kinetic Theory. Phys. Rev. Lett. 2018. vol. 121, 130602
29. Bhatnagar P. L., Gross E. P., Krook M. A Model for Collision Processes in Gases. I. Small Amplitude Processes in Charged and Neutral One-Component Systems. Phys. Rev. 1954, pp. 511-525
30. He X., Chen S., Zhang R. A Lattice Boltzmann Scheme for Incompressible Multiphase Flow and Its Application in Simulation of Rayleigh-Taylor Instability. J.Comp. Phys. 1999, vol. 152, pp. 642-663.
31. Huang H., Sukop M. C., Lu X. Multiphase lattice Boltzmann methods: theory and application. New Jersey: Wiley-Blackwell, 2015. 373 p.
32. Kull H. J. Theory of the Rayleigh-Taylor instability. Phys. Reports 1991, vol. 206, pp. 197-325.
33. Landau L. D., Lifshitz E. M. Fluid Mechanics: Volume 6 (Course of Theoretical Physics S). Butterworth-Heinemann, 560 p.
34. Lee H.G., Kim J. Numerical simulation of the three-dimensional Rayleigh-Taylor instability. Computers & Mathematics with Applications. 2013, vol. 66, no. 8, pp. 1466-1474
Просьба ссылаться на эту статью в русскоязычных источниках следующим образом:
Володин И. В., Алабужев А. А. Моделирование квазистационарного рельефа методом решеточных уравнений Больцмана // Вестник Пермского университета. Физика. 2021. № 1. С. 59-68. doi: 10.17072/19943598-2021-1-59-68
Please cite this article in English as:
Volodin I. V., Alabuzhev A. A. Frozen wave simulation by the lattice Boltzmann method. Bulletin of Perm University. Physics, 2021, no. 1, pp. 59-68. doi: 10.17072/1994-3598-2021-1-59-68
Сведения об авторах
1. Иван Валерьевич Володин, аспирант, Пермский государственный национальный исследовательский университет, ул. Букирева, д. 15, Пермь, 614990
2. Алексей Анатольевич Алабужев, канд. физ.-мат. наук, научный сотрудник, Институт механики сплошных сред УрО РАН, ул. Акад. Королева, д. 1., Пермь, 614013
Author information
1. Ivan V. Volodin, PhD student, Perm State University, Bukirev str. 15, 614990, Perm, Russia
2. Alexey A. Alabuzhev, Candidate of Physics and Mathematics, Researcher, Institute of Continuous Media Mechanics UB RAS, Acad. Koroleva str. 1, 614013, Perm, Russia
