Моделирование движения капли в ограниченной геометрии тема диссертации и автореферата по ВАК РФ 00.00.00, кандидат наук Гуськова Мария Сергеевна

• Гуськова Мария Сергеевна
• кандидат наук
• 2023

Оглавление диссертации кандидат наук Гуськова Мария Сергеевна

1.1 Алгоритм вычислений

1.2 Модель Шана-Чена

1.3 Метод погруженной границы

1.3.1 Сила и момент силы

2 Моделирование колебаний капель

2.1 Колебания одной капли

2.1.1 Параметры моделирования

2.1.2 Результаты моделирования

2.2 Колебания трех капель

2.2.1 Параметры моделирования

2.2.2 Результаты моделирования

2.3 Ограничения подхода

2.3.1 Влияние симметрии

2.3.2 Растворение капли

2.4 Влияние модели 1)2(^19 или 1)2(^27

2.5 Технические детали

2.6 Выводы

3 Некоторые аналитические модели кристаллов

3.1 Цепочка пар точечных вихрей

3.1.1 Математическая модель

3.1.2 Движение пар точечных вихрей

3.2 Интегрирование ОДУ для диполей

3.2.1 Моделирование

3.2.2 Сравнение спектров

3.3 Цепочка классических осцилляторов

3.3.1 Атомы одинаковых масс

3.3.2 Атомы двух масс

3.3.3 Моделирование системы из атомов разных масс

3.4 Формулы для расчета спектра

3.5 Технические детали

3.6 Выводы

4 Моделирование микрофлюидных кристаллов

4.1 Постановка задачи

4.1.1 Внешняя сила

4.1.2 Уравнения движения частицы

4.2 Коллективные колебания около центральной оси

4.3 Коллективные колебания при 7 < 1/2

4.3.1 Оптическая ветвь в продольных колебаниях

4.3.2 Оптическая ветвь в поперечных колебаниях

4.4 Технические детали

4.5 Выводы

Заключение

Приложение А Влияние симметрии и границ на одну частицу

А.1 Параметры моделирования

А.2 Результаты моделирования

А.2.1 Влияние несимметричности

A.2.2 Влияние размера решетки

A.З Заключение

Приложение В Выбор граничных условия для моделирования

потока Пуазейля. Оценка их точности

B.1 Граничные условия на скорость

B. 1.1 Справа и слева устанавливается профиль Пуазейля

(Palabos-coursera)

В. 1.2 Алгоритм

В.2 Справа и слева устанавливается профиль Пуазейля

[Zou and Не, 1997]

В.2.1 Левый нижний угол

В.2.2 Левый верхний угол

В.2.3 Outlet

В.2.4 Правый нижний угол

В.2.5 Правый верхний угол

В.2.6 Алгоритм

В.2.7 Условия на давление [Zou and Не, 1997]

В.2.8 Внешняя сила

В.2.9 Алгоритм

В.З Граничные условия Bounce-back

В.4 Bounce-back

В.4.1 Bounce-back

В.4.2 Операторы столкновений

В.4.3 BGK

В.4.4 IncBGK [Не and Luo, 1997]

В.4.5 Моделирование

В.4.6 Скорость-Скорость

В.4.7 Давление-Давление

В.4.8 Скорость-Давление

В.4.9 Внешняя сила

В.4.10 Сравнение четырех способов получения потока Пуазейля

В.5 Влияние алгоритма bounce-back

В.5.1 Влияние оператора столкновений

В.5.2 Local boundary или Zou, Не

B.5.3 Движение одной частицы

B.6 Заключение

Приложение С Положение равновесия в потоке Пуазейля

C.1 Одна частица в потоке Пуазейля

C. 1.1 Зависимость относительной скорости от

С.2 Равновесное положение цепочки частиц

С.2.1 Расчет относительной скорости

Приложение D Спектры

Литература

Введение

Введение диссертации (часть автореферата) на тему «Моделирование движения капли в ограниченной геометрии»

Актуальность

Микрофлюидика является разделом гидродинамики, в котором изучается течение жидкости в малых объемах [Bazaz et al., 2020]. Микрофлюидные устройства часто решают задачи медицины и микробиологии: сортировка клеток крови по размерам [Li et al., 2007], упорядочивание клеток в течении [Hur et al., 2010, Edd et al., 2008], обнаружение патогенов [Schaaf and Stark, 2020], выполнение операций добавления меток к клеткам для их дальнейшего секвенирования [Edd et al., 2008] и другие.

Частицы или капли другой жидкости могут формировать в течении упорядоченную структуру, в которой частицы находятся примерно на одинаковом расстоянии друг от друга. Такие структуры называются "плавающими" кристаллами из-за их свойства периодичности [Beatus et al., 2006]. Кристаллы так же могут формироваться из-за внешних сил (электрического или магнитного поля) или из-за гидродинамического взаимодействия частиц друг с другом [Del Giudice et al., 2021]. По типу частиц кристаллы можно классифицировать на кристаллы из капель, кристаллы из инерциальных частиц, кристаллы из вязкоупругих частиц и обобщенные (compartmentalised). Микрофлюидные кристаллы интересны из-за фундаментальных задач и из-за практического применения в синтезе коллоидных структур, при выполнении цифровой полимеразной цепной реакции (ПЦР) и анализа одиночных клеток (вид биологических экспериментов, когда нужно качественно и количественно описать найденные молекулы РНК).

Кристаллы из капель исследуются в физических экспериментах [Beatus et al., 2006, Schiller et al., 2015] и в компьютерном моделировании [Uspal and Doyle, 2012, Schaaf and Stark, 2020, Janssen et al., 2012]. В работе [Beatus et al., 2006] 2006 года описан эксперимент, где в одномерном кристалле из сжатых капель воды, движущихся в масле, обнаружены поперечные и продольные коллективные колебания. Эти колебания похожи на колебания фононов в твердых телах. Причина возникновения спектра объясняется дипольным взаимодействием микрочастиц за счет большой разницы скоростей микрочастиц и окружающего потока. Представляется важным исследовать

вопрос: является ли этот механизм единственным для возникновения фононного спектра?

В диссертации исследуются микрофлюидные кристаллы из капель и частиц, которые формируются без воздействия внешнего поля. Компьютерное моделирование проводится с использованием метода решеточного уравнения Больцмана (LBM) [Krüger et al., 2013].

Постановка проблемы

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

Исследование проведено с использованием компьютерного моделирования: численное интегрирование обыкновенных дифференциальных уравнений, моделирование потока и помещенных в него капель другой жидкости с помощью трехмерной модели метода решеточного уравнения Больцмана (LBM) и модели Шана-Чена и двумерное моделирование потока и помещенных в него твердых частиц с помощью LBM и метода погруженной границы. Основными задачами работы являются:

1. Исследовать влияние колебаний капли на движение цепочки капель.

2. Проверить, действительно ли дипольное взаимодействие может генерировать наблюдаемый в экспериментах спектр.

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

Дополнительные задачи:

1. Исследование влияния начальных условий (симметрия, угловая скорость и линейная скорость) на движение одной частицы в ограниченном пространстве.

2. Выбор граничных условий для моделирования потока Пуазейля, оценка их точности.

3. Проверка эффекта Сегре-Зильберберга и расчет равновесных положений для цепочек частиц.

Степень разработанности темы исследования

В работе [Beatus et al., 2012] представлены результаты серии экспериментов над микрофлюидными кристаллами. Кристалл состоит из капель воды, которые движутся в масле. Капли впрыскиваются в поток масла через примерно одинаковые расстояния и увлекаются потоком масла. Из-за трения капель о верхнюю и нижнюю стенки канала, скорость капель примерно в 5 раз меньше скорости основного потока. Такая существенная разница в скоростях приводит к возникновению дипольного взаимодействия между каплями. Капли в кристалле расположены не строго на одинаковых расстояниях и на центральной оси, а с некоторым случайным смещением, вызванным несовершенством канала (неровность стенок) и асимметрией генератора капель. В одномерной цепочке капель возникают продольные и поперечные колебания, спектр которых похож на фононный спектр. Авторы [Beatus et al., 2012] объясняют возникновение фононного спектра именно дипольным взаимодействием. Для проверки они интегрируют обыкновенные дифференциальные уравнения для диполей и полученный численно спектр похож на реальный. Аналогично в работе [Liu et al., 2012] представлены результаты численного интегрирования уравнений движения капель, полученных в [Beatus et al., 2006].

На движение частиц в канале при низких числах Рейнольдса влияет близость стенок, которые тормозят частицу, а также вязкость жидкости, скорость частиц относительно скорости потока и вращение частиц. В экспериментальной работе [Segre and Silberberg, 1961] обнаружено, что частицы в цилиндрическом течении Пуазейля концентрируются в кольце на некотором расстоянии от стенок и центральной оси. В другой работе [Oliver, 1962] показано, что при отсутствии вращения частицы сходятся на центральную ось. В работе [Saffman, 1965] получена формула для подъемной силы, которая действует на частицу, 81.2Va2K1/2/i'1/2 + о(ъ>—1/2), где V — скорость потока,

а—радиус частицы, к— градиент скорости потока, и — кинематическая вязкость. Ограничения аналитических работ: совсем небольшие числа Рейнольдса [Saffman, 1965, Rubinow and Keller, 1961] или умеренные числа Рейнольдса [Asmolov et al., 2018], небольшое значение отношения диаметра частицы к диаметру канала [Saffman, 1965, Rubinow and Keller, 1961], невозможность оценить влияние колебаний капли и отраженных волн от стенок канала. В работах компьютерного моделирования, напротив, сложнее моделировать небольшие числа Рейнольдса из-за необходимости использовать совсем низкую скорость или большую вязкость. В работах [Yang et al., 2005, Esipov et al., 2020] исследовано влияние числа Рейнольдса на положение равновесия частиц, при этом увеличение числа Рейнольдса приводит к тому, что равновесное положение сдвигается к стенке канала. Другим параметром задачи является отношение диаметра частицы к диаметру канала 7. Увеличение 7 приводит к тому, что равновесное положение частиц сдвигается к центру канала.

Для имитационного моделирования инерциальной микрофлюидики используются методы на основе решения уравнений Навье-Стокса или методы LBM [Bazaz et al., 2020]. Инерциальной микрофлюидикой называется миграция случайно распределенных частиц к некоторым положениям равновесия внутри микроканала в зависимости от размера частицы. Методы, основанные на решении уравнений Навье-Стокса, испытывают трудности в определении интерфейса между каплей (или частицей) и потоком, также эти методы затратны по времени. LBM имеет преимущество из-за относительной простоты параллелизации алгоритма и множества моделей для различных физических эффектов, в том числе и для моделирования капель одной жидкости внутри другой.

Описание методологии исследования

Диссертационное исследование опирается на результаты имитационного моделирования. Для моделирования капель жидкости, помещенных в другую жидкость, используется LBM и модель Шана-Чена для многокомпонентных жидкостей [Shan and Chen, 1993]. Для оценки точности моделирования проверяется аналитическая формула Рэлея для частоты колебаний капли [Rayleigh, 1879]. Для моделирования движущихся твердых частиц в течении

используется метод погруженной границы [Inamuro, 2012] вместе с LBM. Для оценки правомерности использования этих методов полученные моделированием спектры колебаний микрофлюидного кристалла сравнивались с результатами физического эксперимента [Beatus et al., 2006].

Основные результаты исследования Влияние колебаний капли на цепочку капель

В одномерной цепочке капель иногда возникают неустойчивости, которые приводят к разрушению кристалла [Beatus et al., 2012]. В диссертационной работе было выполнено моделирование колебаний капли для оценки влияния на соседние капли. Получилось, что колебания капли могут повлиять на устойчивость цепочки капель, если скорость впрыскиваемых капель v¿7 умноженная на расстояние между каплями а, близка к наименьшей частоте Рэлея [Rayleigh, 1879], т.е., v¿ а « Эффект привел бы к деформации

капли, к изменению расстояния между каплями и изменению осевой симметрии цепочки. Возможно, именно этим резонансом можно объяснить нарушение порядка в экспериментах [Beatus et al., 2012].

Используется модель многокомпонентной жидкости Шана-Чена [Shan and Chen, 1993] для двух несмешиваюгцихся жидкостей. В

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

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

Дипольное взаимодействие

Для дополнительной проверки дппольного взаимодействия как причины для возникновения фононного спектра в кристалле капель была реализована цепочка, состоящая из пар точечных вихрей с противоположными по знаку циркуляциями. Спектры для возникших колебаний в цепочке качественно похожи на спектры в [Beatus et al., 2006].

Рассматривается система 2N вихрей с циркуляциями = 1 ,...,N), N = 32 - количество пар вихрей, т.е. диполей. Координаты вихрей Xi,yi(i = 1 ,...,N). Гамильтониан такой системы определяется выражением [Манаков, 1983]

