1. На главную
  2. Электронные версии
  3. Теплоэнергетика
  4. № 1, 2023

Теплоэнергетика, 2023, № 1, стр. 75-86

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

В. И. Мелихов a, b, *, О. И. Мелихов a, b, Г. Ю. Волков a, С. Е. Якуш c, Б. Салех a

a Национальный исследовательский университет “Московский энергетический институт”
111250 Москва, Красноказарменная ул., д. 14, Россия

b Электрогорский научно-исследовательский центр по безопасности АЭС
142530 Московская обл., г. Электрогорск, ул. Святого Константина, д. 6, Россия

c Институт проблем механики им. А.Ю. Ишлинского РАН
119526 Москва, просп. Вернадского, д. 101, корп. 1, Россия

* E-mail: vladimir.melikhov@erec.ru

Поступила в редакцию 25.05.2022
После доработки 30.05.2022
Принята к публикации 28.06.2022

Полный текст (PDF)

Аннотация

Рассмотрено проникновение струи жидкости в пространство, затопленное другой жидкостью, в условиях большой разности плотностей этих жидкостей применительно к фрагментации и перемешиванию струй при различных сценариях аварий на атомных электростанциях (истечение расплава кориума в воду или воды в заполненный жидким свинцом объем). Решение уравнений гидродинамики для нескольких несмешивающихся жидкостей с резкими поверхностями раздела осуществляется методом объема жидкости VOF (Volume of Fluid). Приведены математические модели многофазных течений, реализованные численно в рамках программы OpenFOAM с открытым исходным кодом. Проведена серия верификационных и исследовательских расчетов в двумерной (плоской) постановке, результаты сопоставлены с литературными данными. Выполненные с помощью такого подхода расчеты проникновения струи расплава металла Вуда в глубокий резервуар с водой выявили последовательные стадии эволюции фрагментирующей струи. Получено хорошее качественное и количественное соответствие с опубликованными в литературе результатами расчетов данного процесса с помощью кода ANSYS Fluent. Исследованы закономерности перемешивания струи воды, втекающей в объем, заполненный первоначально неподвижным жидким (расплавленным) свинцом. Показано, что в данном случае фрагментация струи воды происходит практически сразу же после ее входа в жидкий свинец. Исследовано влияние размера сеточной ячейки на характеристики перемешивания фаз. Приведены результаты расчетов характерных размеров образующихся фрагментов, их дисперсного состава и площади межфазной границы при перемешивании. Для обеих исследованных задач построены номограммы распределения по размерам образующихся при дроблении струи фрагментов. Показано, что получающиеся в результате такого процесса фрагменты расплава значительно меньше минимального предельного размера свободно падающих в воде капель расплава.

Ключевые слова: тяжелые аварии, фрагментация струи, метод VOF, численное моделирование, жидкий свинец, гидродинамическая неустойчивость, межфазная поверхность, перемешивание жидкостей, дисперсный состав

Многие задачи обеспечения безопасности атомной энергетики связаны с истечением жидкостей в пространство, заполненное жидкостью с существенно различающимися свойствами (плотность, температура), их взаимодействием, фрагментацией, возникновением многофазной зоны. Так, при тяжелых авариях на реакторах с водяным теплоносителем (ВВЭР, PWR) расплав активной зоны может истекать в виде струи во внутрикорпусное или внекорпусное пространство, заполненное водой [1, 2]. В то же время при аварии, связанной с разрывом теплообменной трубки парогенератора реактора со свинцовым теплоносителем, может происходить истечение воды, находящейся под давлением 20–25 МПа, в пространство, заполненное расплавом свинца при значительно более низком давлении (примерно 1 МПа) [3, 4]. При этом струя истекающей воды быстро дробится, создавая в окружающем пространстве объем, заполненный смесью свинца, пара и воды, в которой потенциально может происходить паровой взрыв [5]. Истечение высокоскоростной струи также может представлять опасность – воздействовать на соседние трубки парогенератора и создавать угрозу развития (эскалации) аварийной ситуации [6].

Формирование многофазной области при струйном истечении жидкости в затопленное пространство – это сложный гидродинамический процесс, в котором проявляются неустойчивости различного типа (Рэлея – Тейлора, Кельвина – Гельмгольца) на границе раздела сред, приводящие к быстрому развитию возмущений, усложнению формы межфазной границы, каскадной фрагментации истекающей жидкости. Еще более сложные процессы связаны со вскипанием жидкости, что добавляет теплофизические аспекты в проблему перемешивания жидкостей. Это ограничивает возможности аналитического описания протекающих процессов, определяя главенствующую роль экспериментальных исследований и численного моделирования.

В настоящее время имеется несколько основных методов моделирования движения несмешивающихся фаз, отделенных одна от другой поверхностями раздела. Базирующиеся на лагранжевом подходе так называемые бессеточные методы SPH (Smoothed Particle Hydrodynamics) [7] и MPS (Moving Particle Semi-implicit) [8] используются для описания движения среды со свободной границей взаимодействующих между собой частиц, которые могут передвигаться в рамках наложенных на них посредством основных уравнений динамики сплошной среды связей. Одним из недостатков подобных методов является необходимость использования большого числа частиц для достижения приемлемой точности моделирования. Повысить точность возможно путем применения лагранжева метода для отслеживания границы раздела с эйлеровым подходом, при котором используется неподвижная сетка, для интегрирования уравнений движения жидкости. Развитие подобных методов началось с метода маркеров и ячеек MAC (Marker-in-Cell) [9]. На сегодняшний день наиболее популярны метод объема жидкости в ячейке VOF (Volume of Fluid) [10] и метод изоповерхностей (Level-set method) [11], в которых межфазная граница описывается неявно как поверхность нулевого уровня функции расстояния, причем эволюция последней описывается уравнением конвективного переноса.

