Введение. Аппараты погружного горения (АПГ) представляют собой теплотехнические устройства для прямого нагрева и выпаривания растворов. Состав нагревается при непосредственном контакте со струей горячих дымовых газов, образующихся в камере сгорания. Преимущество такого способа нагрева — отсутствие теплопередающих поверхностей, на которых выпадают осадки солей. При работе АПГ струя дымовых газов контактирует с жидкостью и разбивается на пузырьки, которые барботируют через слой жидкости и отдают ей тепло. В жидкости возникает турбулентный поток из-за высокоскоростного ввода струи дымовых газов. Значит, течения в аппарате рассчитываются по одной из моделей турбулентности. Известно несколько классов таких моделей [1].

На рис. 1 приводится укрупненная схема классификации турбулентных моделей. Большинство из них относятся к классу RANS — Reynolds averaged Navier – Stokes (осредненные по Рейнольдсу уравнения Навье – Стокса). Суть осреднения заключается в замене истинной скорости на осредненную по ансамблю скорость [2]:

 (1)

Здесь v — истинная скорость; ⟨⟩ — осредненная скорость; а v′ — пульсации скорости. Уравнение Навье – Стокса принимает вид:

 (2)

  (3)

Уравнение (3) представляет собой тензор напряжений Рейнольдса. Уравнение (2) останется незамкнутым до тех пор, пока не определятся компоненты тензора. Модели класса RANS определяют тензор напряжений Рейнольдса.

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

– без ДУ (модели нулевого порядка, описанные в работах Прандтля [4] и Кармана [5]);

– с одним ДУ (Спаларт – Аллмарас) [6];

– двумя и более ДУ (k – ε [7], k – ω [8], k – ε – n² [9]).

Альтернативное моделирование турбулентности в рамках RANS предполагает отход от гипотезы Буссинеска [10]. С математической точки зрения напряжения Рейнольдса представляют собой вторые одноточечные моменты пульсаций скорости, и уравнения переноса этих величин можно строго вывести из выражений Навье – Стокса:

 (4)

Здесь Dijk — диффузионный член, Pij — генерационный член, Fij —перераспределяющий член, eij — диссипативный член. Определение связи между членами уравнения (4) составляет суть проблемы построения моделей переноса напряжений Рейнольдса. Можно использовать набор алгебраических и дифференциальных уравнений.

Вихреразрешающие модели вместо процедуры осреднения напряжений реализуют процедуру фильтрации переменных [11]:

 (5)

Крупномасштабные турбулентные вихри разрешаются напрямую с помощью отфильтрованных переменных. Мелкомасштабные вихри, размер которых меньше размера ячейки расчетной сетки, разрешаются с применением подсеточных моделей. Наиболее известная из них — симуляция крупных вихрей (large eddy simulation, LES). В пределе, при стремящемся к нулю размере ячейки LES переходит в прямое численное моделирование турбулентности (direct numerical modeling, DNS), применяемое в настоящий момент лишь для фундаментальных исследований.

Комбинация RANS и разрешения вихрей — это база для класса гибридных моделей. Среди них отметим модель оторванных вихрей (detached eddy simulation, DES) [12], а также группу масштабно-адаптивных моделей (scale-adaptive simulation, SAS) [13]. Многообразие, с одной стороны, дает возможность гибкого выбора турбулентной модели, подходящей для конкретного расчетного случая, а с другой стороны, требует тщательного анализа и проверки моделей.

Для численных экспериментов используется абстрактное представление АПГ — ячейка с затопленной струей и свободной поверхностью [14]. Ее схема представлена на рис. 2.

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

Рекомендации по описанию этих факторов противоречат друг другу1. Предлагается многоскоростной континуум для моделирования пузырькового потока и односкоростной для моделирования свободной поверхности. Локально нестабильные потоки требуют применения вихреразрешающих моделей турбулентности, а ламинарные и турбулентные области — разрешения ламинарно-турбулентного перехода. Исследователи газлифтных [15] и барботажных [16] систем применяют модель k – ε, однако в таких системах отсутствует скоростной ввод газожидкостной струи. Таким образом, весьма затруднительно сделать окончательный выбор турбулентной модели для рассматриваемой ячейки, основываясь лишь на предыдущих публикациях. В то же время национальная экономика нуждается в надежной технике для выпаривания концентрированных рассолов. Такие рассолы — это отходы производства калийных удобрений либо промежуточный продукт при подземном выщелачивании солей.