где скобка Пуассона определяется формулой {А, В} = к—1(firff" — f^Hr)•

Используются периодические граничные условия вдоль оси X. Вихри в паре с индексом i на каждом шаге интегрирования смещаются на один и тот же вектор (xi,yi) = (O.5(i0 + х1), 0.5(i¡¡ + у1)), который является средним смещением для вихрей в этой паре. Среднее считается для того, чтобы не менялся размер диполя, иначе вихри будут разбегаться.

На Рис.1 показаны спектры колебаний для координат где £ это

расстояние по X между соседними парами вихрей и у это у—координата верхнего из вихрей в диполе. Спектры были получены с помощью двумерного быстрого преобразования Фурье по пространству и времени. Спектры качественно совпадают с результатами [Beatus et al., 2012].

Фононный спектр

В диссертационной работе с использованием имитационного моделирования был получен фононный спектр для кристалла из жестких частиц, которые движутся примерно со скоростью потока. Т.к. относительная скорость частиц относительно воды невелика (< 5%), то дипольное взаимодействие между

Уравнения движения записываются в следующем виде:

{хг,Н },уг = {уг,Н},

(1)

(а) Спектр колебаний но х

3 0.0

-0.1

(Ь) Спектр колебаний но у

Рис. 1: Спектры колебаний для пар точечных вихрей.

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

С использованием метода погруженной границы вместе с LBM для расчета течения [Inamuro, 2012] было проведено моделирование движения жестких круглых частиц в течении Пуазейля. Рассматривалась плоская задача. В работе [Beatus et al., 2012] рассматривается квазидвумерная задача: высота канала много меньше его длины и ширины, из-за этого возникает трение капли о стенки и капли движутся медленнее основного потока.

Расчетная область является длинным каналом, по оси Y сверху и снизу канала граничные условия без прилипания, вдоль оси X периодические граничные условия. В начальный момент времени частицы радиуса R = 10 помещаются на центральную ось канала шириной W = 40 с расстоянием между центрами частиц а = 4R. К координатам частиц добавляется небольшое случайное возмущение. Число Рейнольдса выбиралось для разных параметров моделирования так, чтобы оно было равно 1 или меньше. Число Рейнольдса Re ~ 1.

Для таких параметров канала равновесное положение частиц находится на центральной оси, частицы пытаются сойтись к этому положению, но гидродинамическое взаимодействие с остальными приводит к колебаниям поперек оси. Поперечные колебания более выражены, чем продольные, и их спектр схож с результатами эксперимента [Beatus et al., 2012]. Продольные колебания в нашей постановке имеют больший период, чем поперечные и требуют длительного времени моделирования для получения

20 10 0 -10 -20 -30

I

ТТ]^ ./ |