Принципиально иной подход к моделированию многофазных течений с поверхностью раздела используется в методе решеточных уравнений Больцмана [12], который рассматривает течение жидкости как движение ансамбля псевдочастиц, имеющих некоторую функцию распределения по дискретным скоростям. В этом случае в двухфазном течении процесс разделения фаз и поверхностные эффекты на границе раздела определяются силами взаимодействия между молекулами, составляющими фазы, для описания которых предложено несколько разных способов. Поэтому в данном методе не требуется специальное выделение поверхности раздела, которая рассчитывается сквозным образом на размерах нескольких узлов решетки (аналогично методам сквозного счета ударных волн в газодинамике).

Следует сказать, что все упомянутые ранее методы расчета течений с поверхностями раздела имеют как свои достоинства, так и недостатки [13], поэтому их применимость к конкретным задачам требует тщательного предварительного анализа для выяснения адекватности используемого метода. Данная работа посвящена изучению возможностей метода объема жидкости в ячейке VOF, реализованного в числе входящих в состав кода OpenFOAM [14] расчетных модулей, для описания характеристик многофазной смеси, возникающей при дроблении струи жидкости в затопленном пространстве.

Метод VOF неоднократно применялся для анализа различных аспектов процесса дробления струи. В [15] с его помощью была решена задача о стационарном истечении струи жидкости в бассейн, заполненный другой жидкостью. Метод VOF, реализованный в коммерческом CFD-коде ANSYS Fluent, использовался для определения длины распада струи расплава в воде в зависимости от различных параметров (относительной скорости движения расплава, его плотности и вязкости и поверхностного натяжения) [16]. Изучался также спектр размеров образующихся капель. При проведенном в [17] численном моделировании проникновения струи легкой жидкости в более тяжелую жидкость с помощью кода ANSYS Fluent с использованием метода VOF были успешно предсказаны наблюдавшиеся в экспериментах закономерности: начальное движение струи с практически постоянной скоростью и последующее зависание переднего фронта струи на некоторой высоте.

Анализ фрагментации струи расплава в воде с учетом его затвердевания, проведенный в [18] с помощью метода VOF, показал хорошее совпадение длины распада струи, скорости продвижения ее переднего фронта с экспериментальными данными. Результаты применения кода OpenFOAM с использованием метода VOF для задач фрагментации расплава под воздействием внешних факторов (импульс давления или удар водяной струи) представлены в [19, 20].

В настоящей работе приведены математические модели многофазных течений, рассмотрена их численная реализация в рамках программного кода OpenFOAM с открытым исходным кодом [14]. Проведена серия верификационных расчетов, в ходе которых рассматривалось проникновение струи одной жидкости (металла Вуда) в пространство, заполненное другой жидкостью (водой). Использовалась двумерная (плоская) постановка задачи, результаты, полученные на различных сетках, сопоставлялись с литературными данными. Для системы жидкий свинец – вода рассмотрены гидродинамические аспекты формирования многофазной струи за счет развития неустойчивостей Рэлея – Тейлора и Кельвина – Гельмгольца на межфазной границе. Исследовано влияние размера сеточной ячейки на характеристики перемешивания фаз. Приведены результаты по площади межфазной границы и размерам фрагментов при перемешивании.

ПОСТАНОВКА ЗАДАЧИ

Рассматривается совместное течение двух несмешивающихся жидкостей, обладающих различными свойствами (плотность, вязкость) при струйном истечении одной из них (“активная” жидкость) в полупространство, заполненное второй жидкостью. Гидродинамическое взаимодействие жидкостей происходит на подвижной межфазной поверхности, которая может приобретать весьма сложную форму в процессе развития неустойчивостей по механизмам Рэлея – Тейлора и Кельвина – Гельмгольца. Для моделирования таких течений используется метод объема (VOF) [10], который будет рассмотрен далее.

Метод VOF относится к классу эйлеровых методов. Его ключевым элементом являются пространственные распределения для объемных долей каждой фазы, которые эволюционируют в едином поле скорости U. На резких межфазных границах между несмешивающимися жидкостями объемные доли фаз изменяются скачком от 0 до 1. В численной реализации, однако, межфазные границы размываются на толщину нескольких ячеек сетки, что позволяет рассматривать объемные доли как непрерывные функции и применять для их расчета соответствующие численные методы. Фактически в методе VOF совокупность несмешивающихся жидкостей заменяется на одну эффективную жидкость с переменными свойствами, зависящими от свойств и объемных долей фаз. Для объемных долей фаз при этом решаются соответствующие уравнения переноса с помощью специальных численных методов, не допускающих размытие межфазных границ.