Материалы и методы. Для анализа рассматривались следующие модели турбулентности:

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

Рассматриваемые модели реализуют разные подходы к моделированию турбулентности:

– осреднение Рейнольдсовых напряжений с использованием гипотезы Буссинеска (k – ε, k – ω, SST) и без нее (EARSM);

– фильтрация коротковолновых колебаний скорости (LES и DES).

Кроме того, эти модели внедрены в вычислительные пакеты: CFX, Fluent, FlowVision, «Логос».

Уравнения рассматриваемых турбулентных моделей включались в общую физико-математическую модель исследуемой гидродинамической системы. Общая модель описывает многофазное течение, в котором несущая фаза — жидкость, а раздробленная — пузырьки газа и твердые частицы. Для расчета течения газовой фазы использовался подход «Эйлер – Эйлер» [18], а для твердой — «Эйлер – Лагранж» [19]. Межфазное взаимодействие вычисляется с помощью эмпирических соотношений:

На срезе сопла камеры сгорания вводятся пузырьки газа, для которых задается постоянный диаметр 3 мм. Плотность газа — переменная величина, что является следствием сильных изменений температуры газовой фазы вдоль линии тока.

Коалесценция пузырьков не рассчитывалась. В качестве газа задан воздух, так как теплофизические свойства дымовых газов близки к его свойствам при коэффициенте избытка воздуха 1,7. Место ввода твердой фазы расположено в гипоцентре камеры сгорания на некотором удалении от среза сопла. Удаление подобрано так, что область ввода находится внутри газожидкостной струи, в зоне максимальной теплонапряженности. Именно в этом месте в реальном процессе испаряется жидкость и, следовательно, выпадет осадок. Открытая граница позволяет жидкости выходить из расчетной области, а газу проникать в нее и частично выходить. На данном участке границы фактически задаются параметры газовой фазы за расчетной областью: все стенки — непроницаемые и адиабатические плюс условия прилипания. Дополнительно для нижней стенки (рис. 1) вводится условие полной абсорбции твердой фазы. Таким образом, имитируется осаждение твердых частиц на дно аппарата.

Уравнения физико-математической модели решались методом конечных объемов с помощью программного пакета для моделирования задач гидродинамики Ansys CFX 2019 R3.

Описание математической модели. Математическая модель ячейки основывается на уравнениях «Эйлер – Эйлер» и «Эйлер – Лагранж». Ниже представлены уравнения неразрывности, сохранения импульса и энергии:

 (6)

 (7)

  (8)

Тензор вязких напряжений из (7) и (8) вычисляется по соответствующим уравнениям упомянутых выше моделей турбулентности.

Закон движения для Лагранжевых частиц:

  (9)

Силы в (7) и (9):

В представленных уравнениях α — объемная доля; ρ — плотность, кг/м3; v — скорость, м/с; t — время, с; p — давление, Па; H — удельная энтальпия, Дж/кг; λ — теплопроводность, Вт/м∙К; T — температура, К; Q — количество тепла, поступившее из фазы j, Дж/м3; cp — теплоемкость при постоянном давлении, Дж/кг∙К; d — диаметр, м; m — масса, кг. Векторы , , , представляют собой ускорение свободного падения, м/с², удельную силу межфазного взаимодействия, Н/м3, удельную силу присоединенной массы, Н/м3 и силу тяжести, Н; τ — тензор вязких напряжений. Индексы i, j указывают на произвольные фазы. Индексы g, l, s отвечают газовой, жидкой и твердой фазам соответственно. CVM — коэффициент присоединенной массы, задан 0,5. CD — коэффициент сопротивления, вычисляемый по корреляции Шиллера – Науманна для твердых частиц или Ишии – Зубера для пузырьков газа.

К базовым уравнениям добавлены описанные ниже граничные и начальные условия.

Первое — срез сопла камеры сгорания (условия соответствуют опытной установке [23]):

Второе — открытая граница:

Третье — вход жидкости:

Четвертое — адиабатические стенки с учетом условия прилипания:

Пятое — начальные условия:

Свободная поверхность жидкости ограничена плоскостью z = 0,8 м.

Приняты следующие размеры расчетной области: длина и ширина — 1,3 м, высота — 1 м. Количество узлов расчетной сетки по трем измерениям — 130:130:100. Использовалась равномерная тетраэдрическая сетка. Для проверки независимости результатов расчетов от размера ячейки сетки проводились численные эксперименты с использованием сеток различной крупности. Результаты этих экспериментов представлены в таблице 1.

Как видно из таблицы 1, средняя энергия турбулентных пульсаций (ключевой параметр с точки зрения настоящего исследования) для сеток с ячейками 15 и 10 мм отличается на 0,001 Дж/кг, что составляет 2 % от полученных значений. Таким образом, для дальнейших расчетов можно использовать как сетку с ячейкой 15 мм, так и сетку с ячейкой 10 мм. Однако сетка с ячейкой 10 мм продемонстрировала бо́льшую устойчивость счета, поэтому именно ее приняли для дальнейшей работы.

Численные эксперименты проводились в стационарном и нестационарном режимах. Шаг по времени в нестационарном режиме — 0,1 с. В обоих режимах поочередно применялись одноконтинуальная (гомогенная) и многоконтинуальная (гетерогенная) постановки при моделировании турбулентности. Остальные физические величины вычислялись в гетерогенной постановке.

Результаты исследования. В результате вычислительных экспериментов получены следующие данные:

В таблицах тип данных, которые нельзя получить в ходе эксперимента, обозначили n/a. Символ х указывает на случаи, когда расчет не завершен из-за ошибки решателя.

На рис. 3 и 4 представлены поля скоростей жидкости с наложенными на них линиями тока пузырьков газа.

Как можно видеть, масса осевшей твердой фазы находится в диапазоне 1,18–1,38 кг, что говорит о низком разбросе значений. Разница между двумя постановками модели SST оказалась несущественной (6 %). Наименьшую разницу результатов дали гомогенная SST и DES модели (3 %), несмотря на различия в характере течения.

Обсуждение. Сравнение величин удельных турбулентных пульсаций (таблица 2) показывает, что все они имеют одинаковый порядок — десятые и сотые доли Джоулей на килограмм. При этом гомогенная постановка дает более высокие значения энергий пульсаций, чем гетерогенная.

Величины средних скоростей жидкости (таблица 3) также принадлежат к одному порядку — десятым долям метров в секунду. Разброс значений энергий пульсации и скоростей невелик.

Значения объема жидкости (таблица 4) в расчетной области принимают близкие значения для стационарного режима. Исключение — гетерогенная SST модель со значением объема на 0,04 м3 ниже среднего для остальных моделей. В нестационарном же режиме модели DES и гомогенная SST дают значения объема около 1,19 м3, тогда как остальные две — около 1,3 м3. Это указывает на низкую консервативность первых двух моделей, т. е. они могут давать неприемлемые результаты на длинных временах.

Распределение значений средней температуры (таблица 5) имеет иную тенденцию. Так, в нестационарных расчетах разница между двумя постановками SST модели минимальна (0,7 °С), тогда как разница между SST и вихреразрешающими моделями составляет несколько градусов. Для стационарного режима характерен разброс значений средней температуры. Минимум (28,2 °С) дает модель k – ω, максимум (58,7 °С) — гетерогенная RNG k – ε.

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

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

Стационарный режим расчетов оказался устойчивее нестационарного. Это можно объяснить накоплением численных ошибок в нестационарном решении [24]. Данный эффект особенно заметен для многофазных задач, где решается большее число уравнений (кратно количеству фаз). Таким образом, стационарный режим можно рекомендовать для моделирования установившихся режимов, когда зависимость переменных от времени не имеет значения.

Рассматривая поля скоростей, полученные в нестационарном режиме (рис. 3), мы можем обнаружить, что модели DES и LES выводят структуру потока с бо́льшим числом вихрей, чем SST. Свободная поверхность жидкости оказалась сильно деформированной в расчетах с моделями DES и SST в гомогенной постановке.

