Моделирование одно- и двухфазных течений бинарных и трехкомпонентных жидких сред

Во введении обоснована актуальность диссертационной работы,
определены цели и задачи исследования, продемонстрированы научная новизна
и достоверность результатов, описана их практическая и теоретическая
значимость.
Глава 1 посвящена исследованию динамики гетерогенных двухфазных
систем методом фазового поля.
В разделе 1.1 дается литературный обзор, в котором описывается история
развития и текущее состояние идей по изучению двухфазных систем с
межфазной границей. Приводится история развития теории фазового поля,
дается детальное описание данной теории, обосновываются ее преимущества по
сравнению с другими подходами.
Теория фазового поля – математическая модель представления
гетерогенной системы с межфазной границей, сама граница раздела при этом
имеет малую, но конечную толщину и определяется с помощью поля
концентрации. Уравнения теории фазового поля в безразмерном виде
записываются следующим образом1,2:

div = 0,(1)

11
+ ( ∙ ∇) = −∇ + ∇2 − ∇ ,(2)
Re
1
+ ( ∙ ∇) = ∇2 ,(3)
Pe
0
==− Cn∇2 .(4)

Здесь , , , – скорость, давление, концентрация и химический потенциал
соответственно, – полная функция свободной энергии, 0 – «классическая»
часть функции свободной энергии без учета капиллярных сил. Уравнения (1)-(4)
содержат следующие безразмерные параметры – числа Рейнольдса (Re), Пекле
(Pe), Маха (M), Кана (Cn):
∗ ∗
Re =, Pe =, =, Cn =,
∗ ∗ ∗ ∗ 2

где ∗ , ∗ , ∗ – характерные значения плотности, динамической вязкости и
химического потенциала, – характерный размер системы, – характерная
(максимальная) скорость, , – коэффициенты мобильности и капиллярности.

Система (1)-(4) включает в себя уравнения Навье-Стокса ((1)-(2)) для
несжимаемой изотермической жидкости, где в уравнении движения
присутствует дополнительный член, −1 ∇ , ответственный за вклад в скорость
капиллярных эффектов. При записи уравнений (3) учтено, что движущей силой
процесса диффузии является градиент химического потенциала, который
определяется функцией свободной энергии . Одной из основных особенностей
метода фазового поля является предположение, что функция свободной энергии
зависит не только от концентрации (состава системы), но и от градиентов
Cn
концентрации. Данная функция имеет следующий вид: = 0 ( ) + (∇ )2 .
Lowengrub J. and Truskinovsky L. Quasi-incompressible Cahn-Hilliard fluids and topological
transitions // Proc. R. Soc. Lond. A. – 1998. – Vol. 454. – Pp. 2617-2654.
Jacqmin D. Calculation of Two-Phase Navier–Stokes Flows Using Phase-Field Modeling //
Journal of Computational Physics. – 1999. – Vol. 155, no.1. – Pp. 96-127.
Cn
Градиентное слагаемое(∇ )2 отвечает за капиллярные эффекты, 0 ( ) –
«классическая» функция свободной энергии, зависящая только от концентрации.
Функция 0 ( ) обычно имеет вид двухъямного потенциала (см. рис. 1): в
рассматриваемой системе физические значения концентрации изменяются от –
0.5 до 0.5. Предельные значения ± 0.5 соответствуют чистым первой и второй
фазам, при этом данные значения могут и не являться равновесными.
Равновесные значения определяются минимумами функции 0 – подбирая вид
этой функции, можно задать термодинамику смешения. В диссертации в
большинстве случаев используется функция 0 = 2 + 4 (рис. 1 штриховая
кривая), где – параметр, пропорциональный температуре системы, он также
задает равновесные значения концентрации = √− /2. Значения > 0
соответствуют однофазной гомогенной системе (компоненты полностью
растворимы).

Рисунок 1. Разные варианты вида функции 0 ( ).