Для рассматриваемой в настоящей работе задачи о перемешивании двух несжимаемых жидкостей система определяющих уравнений формулируется следующим образом. Объемную долю “активной” жидкости можно обозначить как ${{\alpha }_{j}},$ объемную долю жидкости в затопленном пространстве – ${{\alpha }_{a}}$ (далее те же индексы используются и для других характеристик), причем ${{\alpha }_{j}} + {{\alpha }_{a}} = 1.$ Плотность $\rho $ и коэффициент динамической вязкости $\mu $ эффективной жидкости в методе VOF определяются как

$\rho = {{\alpha }_{j}}{{\rho }_{j}} + {{\alpha }_{a}}{{\rho }_{a}};\,\,\,\,\mu = {{\alpha }_{j}}{{\mu }_{j}} + {{\alpha }_{a}}{{\mu }_{a}}.$

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

$\begin{gathered} \frac{{\partial {{\rho }_{k}}{{\alpha }_{k}}}}{{\partial t}} + \nabla \left( {{{\rho }_{k}}{{\alpha }_{k}}U} \right) = 0,\,\,\,\,k = j,\,\,\,\,a; \\ \frac{{\partial \rho U}}{{\partial t}} + \nabla \left( {\rho U \otimes U} \right) = - \nabla p + \nabla {{\tau }} + \rho g + {{F}_{\sigma }}, \\ \end{gathered} $
где t – время; $p$ – давление; τ = = $\mu \left[ {\nabla U + \nabla {{U}^{T}} - \tfrac{2}{3}{\mathbf{I}}\left( {\nabla U} \right)} \right]$ – тензор вязких напряжений; ${\mathbf{I}}$ – единичный тензор; $g$ – ускорение свободного падения (силы тяжести); ${{F}_{\sigma }}$ – сила поверхностного натяжения на границе раздела фаз.

В настоящей работе не учитываются тепловые эффекты при взаимодействии жидкостей, поэтому уравнение энергии не рассматривается, а свойства фаз (плотность и вязкость) считаются постоянными. В этом случае из уравнений неразрывности фаз и условия ${{\alpha }_{j}} + {{\alpha }_{a}} = 1$ следуют условие бездивергентности поля скорости $\nabla {\mathbf{U}} = 0$ и уравнение переноса для объемной доли активной жидкости

$\frac{{\partial {{\alpha }_{j}}}}{{\partial t}} + \left( {\nabla {\mathbf{U}}} \right){{\alpha }_{j}} = 0.$

Сила поверхностного натяжения, действующая на границе раздела фаз, в методе VOF преобразуется в эквивалентную объемную силу, сосредоточенную в пределах распределенного переходного слоя между фазами. Используется широко распространенная модель CSF (Continuum Surface Force) [21], согласно которой эта сила вычисляется как

${{F}_{\sigma }} = \sigma \kappa \nabla {{\alpha }_{j}};\,\,\,\,\kappa = - \nabla \left( {\frac{{\nabla {{\alpha }_{j}}}}{{\left| {\nabla {{\alpha }_{j}}} \right|}}} \right),$
где $\sigma $ – коэффициент поверхностного натяжения; $\kappa $ – кривизна межфазной поверхности, определенная через дивергенцию единичной нормали к поверхности, которая, в свою очередь, определяется через градиент объемной доли.

В настоящей работе рассмотрены задачи о струйном проникновении “активной” жидкости в затопленное другой жидкостью пространство. Использована плоская постановка задачи, что позволяет моделировать асимметричные течения и крупные вихри, возникающие из-за развития неустойчивости, даже при симметричных начальных и граничных условиях. В работе [22] показано, что результаты анализа задач о гидродинамическом взаимодействии двух жидкостей в плоской постановке позволяют делать качественные выводы об этом процессе и в трехмерном варианте.

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

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

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

РЕАЛИЗАЦИЯ МЕТОДА VOF В ПАКЕТЕ OpenFOAM

Для реализации описанной математической модели использовали открытую интегрируемую платформу для численного моделирования задач механики сплошных сред – OpenFOAM [14]. На момент выполнения работы использовали наиболее современную версию программы OpenFOAM-v2106.

Программный комплекс OpenFOAM содержит набор расчетных модулей, в которых реализованы различные математические модели механики сплошных сред. Язык реализации – C++, кроме того, в пакете OpenFOAM широко используются новые типы данных (классы) и макросы для упрощения операций с ними.

Основным расчетным модулем для двухфазных задач является interFoam, позволяющий проводить гидродинамический анализ методом VOF для двух несжимаемых жидкостей без учета фазовых переходов. Совместное решение уравнения неразрывности и нестационарного уравнения импульса для эффективной жидкости осуществляется методом PISO [23], основанным на схеме типа предиктор – корректор. В этом методе на этапе предиктора находится предварительное поле скорости, которое затем уточняется вместе с полем давления в ходе нескольких коррекций. Особенностью реализации метода VOF в пакете OpenFOAM является переход от статического давления к давлению за вычетом гидростатического столба, что позволяет “перенести” влияние сил плавучести только в область границы раздела и объединить силы плавучести и объемные силы поверхностного натяжения в один источник, содержащий градиент объемной доли.

Ключевую роль в реализации метода VOF играет выбор численного метода, обеспечивающего сохранение резких градиентов объемной доли вблизи поверхности раздела фаз. В расчетных модулях OpenFOAM применяется метод MULES, в котором вводится искусственная “сжимающая” скорость, создающая дополнительную “антидиффузию” вблизи границы раздела фаз, препятствующую размыванию границы.