Шш\

-10 -5 0 5 10

20 18 16 14 12 10

20 10 0 -10 -20 -30

уу я гг -

16 14

1-12 10

8

-10

-5

10

(а) Спектр колебаний но х

(Ь) Спектр колебаний но у

Рис. 2: Спектры колебаний кристалле, расположенном на оси^/Ж = 0.748 в канале шириной ]У = 80

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

Фононный спектр для разных положений равновесия

В диссертационной работе был численно получен спектр продольных и поперечных колебаний для кристалла, расположенного близко к стенке.

Для широких каналов положение равновесия на центральной оси является неустойчивым, а положение на некоторой оси между центром и стенкой устойчивым. Например, для ширины канала ЦТ = 80 и радиуса частицы Я = 10, такой осью является у/]У ~ 0.748 и симметричная ей относительно оси канала у/}¥ ~ 0.252. Написать про одну ветку и что совпадает с физическим

Цепочка частиц, размещенная в равновесном положении около стенки, с некоторым небольшим начальным возмущением координат будет колебаться в продольном и поперечном направлении. В отличии от случая симметричного расположения цепочки и колебаний вокруг центральной оси в спектре (Рис.2) возникает вторая ветвь колебаний.

Положения, выносимые на защиту

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

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

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

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

Научная новизна

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