В разделе 1.2 с помощью метода фазового поля изучается динамика
вытеснения одной жидкости другой в капиллярных трубках. Рассматривается
несколько задач: 1) вытеснение несмешивающихся жидкостей в одиночном
капилляре, 2) вытеснение смешивающихся жидкостей в одиночном капилляре с
неравновесным капиллярным давлением, 3) вытеснение в матрице однородно
распределенных капилляров.
Капилляр моделируется двухмерным слоем, в начальный момент времени
одна жидкость насыщает капилляр, вторая прокачивается под постоянным
перепадом давления с левого конца капилляра.
Обнаружено, что возможны два качественно различных сценария
вытеснения (см. рис. 2). В случае малых капиллярных сил ( = 0.01) фронт
вытеснения имеет пальцеобразную форму, профиль скорости течения лишь
незначительно отклоняется от течения Пуазейля вблизи границы раздела. Во
втором случае ( = 0.0001) капиллярные эффекты более существенны, и
вытеснение происходит поршнеобразным способом.

Рисунок 2. Вытеснение одной жидкости другой в капилляре для
несмешивающихся систем, Cn = 4∙10-4 , Pe = 104 , Re = 1, = – 0.5 в момент
времени = 5, длина капилляра = 6; слева: = 0.01, справа: = 0.0001.

Ряд полученных в расчетах характеристик течения согласуется с
известными аналитическими выражениями. Расход жидкости согласуется с
формулой Уошборна: = Re ∙ (12 )−1 (Δ − ), где Δ = 8 /Re – заданный
перепад давления, – капиллярное давление, – длина капилляра, – расход
жидкости. Величина Δ и коэффициент «12» подобраны для удобства, чтобы в
предельном случае течения Пуазейля максимум скорости был равен 1.
Хорошее согласие с аналитическими формулами получается и для
величины капиллярного давления. Например, для случая на рис. 2 при = 0.0001
безразмерное капиллярное давление ≈ 28, полученное численно из поля
давления, близко к значению капиллярного давления ≈ 26, полученному из
закона Лапласа (при этом значения коэффициента поверхностного натяжения и
кривизны находятся численно с использованием отдельного алгоритма, который
не зависит от давления). Таким образом, теория фазового поля дает хороший
метод оценки капиллярного давления.
Результаты моделирования вытеснения одной жидкости другой в случае
смешивающихся систем (рис. 3) показали, что неравновесное капиллярное
давление для таких систем быстрее всего убывает на начальной стадии, однако
при этом не падает до нуля. Таким образом, при рассмотрении динамики
вытеснения в капиллярах явление неравновесного капиллярного давления
является важным и существенным, и данным эффектом нельзя пренебрегать.
Капиллярное давление является сдерживающим фактором для скорости
течения в капилляре, максимум скорости достигается в пределе однофазного
течения с нулевым капиллярным давлением (течение Пуазейля), что реализуется
при больших значениях параметра в функции свободной энергии 0 = 2 +
4 ( задает степень смешения компонентов, т.е. равновесное значение
концентрации). Минимальное физически релевантное значение = – 0.5, что
соответствует полностью несмешивающимся системам. Было обнаружено, что
максимальное значение капиллярного давления (минимальная скорость
вытеснения) достигается не при = – 0.5, а при значении , немного
превышающем данное значение (см., рис. 3: = – 0.3, бирюзовые линии).
Данный эффект объясняется следующим образом. Как следует из закона
Лапласа, капиллярное давление пропорционально произведению коэффициента
поверхностного натяжения на кривизну границы: = ∙ κ. С увеличением
величина коэффициента поверхностного натяжения уменьшается, мениск
сильнее вытягивается, что приводит к его заострению – кривизна поверхности
при этом возрастает. При определенных параметрах кривизна увеличивается
на бо́льшую величину, чем уменьшается коэффициент поверхностного
натяжения , и капиллярное давление, таким образом, уменьшается, что
приводит к более медленному вытеснению. При дальнейшем увеличении
значение начинает уменьшаться более существенно, и капиллярное давление
также уменьшается.