В настоящей работе расчетный модуль interFoam использовали без изменений, однако его функциональность была расширена для того, чтобы можно было вычислять и текущую площадь межфазной поверхности (модифицированный расчетный модуль получил название interAreaFoam). Для этого сначала в каждой ячейке расчетной области вычисляли удельную площадь межфазной поверхности ${{s}_{V}} = {{{\text{d}}{{S}_{i}}} \mathord{\left/ {\vphantom {{{\text{d}}{{S}_{i}}} {{\text{d}}{{V}_{i}}}}} \right. \kern-0em} {{\text{d}}{{V}_{i}}}},$ равную отношению площади сечения i-й ячейки ${\text{d}}{{S}_{i}}$ изоповерхностью, на которой объемная доля активной жидкости равна ${{\alpha }_{j}} = 0.5,$ к объему этой ячейки ${\text{d}}{{V}_{i}}$ (при отсутствии пересечения ${{s}_{V}} = 0$), после чего полную площадь находили как объемный интеграл $S = \int_V {{{s}_{V}}{\text{d}}V} .$ Кроме того, в процессе расчета вычисляли текущий полный объем втекающей жидкости в расчетной области ${{V}_{j}} = \int_V {{{\alpha }_{j}}{\text{d}}V} .$ Две указанные характеристики позволяют оценить средний эквивалентный диаметр d фрагментов втекающей жидкости

(1)
$d = \frac{{4{{V}_{j}}}}{S}.$

В OpenFOAM все сетки считаются трехмерными, так что при плоской постановке задачи абсолютные значения $S$ и ${{V}_{j}}$ зависят от размера ячеек по третьему (неиспользуемому) координатному направлению, однако средний размер фрагментов (1) от этой величины не зависит.

Расчеты проводили на рабочей станции класса Intel Core i7 с использованием распараллеливания на восемь потоков, реализованного в стандартном варианте OpenFOAM.

РАСЧЕТЫ ДРОБЛЕНИЯ СТРУИ РАСПЛАВА В ВОДЕ

Для верификации описанных ранее моделей и численных схем в коде OpenFOAM-v2106 были проведены расчеты проникновения струи расплава металла Вуда в глубокий резервуар с водой. Расчеты проводили для тех же условий, что и в работе [16], где моделирование распада струи при проникновении в бассейн с водой выполняли коммерческим кодом ANSYS Fluent. Металл Вуда имеет низкую температуру плавления (68.5°С), поэтому взаимодействие его расплава с водой исследовали без учета вскипания воды.

Расчетная область имела высоту 0.5 м и ширину 0.075 м, ширина струи на входе в область составляла 0.005 м. Параметры веществ, применявшиеся в расчетах, представлены в табл. 1.

Таблица 1.

Параметры расчетов взаимодействия струи расплава с водой [16]

Параметр Вода Расплав (металл Вуда)
ρ, кг/м3 998.2 9700
μ, Па · с 0.001 0.00194
σ, Н/м 1.0

Различия в возможностях использованных кодов заключаются в том, что расчеты [16] проводили в ANSYS Fluent с адаптивным локальным измельчением сетки (минимальный размер ячейки составлял 0.125 мм), тогда как в данной работе расчеты в OpenFOAM проводили на фиксированных сетках с квадратными ячейками со сторонами 0.5, 0.25 и 0.125 мм (число ячеек составляло 150 000, 600 000 и 2 400 000 соответственно). Кроме того, полного совпадения результатов нельзя ожидать также вследствие стохастичности течений с неустойчивой межфазной границей. Поэтому речь может идти, в основном, о сравнении формы струи, проникающей в теплоноситель на различных этапах процесса.

Трудность количественного сопоставления также заключается в том, что в работе [16] не приведены точные моменты времени, в которые показаны пространственные распределения, демонстрирующие форму струи и распространение сорвавшихся с ее поверхности фрагментов. Поэтому для сопоставления были выбраны несколько моментов, в которые глубина проникновения лидирующего фронта струи h в обоих расчетах совпадает. При этом линейный масштаб рисунков из работы [16] восстанавливали, исходя из известной ширины струи (5 мм) и ширины расчетной области 75 мм. Для сопоставления был выбран вариант со скоростью струи на входе ${{U}_{j}} = 16$ м/с (число Вебера, построенное по плотности жидкости в резервуаре, ${\text{W}}{{{\text{e}}}_{a}} = {{{{\rho }_{a}}U_{j}^{2}{{D}_{j}}} \mathord{\left/ {\vphantom {{{{\rho }_{a}}U_{j}^{2}{{D}_{j}}} \sigma }} \right. \kern-0em} \sigma } = 1280$).

На рис. 1 представлены распределения объемной доли расплава, демонстрирующие форму струи при глубинах проникновения 38.5, 80.1 и 114.2 мм. На каждом рисунке показаны результаты расчета программой OpenFOAM на сетках с ячейками размером 0.5, 0.25 и 0.125 мм, а также результаты расчета на адаптивной сетке из работы [16]. Указанные глубины проникновения выбирали по положению переднего фронта струи на рисунках, представленных в работе [16]. Точные моменты времени для соответствующих рисунков в [16] не указаны, поэтому для сравнения выбирали моменты времени 3.42, 7.5 и 11.0 мс, в которые при расчете программой OpenFOAM на самой мелкой сетке достигалась указанная глубина проникновения, этим же моментам соответствуют и расчеты на более грубых сетках.