2. Было показано, что причиной возникновения фононного спектра в микрофлюидном кристалле может быть гидродинамическое взаимодействие со стенками канала. Ранее в работе [Beatus et al, 2006] как причина колебаний было указано дипольное взаимодействие капель. В нашей постановке скорость частиц относительно основного потока невелика (менее 5%), поэтому не возникает заметного дипольного взаимодействия.

3. Продемонстрированао возникновение оптической ветви в спектре продольных и поперечных колебаний в кристалле, расположенном не на оси симметрии.

Апробация результатов исследования

1. "Использование модели со скоростью звука \/2/3 в методе решеточного уравнения Больцмана на примере движения частицы в потоке Пуазейля",

Межвузовская научно-техническая конференция студентов, аспирантов и молодых специалистов имени Е. В. Армейского, 27.02-7.03 2023, Москва

2. "Wave spectrum of the flowing drops", Суперкомпьютерные дни в России, 27-28 сентября 2021 г., Москва

3. "Drop chain simulation with Lattice Boltzmann method", 29th International Conference on Discrete Simulation of Fluid Dynamics, сентябрь 13-17 2021, Viterbo, Italy

4. "LBM simulation of the chain of three drops", Conference on Computer Simulation in Physics and beyond, October 12-16, 2020, Moscow

Публикации

Все публикации входят в международную систему цитирования Scopus.

1. Guskova М. et al. Simulation of the Particle Dynamics in the Two-Dimensional Poiseuille Flow with Low Reynolds Number //Lobachevskii Journal of Mathematics. - 2022. - T. 43. - №. 2. - C. 381-385.

2. Guskova M., Shchur V., Shchur L. Simulation of drop oscillation using the lattice Boltzmann method //Lobachevskii Journal of Mathematics. - 2020. - T. 41. -№. 6. - C. 992-995.

3. Guskova M., Shchur L. Immersed boundary simulation of drop stability //Journal of Physics: Conference Series. - IOP Publishing, 2021. - T. 1740. - №. 1. - C. 012026.

4. Guskova M., Shchur L. Wave Spectrum of Flowing Drops //Russian Supercomputing Days. - Cham : Springer International Publishing, 2021. - C. 283-294.

5. Shchur L., Guskova M. Drop Oscillation Modeling //Russian Supercomputing Days. - Cham : Springer International Publishing, 2020. - C. 198-206.

Содержание диссертации

В главе 1 описаны методы, используемые в диссертационной работе:

метод решеточного уравнения Больцмана (LBM), модель Шана-Чена для

моделирования многокомпонентной жидкости и метод погруженной границы для моделирования движущейся границы в потоке. В разделе 2 представлены результаты моделирования мягких капель и исследовано влияние колебаний капли на соседние капли в цепочке. В главе 3 приведены результаты численного интегрирования известных моделей кристаллов и результаты моделирования диполя, как пары точечных вихрей. Моделирование микрофлюидного кристалла в потоке Пуазейля для разных положений равновесия описано в 4. В приложении А изучено влияние начальных условий и симметрии на движение одной частицы. В приложении В сравниваются различные граничные условия, которые позволяют получить поток Пуазейля. В приложении С моделируется движение одной капли в течении Пуазейля и цепочки для получения равновесных положений (эффект Сэгре-Зильберберга). В приложение Б вынесены некоторые спектры из раздела 4 для удобства.