Рассматривая поля скоростей, полученные в стационарном режиме (рис. 4), мы видим, что ни одна из моделей не дает сильной деформации свободной поверхности, даже SST. С применением последней получается картина течения с чуть более сложным рельефом на свободной поверхности. Стоит отметить результат, полученный с помощью модели k – ω. Здесь свободная поверхность наиболее гладкая, а поле скоростей симметричное, что фактически напоминает эталонную картину движения в АПГ [25]. Остальные модели дали схожие между собой результаты как в гомогенной, так и в гетерогенной постановках (рис. 4 а–е, з): двухвихревой поток с преимущественным течением слева направо расчетной области.

Поля скоростей жидкости соответствуют областям, где объемная доля жидкости превышает 85 %. Это позволяет судить о состоянии свободной поверхности. На всех иллюстрациях мы видим одинаково оформленную область газожидкостной струи. Ее дальнобойность практически равна для всех рассматриваемых моделей и режимов вычислений. Диаметр струи на свободной поверхности варьируется незначительно, за исключением тех случаев, когда свободная поверхность оказалась сильно искривленной (рис. 3 а и 3 в). Векторы скорости жидкости направлены по касательной к линиям тока газовой фазы. Это совпадает с ранее полученными результатами исследования гидродинамической ситуации в АПГ [23].

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

Выше отмечалась зеркальная асимметрия всех полей скоростей, кроме полученного с помощью модели k – ω. Это объясняется несимметричной вихревой структурой потока (три вихря принципиально не могут быть симметричны относительно вертикальной плоскости сечения). Речь идет о применении вихреразрешающих моделей, поэтому рассматриваемая особенность обусловлена входным и выходным течениями, справа и слева соответственно.

Модель SST — это комбинация k – ω и k – ε, и вклад каждой части определяется весовым коэффициентом [26]. Качественное сходство результатов моделирования SST и k – ε объясняется превалирующим весом компоненты k – ε.

Модель EARSM использует иную процедуру осреднения Рейнольдсовых напряжений, а значит, может служить эталоном при сравнении моделей класса RANS. Согласно рис. 4з, структура потока, воспроизводимого моделью EARSM, близка к k – ε и SST, что подтверждает достоверность результатов, полученных с их помощью. Можно сделать вывод, что модель k – ω дает афизичный результат, так как в несимметричной задаче не может быть строгой симметрии поля скоростей.

Ранее отмечалось, что гомогенная модель DES воспроизводит максимальную асимметрию свободной поверхности жидкости. Полученный угол наклона (около 30º) не соответствует реальной гидродинамической ситуации в АПГ, поэтому авторы считают такой результат крайним случаем и полагают нецелесообразным применять гомогенную модель DES.

Говоря о надежности различных моделей, стоит выделить сходство большинства из них как в качественных, так и в количественных результатах. Назовем три крайних случая:

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

Заключение. Выполнено численное моделирование АПГ с использованием различных способов разрешения турбулентных пульсаций. Итоги численных экспериментов позволяют говорить о сходной картине течения всех моделей турбулентности. Выявлены совпадения по четырем признакам:

Все модели, кроме гетерогенной RNG k – ε, k – ω и DES, демонстрируют незначительный разброс таких интегральных параметров, как энергия турбулентных пульсаций, средние скорость, температура и объем жидкости. Те же модели воспроизводят несимметричность вынужденного потока жидкости в аппарате. Модель k – ω оказалась нечувствительной к несимметричности, а DES воспроизвела свободную поверхность с углом наклона 30º, что не соответствует гидродинамике АПГ.

Численные эксперименты показали невысокую устойчивость счета в задаче с одновременным наличием свободной поверхности и плотного пузырькового потока. Половина расчетов не завершились из-за ошибки решателя. При этом обнаружено, что в стационарном режиме устойчивость счета выше. Это позволяет рекомендовать стационарный режим расчета для тех случаев, когда изменение параметров во времени не представляет интереса. Так бывает, например, при исследовании стационарного режима или если не нужно учитывать накопление Лагранжевых частиц. Кроме того, можно рекомендовать SST при расчетах АПГ. Эта модель показала наибольшую устойчивость счета, а также позволяет разрешать ламинарно-турбулентный переход.

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

1. Ansys CFX Reference Guide. Release 2021 R1. Canonsburg: ANSYS Inc; 2021. 450 p. 2. Ansys CFX Tutorials. Release 2024 R2. Canonsburg: Ansys Inc., 2024. 1002 с.