Рис. 1.

Форма струи расплава при глубине проникновения 38.5 мм (а–в), 80.1 мм (д–ж), 114.2 (и–л). (а–в), (д–ж), (и–л) – расчет по OpenFOAM; г, з, м – расчет по ANSYS Fluent [16]; размер ячейки, мм: а, д, и – 0.5 × 0.5; б, е, к – 0.25 × 0.25; в, ж, л – 0.125 × 0.125

Видно, что расчеты по OpenFOAM на сетке с ячейками размером 0.5 × 0.5 мм дают более гладкую форму струи и менее выраженные коротковолновые возмущения на ее поверхности. Расчеты на сетке с ячейками размером 0.25 × 0.25 мм в целом дают весьма схожую картину с решением, полученным в ANSYS Fluent на адаптивной сетке. Результаты расчетов на сетке с ячейками размером 0.125 × 0.125 мм визуально весьма близки к данным, представленным в [16]. Так, хорошее количественное согласие наблюдается для ширины грибообразной структуры, формирующейся у лидирующей кромки струи. Совпадают значения длины волны крупномасштабных возмущений, приводящих к разрыву сплошности ядра струи, особенно хорошо это видно на рис. 1, л в сравнении с рис. 1, м. Расчет на самой мелкой сетке дает выраженные дисперсные капли, отделившиеся от основной струи, а также достаточно протяженные лигаменты (вытянутые фрагменты), периодически отделяющиеся от основной струи.

С точки зрения моделирования взаимодействия расплава с водой применительно к проблеме предварительного перемешивания при паровых взрывах, большой интерес представляет оценка дисперсного состава фрагментов расплава. Такой анализ был осуществлен путем обработки изображений, полученных в расчетах на самой мелкой сетке (размер ячейки 0.125 × 0.125 мм). На рис. 2, а показан фрагмент поля объемной доли расплава, полученный в расчете программой OpenFOAM и соответствующий моменту времени 25 мс, когда фрагментация струи уже приняла квазистационарный характер. Был выбран участок, охватывающий основную часть области справа от когерентного ядра струи вплоть до боковой границы и содержащий свыше 800 фрагментов разного размера. Для повышения контраста изображение было преобразовано в двухцветный (черно-белый) формат. Затем с помощью программы ImageJ [24] были определены площади всех фрагментов, а по ним – эквивалентные диаметры кругов, площади которых равны площадям фрагментов. Весь получившийся спектр диаметров был разбит на восемь равных диапазонов с шагом 0.5 мм, после чего все фрагменты были распределены по этим диапазонам и образовали фракции. Число фрагментов в мелкодисперсных фракциях (${{d}_{{ef}}} < 2.5$ мм) варьировалось от 26 до 356, а в трех крупнодисперсных фракциях – от 5 до 15. На рис. 2, б представлена полученная таким образом функция распределения массовой доли фрагментов m по фракциям.

Рис. 2.

Фрагмент изображения поля объемной доли расплава, использованной для определения дисперсного состава фрагментов (а), и функция распределения массовой доли фрагментов m по их диаметру d (б)

Из рис. 2, б следует, что наибольшую массу фрагментированных частиц составляет фракция, минимальный диаметр фрагментов в которой равен 1.0 мм, а максимальный диаметр 1.5 мм (центральное значение этого диапазона равно 1.25 мм). Вторую по массе фракцию образуют фрагменты диаметром от 0.5 до 1.0 мм. Таким образом, почти половину массы фрагментированного расплава составляют мелкие фрагменты эффективным диаметром от 0.5 до 1.5 мм, что объясняется высокой скоростью струи расплава в воде.

Такую процедуру выполняли неоднократно для разных моментов времени поздней стадии эволюции струи и различных выборок фрагментов (от 500 до 1000). При этом получающиеся распределения фрагментов различаются незначительно, что свидетельствует об их статистической достоверности. Кроме того, важно отметить, что в методе VOF происходит размытие межфазной границы на одну-две ячейки сетки, поэтому для надежной идентификации частицы необходимо, чтобы ее размер был не меньше размера четырех–шести ячеек сетки. Для использовавшейся наиболее мелкой сетки это составляет 0.50–0.75 мм, так что разрешающая способность метода достаточна для получения представленного на рис. 2 распределения по размерам, за исключением наиболее мелких частиц, массовая доля которых невелика.

К сожалению, для рассматриваемого варианта параметров струи расплава в работе [16] не представлены сведения о дисперсном составе образующихся фрагментов, а приведена информация для более низкой входной скорости струи расплава (5 и 10 м/с). При сопоставлении полученного авторами распределения фрагментов по размерам, когда входная скорость струи равна 16 м/с, с упомянутыми ранее распределениями при скоростях 10 и 5 м/с выявляется единый тренд – при увеличении входной скорости струи возрастает доля мелкодисперсных фрагментов.

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