1 Метод решеточных уравнений Больцмана

В кинетической теории состояние молекулы описывается координатами х = (х, у, х) и скоростями £ = (£х, £у). Прямое атомистическое моделирование сложно, т.к. в обычном макроскопическом объеме газа находится ~ 1023 молекул. Другой подход заключается в решении кинетических уравнений: в описании газа через функцию плотности вероятности. функция плотности вероятности скорости / (х, Ь) задает вероятность присутствия молекулы в некотором объеме с заданной скоростью в момент времени

Можно определить следующие макроскопические переменные, как моменты функции распределения скоростей:

• Плотность р (х, £) = т / / (х, где т— это масса одной молекулы

• Скорость потока и (х,1) = т/р§ (х, ^,1) - математическое ожидание скорости

• Температура Т (х,1) = т/ (3Яр) / |£|2 / (х, £,1) д3^ - математическое ожидание квадрата модуля скорости |£|2. Я- универсальная газовая постоянная

• Тензор напряжений Раф = т / (£а — иа) — ир) / (х, с13^

Динамика в пространстве-времени для функции распределения скоростей обычно определяется уравнением Больцмана. Уравнение Больцмана - это уравнение баланса для плотностей частиц.

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

= ..

\ГгапвротЪ ^^ ^оШвгоп

9/ + , д/ г ^ Ж" + ^ = °]

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

Основываясь на теории 13 моментов Града [Grad, 1949, Grad, 1958], функции распределения в пространстве скоростей раскладываются в ряд по ортонормальным многочленам Эрмита.

1

/ (г, г) « Г (г, = ш (^ -а(п) (г, *): н(п)

п=0

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

гр

а(0) = 1 ¡й3^ = р,

а(1) = I = ри.

Интеграл по пространству скоростей должен вычисляться для каждого из коэффициентов разложения по многочленам Эрмита или для вычисления макроскопической переменной. Чтобы аппроксимировать интегралы применяются квадратуры Гаусса, которые дают точные значения для полиномов вплоть до соответствующей степени.

«(п) (г, = { / (г, Ь г) н^ (£) ^ = ^^ /" (Г, ь г) н^ (е)

Функция плотности вероятности вычисляется в квадратурных точках: это и есть переменные, которые моделируются в решеточном методе Больцмана.

ь (Г, г) = ^/ (Г, г) /ш (6)

Функции также называются популяциям и частиц, % пробегает значения от 1 до числа дискретизированных частиц.

1.1 Алгоритм вычислений

Дискретизируя уравнение Больцмана в пространстве скоростей, в пространстве координат и по времени, получается решеточное уравнение Больцмана:

ft (x + CiAt,t + At) = fi(x,t) + üi([x\,t).

Согласно этому выражению частицы f¡(x,t) передвигаются со скоростью c¡ в соседнюю точку х + c¡ At за врем я At. В то же время на частицы воздействует оператор столкновений Этот оператор моделирует столкновения частиц, перераспределяя частицы по популяциям f¡ в каждом узле.

Существует много разных приближений операторов столкновения простейший оператор, который также может использоваться для моделирования Навье-Стокса, это Bhatnagar-Gross-Krook (BGK) оператор [Bhatnagar et al., 1954]:

f-- feq

Q (f) = - At.

T

Популяции релаксируют к равновесному состоянию flq со скоростью, определяемой временем релаксации т.

Равновесное состояние задается как:

оеп, (, и Ci (и • Ci)2 и • и\

fft) = «>-Р{! + -¿Г + ЧТ^ - ) ' (2>

где веса w¡ соответствуют табору ско ростей, cs скорость звука. Равновесное состояние такое, что его моменты равны моментам f¡7 т.е. Vif? = £ifi = р и = ри. Равновесное состояние f¡q зависит только

от локальной плотности р и ^^^^^стп по тока и. Формулы для вычислений: P(x,t) = Ег f¡(x,t)) pu(x,t) = Cifi(x,t).

Рассмотрим модель D3Q27: В каждом узле решетки определен набор скоростей аИ наборы, использованные в работе, показаны на Рис. 3. В представлении D3Q19 [Krüger et al., 2017], 19 скоростей c¡ указывают из центра куба в середины 12 ребер, в центры 6 граней и одна скорость в центр куба (нулевая), i = 0,1,2,..., 18. В представлении D3Q27 к 19 скоростям из D3Q19 добавляются 8 скоростей, указывающих из центра в вершины, i = 0,1,2,..., 26, соответственно. На Рис. 3 цветом показаны скорости, сгруппированные по длинам. Популяциями называются переменные, представляющие плотность вероятности функции распределения частиц.

Каждый шаг моделирования делится на две части (как и в решеточном газе):

Рис. 3: Дискретные направления скоростей в трехмерной модели. Левая колонка: модель ЭЗдЮ. Правая колонка: модель 03(^27.

• Столкновение

• Распространение.

Решеточное уравнение ВСК (полностью дискретизировашюе уравнение Больцмана с оператором столкновений ВСК) выглядит следующим образом:

Мх + с—г, г + —) = Мх, I) — — (Мх, I) — ^ (X, г))

т

Можно разложить это уравнение на две раздельные части: 1. Первая часть - это столкновение (или релаксация),

¡:(х,г) = Мх,1) — — (Мх,г) — ^ (х,г)),

т

где /* представляет функцию плотности вероятности поеле столкновений и может быть найдена из по формуле (2). Обычно и эффективно реализовать столкновение можно в форме:

П(х,г) = Мх,г)(1 — + /Г (х,г)—

Выражение упрощается для т/At = 1, где ff(x,t) = fftq(x,t) 2. Шаг распространения

f(x + Ci At, t + At) = f*(x,t)

Т.о. идея метода заключается в реализации двух частей: в столкновении и распространении. Столкновение - это просто алгебраические локальные операции. Сначала вычисляется плотность р и макроскопическая скорость и дЛЯ того? чтобы найти равновесные функции плотности вероятности f^q по (2) и функцию плотности вероятности после столкновения /*. После столкновения, полученная функция плотности вероятности f- ^ распространяется в соседние узлы. После того, как оба этапа сделаны, прошел один шаг по времени, операции повторяются.

1.2 Модель Шана-Чена

Модель Шана-Чена применяется для моделирования многофазных или многокомпонентных потоков [Succi, 2018]. Модель основана на псевдопотенциале и предложена в [Shan and Chen, 1993].

Для каждой фазы или для каждой жидкости вводятся свои функции плостности вероятности f¡R, fs, условно можно обозначить разные жидкости разными цветами, например, красным и синим. Каждая жидкость должна удовлетворять собственному LBGK уравнению:

Источник Sis отвечает за мезоскопическое взаимодействие жидкостей и описывает разделение из-за эффектов поверхностного натяжения. В модели Шана-Чена

где сила Ра для каждого из в видов получается из парного потенциала взаимодействия:

Ai fis = -us(fis - f¡s ) + Sis,i = 1,...,b,s = {Red, Blue]

s' i

Физический смысл следующий: в каждом узле решетки х

равнодействующая сила действующая на частицу вида в является суммой обменов импульсов с частицами всех других видов в соседних узлах Х{ = х + Матрица псевдопотенциалов У88'(х,у) выбирается в виде пропагатора:

V,,' (х,у) = Ъа{х)ва8' (х,у)Ъ8' (у),

где С88'(х,у) — функция Грина, выражающая парное взаимодействие между видами в и в' в узлах х и у соответственно, и ^8(х) = локальный

функционал плотности.

Для большинства практических целей это взаимодействие принимает простую форму:

^ I телн |ж - у\ < стах

С88' (х,у) = <

I 0 иначе

Стах— определяет максимальную амплитуду дискретных скоростей, например,

Список литературы диссертационного исследования кандидат наук Гуськова Мария Сергеевна, 2023 год

Литература

[Манаков, 1983] Манаков, C.B. и Щур, . (1983). Стохастичность в двухчастичном рассеянии. Письма в ЖЭТФ, 37(1):45-48.

[Asmolov et al., 2018] Asmolov, E. S., Dubov, A. L., Nizkaya, T. V., Harting, J., and Vinogradova, О. I. (2018). Inertial focusing of finite-size particles in microchannels. Journal of fluid mechanics, 840:613-630.

[Bazaz et al., 2020] Bazaz, S. R., Mashhadian, A., Ehsani, A., Saha, S. C., Kriiger, T., and Warkiani, M. E. (2020). Computational inertial microfluidics: A review. Lab on a Chip, 20(6):1023-1048.

[Beatus et al., 2012] Beatus, T., Bar-Ziv, R. H., and Tlusty, T. (2012). The physics of 2D microfluidic droplet ensembles. Physics Reports, 516(3):103—145.

[Beatus et al., 2006] Beatus, T., Tlusty, T., and Bar-Ziv, R. (2006). Phonons in a one-dimensional microfluidic crystal. Nature Physics, 2(ll):743-748.

[Bhatnagar et al., 1954] Bhatnagar, P. L., Gross, E. P., and Krook, M. (1954). A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component systems. Physical review, 94(3):511.

[Del Giudice et al., 2021] Del Giudice, F., D'Avino, G., and Maffettone, P. L. (2021). Microfluidic formation of crystal-like structures. Lab on a Chip, 21(ll):2069-2094.

[Di Ilio et al., 2019] Di Ilio, G., Chiappini, D., Ubertini, S., Bella, G., and Succi, S. (2019). A moving-grid approach for fluid-structure interaction problems with hybrid lattice Boltzmann method. Computer Physics Communications, 234:137 145.

[Edd et al., 2008] Edd, J. F., Di Carlo, D., Humphry, K. J., Kôster, S., Irimia, D., Weitz, D. A., and Toner, M. (2008). Controlled encapsulation of single-cells into monodisperse picolitre drops. Lab on a Chip, 8(8):1262-1264.

[Esipov et al., 2020] Esipov, D. V., Chirkov, D. V., Kuranakov, D. S., and Lapin, V. N. (2020). Direct numerical simulation of the segre-silberberg effect using immersed boundary method. Journal of Fluids Engineering, 142(11).

[Freitas et al., 2018] Freitas, V. M., Hilfenhaus, G., and Iruela-Arispe, M. L. (2018). Metastasis of Circulating Tumor Cells: Speed Matters. Developmental Cell, 45(l):3-5.

[Garstecki and Whitesides, 2006] Garstecki, P. and Whitesides, G. M. (2006). Flowing crystals: Nonequilibrium structure of foam. Physical Review Letters, 97(2):1-4.

[Grad, 1949] Grad, H. (1949). On the kinetic theory of rarefied gases. Communications on pure and applied mathematics, 2(4):331-407.

[Grad, 1958] Grad, H. (1958). Principles of the kinetic theory of gases. In Thermodynamik der Gase/Thermodynamics of Gases, pages 205-294. Springer.

[Guo et al., 2002] Guo, Z., Zheng, C., and Shi, B. (2002). Discrete lattice effects on the forcing term in the lattice boltzmann method. Physical review E, 65(4):046308.

[Guskova, 2023] Guskova, M. (2023). Lbm simulation of particles in poiseuille flow, https://github.com/msguskova/lbm_chain_of.particles.

[Guskova et al., 2022] Guskova, M., Burovski, E., Shchur, V., and Shchur, L. (2022). Simulation of the particle dynamics in the two-dimensional poiseuille flow with low reynolds number. Lobachevskii Journal of Mathematics, 43(2):381-385.

[Guskova and Shchur, 2021a] Guskova, M. and Shchur, L. (2021a). Immersed boundary simulation of drop stability. In Journal of Physics: Conference Series, volume 1740, page 012026. IOP Publishing.

[Guskova and Shchur, 2021b] Guskova, M. and Shchur, L. (2021b). Wave spectrum of flowing drops. In Russian Supercomputing Days, pages 283-294. Springer.

[Guskova et al., 2020] Guskova, M., Shchur, V., and Shchur, L. (2020). Simulation of drop oscillation using the lattice boltzmann method. Lobachevskii Journal of Mathematics, 41(6):992-995.

[He and Luo, 1997] He, X. and Luo, L.-S. (1997). Lattice boltzmann model for the incompressible navier-stokes equation. Journal of statistical Physics, 88(3):927-944.

[HeetaL, 1997] He, X., Zou, Q., Luo, L.-S., and Dembo, M. (1997). Analytic solutions of simple flows and analysis of nonslip boundary conditions for the lattice boltzmann bgk model. Journal of Statistical Physics, 87(1) :115—136.

[Hur et al., 2010] Hur, S. C., Tse, H. T. K., and Di Carlo, D. (2010). Sheathless inertial cell ordering for extreme throughput flow cytometry. Lab on a Chip, 10(3):274-280.

[Inamuro, 2012] Inamuro, T. (2012). Lattice Boltzmann methods for moving boundary flows. Fluid Dynamics Research, 44(2).

[Inamuro et al., 2000] Inamuro, T., Maeba, K., and Ogino, F. (2000). Flow between parallel walls containing the lines of neutrally buoyant circular cylinders. International journal of multiphase flow, 26(12):1981-2004.

[Janssen et al., 2012] Janssen, P. J., Baron, M. D., Anderson, P. D., Blawzdziewicz, J., Loewenberg, M., and Wajnryb, E. (2012). Collective dynamics of confined rigid spheres and deformable drops. Soft Matter, 8(28):7495-7506.

[Kostenetskiy et al., 2021] Kostenetskiy, P., Chulkevich, R., and Kozyrev, V. (2021). Hpc resources of the higher school of economics. In Journal of Physics: Conference Series, volume 1740, page 012050. IOP Publishing.

[Krüger et al., 2013] Krüger, T., Gross, M., Raabe, D., and Varnik, F. (2013). Crossover from tumbling to tank-treading-like motion in dense simulated suspensions of red blood cells. Soft Matter, 9(37):9008-9015.

[Krüger et al., 2017] Krüger, T., Kusumaatmaja, II., Kuzmin, A., Shardt, O., Silva, G., and Viggen, E. M. (2017). The lattice boltzmann method. Springer International Publishing, 10(978-3) :4—15.

[Lai and Peskin, 2000] Lai, M.-C. and Peskin, C. S. (2000). An immersed boundary method with formal second-order accuracy and reduced numerical viscosity. Journal of computational Physics, 160(2):705-719.

[Lamb, 1932] Lamb, H. (1932). Hydrodynamics dover publications inc. New York, 738.

[Landau and Lifshitz, 1960] Landau, L. and Lifshitz, E. (1960). Fluid mechanics (course of theoretical physics, v. vi).

[Latt et al., 2020] Latt, J., Malaspinas, O., Kontaxakis, D., Parmigiani, A., Lagrava, D., Brogi, F., Belgacem, M. B., Thorimbert, Y., Leclaire, S., Li, S., et al. (2020). Palabos: Parallel lattice boltzmann solver. Computers & Mathematics with Applications.

[Levitov and Shitov, 2003] Levitov, L. and Shitov, A. (2003). Green's functions, problems and solutions.

[Li et al., 2007] Li, N., Kamei, D. T., and Ho, C.-M. (2007). On-chip continuous blood cell subtype separation by deterministic lateral displacement. In 2007 2nd IEEE International Conference on Nano/Micro Engineered and Molecular Systems, pages 932-936. IEEE.

[Liu et al, 2012] Liu, B., Goree, J., and Feng, Y. (2012). Waves and instability in a one-dimensional microfluidic array. Physical Review E, 86(4):046309.

[Malet-Engra et al., 2015] Malet-Engra, G., Yu, W., Oldani, A., Rey-Barroso, J., Gov, N. S., Scita, G., and Dupre, L. (2015). Collective cell motility promotes chemotactic prowess and resistance to chemorepulsion. Current Biology, 25(2):242-250.

[of Geneva, 2019] of Geneva, U. (2019). Palabos user guide.

[Oliver, 1962] Oliver, D. (1962). Influence of particle rotation on radial migration in the poiseuille flow of suspensions. Nature, 194(4835):1269—1271.

[Raven and Marmottant, 2009] Raven, J. P. and Marmottant, P. (2009). Microfluidic crystals: Dynamic interplay between rearrangement waves and flow. Physical Review Letters, 102(8): 1-4.

[Rayleigh, 1879] Rayleigh, L. (1879). On the capillary phenomena of jets. Proc. R. Soc. London, 29(196-199):71—97.

[Rubinow and Keller, 1961] Rubinow, S. I. and Keller, J. B. (1961). The transverse force on a spinning sphere moving in a viscous fluid. Journal of Fluid Mechanics, 11(3):447—459.

[Saffman, 1965] Saffman, P. G. (1965). The lift on a small sphere in a slow shear flow. Journal of fluid mechanics, 22(2):385-400.

[Schaaf and Stark, 2020] Schaaf, C. and Stark, H. (2020). Particle pairs and trains in inertial microfluidics. The European Physical Journal E, 43:1-13.

[Schiller et al., 2015] Schiller, U. D., Fleury, J.-B., Seemann, R., and Gompper, G. (2015). Collective waves in dense and confined microfluidic droplet arrays. Soft Matter, 11(29):5850—5861.

[Segre and Silberberg, 1961] Segre, G. and Silberberg, A. (1961). Radial particle displacements in poiseuille flow of suspensions. Nature, 189(4760) :209-210.

[Shan and Chen, 1993] Shan, X. and Chen, H. (1993). Lattice Boltzmann model for simulating flows with multi phases and components. Physical Review E, 47(3):1815—1819.

[Shchur and Guskova, 2020] Shchur, L. and Guskova, M. (2020). Drop oscillation modeling. Communications in Computer and Information Science, 1331:198-206.

[Succi, 2018] Succi, S. (2018). The lattice Boltzmann equation: for complex states of flowing matter. Oxford University Press.

[Suzuki and Inamuro, 2011] Suzuki, K. and Inamuro, T. (2011). Effect of internal mass in the simulation of a moving body by the immersed boundary method. Computers & Fluids, 49(1):173-187.

[Trinh and Wang, 1982] Trinh, E. and Wang, T. (1982). Large-amplitude free and driven drop-shape oscillations: experimental observations. Journal of Fluid Mechanics, 122:315-338.

[Uspal and Doyle, 2012] Uspal, W. E. and Doyle, P. S. (2012). Collective dynamics of small clusters of particles flowing in a quasi-two-dimensional microchannel. Soft Matter, 8(41) :10676-10686.

[Virtanen et al., 2020] Virtanen, P., Gommers, R., Oliphant, T. E., Haberland, M., Reddy, T., Cournapeau, D., Burovski, E., Peterson, P., Weckesser, W., Bright, J., et al. (2020). Scipy 1.0: fundamental algorithms for scientific computing in python. Nature methods, 17(3) :261—272.

[Wang et al., 2008] Wang, Z., Fan, J., and Luo, K. (2008). Combined multi-direct forcing and immersed boundary method for simulating flows with moving particles. International Journal of Multiphase Flow, 34(3):283-302.

[Yang et al., 2005] Yang, B. H., Wang, J., Joseph, D. D., Hu, H. H., Pan, T.-W., and Glowinski, R. (2005). Migration of a sphere in tube flow. Journal of Fluid mechanics, 540:109-131.

[Zou and He, 1997] Zou, Q. and He, X. (1997). On pressure and velocity boundary conditions for the lattice boltzmann bgk model. Physics of fluids, 9(6):1591—1598.