Рисунок 3, Зависимость от времени расхода (a), отношения капиллярного
давления к перепаду давлений между концами (b), коэффициента
поверхностного натяжения (c). Параметры: Cn = 10-4 , Re = 1, =
6.25∙10-5 , Pe = 300, = 40.

Рисунок 4, насыщенность среды,
размер матрицы капилляров:
2×2 (штриховая линия), 3×3
(штрихпунктирная), 4×4
(штрихпунктирная с двумя точками),
5×5 (пунктирная), 6×6 (сплошная).
Параметры: = −0.5, Cn = 10-4 ,
Re = 1, = 2.5∙10-4 , Pe = 104.

Моделирование для сети капилляров показало, что вытеснение в системе
из равномерно распределенных капиллярных трубок качественно повторяет все
особенности вытеснения в одиночном капилляре. В частности, обнаружено, что
уравнение Уошборна также хорошо согласуется с численными расчетами и в
случае сети капилляров. При увеличении размеров матрицы характеристики
течения сходятся к некоторым предельным значениям (рис. 4).

В разделе 1.3 посвящен неустойчивости Релея-Тейлора. Приводится обзор
состояния проблемы. В диссертации численно исследуется задача о
неустойчивости Релея-Тейлора при наличии смешивающейся границы раздела,
при этом система в начальный момент времени не находится в состоянии
термодинамического равновесия (начальное значение концентрации не является
равновесным). Изучается как влияние диффузии на развитие неустойчивости
Релея-Тейлора, так и влияние динамики развития неустойчивости на диффузию
(перемешивание компонентов).
На рис. 5 приведены зависимости от времени средней концентрации одной
из фаз (что характеризует степень смешения компонентов). Сплошной линией
показан случай, когда более тяжелая жидкость располагается поверх легкой,
после чего наблюдается развитие неустойчивости, и в конечном счете две фазы
меняются местами: легкая жидкость – вверху, тяжелая – внизу. Для сравнения
пунктирной линией показан чисто диффузионный процесс в случае, когда в
начальный момент легкая жидкость находится поверх тяжелой. Видно, что
развитие неустойчивости Релея-Тейлора способствует существенной
интенсификации процесса диффузии.

Рисунок 5. Зависимость от
времени средней концентрации в
одной из фаз. Черные линии –
Pe = 107 , синие линии – Pe =
105 . Сплошные линии
соответствуют задаче о
неустойчивости Релея-Тейлора,
пунктирные линии – чисто
диффузионному процессу.

Глава 2 посвящена исследованию конвекции трехкомпонентных смесей с
эффектом Соре.
В разделе 2.1 приводится литературный обзор, в котором описывается
состояние исследований конвекции однокомпонентных и многокомпонентных
смесей. Дается обоснование и актуальность задач, рассматриваемых в разделах
2.2 и 2.3 диссертации. Приводится описание модели и методов, используемых
при численном решении исследуемых задач.
В разделе 2.2 описываются результаты исследования устойчивости и
надкритических режимов конвекции трехкомпонентных смесей с эффектом
Соре в горизонтальном слое с заданным тепловым потока на верхней и нижней
границах. Рассматривалась характерная жидкая смесь со следующими
параметрами: число Прандтля Pr = 10, числа Шмидта 1 = 100, 2 = 1000 при
трех значениях отношения разделения одного из компонентов 1 : 1 = 0, – 0.4,
0.4.
Линейная задача устойчивости решалась методом пристрелки. Получены
карты устойчивости на плоскости параметров (Ra, Ψ), где Ra – число Релея,
которое, в частности, отвечает за величину нагрева слоя; Ψ – суммарное
отношение разделения, которое отвечает за эффект термодиффузии. Карты
устойчивости представлены на рис. 6: имеются как монотонные моды
неустойчивости (сплошные линии), так и колебательные моды (штриховые
линии).