На рис. 3 для представленных на рис. 1 результатов расчетов программой OpenFOAM на сетке с ячейками размером 0.125 × 0.125 мм показана зависимость среднего диаметра фрагментов, вычисленного по формуле (1), от продолжительности процесса. Здесь же показана зависимость площади межфазной поверхности S, рассчитанной на единицу длины в поперечном к расчетной области направлении, от времени. Видно, что на начальном этапе проникновения струи с формированием грибообразной структуры у лидирующей кромки, когда фрагменты еще не начали формироваться, формально определенный по (1) диаметр фрагментов сначала растет и достигает максимума примерно 5 мм, что соответствует диаметру струи. После этого, когда уже началась интенсивная фрагментация струи расплава, диаметр фрагментов резко уменьшается и выходит на стационар через 10–15 мс. Средний размер фрагментов при этом составляет 1.2 мм, что хорошо согласуется с функцией распределения, представленной на рис. 2, б.

Рис. 3.

Диаметр фрагментов d (1) и площадь межфазной поверхности на единицу поперечной длины S (2) при расчете на сетке с ячейками размером 0.125 × × 0.125 мм во времени t при истечении металла Вуда в резервуар с водой

Интересно сопоставить размеры фрагментов, получающихся при втекании высокоскоростной струи расплава в воду, с предельным размером его фрагментов при свободном падении в воде. Установившееся падение такой массы определяется балансом силы тяжести с одной стороны и силы сопротивления и архимедовой силы – с другой. Для рассматриваемой плоской постановки задачи это условие запишется в виде

(2)
${{C}_{d}}\left( {\operatorname{Re} } \right)d\frac{{{{\rho }_{a}}{{u}^{2}}}}{2} = \frac{{\pi {{d}^{2}}}}{4}\left( {{{\rho }_{j}} - {{\rho }_{a}}} \right)g,$
где $u$ – скорость падающего цилиндрического фрагмента расплава; ${{C}_{d}}$ – коэффициент сопротивления; $\operatorname{Re} = {{ud{{\rho }_{a}}} \mathord{\left/ {\vphantom {{ud{{\rho }_{a}}} {{{\mu }_{a}}}}} \right. \kern-0em} {{{\mu }_{a}}}}$ – число Рейнольдса.

Если скорость такого фрагмента будет достаточно большой, то фрагментация продолжится. Критерий прекращения фрагментации определяется числом Вебера [25, 26]

(3)
${\text{We}} = \frac{{d{{\rho }_{a}}{{u}^{2}}}}{\sigma } = {\text{W}}{{{\text{e}}}_{*}},$
где ${\text{W}}{{{\text{e}}}_{*}}$ – критическое число Вебера для колебательного режима распада капли, составляющее по разным оценкам от $2\pi $ [25] до 8–12 [26].

При применении соотношений (2), (3) для рассмотренных значений параметров задачи минимальный предельный размер свободно падающих фрагментов расплава составляет ${{d}_{*}} \approx 6.5$ мм (в расчетах для определения Cd использовали диаграмму [27], аппроксимированную на разных участках полиномами второй и третьей степеней). Эта величина значительно больше, чем при исследованной ранее фрагментации струи расплава, что, в первую очередь, объясняется высокой скоростью струи. В то же время при подстановке в (3) входной скорости струи $u = {{U}_{j}}$ верхняя граница диаметра фрагментов равняется ${{d}_{*}} = {{12\sigma } \mathord{\left/ {\vphantom {{12\sigma } {\left( {{{\rho }_{a}}U_{j}^{2}} \right)}}} \right. \kern-0em} {\left( {{{\rho }_{a}}U_{j}^{2}} \right)}} \approx 0.05$ мм, что также весьма далеко от полученного в расчетах среднего диаметра фрагментов. Это свидетельствует о необходимости использования детальных моделей взаимодействия жидкостей для описания процессов фрагментации.

Таким образом, проведенные верификационные расчеты подтвердили правильность работы расчетного модуля interAreaFoam программы OpenFOAM применительно к моделированию распада струи расплава при взаимодействии с водой. Далее этот расчетный модуль используется для моделирования струйного течения в обратной конфигурации – попадании легкой жидкости в пространство, заполненное тяжелым расплавом.

ИСТЕЧЕНИЕ СТРУИ ВОДЫ В ЖИДКИЙ СВИНЕЦ

В этом разделе рассмотрены результаты расчетов с помощью модуля interAreaFoam истечения жидкости меньшей плотности (воды) в виде струи в пространство, заполненное тяжелой жидкостью (расплавом свинца). Параметры жидкостей представлены в табл. 2 (свойства воды соответствуют ее состоянию при более высокой температуре, чем в табл. 1). Цель расчетов состояла в установлении закономерностей перемешивания, связанных с его гидродинамическими механизмами, обусловленными развитием неустойчивости Рэлея – Тейлора и Кельвина – Гельмгольца на границах раздела фаз. Задачу рассматривали в чисто гидродинамической постановке, теплообмен между жидкостями и последующее вскипание воды на поверхности расплава на данном этапе исследований не моделировали.

Таблица 2.

Параметры расчетов взаимодействия струи воды с расплавом свинца

Параметр Вода Расплав (свинец)
ρ, кг/м3 897 10 417
μ, Па · с 0.00016 0.00185
σ, Н/м 0.5

Расчетная область имела длину 1.0 м, ширину 0.6 м, ширина струи составляла 12 мм. Струя входила в область из середины боковой стенки, противоположная боковая стенка была твердой. Верхняя и нижняя границы области были открытыми. Силу тяжести не учитывали. Расчеты проводили на сетке размером 2000 × 1200 мм (2.4 × × 106 ячеек), размер ячеек составлял 0.5 × 0.5 мм.

На рис. 4 показано истечение струи воды со скоростью 20 м/с. На протяжении первых 10 мс заполненный водой объем является практически симметричным, однако вследствие развития неустойчивости Кельвина – Гельмгольца на границе раздела фаз происходит сильное возмущение межфазной поверхности, хорошо видное на рис. 4, а (момент времени 30 мс). Дальнейшее развитие этой неустойчивости приводит к разрушению ядра струи на отдельные фрагменты (см. рис. 4, б). Истекающая вода достаточно быстро теряет свой импульс, передавая его значительно более тяжелому свинцу, вследствие чего остронаправленная струя воды в свинце не возникает. На более поздних стадиях процесса перемешивания форма струи становится несимметричной, от нее отрываются крупные фрагменты (см. рис. 4, г–е), течение приобретает нерегулярный (стохастический) характер, в результате чего форма межфазной поверхности значительно усложняется. В двухфазной области возникает значительное перемешивание фаз, в том числе крупномасштабными вихрями, однако приведенные на рис. 4 распределения объемной доли жидкости показывают, что большая часть воды оказывается сосредоточенной в довольно крупных образованиях.

Рис. 4.

Проникновение струи воды в жидкий свинец. t, мс: а – 30; б – 50; в – 100; г – 150; д – 200; е – 300

На рис. 5 показано распределение образовавшихся фрагментов по размерам (массовые доли дисперсных фракций с интервалом 1 мм) по информации, представленной на рис. 4, е, построенное таким же способом, что и распределение на рис. 2, б. Видно, что, в отличие от фрагментации струи расплава в воде, дробление струи воды в расплаве приводит к образованию более крупных фрагментов при сопоставимых скоростях струй. Самая большая по массе фракция имеет средний размер 4.5 мм, в то время как для струи расплава диаметр фракции составляет 1.25 мм. Другое отличие распределения, представленного на рис. 5, от распределения на рис. 2, б заключается в наличии нерегулярного “хвоста” в области крупных фрагментов. Такая картина связана с малым количеством фрагментов в таких фракциях (от одного до пяти). Это приводит к тому, что среднемассовый размер, вычисленный по всем фрагментам, равный 7.96 мм, почти в 2 раза превышает средний размер фрагментов самой большой фракции (4.5 мм).

Рис. 5.

Функция распределения по диаметру d массовой доли дисперсных фрагментов воды m, образовавшихся при дроблении струи воды в расплавленном свинце

На рис. 6 приведены зависимости от времени эффективного диаметра $d$ и площади межфазной поверхности S, полученные в данном расчете программой OpenFOAM. Видно, что на начальном этапе, когда поверхность раздела фаз близка к цилиндрической (с учетом плоской постановки задачи), эффективный диаметр растет со временем. После потери устойчивости границы вода разбивается на фрагменты, эффективный диаметр которых резко уменьшается, а через 0.05 с изменение диаметра становится более плавным. К концу расчета средний диаметр водяных капель достигает примерно 8 мм, что соответствует среднемассовому размеру, вычисленному по распределению на рис. 5.

Рис. 6.

Средний размер фрагментов $d$ (1) и площадь межфазной поверхности на единицу поперечной длины $S$ (2) при проникновении струи воды в жидкий свинец

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

В работе [28] был выполнен линейный анализ устойчивости цилиндрической струи воды, окруженной паровой пленкой, в расплавленном свинце. При этом значения параметров задачи в [28] были близки к параметрам, рассматриваемым в настоящей работе. В [28] было получено, что длина распада струи воды при малых толщинах паровой пленки составляет примерно 1 см, что соответствует численным результатам работы авторов. Однако характерный размер образующихся капель воды в [28] был около 1 мм, что на порядок меньше значений, наблюдающихся при численном моделировании. Скорее всего, различие не связано с размером ячеек, используемых в расчетах, поскольку возникающие возмущения и искривления поверхностей раздела достаточно хорошо разрешаются на сеточном уровне (пять – семь ячеек на характерный масштаб длины). Возможно, причина состоит в том, что сама постановка задачи для линейного анализа не полностью соответствует процессу, который исследовался численно, а именно исходное предположение о стационарной струе воды, движущейся в окружающем свинце, при численном анализе авторов не подтверждается – “наконечник” струи сразу же вступает во взаимодействие с тяжелым свинцом, что существенно деформирует гидродинамическую картину, предполагаемую при линейном анализе устойчивости.

ВЫВОДЫ

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

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

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

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

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