1 = 0
1 = 0.4

Рисунок 6. Карта устойчивости
трехкомпонентной смеси на
плоскости параметров число Релея
(Ra) – суммарное отношение
1 = – 0.4разделения (Ψ). Сплошные линии –
монотонная длинноволновая
неустойчивость, штриховые линии –
колебательная длинноволновая
неустойчивость. Областям
неустойчивости соответствуют
серые зоны.

Т.П. Любимовой и Е.С. Садиловым получены аналитические формулы для
границ длинноволновых монотонной и колебательной неустойчивостей [5].
Результаты длинноволнового анализа и численных расчетов совпадают.
Численно обнаружено, что во всех случаях наиболее опасной является именно
длинноволновая неустойчивость. Случай 1 = 0 соответствует частному случаю
бинарной жидкости. Карты устойчивости трехкомпонентных смесей ( 1 ≠ 0)
подобны случаю бинарной системы, но наблюдается сдвиг полученных кривых
по оси Ψ на величину, приблизительно равную 1 .
В нелинейных расчетах рассматривалась вытянутая ячейка слоя с
отношением сторон 5:1. Найдено, что длинноволновая неустойчивость в
действительности реализуется лишь вблизи порога возникновения конвекции.
На рис. 7 представлена зависимость модуля максимального значения функции
тока от числа Релея для монотонной моды неустойчивости. Видно, что по мере
уменьшения нагрева области (уменьшение |Ra|) наблюдается серия перестроек
к режимам с меньшим числом вихрей.
Для колебательной неустойчивости обнаружено, что при небольших |Ψ|
колебания существуют лишь при значениях Ra, немногим превышающим
пороговое, тогда как при больших |Ψ| колебания наблюдались во всей области
выше порога возникновения конвекции.