Список литературы

  1. Мелихов В.И., Мелихов О.И., Якуш С.Е. Термическое взаимодействие высокотемпературных расплавов с жидкостями // ТВТ. 2022. Т. 60. № 2. С. 1–39.

  2. Мелихов В.И., Мелихов О.И., Якуш С.Е. Гидродинамика и теплофизика паровых взрывов. М.: Изд-во ИПМех РАН, 2020.

  3. Безносов А.В., Бокова Т.А. Оборудование энергетических контуров с тяжелыми жидкометаллическими теплоносителями в атомной энергетике. Н. Новгород: Нижегород. гос. техн. ун-т им. Р.Е. Алексеева, 2012.

  4. Review of thermal-hydraulic issues and studies of lead-based fast reactors / Y. Zhang, C. Wang, Z. Lan, S. Wei, R. Chen, W. Tian // Renewable Substainable Energy Rev. 2020. V. 120. P. 109625. https://doi.org/10.1016/j.rser.2019.109625

  5. Iskhakov A.S., Melikhov V.I., Melikhov O.I. Hugoniot analysis of energetic molten lead-water interaction // Annals Nucl. Energy. 2019. V. 129. P. 437–449. https://doi.org/10.1016/j.anucene.2019.02.018

  6. Steam generator tube rupture in lead-cooled fast reactors: Estimation of impact on neighboring tubes / A.S. Iskhakov, V.I. Melikhov, O.I. Melikhov, S.E. Yakush // Nucl. Eng. Des. 2019. V. 341. P. 198–208. https://doi.org/10.1016/j.nucengdes.2018.11.001

  7. Monaghan J.J. Simulating free-surface flows with SPH // J. Comput. Phys. 1994. V. 110. P. 399–406.

  8. Koshizuka S., Oka Y. Moving-particle semi-implicit method for fragmentation of incompressible fluid // Nucl. Sci. Eng. 1996. V. 123. P. 421–434.

  9. Harlow F., Welch J.E. Numerical calculation of time-dependent viscous incompressible flow of fluid with a free surface // Phys. Fluids. 1965. V. 8. P. 2182–2189.

  10. Hirt C.W., Nichols B.D. Volume of Fluid (VOF) method for the dynamics of free boundaries // J. Comput. Phys. 1981. V. 39. P. 201–225.

  11. Olsson E., Kreiss G. A conservative level set method for two phase flow // J. Comput. Phys. 2005. V. 210. P. 225–246.

  12. Saito S., Abe Y., Koyama K. Lattice Boltzmann modeling and simulation of liquid jet breakup // Phys. Rev. 2017. V. 96. № 1. P. 013317.

  13. Soligo G., Roccon A., Soldati A. Turbulent flows with drops and bubbles: What numerical simulations can tell us – Freeman Scholar Lecture // J. Fluids Eng. 2021. V. 143. № 8. P. 080801.

  14. OpenFOAM: The open source CFD toolbox. https://www.openfoam.com/

  15. Richards J.R., Beris A.N., Lenhoff A.M. Steady laminar flow of liquid-liquid jets at high Reynolds numbers // Phys. Fluids A 5. 1993. V. 7. P. 1703–1717.

  16. Thakre S., Manickam L., Ma W. A numerical simulation of jet breakup in melt coolant interactions // Ann. Nucl. Energy. 2015. V. 80. P. 467–475.

  17. Experimental and computational study on jet final depth under coolant jet-melt interaction / J. Chen, Y. Zhou, M. Zhong, J. Wang, H. Gong, Q. Huang // J. Nucl. Sci. Technol. 2017. V. 54. No. 9. P. 1002–1010.

  18. Numerical simulation of metal jet breakup, cooling and solidification in water / Y. Zhou, J. Chen, M. Zhong, J. Wang, M. Lv // Int. J. Heat Mass Transfer. 2017. V. 109. P. 1100–1109.

  19. Numerical modelling of melt droplet interaction with water / S.E. Yakush, N.S. Sivakov, V.I. Melikhov, O.I. Melikhov // J. Phys. Conf. Ser. 2021. V. 2057. P. 012057.

  20. Modelling of water jet impact on molten metal / S.E. Yakush, N.S. Sivakov, V.I. Melikhov, O.I. Melikhov // J. Phys. Conf. Ser. 2021. V. 2119. P. 012073.

  21. Brackbill J., Kothe D., Zemach C. A continuum method for modeling surface tension // J. Comput. Phys. 1992. V. 100. P. 335–354.

  22. Duan R.-Q., Koshizuka S., Oka Y. Two-dimensional simulation of drop deformation and breakup at around the critical Weber number // Nuc. Eng. Des. 2003. V. 225. P. 37–48.

  23. Issa R.I. Solution of the implicitly discretised fluid flow equations by operator-splitting // J. Comput. Phys. 1986. V. 62. P. 40–65.

  24. ImageJ User Guide. https://imagej.nih.gov/ij/docs/ guide/user-guide.pdf

  25. Нигматулин Р.И. Динамика многофазных сред. Ч. I. М.: Наука, 1987.

  26. Pilch M., Erdman C.A. Use of breakup time data and velocity history data to predict the maximum size of stable fragments for acceleration-induced breakup of a liquid drop // Int. J. Multiphys. Flow. 1987. V. 13. P. 741–757.

  27. Идельчик И.Е. Справочник по гидравлическим сопротивлениям / под ред. М.О. Штейнберга. 3-е изд., перераб. и доп. М.: Машиностроение, 1992.

  28. Finoshkina D.V., Melikhov O.I., Melikhov V.I. Evaluation of melt-water premixture formation due to hydrodynamic instabilities // Proc. of the 30th Intern. Conf. “Nuclear Energy for New Europe”. Bled, Slovenia, 6–9 Sept. 2021. Društvo jedrskih strokovnjakov Slovenije (Nuclear Society of Slovenia), Ljubljana, Slovenia, 2021. P. 412.1–412.8.

Дополнительные материалы отсутствуют.

Инструменты

следующая статья выпуска предыдущая статья выпуска содержание выпуска

Теплоэнергетика

Архивы выпусков Информация о журнале Отправить рукопись в журнал