Рисунок 7. Зависимости модуля максимального значения функции
тока от числа Релея, цифры в кружках соответствуют числу
вихрей. 1 = 0, Ψ = – 0.2, подогрев сверху (Ra < 0). В разделе 2.3 приводятся результаты исследования конвекции трехкомпонентной смеси толуол-метанол-циклогексан с массовыми долями компонентов 0.62/0.31/0.07. Работа выполнялась в рамках анализа и интерпретации эксперимента DCMIX-2 на Международной космической станции. Рассматривалась замкнутая область с отношением сторон 3:1 при нагреве сверху с приложенной разницей температур 6 ℃. В работе3 представлены результатыобработкиэкспериментальныхданныхпоизмерению коэффициентов Соре. Двумя различными исследовательскими группами получены два варианта коэффициентов Соре. В обоих случаях оказалось, что суммарное отношение разделения Ψ близко к нулю, однако в первом случае оно положительно Ψ = 0.0154, а во втором – отрицательно Ψ = – 0.0006. В диссертации исследованы оба варианта. В первом случае (Ψ = 0.0154) найдено, что в итоге система приходит к состоянию механического равновесия с линейными распределениями температуры и концентрации компонентов. При этом, однако, для концентрации 1 (толуол) наблюдается нетривиальное поведение, обусловленное конкуренцией процессов диффузии, перекрестной диффузии и термодиффузии: на начальном этапе (до ~ 20) толуол движется к верхней горячей стенке, затем «разворачивается» и движется вниз, после чего достигается состояние равновесия (см. рис. 8.) Рисунок 8. Распределения 1 (a) и 2 (b) с высотой концентраций в различные моменты времени, = 1.5. Во втором случае (Ψ = – 0.0006) обнаружено, что поведение системы зависит от амплитуды начальных возмущений, которая задавалась в виде значения завихренности в некоторой точке области. При небольших значениях амплитуды возмущений в системе реализуется сценарий, качественно близкий к тому, что и для случая Ψ = 0.0154. Однако, когда амплитуда возмущений превышает некоторое критическое значение, в системе развивается неустойчивость, в результате чего формируется одновихревое конвективное течение. В приложении к диссертации описаны параллельный алгоритм расчета на графических процессорах и его программная реализации на CUDA C. Mialdun A., Ryzhkov I., Khlybov O., Lyubimova T., Shevtsova V. Measurement of Soret coefficients in a ternary mixture of toluene–methanol–cyclohexane in convection-free environment // Journal of Chemical Physics. – 2018. – Vol. 148. – P. 044506. ЗАКЛЮЧЕНИЕ Проведено прямое численное моделирование вытеснения одной жидкости другой в капилляре под действием перепада давления. Найдено, что в зависимости от значений параметров могут реализоваться сценарии пальцеобразного или поршневого вытеснения. Получен эффективный способ оценки капиллярного давления. Рассчитано динамически меняющееся капиллярное давление при вытеснении одной жидкости другой, смешивающейся с первой, показывающее не изученный ранее эффект неравновесного поверхностного натяжения. Найдено, что на первых стадиях капиллярное давление стремительно падает, а в дальнейшем изменяется значительно медленнее. Неравновесное капиллярное давление – эффект, который может оказывать значительное влияние на интенсивность вытеснения. Понимание данного эффекта поможет улучшить управление динамикой вытеснения многофазных сред. Проведено моделирование вытеснения одной жидкости другой, несмешивающейся с первой, в матрице капилляров. Впервые проведен расчет для достаточно больших матриц (10×10), при которых характеристики течения сходятся к некоторым предельным значениям. Исследованы возникновение и развитие неустойчивости Релея-Тейлора в изотермической гетерогенной системе двух смешивающихся жидкостей. Найдено, что эволюция системы состоит из двух этапов: конвективного и диффузионного. Во время обоих этапов и диффузия, и конвекция оказывают существенное влияние на развитие неустойчивости. Исследованы возникновение и нелинейные режимы конвекции трехкомпонентных смесей с эффектом Соре в плоском горизонтальном слое с заданным тепловым потоком на границах. Обнаружено, что при всех рассмотренных значениях параметровнаиболее опаснойявляется длинноволновая монотонная или колебательная неустойчивости. Нелинейные расчеты показали, что длинноволновый характер неустойчивости, обнаруженный с помощью линейного анализа, реализуется при малых надкритичностях и сохраняется при дальнейшем увеличении нагрева. Для обратной ситуации, когда при проведении моделирования интенсивность нагрева постепенно уменьшается, обнаружено, что длинноволновый режим реализуется лишь вблизи порога возникновения конвекции. Исследована индуцированная эффектом Соре конвекция смеси толуол- метанол-циклогексан при значениях суммарного отношения разделения (Ψ), близких к нулю. Обнаружено, что при положительных значениях Ψ на начальных этапах поведение системы может сильно осложняться взаимодействием диффузии, перекрестной диффузии и термодиффузии; при этом общий концентрационный вклад в градиент плотности может приводить к возникновению конвективного течения. Устойчивость системы при отрицательных значениях суммарного отношения разделения зависит от амплитуды начального возмущения. Численные расчеты позволили убедиться в корректности методики обработки данных космического эксперимента DCMIX- 2. Рекомендации, перспективы дальнейшей разработки темы Следующим шагом при изучении конвекции многокомпонентных гомогенных смесей с заданным тепловым потоком на границах будет являться исследование трехмерных конвективных режимов. При исследовании гетерогенных систем методом фазового поля следующим этапом может стать рассмотрение неизотермических задач. Функция свободной энергии в этом случае будет иметь более сложный вид и будет зависеть не только он концентрации, но и от температуры среды, что, по-видимому, повлечет за собой расширение математической модели теории фазового поля. Создание такой теории позволит в рамках единой математической модели решать широкий класс задач.